Skip to content

GitLab

Commit 58448660 authored by Sabine's avatar Sabine
Browse files

first working version with backward deposition

parent 0581cac7
Hide whitespace changes
Inline Side-by-side
src/advance.f90
View file @ 58448660
......@@ -568,6 +568,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
! where f(n) is the exponential term
prob(ks)=1.+(prob(ks)-1.)* &
exp(-vdepo(ks)*abs(dt)/(2.*href))
! write(*,*) 'Prob calc: ',zt,href,prob(ks)
endif
end do
endif
......
src/advance_rec.f90
View file @ 58448660
......@@ -19,10 +19,8 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
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
!*****************************************************************************
! *
! Calculation of turbulent particle trajectories utilizing a *
......@@ -65,32 +63,14 @@ subroutine advance_rec(itime,nrelpoint,ldt,up,vp,wp, &
!*****************************************************************************
! *
! 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 *
! *
......@@ -100,42 +80,15 @@ subroutine advance_rec(itime,nrelpoint,ldt,up,vp,wp, &
use par_mod
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
do ks=1,nspec
......@@ -144,15 +97,6 @@ subroutine advance_rec(itime,nrelpoint,ldt,up,vp,wp, &
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
! Determine whether lat/long grid or polarstereographic projection
! is to be used
......@@ -206,269 +150,6 @@ subroutine advance_rec(itime,nrelpoint,ldt,up,vp,wp, &
jyp=jy+1
! 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
!************************************
......@@ -486,12 +167,11 @@ subroutine advance_rec(itime,nrelpoint,ldt,up,vp,wp, &
! this formulation means that prob(n) = 1 - f(0)*...*f(n)
! 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
View file @ 58448660
......@@ -53,7 +53,7 @@ subroutine interpol_vdep(level,vdepo)
real :: y(2),vdepo
! a) Bilinear horizontal interpolation
! write(*,*) 'interpol: ',dt1,dt2,dtt,lsynctime,ix,jy
do m=1,2
indexh=memind(m)
......
src/par_mod.f90
View file @ 58448660
......@@ -185,7 +185,7 @@ module par_mod
! Maximum number of particles, species, and similar
!**************************************************
integer,parameter :: maxpart=20000000
integer,parameter :: maxpart=40000000
integer,parameter :: maxspec=2
real,parameter :: minmass=0.0001
......
src/releaseparticles.f90
View file @ 58448660
......@@ -208,7 +208,6 @@ subroutine releaseparticles(itime)
!*************************************
ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux
! Interpolation of topography and density
!****************************************
......
src/timemanager.f90
View file @ 58448660
......@@ -106,7 +106,7 @@ subroutine timemanager
integer :: loutnext,loutstart,loutend
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)
! integer(kind=2) :: cbt(maxpart)
......@@ -551,18 +551,17 @@ subroutine timemanager
if (DRYBKDEP) then
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))
else
decfact=1.
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)* &
exp(real(abs(ldeltat))*decay(ks))
......@@ -571,6 +570,7 @@ subroutine timemanager
xmass1(j,ks)=0
xscav_frac1(j,ks)=0.
endif
! write (*,*) 'xscav: ',j,ks,xscav_frac1(j,ks)
endif
enddo
endif
......@@ -599,6 +599,7 @@ subroutine timemanager
call advance(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))
! write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
! Calculate the gross fluxes across layer interfaces
!***************************************************
......
src/wetdepo.f90
View file @ 58448660
......@@ -19,8 +19,8 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine wetdepo(itime,ltsample,loutnext,forreceptor)
! i i i i
subroutine wetdepo(itime,ltsample,loutnext)
! i i i
!*****************************************************************************
! *
! Calculation of wet deposition using the concept of scavenging coefficients.*
......@@ -29,8 +29,6 @@ subroutine wetdepo(itime,ltsample,loutnext,forreceptor)
! grid cell, but that only a fraction of the grid cell experiences rainfall. *
! 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 *
! *
......@@ -93,7 +91,7 @@ subroutine wetdepo(itime,ltsample,loutnext,forreceptor)
integer :: blc_count, inc_count
real :: Si_dummy, wetscav_dummy
logical :: readclouds_this_nest,forreceptor
logical :: readclouds_this_nest
! Compute interval since radioactive decay of deposited mass was computed
......@@ -174,6 +172,22 @@ subroutine wetdepo(itime,ltsample,loutnext,forreceptor)
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
do il=2,nz
if (height(il).gt.ztra1(jpart)) then
......@@ -398,31 +412,37 @@ subroutine wetdepo(itime,ltsample,loutnext,forreceptor)
else
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
if (decay(ks).gt.0.) then
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
! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment