Changeset 5844866 in flexpart.git

Timestamp:

Oct 12, 2016, 1:14:19 PM (8 years ago)

Author:

Sabine <sabine.eckhardt@…>

Branches:

master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie

Children:

df4a68e

Parents:

0581cac

Message:

first working version with backward deposition

Location:

src

Files:

• advance.f90 (modified) (1 diff)
• advance_rec.f90 (modified) (6 diffs)
• interpol_vdep.f90 (modified) (1 diff)
• par_mod.f90 (modified) (1 diff)
• releaseparticles.f90 (modified) (1 diff)
• timemanager.f90 (modified) (5 diffs)
• wetdepo.f90 (modified) (6 diffs)

src/advance.f90

@@ -569 +569 @@
                prob(ks)=1.+(prob(ks)-1.)* &
                     exp(-vdepo(ks)*abs(dt)/(2.*href))
+!              write(*,*) 'Prob calc: ',zt,href,prob(ks)
           endif
         end do

src/advance_rec.f90

@@ -20 +20 @@
 !**********************************************************************
 
-subroutine advance_rec(itime,nrelpoint,ldt,up,vp,wp, &
-       usigold,vsigold,wsigold,nstop,xt,yt,zt,prob,icbt)
-  !                     i    i  i/oi/oi/o
-  !  i/o     i/o     i/o     o  i/oi/oi/o i/o  i/o
+subroutine advance_rec(itime,xt,yt,zt,prob)
+  !                    i     i  i  i  o
   !*****************************************************************************
   !                                                                            *
@@ -66 +64 @@
   !                                                                            *
   ! Variables:                                                                 *
-  ! icbt               1 if particle not transferred to forbidden state,       *
-  !                    else -1                                                 *
-  ! dawsave            accumulated displacement in along-wind direction        *
-  ! dcwsave            accumulated displacement in cross-wind direction        *
-  ! dxsave             accumulated displacement in longitude                   *
-  ! dysave             accumulated displacement in latitude                    *
-  ! h [m]              Mixing height                                           *
-  ! lwindinterv [s]    time interval between two wind fields                   *
   ! itime [s]          time at which this subroutine is entered                *
   ! itimec [s]         actual time, which is incremented in this subroutine    *
   ! href [m]           height for which dry deposition velocity is calculated  *
-  ! ladvance [s]       Total integration time period                           *
   ! ldirect            1 forward, -1 backward                                  *
   ! ldt [s]            Time step for the next integration                      *
   ! lsynctime [s]      Synchronisation interval of FLEXPART                    *
   ! ngrid              index which grid is to be used                          *
-  ! nrand              index for a variable to be picked from rannumb          *
-  ! nstop              if > 1 particle has left domain and must be stopped     *
   ! prob               probability of absorption due to dry deposition         *
-  ! rannumb(maxrand)   normally distributed random variables                   *
-  ! rhoa               air density                                             *
-  ! rhograd            vertical gradient of the air density                    *
-  ! up,vp,wp           random velocities due to turbulence (along wind, cross  *
-  !                    wind, vertical wind                                     *
-  ! usig,vsig,wsig     mesoscale wind fluctuations                             *
-  ! usigold,vsigold,wsigold  like usig, etc., but for the last time step       *
   ! vdepo              Deposition velocities for all species                   *
   ! xt,yt,zt           Particle position                                       *
@@ -101 +81 @@
   use com_mod
   use interpol_mod
-  use hanna_mod
-  use cmapf_mod
-  use random_mod, only: ran3
 
   implicit none
 
   real(kind=dp) :: xt,yt
-  real :: zt,xts,yts,weight
-  integer :: itime,itimec,nstop,ldt,i,j,k,nrand,loop,memindnext,mind
-  integer :: ngr,nix,njy,ks,nsp,nrelpoint
-  real :: dz,dz1,dz2,xlon,ylat,xpol,ypol,gridsize
-  real :: ru,rv,rw,dt,ux,vy,cosfact,xtn,ytn,tropop
-  real :: prob(maxspec),up,vp,wp,dxsave,dysave,dawsave
-  real :: dcwsave
-  real :: usigold,vsigold,wsigold,r,rs
-  real :: uold,vold,wold,vdepo(maxspec)
-  !real uprof(nzmax),vprof(nzmax),wprof(nzmax)
-  !real usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
-  !real rhoprof(nzmax),rhogradprof(nzmax)
-  real :: rhoa,rhograd,delz,dtf,rhoaux,dtftlw,uxscale,wpscale
-  integer(kind=2) :: icbt
-  real,parameter :: eps=nxmax/3.e5,eps2=1.e-9
-  real :: ptot_lhh,Q_lhh,phi_lhh,ath,bth !modified by mc 
-  real :: old_wp_buf,dcas,dcas1,del_test !added by mc
-  integer :: i_well,jj,flagrein !test well mixed: modified by mc
-
-
-  integer :: idummy = -7
-  real    :: settling = 0.
-
-
-  nstop=0
-  do i=1,nmixz
-    indzindicator(i)=.true.
-  end do
-
+  real :: zt,xtn,ytn
+  integer :: itime,i,j,k,memindnext
+  integer :: nix,njy,ks
+  real :: prob(maxspec),vdepo(maxspec)
+  real,parameter :: eps=nxmax/3.e5
 
   if (DRYDEP) then    ! reset probability for deposition
@@ -144 +97 @@
     end do
   endif
-
-  dxsave=0.           ! reset position displacements
-  dysave=0.           ! due to mean wind
-  dawsave=0.          ! and turbulent wind
-  dcwsave=0.
-
-  itimec=itime
-
-  nrand=int(ran3(idummy)*real(maxrand-1))+1
 
 
@@ -207 +151 @@
 
 
-  ! Compute maximum mixing height around particle position
-  !*******************************************************
-
-  h=0.
-  if (ngrid.le.0) then
-    do k=1,2
-      mind=memind(k) ! eso: compatibility with 3-field version
-      do j=jy,jyp
-        do i=ix,ixp
-          if (hmix(i,j,1,mind).gt.h) h=hmix(i,j,1,mind)
-        end do
-      end do
-    end do
-    tropop=tropopause(nix,njy,1,1)
-  else
-    do k=1,2
-      mind=memind(k)
-      do j=jy,jyp
-        do i=ix,ixp
-          if (hmixn(i,j,1,mind,ngrid).gt.h) h=hmixn(i,j,1,mind,ngrid)
-        end do
-      end do
-    end do
-    tropop=tropopausen(nix,njy,1,1,ngrid)
-  endif
-
-  zeta=zt/h
-
-
-
-  !*************************************************************
-  ! If particle is in the PBL, interpolate once and then make a
-  ! time loop until end of interval is reached
-  !*************************************************************
-
-  if (zeta.le.1.) then
-
-  ! BEGIN TIME LOOP
-  !================
-
-    loop=0
-100   loop=loop+1
-      if (method.eq.1) then
-        ldt=min(ldt,abs(lsynctime-itimec+itime))
-        itimec=itimec+ldt*ldirect
-      else
-        ldt=abs(lsynctime)
-        itimec=itime+lsynctime
-      endif
-      dt=real(ldt)
-
-      zeta=zt/h
-
-
-      if (loop.eq.1) then
-        if (ngrid.le.0) then
-          xts=real(xt)
-          yts=real(yt)
-          call interpol_all(itime,xts,yts,zt)
-        else
-          call interpol_all_nests(itime,xtn,ytn,zt)
-        endif
-
-      else
-
-
-  ! Determine the level below the current position for u,v,rho
-  !***********************************************************
-
-        do i=2,nz
-          if (height(i).gt.zt) then
-            indz=i-1
-            indzp=i
-            goto 6
-          endif
-        end do
-6       continue
-
-  ! If one of the levels necessary is not yet available,
-  ! calculate it
-  !*****************************************************
-
-        do i=indz,indzp
-          if (indzindicator(i)) then
-            if (ngrid.le.0) then
-              call interpol_misslev(i)
-            else
-              call interpol_misslev_nests(i)
-            endif
-          endif
-        end do
-      endif
-
-
-  ! Vertical interpolation of u,v,w,rho and drhodz
-  !***********************************************
-
-  ! Vertical distance to the level below and above current position
-  ! both in terms of (u,v) and (w) fields
-  !****************************************************************
-
-      dz=1./(height(indzp)-height(indz))
-      dz1=(zt-height(indz))*dz
-      dz2=(height(indzp)-zt)*dz
-
-      u=dz1*uprof(indzp)+dz2*uprof(indz)
-      v=dz1*vprof(indzp)+dz2*vprof(indz)
-      w=dz1*wprof(indzp)+dz2*wprof(indz)
-      rhoa=dz1*rhoprof(indzp)+dz2*rhoprof(indz)
-      rhograd=dz1*rhogradprof(indzp)+dz2*rhogradprof(indz)
-
-
-  ! Compute the turbulent disturbances
-  ! Determine the sigmas and the timescales
-  !****************************************
-
-      if (turbswitch) then
-        call hanna(zt)
-      else
-        call hanna1(zt)
-      endif
-
-
-  !*****************************************
-  ! Determine the new diffusivity velocities
-  !*****************************************
-
-  ! Horizontal components
-  !**********************
-
-      if (nrand+1.gt.maxrand) nrand=1
-      if (dt/tlu.lt..5) then
-        up=(1.-dt/tlu)*up+rannumb(nrand)*sigu*sqrt(2.*dt/tlu)
-      else
-        ru=exp(-dt/tlu)
-        up=ru*up+rannumb(nrand)*sigu*sqrt(1.-ru**2)
-      endif
-      if (dt/tlv.lt..5) then
-        vp=(1.-dt/tlv)*vp+rannumb(nrand+1)*sigv*sqrt(2.*dt/tlv)
-      else
-        rv=exp(-dt/tlv)
-        vp=rv*vp+rannumb(nrand+1)*sigv*sqrt(1.-rv**2)
-      endif
-      nrand=nrand+2
-
-
-      if (nrand+ifine.gt.maxrand) nrand=1
-      rhoaux=rhograd/rhoa
-      dtf=dt*fine
-
-      dtftlw=dtf/tlw
-
-  ! Loop over ifine short time steps for vertical component
-  !********************************************************
-
-      do i=1,ifine
-
-  ! Determine the drift velocity and density correction velocity
-  !*************************************************************
-
-        if (turbswitch) then
-          if (dtftlw.lt..5) then
-  !*************************************************************
-  !************** CBL options added by mc see routine cblf90 ***
-            if (cblflag.eq.1) then  !modified by mc
-              if (-h/ol.gt.5) then  !modified by mc
-              !if (ol.lt.0.) then   !modified by mc  
-              !if (ol.gt.0.) then   !modified by mc : for test
-                  !print  *,zt,wp,ath,bth,tlw,dtf,'prima'
-                  flagrein=0
-                  nrand=nrand+1
-                  old_wp_buf=wp
-                  call cbl(wp,zt,ust,wst,h,rhoa,rhograd,sigw,dsigwdz,tlw,ptot_lhh,Q_lhh,phi_lhh,ath,bth,ol,flagrein) !inside the routine for inverse time
-                  wp=(wp+ath*dtf+bth*rannumb(nrand)*sqrt(dtf))*real(icbt) 
-                  ! wp=(wp+ath*dtf+bth*gasdev2(mydum)*sqrt(dtf))*real(icbt) 
-                  delz=wp*dtf
-                  if (flagrein.eq.1) then
-                      call re_initialize_particle(zt,ust,wst,h,sigw,old_wp_buf,nrand,ol)
-                      wp=old_wp_buf
-                      delz=wp*dtf
-                      nan_count=nan_count+1
-                  end if
-                  !print  *,zt,wp,ath,bth,tlw,dtf,rannumb(nrand+i),icbt
-                  !pause                  
-              else 
-                  nrand=nrand+1
-                  old_wp_buf=wp
-                  ath=-wp/tlw+sigw*dsigwdz+wp*wp/sigw*dsigwdz+sigw*sigw/rhoa*rhograd  !1-note for inverse time should be -wp/tlw*ldirect+... calculated for wp=-wp
-                                                                                      !2-but since ldirect =-1 for inverse time and this must be calculated for (-wp) and
-                                                                                      !3-the gaussian pdf is symmetric (i.e. pdf(w)=pdf(-w) ldirect can be discarded
-                  bth=sigw*rannumb(nrand)*sqrt(2.*dtftlw)
-                  wp=(wp+ath*dtf+bth)*real(icbt)  
-                  delz=wp*dtf
-                  del_test=(1.-wp)/wp !catch infinity value
-                  if (isnan(wp).or.isnan(del_test)) then 
-                      nrand=nrand+1                      
-                      wp=sigw*rannumb(nrand)
-                      delz=wp*dtf
-                      nan_count2=nan_count2+1
-                      !print *,'NaN coutner equal to:', nan_count,'reduce ifine if this number became a non-negligible fraction of the particle number'
-                  end if  
-              end if
-  !******************** END CBL option *******************************            
-  !*******************************************************************            
-            else
-                 wp=((1.-dtftlw)*wp+rannumb(nrand+i)*sqrt(2.*dtftlw) &
-                 +dtf*(dsigwdz+rhoaux*sigw))*real(icbt) 
-                 delz=wp*sigw*dtf
-            end if
-          else
-            rw=exp(-dtftlw)
-            wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2) &
-                 +tlw*(1.-rw)*(dsigwdz+rhoaux*sigw))*real(icbt)
-            delz=wp*sigw*dtf
-          endif
-          
-        else
-          rw=exp(-dtftlw)
-          wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2)*sigw &
-               +tlw*(1.-rw)*(dsigw2dz+rhoaux*sigw**2))*real(icbt)
-          delz=wp*dtf
-        endif
-
-  !****************************************************
-  ! Compute turbulent vertical displacement of particle
-  !****************************************************
-
-        if (abs(delz).gt.h) delz=mod(delz,h)
-
-  ! Determine if particle transfers to a "forbidden state" below the ground
-  ! or above the mixing height
-  !************************************************************************
-
-        if (delz.lt.-zt) then         ! reflection at ground
-          icbt=-1
-          zt=-zt-delz
-        else if (delz.gt.(h-zt)) then ! reflection at h
-          icbt=-1
-          zt=-zt-delz+2.*h
-        else                         ! no reflection
-          icbt=1
-          zt=zt+delz
-        endif
-
-        if (i.ne.ifine) then
-          zeta=zt/h
-          call hanna_short(zt)
-        endif
-
-      end do
-      if (cblflag.ne.1) nrand=nrand+i
-
-  ! Determine time step for next integration
-  !*****************************************
-
-      if (turbswitch) then
-        ldt=int(min(tlw,h/max(2.*abs(wp*sigw),1.e-5), &
-             0.5/abs(dsigwdz))*ctl)
-      else
-        ldt=int(min(tlw,h/max(2.*abs(wp),1.e-5))*ctl)
-      endif
-      ldt=max(ldt,mintime)
-
   ! Determine probability of deposition
   !************************************
@@ -487 +168 @@
   !   where f(n) is the exponential term
                prob(ks)=1.+(prob(ks)-1.)* &
-                    exp(-vdepo(ks)*abs(dt)/(2.*href))
+                    exp(-vdepo(ks)*abs(ideltas)/(2.*href))
           endif
         end do
       endif
 
-  endif
 
 end subroutine advance_rec

src/interpol_vdep.f90

@@ -54 +54 @@
 
   ! a) Bilinear horizontal interpolation
-
+! write(*,*) 'interpol: ',dt1,dt2,dtt,lsynctime,ix,jy
   do m=1,2
     indexh=memind(m)

src/par_mod.f90

@@ -186 +186 @@
   !**************************************************
 
-  integer,parameter :: maxpart=20000000
+  integer,parameter :: maxpart=40000000
   integer,parameter :: maxspec=2
   real,parameter :: minmass=0.0001

src/releaseparticles.f90

@@ -209 +209 @@
 
             ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux
-
   ! Interpolation of topography and density
   !****************************************

src/timemanager.f90

@@ -107 +107 @@
   integer :: ix,jy,ldeltat,itage,nage
   integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme 
-  real :: outnum,weight,prob(maxspec),decfact
+  real :: outnum,weight,prob_rec(maxspec),prob(maxspec),decfact
   ! real :: uap(maxpart),ucp(maxpart),uzp(maxpart)
   ! real :: us(maxpart),vs(maxpart),ws(maxpart)
@@ -552 +552 @@
       do ks=1,nspec
          if  ((xscav_frac1(j,ks).lt.0)) then
-         call advance_rec(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
-            us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
-            cbt(j))
+         call advance_rec(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec)
             if (decay(ks).gt.0.) then             ! radioactive decay
                 decfact=exp(-real(abs(lsynctime))*decay(ks))
@@ -561 +559 @@
             endif
             if (DRYDEPSPEC(ks)) then        ! dry deposition
-               drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
+               drydeposit(ks)=xmass1(j,ks)*prob_rec(ks)*decfact
                xscav_frac1(j,ks)=xscav_frac1(j,ks)*(-1.)* &
                drydeposit(ks)/xmass1(j,ks)
+!               write (*,*) 'notance: ',prob(ks),xmass1(j,ks),ztra1(j)
                if (decay(ks).gt.0.) then   ! correct for decay (see wetdepo)
                   drydeposit(ks)=drydeposit(ks)* &
@@ -572 +571 @@
                 xscav_frac1(j,ks)=0.
              endif
+!         write (*,*) 'xscav: ',j,ks,xscav_frac1(j,ks)
          endif
        enddo
@@ -600 +600 @@
             us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
             cbt(j))
+!        write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
 
   ! Calculate the gross fluxes across layer interfaces

src/wetdepo.f90

@@ -20 +20 @@
 !**********************************************************************
 
-subroutine wetdepo(itime,ltsample,loutnext,forreceptor)
-!                  i      i        i            i
+subroutine wetdepo(itime,ltsample,loutnext)
+!                  i      i        i
 !*****************************************************************************
 !                                                                            *
@@ -30 +30 @@
 ! This fraction is parameterized from total cloud cover and rates of large   *
 ! scale and convective precipitation.                                        *
-! SEC: if forrecptor is true then the wetdeposition fraction is only applied *
-! on the xscav_frac and not on the xmass                                     *
 !                                                                            *
 !    Author: A. Stohl                                                        *
@@ -94 +92 @@
   integer :: blc_count, inc_count
   real    :: Si_dummy, wetscav_dummy
-  logical :: readclouds_this_nest,forreceptor
+  logical :: readclouds_this_nest
 
 
@@ -174 +172 @@
            memtime(1),memtime(2),interp_time,lsp,convp,cc)
     endif
+
+!  If total precipitation is less than 0.01 mm/h - no scavenging occurs
+!  sec this is just valid if release is over a point
+    if ((lsp.lt.0.01).and.(convp.lt.0.01)) then
+          if (WETBKDEP) then
+             do ks=1,nspec
+                if (xscav_frac1(jpart,ks).lt.0) then ! first timestep no scavenging
+                   xmass1(jpart,ks)=0.
+                   xscav_frac1(jpart,ks)=0.
+!                  write (*,*) 'paricle removed - no scavenging: ',jpart,ks
+                endif
+             end do
+          endif
+          goto 20
+    endif
+
 
 ! get the level were the actual particle is in
@@ -399 +413 @@
         kp=1
       endif
-      if (forreceptor .eqv. .false.) then
-         if (restmass .gt. smallnum) then
-           xmass1(jpart,ks)=restmass
+      if (restmass .gt. smallnum) then
+        xmass1(jpart,ks)=restmass
 !   depostatistic
 !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
 !   depostatistic
-         else
-           xmass1(jpart,ks)=0.
-        endif
-     else ! for the backward deposition calculation
-         if (wetdeposit(ks).gt.0) then ! deposition occured
-                xscav_frac1(jpart,ks)=xscav_frac1(jpart,ks)*(-1.)* &
-                   wetdeposit(ks)/xmass1(jpart,ks)
-!                write (*,*) 'paricle kept: ',jpart,ks,wetdeposit(ks),xscav_frac1(jpart,ks)
-         else
-                xmass1(jpart,ks)=0.
-                xscav_frac1(jpart,ks)=0.
-         endif
-     endif
+      else
+        xmass1(jpart,ks)=0.
+      endif
 !   Correct deposited mass to the last time step when radioactive decay of
 !   gridded deposited mass was calculated
@@ -423 +426 @@
         wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
       endif
+
+      if (WETBKDEP) then
+! the calculation of the scavenged mass shall only be done once after release
+! xscav_frac1 was initialised with a negative value
+          if (xscav_frac1(jpart,ks).lt.0) then
+             if (wetdeposit(ks).eq.0) then
+! terminate particle
+                xmass1(jpart,ks)=0.
+                xscav_frac1(jpart,ks)=0.
+             else
+                xscav_frac1(jpart,ks)=xscav_frac1(jpart,ks)*(-1.)* &
+                   wetdeposit(ks)/xmass1(jpart,ks)
+!                write (*,*) 'paricle kept: ',jpart,ks,wetdeposit(ks),xscav_frac1(jpart,ks)
+             endif
+          endif
+      endif
+
 
     end do !all species