Skip to content
GitLab
Commit 54cbd6cb authored by Sabine's avatar Sabine
Browse files

all f90 files for dry/wet backward mode

parent dced13c8
Hide whitespace changes
Inline Side-by-side
src/com_mod.f90
View file @ 54cbd6cb
......@@ -575,7 +575,8 @@ module com_mod
integer :: numxgridn,numygridn
real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
!real outheight(maxzgrid),outheighthalf(maxzgrid)
logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC,SCAVDEP
logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC
logical :: DRYBKDEP,WETBKDEP
! numxgrid,numygrid number of grid points in x,y-direction
! numxgridn,numygridn number of grid points in x,y-direction for nested output grid
......@@ -594,7 +595,7 @@ module com_mod
! WETDEP .true., if wet deposition is switched on
! OHREA .true., if OH reaction is switched on
! ASSSPEC .true., if there are two species asscoiated
! SCAVDEP .true., for bkwd runs, where mass deposited and source regions is calculated
! DRYBKDEP,WETBKDEP .true., for bkwd runs, where mass deposited and source regions is calculated - either for dry or for wet deposition
! (i.e. transfer of mass between these two occurs
......
src/conccalc.f90
View file @ 54cbd6cb
......@@ -191,7 +191,7 @@ subroutine conccalc(itime,weight)
(yl.gt.real(numygrid-1)-0.5)).or.(usekernel.eq.0)) then ! no kernel, direct attribution to grid cell
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
......@@ -233,7 +233,7 @@ subroutine conccalc(itime,weight)
if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
if ((jy.ge.0).and.(jy.le.numygrid-1)) then
w=wx*wy
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
......@@ -250,7 +250,7 @@ subroutine conccalc(itime,weight)
if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
w=wx*(1.-wy)
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
......@@ -270,7 +270,7 @@ subroutine conccalc(itime,weight)
if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
w=(1.-wx)*(1.-wy)
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
......@@ -287,7 +287,7 @@ subroutine conccalc(itime,weight)
if ((jy.ge.0).and.(jy.le.numygrid-1)) then
w=(1.-wx)*wy
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
......@@ -328,7 +328,7 @@ subroutine conccalc(itime,weight)
(yl.gt.real(numygridn-1)-0.5)) then ! no kernel, direct attribution to grid cell
if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
(jy.le.numygridn-1)) then
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
......@@ -370,7 +370,7 @@ subroutine conccalc(itime,weight)
if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
if ((jy.ge.0).and.(jy.le.numygridn-1)) then
w=wx*wy
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
......@@ -387,7 +387,7 @@ subroutine conccalc(itime,weight)
if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
w=wx*(1.-wy)
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
......@@ -407,7 +407,7 @@ subroutine conccalc(itime,weight)
if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
w=(1.-wx)*(1.-wy)
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
......@@ -424,7 +424,7 @@ subroutine conccalc(itime,weight)
if ((jy.ge.0).and.(jy.le.numygridn-1)) then
w=(1.-wx)*wy
if (SCAVDEP) then
if (DRYBKDEP.or.WETBKDEP) then
do ks=1,nspec
griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
......
src/concoutput.f90
View file @ 54cbd6cb
......@@ -245,29 +245,32 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
do ks=1,nspec
write(anspec,'(i3.3)') ks
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if (ldirect.eq.1) then
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
atime//'_'//anspec,form='unformatted')
else
open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
atime//'_'//anspec,form='unformatted')
endif
write(unitoutgrid) itime
endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
if (DRYBKDEP.or.WETBKDEP) then !scavdep output
if (DRYBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgridppt) itime
endif
if (iout.eq.6) then !scavdep output
open(unitoutgrid,file=path(2)(1:length(2))//'grid_scav_'//adate// &
if (WETBKDEP) &
open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgrid) itime
endif
else
if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
if (ldirect.eq.1) then
open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
atime//'_'//anspec,form='unformatted')
else
open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
atime//'_'//anspec,form='unformatted')
endif
write(unitoutgrid) itime
endif
if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
atime//'_'//anspec,form='unformatted')
write(unitoutgridppt) itime
endif
endif ! if deposition output
do kp=1,maxpointspec_act
do nage=1,nageclass
......
src/readcommand.f90
View file @ 54cbd6cb
......@@ -295,6 +295,8 @@ subroutine readcommand
!Af IND_RECEPTOR switches between different units for concentrations at the receptor
!Af 1 = mass units
!Af 2 = mass mixing ratio units
!Af 3 = wet deposition in outputfield
!Af 4 = dry deposition in outputfield
if ( ldirect .eq. 1 ) then ! FWD-Run
!Af set release-switch
......@@ -317,11 +319,24 @@ subroutine readcommand
ind_samp = 0
endif
!Af set release-switch
if (ind_receptor .eq. 1) then !mass
WETBKDEP=.false.
DRYBKDEP=.false.
select case (ind_receptor)
case (1) ! 1 .. concentration at receptor
ind_rel = 1
else ! mass mix
case (2) ! 2 .. mixing ratio at receptor
ind_rel = 0
endif
case (3) ! 3 .. wet deposition in outputfield
ind_rel = 3
write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE #### '
WETBKDEP=.true.
allocate(xscav_frac1(maxpart,maxspec))
case (4) ! 4 .. dry deposition in outputfield
ind_rel = 4
write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE #### '
DRYBKDEP=.true.
allocate(xscav_frac1(maxpart,maxspec))
end select
endif
!*************************************************************
......@@ -338,21 +353,6 @@ subroutine readcommand
stop
endif
if ((ldirect.eq.-1).and.(iout.eq.6)) then
if ((ind_receptor .eq. 1) .and. (ind_source .eq. 1)) then
write(*,*) ' #### FLEXPART SCAVENGING DEPOSIT BACKWARD MODE #### '
SCAVDEP=.true.
allocate(xscav_frac1(maxpart,maxspec))
else
write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND: ####'
write(*,*) '#### FOR SCAVDEP MODE ind_source and ####'
write(*,*) '#### ind_receptor have to be 1 ! ####'
stop
endif
else
SCAVDEP=.false.
endif
! Check input dates
!******************
......@@ -400,7 +400,7 @@ subroutine readcommand
if ((iout.lt.1).or.(iout.gt.6)) then
write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, 5 OR 6 FOR #### '
write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR #### '
write(*,*) ' #### STANDARD FLEXPART OUTPUT OR 9 - 13 #### '
write(*,*) ' #### FOR NETCDF OUTPUT #### '
stop
......
src/releaseparticles.f90
View file @ 54cbd6cb
......@@ -56,7 +56,6 @@ subroutine releaseparticles(itime)
!real xaux,yaux,zaux,ran1,rfraction,xmasssave(maxpoint)
real :: xaux,yaux,zaux,rfraction
real :: topo,rhoaux(2),r,t,rhoout,ddx,ddy,rddx,rddy,p1,p2,p3,p4
real :: rhosum(nspec)
real :: dz1,dz2,dz,xtn,ytn,xlonav,timecorrect(maxspec),press,pressold
real :: presspart,average_timecorrect
integer :: itime,numrel,i,j,k,n,ix,jy,ixp,jyp,ipart,minpart,ii
......@@ -86,9 +85,6 @@ subroutine releaseparticles(itime)
minpart=1
do i=1,numpoint
do k=1,nspec
rhosum(k)=0
end do
if ((itime.ge.ireleasestart(i)).and. &! are we within release interval?
(itime.le.ireleaseend(i))) then
......@@ -185,17 +181,11 @@ subroutine releaseparticles(itime)
do k=1,nspec
xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
*timecorrect(k)/average_timecorrect
if (SCAVDEP) then ! if there is no scavenging in wetdepo it will be set to 0
if (DRYBKDEP.or.WETBKDEP) then ! if there is no scavenging in wetdepo it will be set to 0
! if ( henry(k).gt.0 .or. &
! crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. &
! ccn_aero(k).gt.0. .or. in_aero(k).gt.0. ) then
xscav_frac1(ipart,k)=-1.
! write(*,*) '190: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i
! xscav_frac1(ipart,k)=(-1.)/real(npart(i)) &
! *timecorrect(k)/average_timecorrect
! else
! xscav_frac1(ipart,k)=0
! endif
endif
! Assign certain properties to particle
!**************************************
......@@ -388,16 +378,9 @@ subroutine releaseparticles(itime)
do k=1,nspec
xmass1(ipart,k)=xmass1(ipart,k)*rhoout
if (SCAVDEP) then
xscav_frac1(ipart,k)=xscav_frac1(ipart,k)
!mctest xscav_frac1(ipart,k)=xscav_frac1(ipart,k)*rhoout
rhosum(k)=rhosum(k)+rhoout
! write(*,*) '391: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i
endif
end do
endif
numpart=max(numpart,ipart)
goto 34 ! Storage space has been found, stop searching
endif
......@@ -406,16 +389,6 @@ subroutine releaseparticles(itime)
34 minpart=ipart+1
end do ! ipart=minpart,maxpart
if (SCAVDEP) then
do ipart=minpart,maxpart
do k=1,nspec
if (xscav_frac1(ipart,k).lt.0) then
!mctest xscav_frac1(ipart,k)=xscav_frac1(ipart,k)/rhosum(k)
! write(*,*) '409: ',xscav_frac1(ipart,k),k,ipart,rhosum(k),rhoout,i
endif
end do
end do
endif
endif ! j=1,numrel
end do
......
src/timemanager.f90
View file @ 54cbd6cb
......@@ -114,6 +114,7 @@ subroutine timemanager
real(dep_prec) :: drydeposit(maxspec),wetgridtotalunc,drygridtotalunc
real :: xold,yold,zold,xmassfract
real, parameter :: e_inv = 1.0/exp(1.0)
logical :: firstdepocalc
!double precision xm(maxspec,maxpointspec_act),
! + xm_depw(maxspec,maxpointspec_act),
! + xm_depd(maxspec,maxpointspec_act)
......@@ -170,7 +171,7 @@ subroutine timemanager
if (verbosity.gt.0) then
write (*,*) 'timemanager> call wetdepo'
endif
call wetdepo(itime,lsynctime,loutnext)
call wetdepo(itime,lsynctime,loutnext,.false.)
endif
if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
......@@ -540,12 +541,60 @@ subroutine timemanager
yold=ytra1(j)
zold=ztra1(j)
! RECEPTOR: dry/wet depovel
!****************************
! Before the particle is moved
! the calculation of the scavenged mass shall only be done once after release
! xscav_frac1 was initialised with a negative value
do ks=1,nspec
if (DRYBKDEP.and.(xscav_frac1(j,ks).lt.0)) then
if (ks.eq.1) 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))
endif
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
xscav_frac1(j,ks)=xscav_frac1(j,ks)*(-1.)* &
drydeposit(ks)/xmass1(j,ks)
if (decay(ks).gt.0.) then ! correct for decay (see wetdepo)
drydeposit(ks)=drydeposit(ks)* &
exp(real(abs(ldeltat))*decay(ks))
endif
else
xmass1(j,ks)=0
xscav_frac1(j,ks)=0.
endif
endif
enddo
firstdepocalc=.false.
do ks=1,nspec
if ((WETBKDEP).and.(xscav_frac1(j,ks).lt.0) &
.and.firstdepocalc.eqv..false.) then
! Backward wetdeposition and first timestep after release
call wetdepo(itime,lsynctime,loutnext,.true.)
firstdepocalc=.true.
endif
enddo
! Integrate Lagevin equation for lsynctime seconds
!*************************************************
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))
if (verbosity.gt.0) then
if (j.eq.1) then
write (*,*) 'timemanager> call advance'
endif
endif
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))
! Calculate the gross fluxes across layer interfaces
!***************************************************
......
src/wetdepo.f90
View file @ 54cbd6cb
......@@ -19,8 +19,8 @@
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine wetdepo(itime,ltsample,loutnext)
! i i i
subroutine wetdepo(itime,ltsample,loutnext,forreceptor)
! i i i i
!*****************************************************************************
! *
! Calculation of wet deposition using the concept of scavenging coefficients.*
......@@ -29,6 +29,8 @@ subroutine wetdepo(itime,ltsample,loutnext)
! 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 *
! *
......@@ -91,7 +93,7 @@ subroutine wetdepo(itime,ltsample,loutnext)
integer :: blc_count, inc_count
real :: Si_dummy, wetscav_dummy
logical :: readclouds_this_nest
logical :: readclouds_this_nest,forreceptor
! Compute interval since radioactive decay of deposited mass was computed
......@@ -172,22 +174,6 @@ subroutine wetdepo(itime,ltsample,loutnext)
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 (SCAVDEP) 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
......@@ -412,37 +398,31 @@ subroutine wetdepo(itime,ltsample,loutnext)
else
kp=1
endif
if (restmass .gt. smallnum) then
xmass1(jpart,ks)=restmass
if (forreceptor .eqv. .false.) then
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
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
! 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 (SCAVDEP) 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