diff --git a/options/SPECIES/spec_overview b/options/SPECIES/spec_overview
old mode 100644
new mode 100755
diff --git a/src/advance.f90 b/src/advance.f90
index ed57db1be7be8870f18ae389c9bc9d7438dcdd28..6a903cff07f5fd7c68143b47849ceb10e02d22ed 100644
--- a/src/advance.f90
+++ b/src/advance.f90
@@ -116,6 +116,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
real :: dcwsave
real :: usigold,vsigold,wsigold,r,rs
real :: uold,vold,wold,vdepo(maxspec)
+ real :: h1(2)
!real uprof(nzmax),vprof(nzmax),wprof(nzmax)
!real usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
!real rhoprof(nzmax),rhogradprof(nzmax)
@@ -222,18 +223,44 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
jyp=jy+1
+ ! Determine the lower left corner and its distance to the current position
+ !*************************************************************************
+
+ ddx=xt-real(ix)
+ ddy=yt-real(jy)
+ rddx=1.-ddx
+ rddy=1.-ddy
+ p1=rddx*rddy
+ p2=ddx*rddy
+ p3=rddx*ddy
+ p4=ddx*ddy
+
+ ! Calculate variables for time interpolation
+ !*******************************************
+
+ dt1=real(itime-memtime(1))
+ dt2=real(memtime(2)-itime)
+ dtt=1./(dt1+dt2)
+
! 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
+ mind=memind(k) ! eso: compatibility with 3-field version
+ if (interpolhmix) then
+ h1(k)=p1*hmix(ix ,jy ,1,mind) &
+ + p2*hmix(ixp,jy ,1,mind) &
+ + p3*hmix(ix ,jyp,1,mind) &
+ + p4*hmix(ixp,jyp,1,mind)
+ else
+ 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
+ endif
end do
tropop=tropopause(nix,njy,1,1)
else
@@ -248,6 +275,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
tropop=tropopausen(nix,njy,1,1,ngrid)
endif
+ if (interpolhmix) h=(h1(1)*dt2+h1(2)*dt1)*dtt
zeta=zt/h
@@ -445,6 +473,14 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
delz=wp*dtf
endif
+ if (turboff) then
+!sec switch off turbulence
+ up=0.0
+ vp=0.0
+ wp=0.0
+ delz=0.
+ endif
+
!****************************************************
! Compute turbulent vertical displacement of particle
!****************************************************
@@ -646,6 +682,12 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
nrand=nrand+1
endif
+ if (turboff) then
+!sec switch off turbulence
+ ux=0.0
+ vy=0.0
+ wp=0.0
+ endif
! If particle represents only a single species, add gravitational settling
! velocity. The settling velocity is zero for gases
diff --git a/src/com_mod.f90 b/src/com_mod.f90
index 019db1a3d6b67d886a4ba0a01959d8d60be4c5a7..9ac45f1e75cf5aa0fb57bc54b52513d3389b118e 100644
--- a/src/com_mod.f90
+++ b/src/com_mod.f90
@@ -576,8 +576,10 @@ module com_mod
integer :: numxgridn,numygridn
real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
!real outheight(maxzgrid),outheighthalf(maxzgrid)
+
logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),&
& 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
@@ -597,6 +599,7 @@ module com_mod
! WETDEPSPEC .true., if wet deposition is switched on for that species
! OHREA .true., if OH reaction is switched on
! ASSSPEC .true., if there are two species asscoiated
+ ! 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
@@ -667,6 +670,7 @@ module com_mod
real(kind=dp), allocatable, dimension(:) :: xtra1, ytra1
real, allocatable, dimension(:) :: ztra1
real, allocatable, dimension(:,:) :: xmass1
+ real, allocatable, dimension(:,:) :: xscav_frac1
! eso: Moved from timemanager
real, allocatable, dimension(:) :: uap,ucp,uzp,us,vs,ws
@@ -687,7 +691,8 @@ module com_mod
! numparticlecount counts the total number of particles that have been released
! xtra1,ytra1,ztra1 spatial positions of the particles
! xmass1 [kg] particle masses
-
+ ! xscav_frac1 fraction of particle masse which has been scavenged at receptor
+
!*******************************************************
@@ -750,6 +755,11 @@ module com_mod
integer :: mpi_mode=0 ! .gt. 0 if running MPI version
logical :: lroot=.true. ! true if serial version, or if MPI .and. root process
+ logical :: usekernel=.false. ! true if the output kernel shall be switched on
+ logical :: interpolhmix=.false. ! true if the hmix shall be interpolated
+ logical :: turboff=.false. ! true if the turbulence shall be switched off
+
+
contains
subroutine com_mod_allocate_part(nmpart)
!*******************************************************************************
diff --git a/src/conccalc.f90 b/src/conccalc.f90
index 36ad826f618c31cdb20ee0ec35b978c4eae2262a..9d15983006d85ed114d12463cb80e1fa5b0b279e 100644
--- a/src/conccalc.f90
+++ b/src/conccalc.f90
@@ -20,7 +20,7 @@
!**********************************************************************
subroutine conccalc(itime,weight)
- ! i i
+ ! i i
!*****************************************************************************
! *
! Calculation of the concentrations on a regular grid using volume *
@@ -58,13 +58,13 @@ subroutine conccalc(itime,weight)
real :: rhoprof(2),rhoi
real :: xl,yl,wx,wy,w
real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
-
+! integer xscav_count
! For forward simulations, make a loop over the number of species;
! for backward simulations, make an additional loop over the
! releasepoints
!***************************************************************************
-
+! xscav_count=0
do i=1,numpart
if (itra1(i).ne.itime) goto 20
@@ -75,7 +75,8 @@ subroutine conccalc(itime,weight)
end do
33 continue
-
+! if (xscav_frac1(i,1).lt.0) xscav_count=xscav_count+1
+
! For special runs, interpolate the air density to the particle position
!************************************************************************
!***********************************************************************
@@ -171,7 +172,6 @@ subroutine conccalc(itime,weight)
jy=int(yl)
if (yl.lt.0.) jy=jy-1
- ! if (i.eq.10000) write(*,*) itime,xtra1(i),ytra1(i),ztra1(i),xl,yl
! For particles aged less than 3 hours, attribute particle mass to grid cell
@@ -182,17 +182,25 @@ subroutine conccalc(itime,weight)
if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
(xl.gt.real(numxgrid-1)-0.5).or. &
- (yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell
+ (yl.gt.real(numygrid-1)-0.5))) 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
- do ks=1,nspec
- gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
- end do
+ end do
+ endif
endif
- else ! attribution via uniform kernel
+ else ! attribution via uniform kernel
ddx=xl-real(ix) ! distance to left cell border
ddy=yl-real(jy) ! distance to lower cell border
@@ -219,46 +227,76 @@ 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
- do ks=1,nspec
- gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
- end do
+ end do
+ endif
endif
if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
w=wx*(1.-wy)
- do ks=1,nspec
- gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
- end do
+ end do
+ endif
endif
- endif
+ endif !ix ge 0
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)
- do ks=1,nspec
- gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
- end do
+ end do
+ endif
endif
if ((jy.ge.0).and.(jy.le.numygrid-1)) then
w=(1.-wx)*wy
- do ks=1,nspec
- gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
- end do
+ end do
+ endif
endif
- endif
- endif
-
-
+ endif !ixp ge 0
+ endif
!************************************
! Do everything for the nested domain
@@ -281,14 +319,22 @@ subroutine conccalc(itime,weight)
if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
(xl.gt.real(numxgridn-1)-0.5).or. &
- (yl.gt.real(numygridn-1)-0.5)) then ! no kernel, direct attribution to grid cell
+ (yl.gt.real(numygridn-1)-0.5).or.(.not.usekernel)) 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
- do ks=1,nspec
- griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight
- end do
+ end do
+ endif
endif
else ! attribution via uniform kernel
@@ -318,20 +364,36 @@ 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
- do ks=1,nspec
- griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
- end do
+ end do
+ endif
endif
if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
w=wx*(1.-wy)
- do ks=1,nspec
- griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
- end do
+ end do
+ endif
endif
endif
@@ -339,28 +401,44 @@ 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)
- do ks=1,nspec
- griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
- end do
+ end do
+ endif
endif
if ((jy.ge.0).and.(jy.le.numygridn-1)) then
w=(1.-wx)*wy
- do ks=1,nspec
- griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+ 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)+ &
+ xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+ end do
+ else
+ do ks=1,nspec
+ griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
xmass1(i,ks)/rhoi*weight*w
- end do
+ end do
+ endif
endif
endif
endif
-
endif
endif
20 continue
end do
+! write(*,*) 'xscav count:',xscav_count
!***********************************************************************
! 2. Evaluate concentrations at receptor points, using the kernel method
diff --git a/src/concoutput.f90 b/src/concoutput.f90
index 4a8bb0e553b8d5378ef81b8b1efd4c3b29aec19d..a237705a33eecff7b7241ef1809ec03a1375426a 100644
--- a/src/concoutput.f90
+++ b/src/concoutput.f90
@@ -245,23 +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// &
+
+ if (DRYBKDEP.or.WETBKDEP) then !scavdep output
+ if (DRYBKDEP) &
+ open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
+ atime//'_'//anspec,form='unformatted')
+ if (WETBKDEP) &
+ open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
+ atime//'_'//anspec,form='unformatted')
+ write(unitoutgrid) itime
+ 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// &
+ else
+ open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
atime//'_'//anspec,form='unformatted')
+ endif
+ write(unitoutgrid) itime
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 ((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
+ write(unitoutgridppt) itime
+ endif
+ endif ! if deposition output
do kp=1,maxpointspec_act
do nage=1,nageclass
@@ -341,7 +350,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Concentration output
!*********************
- if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+ if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
! Wet deposition
sp_count_i=0
diff --git a/src/drydepokernel.f90 b/src/drydepokernel.f90
index ce7c3c8a3e870df194b3039f604941b585e3635f..6ba8d2dec2955a308b0c657dd513b145af02017f 100644
--- a/src/drydepokernel.f90
+++ b/src/drydepokernel.f90
@@ -101,36 +101,40 @@ subroutine drydepokernel(nunc,deposit,x,y,nage,kp)
if ((abs(deposit(ks)).gt.0).and.DRYDEPSPEC(ks)) then
- if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
+ if (.not.usekernel) then
+ drygridunc(ix,jy,ks,kp,nunc,nage)= &
+ drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)
+ else
+ if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
- w=wx*wy
- drygridunc(ix,jy,ks,kp,nunc,nage)= &
- drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
- continue
- endif
+ w=wx*wy
+ drygridunc(ix,jy,ks,kp,nunc,nage)= &
+ drygridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
+ endif
- if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
+ if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
(jyp.le.numygrid-1)) then
w=(1.-wx)*(1.-wy)
drygridunc(ixp,jyp,ks,kp,nunc,nage)= &
drygridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
- endif
+ endif
- if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
+ if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
- w=(1.-wx)*wy
+ w=(1.-wx)*wy
drygridunc(ixp,jy,ks,kp,nunc,nage)= &
drygridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
- endif
+ endif
- if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
+ if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
(jyp.le.numygrid-1)) then
- w=wx*(1.-wy)
+ w=wx*(1.-wy)
drygridunc(ix,jyp,ks,kp,nunc,nage)= &
drygridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
- endif
+ endif
- endif
+ endif ! kernel
+ endif ! deposit>0
end do
end if
diff --git a/src/ecmwf_mod.f90 b/src/ecmwf_mod.f90
index 063cb9542e665d8cf5b9e37ba1a8636847a1b7bd..a22cc0c1fde04e84ae892ff150c26a6e94617463 100644
--- a/src/ecmwf_mod.f90
+++ b/src/ecmwf_mod.f90
@@ -42,8 +42,10 @@ module wind_mod
!*********************************************
! integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new
- integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
- integer,parameter :: nxshift=359
+ integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new
+! integer,parameter :: nxmax=721,nymax=361,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 0.5
+ integer,parameter :: nxshift=359
+! integer,parameter :: nxshift=718
! integer,parameter :: nxshift=0
!*********************************************
diff --git a/src/get_vdep_prob.f90 b/src/get_vdep_prob.f90
new file mode 100644
index 0000000000000000000000000000000000000000..4b0ac14f9e9f1a888f3fd778aa241a92f9755eee
--- /dev/null
+++ b/src/get_vdep_prob.f90
@@ -0,0 +1,145 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
+! *
+! This file is part of FLEXPART. *
+! *
+! FLEXPART is free software: you can redistribute it and/or modify *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or *
+! (at your option) any later version. *
+! *
+! FLEXPART is distributed in the hope that it will be useful, *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+! GNU General Public License for more details. *
+! *
+! You should have received a copy of the GNU General Public License *
+! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
+!**********************************************************************
+
+subroutine get_vdep_prob(itime,xt,yt,zt,prob)
+ ! i i i i o
+ !*****************************************************************************
+ ! *
+ ! Calculation of the probability for dyr deposition *
+ ! *
+ ! Particle positions are read in - prob returned *
+ ! *
+ !*****************************************************************************
+ ! *
+ ! Variables: *
+ ! 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 *
+ ! 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 *
+ ! prob probability of absorption due to dry deposition *
+ ! vdepo Deposition velocities for all species *
+ ! xt,yt,zt Particle position *
+ ! *
+ !*****************************************************************************
+
+ use point_mod
+ use par_mod
+ use com_mod
+ use interpol_mod
+
+ implicit none
+
+ real(kind=dp) :: xt,yt
+ 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
+ depoindicator(ks)=.true.
+ prob(ks)=0.
+ end do
+ endif
+
+
+ ! Determine whether lat/long grid or polarstereographic projection
+ ! is to be used
+ ! Furthermore, determine which nesting level to be used
+ !*****************************************************************
+
+ if (nglobal.and.(yt.gt.switchnorthg)) then
+ ngrid=-1
+ else if (sglobal.and.(yt.lt.switchsouthg)) then
+ ngrid=-2
+ else
+ ngrid=0
+ do j=numbnests,1,-1
+ if ((xt.gt.xln(j)+eps).and.(xt.lt.xrn(j)-eps).and. &
+ (yt.gt.yln(j)+eps).and.(yt.lt.yrn(j)-eps)) then
+ ngrid=j
+ goto 23
+ endif
+ end do
+23 continue
+ endif
+
+
+ !***************************
+ ! Interpolate necessary data
+ !***************************
+
+ if (abs(itime-memtime(1)).lt.abs(itime-memtime(2))) then
+ memindnext=1
+ else
+ memindnext=2
+ endif
+
+ ! Determine nested grid coordinates
+ !**********************************
+
+ if (ngrid.gt.0) then
+ xtn=(xt-xln(ngrid))*xresoln(ngrid)
+ ytn=(yt-yln(ngrid))*yresoln(ngrid)
+ ix=int(xtn)
+ jy=int(ytn)
+ nix=nint(xtn)
+ njy=nint(ytn)
+ else
+ ix=int(xt)
+ jy=int(yt)
+ nix=nint(xt)
+ njy=nint(yt)
+ endif
+ ixp=ix+1
+ jyp=jy+1
+
+
+ ! Determine probability of deposition
+ !************************************
+
+ if ((DRYDEP).and.(zt.lt.2.*href)) then
+ do ks=1,nspec
+ if (DRYDEPSPEC(ks)) then
+ if (depoindicator(ks)) then
+ if (ngrid.le.0) then
+ call interpol_vdep(ks,vdepo(ks))
+ else
+ call interpol_vdep_nests(ks,vdepo(ks))
+ endif
+ endif
+ ! correction by Petra Seibert, 10 April 2001
+ ! this formulation means that prob(n) = 1 - f(0)*...*f(n)
+ ! where f(n) is the exponential term
+ prob(ks)=vdepo(ks)
+! prob(ks)=vdepo(ks)/2./href
+! instead of prob - return vdepo -> result kg/m2/s
+ endif
+ end do
+ endif
+
+
+end subroutine get_vdep_prob
+
diff --git a/src/get_wetscav.f90 b/src/get_wetscav.f90
new file mode 100644
index 0000000000000000000000000000000000000000..57574d6322aa2c54c2259fe55fca741b873f38ed
--- /dev/null
+++ b/src/get_wetscav.f90
@@ -0,0 +1,344 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
+! *
+! This file is part of FLEXPART. *
+! *
+! FLEXPART is free software: you can redistribute it and/or modify *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or *
+! (at your option) any later version. *
+! *
+! FLEXPART is distributed in the hope that it will be useful, *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+! GNU General Public License for more details. *
+! *
+! You should have received a copy of the GNU General Public License *
+! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
+!**********************************************************************
+
+subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc_count,wetscav)
+! i i i i i o o o o
+!*****************************************************************************
+! *
+! Calculation of wet deposition using the concept of scavenging coefficients.*
+! For lack of detailed information, washout and rainout are jointly treated. *
+! It is assumed that precipitation does not occur uniformly within the whole *
+! 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. *
+! *
+! Author: A. Stohl *
+! *
+! 1 December 1996 *
+! *
+! Correction by Petra Seibert, Sept 2002: *
+! use centred precipitation data for integration *
+! Code may not be correct for decay of deposition! *
+! *
+!*****************************************************************************
+! *
+! Variables: *
+! cc [0-1] total cloud cover *
+! convp [mm/h] convective precipitation rate *
+! grfraction [0-1] fraction of grid, for which precipitation occurs *
+! ix,jy indices of output grid cell for each particle *
+! itime [s] actual simulation time [s] *
+! jpart particle index *
+! lfr, cfr area fraction covered by precipitation for large scale *
+! and convective precipitation (dependent on prec. rate) *
+! loutnext [s] time for which gridded deposition is next output *
+! loutstep [s] interval at which gridded deposition is output *
+! lsp [mm/h] large scale precipitation rate *
+! ltsample [s] interval over which mass is deposited *
+! prec [mm/h] precipitation rate in subgrid, where precipitation occurs*
+! wetgrid accumulated deposited mass on output grid *
+! wetscav scavenging coefficient *
+! *
+! Constants: *
+! *
+!*****************************************************************************
+
+ use point_mod
+ use par_mod
+ use com_mod
+
+ implicit none
+
+ integer :: jpart,itime,ltsample,loutnext,i,j,ix,jy
+ integer :: ngrid,hz,il,interp_time, n
+ integer(kind=1) :: clouds_v
+ integer :: ks, kp
+ integer(selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count
+
+! integer :: n1,n2, icbot,ictop, indcloud !TEST
+ real :: S_i, act_temp, cl, cle ! in cloud scavenging
+ real :: clouds_h ! cloud height for the specific grid point
+ real :: xtn,ytn,lsp,convp,cc,grfraction(3),prec(3),wetscav,totprec
+ real :: restmass
+ real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
+!save lfr,cfr
+
+ real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
+ real, parameter :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
+
+!ZHG aerosol below-cloud scavenging removal polynomial constants for rain and snow
+ real, parameter :: bclr(6) = (/274.35758, 332839.59273, 226656.57259, 58005.91340, 6588.38582, 0.244984/) !rain (Laakso et al 2003)
+ real, parameter :: bcls(6) = (/22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
+ real :: frac_act, liq_frac, dquer_m
+
+ real :: Si_dummy, wetscav_dummy
+ logical :: readclouds_this_nest
+
+
+ wetscav=0.
+
+! Determine which nesting level to be used
+!*****************************************
+ ngrid=0
+ do j=numbnests,1,-1
+ if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
+ (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
+ ngrid=j
+ goto 23
+ endif
+ end do
+23 continue
+
+
+! Determine nested grid coordinates
+!**********************************
+ readclouds_this_nest=.false.
+
+ if (ngrid.gt.0) then
+ xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
+ ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
+ ix=int(xtn)
+ jy=int(ytn)
+ if (readclouds_nest(ngrid)) readclouds_this_nest=.true.
+ else
+ ix=int(xtra1(jpart))
+ jy=int(ytra1(jpart))
+ endif
+
+! Interpolate large scale precipitation, convective precipitation and
+! total cloud cover
+! Note that interpolated time refers to itime-0.5*ltsample [PS]
+!********************************************************************
+ interp_time=nint(itime-0.5*ltsample)
+
+ n=memind(2)
+ if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
+ n=memind(1)
+
+ if (ngrid.eq.0) then
+ call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
+ 1,nx,ny,n,real(xtra1(jpart)),real(ytra1(jpart)),1, &
+ memtime(1),memtime(2),interp_time,lsp,convp,cc)
+ else
+ call interpol_rain_nests(lsprecn,convprecn,tccn, &
+ nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,n,xtn,ytn,1, &
+ memtime(1),memtime(2),interp_time,lsp,convp,cc)
+ endif
+
+! If total precipitation is less than 0.01 mm/h - no scavenging occurs
+ if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+
+! get the level were the actual particle is in
+ do il=2,nz
+ if (height(il).gt.ztra1(jpart)) then
+ hz=il-1
+! goto 26
+ exit
+ endif
+ end do
+!26 continue
+
+
+ if (ngrid.eq.0) then
+ clouds_v=clouds(ix,jy,hz,n)
+ clouds_h=cloudsh(ix,jy,n)
+ else
+ clouds_v=cloudsn(ix,jy,hz,n,ngrid)
+ clouds_h=cloudshn(ix,jy,n,ngrid)
+ endif
+
+! if there is no precipitation or the particle is above the clouds no
+! scavenging is done
+
+ if (clouds_v.le.1) goto 20
+
+! 1) Parameterization of the the area fraction of the grid cell where the
+! precipitation occurs: the absolute limit is the total cloud cover, but
+! for low precipitation rates, an even smaller fraction of the grid cell
+! is used. Large scale precipitation occurs over larger areas than
+! convective precipitation.
+!**************************************************************************
+
+ if (lsp.gt.20.) then
+ i=5
+ else if (lsp.gt.8.) then
+ i=4
+ else if (lsp.gt.3.) then
+ i=3
+ else if (lsp.gt.1.) then
+ i=2
+ else
+ i=1
+ endif
+
+ if (convp.gt.20.) then
+ j=5
+ else if (convp.gt.8.) then
+ j=4
+ else if (convp.gt.3.) then
+ j=3
+ else if (convp.gt.1.) then
+ j=2
+ else
+ j=1
+ endif
+
+
+!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
+! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
+! for now they are treated the same
+ grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
+ grfraction(2)=max(0.05,cc*(lfr(i)))
+ grfraction(3)=max(0.05,cc*(cfr(j)))
+
+
+! 2) Computation of precipitation rate in sub-grid cell
+!******************************************************
+ prec(1)=(lsp+convp)/grfraction(1)
+ prec(2)=(lsp)/grfraction(2)
+ prec(3)=(convp)/grfraction(3)
+
+
+! 3) Computation of scavenging coefficients for all species
+! Computation of wet deposition
+!**********************************************************
+
+
+ if (ngrid.gt.0) then
+ act_temp=ttn(ix,jy,hz,n,ngrid)
+ else
+ act_temp=tt(ix,jy,hz,n)
+ endif
+
+!***********************
+! BELOW CLOUD SCAVENGING
+!***********************
+ if (clouds_v.ge.4) then !below cloud
+
+! For gas: if positive below-cloud parameters (A or B), and dquer<=0
+!******************************************************************
+ if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
+ ! if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then
+ blc_count(ks)=blc_count(ks)+1
+ wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
+
+! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0
+!*********************************************************************************
+ else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then
+ blc_count(ks)=blc_count(ks)+1
+
+!NIK 17.02.2015
+! For the calculation here particle size needs to be in meter and not um as dquer is
+! changed in readreleases
+! For particles larger than 10 um use the largest size defined in the parameterizations (10um)
+ dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
+
+! Rain:
+ if (act_temp .ge. 273. .and. crain_aero(ks).gt.0.) then
+
+! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
+! the below-cloud scavenging (rain efficienty) parameter Crain (=crain_aero) from SPECIES file
+ wetscav=crain_aero(ks)*10**(bclr(1)+(bclr(2)*(log10(dquer_m))**(-4))+ &
+ & (bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)*(log10(dquer_m))**(-2))+&
+ &(bclr(5)*(log10(dquer_m))**(-1))+bclr(6)* (prec(1))**(0.5))
+
+! Snow:
+ elseif (act_temp .lt. 273. .and. csnow_aero(ks).gt.0.) then
+! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
+! the below-cloud scavenging (Snow efficiency) parameter Csnow (=csnow_aero) from SPECIES file
+ wetscav=csnow_aero(ks)*10**(bcls(1)+(bcls(2)*(log10(dquer_m))**(-4))+&
+ &(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)*(log10(dquer_m))**(-2))+&
+ &(bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
+
+ endif
+
+! write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
+
+ endif ! gas or particle
+! endif ! positive below-cloud scavenging parameters given in Species file
+ endif !end BELOW
+
+!********************
+! IN CLOUD SCAVENGING
+!********************
+ if (clouds_v.lt.4) then ! In-cloud
+! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are
+! given in species file, or if gas and positive Henry's constant
+ if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
+ inc_count(ks)=inc_count(ks)+1
+! write(*,*) 'Incloud: ',inc_count
+! if negative coefficients (turned off) set to zero for use in equation
+ if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
+ if (in_aero(ks).lt.0.) in_aero(ks)=0.
+
+!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
+! nested fields
+ if (ngrid.gt.0.and.readclouds_this_nest) then
+ cl=ctwcn(ix,jy,n,ngrid)*(grfraction(1)/cc)
+ else if (ngrid.eq.0.and.readclouds) then
+ cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
+ else !parameterize cloudwater m2/m3
+!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
+ cl=1.6E-6*prec(1)**0.36
+ endif
+
+!ZHG: Calculate the partition between liquid and water phase water.
+ if (act_temp .le. 253.) then
+ liq_frac=0
+ else if (act_temp .ge. 273.) then
+ liq_frac=1
+ else
+ liq_frac =((act_temp-273.)/(273.-253.))**2.
+ end if
+! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
+ frac_act = liq_frac*ccn_aero(ks) +(1-liq_frac)*in_aero(ks)
+
+!ZHG Use the activated fraction and the liqid water to calculate the washout
+
+! AEROSOL
+!********
+ if (dquer(ks).gt.0.) then
+ S_i= frac_act/cl
+
+! GAS
+!****
+ else
+
+ cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+!REPLACE to switch old/ new scheme
+ ! S_i=frac_act/cle
+ S_i=1/cle
+ endif ! gas or particle
+
+! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
+!OLD
+ if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
+ wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
+ else
+ wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
+ endif
+ endif ! positive in-cloud scavenging parameters given in Species file
+ endif !incloud
+
+
+20 continue
+
+end subroutine get_wetscav
diff --git a/src/getfields.f90 b/src/getfields.f90
index fc4061484314c614ce844f37bdcd66372c3f3489..bfc24d25a8b029f958827da19ddc9070806189dc 100644
--- a/src/getfields.f90
+++ b/src/getfields.f90
@@ -137,6 +137,10 @@ subroutine getfields(itime,nstop)
end do
40 indmin=indj
+ if (WETBKDEP) then
+ call writeprecip(itime,memind(1))
+ endif
+
else
! No wind fields, which can be used, are currently in memory
@@ -168,6 +172,10 @@ subroutine getfields(itime,nstop)
end do
60 indmin=indj
+ if (WETBKDEP) then
+ call writeprecip(itime,memind(1))
+ endif
+
endif
lwindinterv=abs(memtime(2)-memtime(1))
diff --git a/src/interpol_rain.f90 b/src/interpol_rain.f90
index 152aeb1a637b87e11c9e9ca0cb33b6970e305ab9..925b8721df13171e0d4fbd36f9cf4375cefc7809 100644
--- a/src/interpol_rain.f90
+++ b/src/interpol_rain.f90
@@ -20,7 +20,7 @@
!**********************************************************************
subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
- ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
+ ny,iwftouse,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
! i i i i i i i
!i i i i i i i i o o o
!****************************************************************************
@@ -76,6 +76,7 @@ subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
real :: yy3(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
+ integer :: iwftouse
@@ -112,36 +113,38 @@ subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
! Loop over 2 time steps
!***********************
- do m=1,2
- indexh=memind(m)
+! do m=1,2
+ indexh=iwftouse
-
- y1(m)=p1*yy1(ix ,jy ,level,indexh) &
+ y1(1)=p1*yy1(ix ,jy ,level,indexh) &
+ p2*yy1(ixp,jy ,level,indexh) &
+ p3*yy1(ix ,jyp,level,indexh) &
+ p4*yy1(ixp,jyp,level,indexh)
- y2(m)=p1*yy2(ix ,jy ,level,indexh) &
+ y2(1)=p1*yy2(ix ,jy ,level,indexh) &
+ p2*yy2(ixp,jy ,level,indexh) &
+ p3*yy2(ix ,jyp,level,indexh) &
+ p4*yy2(ixp,jyp,level,indexh)
- y3(m)=p1*yy3(ix ,jy ,level,indexh) &
+ y3(1)=p1*yy3(ix ,jy ,level,indexh) &
+ p2*yy3(ixp,jy ,level,indexh) &
+ p3*yy3(ix ,jyp,level,indexh) &
+ p4*yy3(ixp,jyp,level,indexh)
- end do
+! end do
!************************************
- ! 2.) Temporal interpolation (linear)
+ ! 2.) Temporal interpolation (linear) - skip to be consistent with clouds
!************************************
- dt1=real(itime-itime1)
- dt2=real(itime2-itime)
- dt=dt1+dt2
+! dt1=real(itime-itime1)
+! dt2=real(itime2-itime)
+! dt=dt1+dt2
- yint1=(y1(1)*dt2+y1(2)*dt1)/dt
- yint2=(y2(1)*dt2+y2(2)*dt1)/dt
- yint3=(y3(1)*dt2+y3(2)*dt1)/dt
+! yint1=(y1(1)*dt2+y1(2)*dt1)/dt
+! yint2=(y2(1)*dt2+y2(2)*dt1)/dt
+! yint3=(y3(1)*dt2+y3(2)*dt1)/dt
+ yint1=y1(1)
+ yint2=y2(1)
+ yint3=y3(1)
end subroutine interpol_rain
diff --git a/src/interpol_vdep.f90 b/src/interpol_vdep.f90
index 8f5caa3248e609155367fabc44b2c8e3361e58ed..eb6b3e1393b68e1f6990ebffb2bd4b5f6882565b 100644
--- a/src/interpol_vdep.f90
+++ b/src/interpol_vdep.f90
@@ -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)
diff --git a/src/makefile b/src/makefile
index e84077a989eff961dbef723b7b0234dcebfc8b19..50005d293b4e30c11069fd92031e34e16f766254 100644
--- a/src/makefile
+++ b/src/makefile
@@ -67,7 +67,7 @@ O_LEV_DBG = g # [0,g]
## LIBRARIES
LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff # -fopenmp
-FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) # -march=native
+FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) #-Warray-bounds -fcheck=all # -march=native
DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
@@ -143,6 +143,7 @@ OBJECTS_GFS = \
OBJECTS = \
advance.o initialize.o \
writeheader.o writeheader_txt.o \
+writeprecip.o \
writeheader_surf.o assignland.o\
part0.o gethourlyOH.o\
caldate.o partdep.o \
@@ -151,7 +152,8 @@ raerod.o readcommand.o \
drydepokernel.o readreceptors.o \
erf.o readavailable.o \
ew.o readreleases.o \
-readdepo.o \
+readdepo.o get_vdep_prob.o \
+get_wetscav.o \
psim.o outgrid_init.o \
outgrid_init_nest.o \
photo_O1D.o readlanduse.o \
@@ -270,6 +272,8 @@ cleanall:
.SUFFIXES = $(SUFFIXES) .f90
## DEPENDENCIES
+get_vdep_prob.o: cmapf_mod.o com_mod.o hanna_mod.o interpol_mod.o par_mod.o \
+ point_mod.o random_mod.o
advance.o: cmapf_mod.o com_mod.o hanna_mod.o interpol_mod.o par_mod.o \
point_mod.o random_mod.o
assignland.o: com_mod.o par_mod.o
@@ -423,10 +427,12 @@ unc_mod.o: par_mod.o
verttransform.o: cmapf_mod.o com_mod.o par_mod.o
verttransform_gfs.o: cmapf_mod.o com_mod.o par_mod.o
verttransform_nests.o: com_mod.o par_mod.o
+get_wetscav.o: com_mod.o par_mod.o point_mod.o
wetdepo.o: com_mod.o par_mod.o point_mod.o
wetdepokernel.o: com_mod.o par_mod.o unc_mod.o
wetdepokernel_nest.o: com_mod.o par_mod.o unc_mod.o
writeheader.o: com_mod.o outg_mod.o par_mod.o point_mod.o
+writeprecip.o: com_mod.o par_mod.o point_mod.o
writeheader_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o
writeheader_nest_surf.o: com_mod.o outg_mod.o par_mod.o point_mod.o
writeheader_surf.o: com_mod.o outg_mod.o par_mod.o point_mod.o
diff --git a/src/par_mod.f90 b/src/par_mod.f90
index e05205deaf4b87d9c9cda26d4dcc5bc89dfd5b5d..301201e8f68f7a7d5a7afbe24713c32fd8a5daba 100644
--- a/src/par_mod.f90
+++ b/src/par_mod.f90
@@ -193,6 +193,7 @@ module par_mod
integer,parameter :: maxpart=10000000
integer,parameter :: maxspec=4
+
real,parameter :: minmass=0.0001
! maxpart Maximum number of particles
@@ -257,7 +258,7 @@ module par_mod
integer,parameter :: unitspecies=1, unitoutrecept=91, unitoutreceptppt=92
integer,parameter :: unitlsm=1, unitsurfdata=1, unitland=1, unitwesely=1
integer,parameter :: unitOH=1
- integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95
+ integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95, unitprecip=101
integer,parameter :: unitboundcond=89
integer,parameter :: unittmp=101
diff --git a/src/readcommand.f90 b/src/readcommand.f90
index 89023976e2e8ebd64995cdad7fdbe8f50f0a1cc3..cc576ce8d5391fb97621db85e77c9c311294ab84 100644
--- a/src/readcommand.f90
+++ b/src/readcommand.f90
@@ -297,6 +297,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
+ ! 3 = wet deposition in outputfield
+ ! 4 = dry deposition in outputfield
if ( ldirect .eq. 1 ) then ! FWD-Run
!Af set release-switch
@@ -319,11 +321,28 @@ 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 #### '
+ write(*,*) ' #### Releaseheight is forced to 0 - 20km #### '
+ write(*,*) ' #### Release is performed above ground lev #### '
+ WETBKDEP=.true.
+ allocate(xscav_frac1(maxpart,maxspec))
+ case (4) ! 4 .. dry deposition in outputfield
+ ind_rel = 4
+ write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE #### '
+ write(*,*) ' #### Releaseheight is forced to 0 - 2*href #### '
+ write(*,*) ' #### Release is performed above ground lev #### '
+ DRYBKDEP=.true.
+ allocate(xscav_frac1(maxpart,maxspec))
+ end select
endif
!*************************************************************
@@ -385,9 +404,9 @@ subroutine readcommand
! Check whether a valid option for gridded model output has been chosen
!**********************************************************************
- if ((iout.lt.1).or.(iout.gt.5)) then
+ if ((iout.lt.1).or.(iout.gt.6)) then
write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND: #### '
- write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5 FOR #### '
+ write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR #### '
write(*,*) ' #### STANDARD FLEXPART OUTPUT OR 9 - 13 #### '
write(*,*) ' #### FOR NETCDF OUTPUT #### '
stop
diff --git a/src/readreleases.f90 b/src/readreleases.f90
index 3c50a20b4b6b3c43938efabb5c23ef0217a4b34c..1f46c27aad78d2560e03ae95d187a297d3910ba8 100644
--- a/src/readreleases.f90
+++ b/src/readreleases.f90
@@ -520,6 +520,20 @@ subroutine readreleases
stop
endif
+ ! If FLEXPART is run for backward deposition force zpoint
+ !*********************************************************************
+ if (WETBKDEP) then
+ zpoint1(numpoint)=0.
+ zpoint2(numpoint)=20000.
+ kindz(numpoint)=1
+ endif
+ if (DRYBKDEP) then
+ zpoint1(numpoint)=0.
+ zpoint2(numpoint)=2.*href
+ kindz(numpoint)=1
+ endif
+
+
! Check whether x coordinates of release point are within model domain
!*********************************************************************
diff --git a/src/releaseparticles.f90 b/src/releaseparticles.f90
index 87a08ee2a60624a7b6e5ec052087d2064a11c1e4..e2f7a35a53343087b6910f0dda7be24c2c91d33e 100644
--- a/src/releaseparticles.f90
+++ b/src/releaseparticles.f90
@@ -175,11 +175,18 @@ subroutine releaseparticles(itime)
! correction factor, by which the number of particles released this time has been
! scaled. Adjust the mass per particle by the species-dependent time correction factor
! divided by the species-average one
+ ! for the scavenging calculation the mass needs to be multiplied with rho of the particle layer and
+ ! divided by the sum of rho of all particles.
!*****************************************************************************
do k=1,nspec
xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
*timecorrect(k)/average_timecorrect
- ! write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
+ 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.
+ endif
! Assign certain properties to particle
!**************************************
end do
@@ -201,7 +208,6 @@ subroutine releaseparticles(itime)
!*************************************
ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux
-
! Interpolation of topography and density
!****************************************
@@ -329,7 +335,7 @@ subroutine releaseparticles(itime)
!Af ind_rel is defined in readcommand.f
- if (ind_rel .eq. 1) then
+ if ((ind_rel .eq. 1).or.(ind_rel .eq. 3).or.(ind_rel .eq. 4)) then
! Interpolate the air density
!****************************
@@ -374,16 +380,15 @@ subroutine releaseparticles(itime)
end do
endif
-
numpart=max(numpart,ipart)
goto 34 ! Storage space has been found, stop searching
endif
- end do
+ end do ! i=1:numpoint
if (ipart.gt.maxpart) goto 996
34 minpart=ipart+1
- end do
- endif
+ end do ! ipart=minpart,maxpart
+ endif ! j=1,numrel
end do
diff --git a/src/timemanager.f90 b/src/timemanager.f90
index 7b40ac91e32d4c5a25a20c6b15074db37d282f11..7ae91dbe0401a95b6f40e80b0657934093994abc 100644
--- a/src/timemanager.f90
+++ b/src/timemanager.f90
@@ -104,16 +104,18 @@ subroutine timemanager
integer :: j,ks,kp,l,n,itime=0,nstop,nstop1
! integer :: ksp
integer :: loutnext,loutstart,loutend
- integer :: ix,jy,ldeltat,itage,nage
+ integer :: ix,jy,ldeltat,itage,nage,idummy
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,wetscav(maxspec)
! real :: uap(maxpart),ucp(maxpart),uzp(maxpart)
! real :: us(maxpart),vs(maxpart),ws(maxpart)
! integer(kind=2) :: cbt(maxpart)
real(sp) :: gridtotalunc
real(dep_prec) :: drydeposit(maxspec),wetgridtotalunc,drygridtotalunc
real :: xold,yold,zold,xmassfract
+ real :: grfraction(3)
real, parameter :: e_inv = 1.0/exp(1.0)
+
!double precision xm(maxspec,maxpointspec_act),
! + xm_depw(maxspec,maxpointspec_act),
! + xm_depd(maxspec,maxpointspec_act)
@@ -144,6 +146,9 @@ subroutine timemanager
! print*, 'Initialized lifetime'
!CGZ-lifetime: set lifetime to 0
+ if (.not.usekernel) write(*,*) 'Not using the kernel'
+ if (turboff) write(*,*) 'Turbulence switched off'
+
write(*,46) float(itime)/3600,itime,numpart
if (verbosity.gt.0) then
@@ -540,12 +545,55 @@ 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
+
+ if (DRYBKDEP) then
+ do ks=1,nspec
+ if ((xscav_frac1(j,ks).lt.0)) then
+ call get_vdep_prob(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec)
+ if (DRYDEPSPEC(ks)) then ! dry deposition
+ xscav_frac1(j,ks)=prob_rec(ks)
+ else
+ xmass1(j,ks)=0.
+ xscav_frac1(j,ks)=0.
+ endif
+ endif
+ enddo
+ endif
+
+ if (WETBKDEP) then
+ do ks=1,nspec
+ if ((xscav_frac1(j,ks).lt.0)) then
+ call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
+ if (wetscav(ks).gt.0) then
+ xscav_frac1(j,ks)=wetscav(ks)* &
+ (zpoint2(npoint(j))-zpoint1(npoint(j)))*grfraction(1)
+ else
+ xmass1(j,ks)=0.
+ xscav_frac1(j,ks)=0.
+ endif
+ endif
+ enddo
+ endif
+
! Integrate Lagevin equation for lsynctime seconds
!*************************************************
+ 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))
+! write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
! Calculate the gross fluxes across layer interfaces
!***************************************************
diff --git a/src/verttransform.f90 b/src/verttransform.f90
index 403c22699c087d53e29adf3803479b9f22d005e5..ad89f2cd6c16642a619c3dc103bfa993c2f6c49e 100644
--- a/src/verttransform.f90
+++ b/src/verttransform.f90
@@ -90,6 +90,9 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
logical :: init = .true.
+ logical :: init_w = .false.
+ logical :: init_r = .false.
+
!ZHG SEP 2014 tests
! integer :: cloud_ver,cloud_min, cloud_max
@@ -102,8 +105,8 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
! character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
! CHARACTER(LEN=3) :: aspec
! integer :: virr=0
- real :: tot_cloud_h
- real :: dbg_height(nzmax)
+ !real :: tot_cloud_h
+ !real :: dbg_height(nzmax)
!ZHG
!*************************************************************************
@@ -122,6 +125,19 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
if (init) then
+
+ if (init_r) then
+
+ open(333,file='heights.txt', &
+ form='formatted')
+ do kz=1,nuvz
+ read(333,*) height(kz)
+ end do
+ close(333)
+ write(*,*) 'height read'
+ else
+
+
! Search for a point with high surface pressure (i.e. not above significant topography)
! Then, use this point to construct a reference z profile, to be used at all times
!*****************************************************************************
@@ -160,6 +176,16 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
pold(ixm,jym)=pint(ixm,jym)
end do
+ if (init_w) then
+ open(333,file='heights.txt', &
+ form='formatted')
+ do kz=1,nuvz
+ write(333,*) height(kz)
+ end do
+ close(333)
+ endif
+
+ endif ! init
! Determine highest levels that can be within PBL
!************************************************
@@ -177,7 +203,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
init=.false.
- dbg_height = height
+! dbg_height = height
endif
@@ -597,13 +623,13 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
lsp=lsprec(ix,jy,1,n)
convp=convprec(ix,jy,1,n)
prec=lsp+convp
- tot_cloud_h=0
+! tot_cloud_h=0
! Find clouds in the vertical
do kz=1, nz-1 !go from top to bottom
if (clwc(ix,jy,kz,n).gt.0) then
! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
- tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
+! tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
! icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n) ! Column cloud water [m3/m3]
ctwc(ix,jy,n) = ctwc(ix,jy,n)+clw(ix,jy,kz,n)
diff --git a/src/wetdepo.f90 b/src/wetdepo.f90
index 0ab6af6bcd70788c333d4e0b6b829e19c9144881..abda6351d8feacbfd201cf515bdc38864d5cc4d7 100644
--- a/src/wetdepo.f90
+++ b/src/wetdepo.f90
@@ -41,20 +41,13 @@ subroutine wetdepo(itime,ltsample,loutnext)
!*****************************************************************************
! *
! Variables: *
-! cc [0-1] total cloud cover *
-! convp [mm/h] convective precipitation rate *
-! grfraction [0-1] fraction of grid, for which precipitation occurs *
! ix,jy indices of output grid cell for each particle *
! itime [s] actual simulation time [s] *
! jpart particle index *
! ldeltat [s] interval since radioactive decay was computed *
-! lfr, cfr area fraction covered by precipitation for large scale *
-! and convective precipitation (dependent on prec. rate) *
! loutnext [s] time for which gridded deposition is next output *
! loutstep [s] interval at which gridded deposition is output *
-! lsp [mm/h] large scale precipitation rate *
! ltsample [s] interval over which mass is deposited *
-! prec [mm/h] precipitation rate in subgrid, where precipitation occurs*
! wetdeposit mass that is wet deposited *
! wetgrid accumulated deposited mass on output grid *
! wetscav scavenging coefficient *
@@ -70,29 +63,12 @@ subroutine wetdepo(itime,ltsample,loutnext)
implicit none
integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
- integer :: ngrid,itage,nage,hz,il,interp_time, n
- integer(kind=1) :: clouds_v
+ integer :: itage,nage,il,interp_time, n
integer :: ks, kp
-! integer :: n1,n2, icbot,ictop, indcloud !TEST
- real :: S_i, act_temp, cl, cle ! in cloud scavenging
- real :: clouds_h ! cloud height for the specific grid point
- real :: xtn,ytn,lsp,convp,cc,grfraction(3),prec(3),wetscav,totprec
+ integer(selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count
+ real :: grfraction(3),wetscav
real :: wetdeposit(maxspec),restmass
real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
-!save lfr,cfr
-
- real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
- real, parameter :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
-
-!ZHG aerosol below-cloud scavenging removal polynomial constants for rain and snow
- real, parameter :: bclr(6) = (/274.35758, 332839.59273, 226656.57259, 58005.91340, 6588.38582, 0.244984/) !rain (Laakso et al 2003)
- real, parameter :: bcls(6) = (/22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
- real :: frac_act, liq_frac, dquer_m
-
- integer(selected_int_kind(16)), dimension(nspec) :: blc_count, inc_count
- real :: Si_dummy, wetscav_dummy
- logical :: readclouds_this_nest
-
! Compute interval since radioactive decay of deposited mass was computed
!************************************************************************
@@ -118,284 +94,39 @@ subroutine wetdepo(itime,ltsample,loutnext)
if (itra1(jpart).lt.itime) goto 20
endif
-! Determine age class of the particle
- itage=abs(itra1(jpart)-itramem(jpart))
- do nage=1,nageclass
- if (itage.lt.lage(nage)) goto 33
- end do
-33 continue
-
-
-! Determine which nesting level to be used
-!*****************************************
-
- ngrid=0
- do j=numbnests,1,-1
- if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
- (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
- ngrid=j
- goto 23
- endif
- end do
-23 continue
-
-
-! Determine nested grid coordinates
-!**********************************
- readclouds_this_nest=.false.
-
- if (ngrid.gt.0) then
- xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
- ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
- ix=int(xtn)
- jy=int(ytn)
- if (readclouds_nest(ngrid)) readclouds_this_nest=.true.
- else
- ix=int(xtra1(jpart))
- jy=int(ytra1(jpart))
- endif
-
-
-! Interpolate large scale precipitation, convective precipitation and
-! total cloud cover
-! Note that interpolated time refers to itime-0.5*ltsample [PS]
-!********************************************************************
- interp_time=nint(itime-0.5*ltsample)
-
- if (ngrid.eq.0) then
- call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
- 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
- memtime(1),memtime(2),interp_time,lsp,convp,cc)
- else
- call interpol_rain_nests(lsprecn,convprecn,tccn, &
- nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
- memtime(1),memtime(2),interp_time,lsp,convp,cc)
- endif
-
-! If total precipitation is less than 0.01 mm/h - no scavenging occurs
- if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
-
-! get the level were the actual particle is in
- do il=2,nz
- if (height(il).gt.ztra1(jpart)) then
- hz=il-1
-! goto 26
- exit
- endif
- end do
-!26 continue
-
- n=memind(2)
- if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
- n=memind(1)
-
- if (ngrid.eq.0) then
- clouds_v=clouds(ix,jy,hz,n)
- clouds_h=cloudsh(ix,jy,n)
- else
- clouds_v=cloudsn(ix,jy,hz,n,ngrid)
- clouds_h=cloudshn(ix,jy,n,ngrid)
- endif
-
-! if there is no precipitation or the particle is above the clouds no
-! scavenging is done
-
- if (clouds_v.le.1) goto 20
-
-! 1) Parameterization of the the area fraction of the grid cell where the
-! precipitation occurs: the absolute limit is the total cloud cover, but
-! for low precipitation rates, an even smaller fraction of the grid cell
-! is used. Large scale precipitation occurs over larger areas than
-! convective precipitation.
-!**************************************************************************
-
- if (lsp.gt.20.) then
- i=5
- else if (lsp.gt.8.) then
- i=4
- else if (lsp.gt.3.) then
- i=3
- else if (lsp.gt.1.) then
- i=2
- else
- i=1
- endif
-
- if (convp.gt.20.) then
- j=5
- else if (convp.gt.8.) then
- j=4
- else if (convp.gt.3.) then
- j=3
- else if (convp.gt.1.) then
- j=2
- else
- j=1
- endif
-
-
-!ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp
-! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
-! for now they are treated the same
- grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
- grfraction(2)=max(0.05,cc*(lfr(i)))
- grfraction(3)=max(0.05,cc*(cfr(j)))
-
-
-! 2) Computation of precipitation rate in sub-grid cell
-!******************************************************
- prec(1)=(lsp+convp)/grfraction(1)
- prec(2)=(lsp)/grfraction(2)
- prec(3)=(convp)/grfraction(3)
-
-
-! 3) Computation of scavenging coefficients for all species
-! Computation of wet deposition
-!**********************************************************
-
- do ks=1,nspec ! loop over species
- wetdeposit(ks)=0.
- wetscav=0.
-
-! Cycle loop if wet deposition for the species is off
- if (WETDEPSPEC(ks).eqv..false.) cycle
-
- if (ngrid.gt.0) then
- act_temp=ttn(ix,jy,hz,n,ngrid)
- else
- act_temp=tt(ix,jy,hz,n)
- endif
-
-
-!***********************
-! BELOW CLOUD SCAVENGING
-!***********************
- if (clouds_v.ge.4) then !below cloud
-
-! For gas: if positive below-cloud parameters (A or B), and dquer<=0
+! Determine age class of the particle - nage is used for the kernel
!******************************************************************
- if ((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) then
- ! if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then
- blc_count(ks)=blc_count(ks)+1
- wetscav=weta_gas(ks)*prec(1)**wetb_gas(ks)
-
-! For aerosols: if positive below-cloud parameters (Crain/Csnow or B), and dquer>0
-!*********************************************************************************
- else if ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.)) then
- blc_count(ks)=blc_count(ks)+1
-
-!NIK 17.02.2015
-! For the calculation here particle size needs to be in meter and not um as dquer is
-! changed in readreleases
-! For particles larger than 10 um use the largest size defined in the parameterizations (10um)
- dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
-
-! Rain:
- if (act_temp .ge. 273. .and. crain_aero(ks).gt.0.) then
-
-! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003,
-! the below-cloud scavenging (rain efficienty) parameter Crain (=crain_aero) from SPECIES file
- wetscav=crain_aero(ks)*10**(bclr(1)+(bclr(2)*(log10(dquer_m))**(-4))+ &
- & (bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)*(log10(dquer_m))**(-2))+&
- &(bclr(5)*(log10(dquer_m))**(-1))+bclr(6)* (prec(1))**(0.5))
-
-! Snow:
- elseif (act_temp .lt. 273. .and. csnow_aero(ks).gt.0.) then
-! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009,
-! the below-cloud scavenging (Snow efficiency) parameter Csnow (=csnow_aero) from SPECIES file
- wetscav=csnow_aero(ks)*10**(bcls(1)+(bcls(2)*(log10(dquer_m))**(-4))+&
- &(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)*(log10(dquer_m))**(-2))+&
- &(bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
-
- endif
-
-! write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
-
- endif ! gas or particle
-! endif ! positive below-cloud scavenging parameters given in Species file
- endif !end BELOW
-
-!********************
-! IN CLOUD SCAVENGING
-!********************
- if (clouds_v.lt.4) then ! In-cloud
-! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are
-! given in species file, or if gas and positive Henry's constant
- if ((ccn_aero(ks).gt.0. .or. in_aero(ks).gt.0.).or.(henry(ks).gt.0.and.dquer(ks).le.0)) then
- inc_count(ks)=inc_count(ks)+1
-! if negative coefficients (turned off) set to zero for use in equation
- if (ccn_aero(ks).lt.0.) ccn_aero(ks)=0.
- if (in_aero(ks).lt.0.) in_aero(ks)=0.
-
-!ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
-! nested fields
- if (ngrid.gt.0.and.readclouds_this_nest) then
- cl=ctwcn(ix,jy,n,ngrid)*(grfraction(1)/cc)
- else if (ngrid.eq.0.and.readclouds) then
- cl=ctwc(ix,jy,n)*(grfraction(1)/cc)
- else !parameterize cloudwater m2/m3
-!ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
- cl=1.6E-6*prec(1)**0.36
- endif
-
-!ZHG: Calculate the partition between liquid and water phase water.
- if (act_temp .le. 253.) then
- liq_frac=0
- else if (act_temp .ge. 273.) then
- liq_frac=1
- else
- liq_frac =((act_temp-273.)/(273.-253.))**2.
- end if
-! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
- frac_act = liq_frac*ccn_aero(ks) +(1-liq_frac)*in_aero(ks)
-
-!ZHG Use the activated fraction and the liqid water to calculate the washout
-
-! AEROSOL
-!********
- if (dquer(ks).gt.0.) then
- S_i= frac_act/cl
-
-! GAS
-!****
- else
-
- cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
-!REPLACE to switch old/ new scheme
- ! S_i=frac_act/cle
- S_i=1/cle
- endif ! gas or particle
-
-! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
-!OLD
- if ((readclouds.and.ngrid.eq.0).or.(readclouds_this_nest.and.ngrid.gt.0)) then
- wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)
- else
- wetscav=incloud_ratio*S_i*(prec(1)/3.6E6)/clouds_h
- endif
+ itage=abs(itra1(jpart)-itramem(jpart))
+ do nage=1,nageclass
+ if (itage.lt.lage(nage)) goto 33
+ end do
+ 33 continue
+ do ks=1,nspec ! loop over species
- endif ! positive in-cloud scavenging parameters given in Species file
- endif !incloud
-!END ZHG TEST
+ if (WETDEPSPEC(ks).eqv..false.) cycle
!**************************************************
! CALCULATE DEPOSITION
!**************************************************
-! if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
-! + ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
+! wetscav=0.
+
+! write(*,*) ks,dquer(ks), crain_aero(ks),csnow_aero(ks)
+! if (((dquer(ks).le.0.).and.(weta_gas(ks).gt.0..or.wetb_gas(ks).gt.0.)) &
+! .or. &
+! ((dquer(ks).gt.0.).and.(crain_aero(ks).gt.0..or.csnow_aero(ks).gt.0.).or. &
+! (ccn_aero(ks).gt0) .or. (in_aero(ks).gt.0) .or. (henry(ks).gt.0))) then
+
+ call get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc_count,wetscav)
+
if (wetscav.gt.0.) then
wetdeposit(ks)=xmass1(jpart,ks)* &
(1.-exp(-wetscav*abs(ltsample)))*grfraction(1) ! wet deposition
-!write(*,*) 'MASS DEPOSITED: PREC, WETSCAV, WETSCAVP', prec(1), wetdeposit(ks), xmass1(jpart,ks)* &
-! (1.-exp(-wetscav_dummy*abs(ltsample)))*grfraction(1), clouds_v
-
-
else ! if no scavenging
wetdeposit(ks)=0.
endif
-
+
restmass = xmass1(jpart,ks)-wetdeposit(ks)
if (ioutputforeachrelease.eq.1) then
kp=npoint(jpart)
@@ -416,8 +147,8 @@ subroutine wetdepo(itime,ltsample,loutnext)
wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
endif
-
- end do !all species
+! endif ! no deposition
+ end do ! loop over species
! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
! Add the wet deposition to accumulated amount on output grid and nested output grid
diff --git a/src/wetdepokernel.f90 b/src/wetdepokernel.f90
index 6188068dad64bb710c1be9f9c7914d6d04e4baa3..65de00ad91f16c500482309417a469a3651d7f8c 100644
--- a/src/wetdepokernel.f90
+++ b/src/wetdepokernel.f90
@@ -94,33 +94,38 @@ subroutine wetdepokernel(nunc,deposit,x,y,nage,kp)
do ks=1,nspec
- if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
+ if (.not.usekernel) then
+ wetgridunc(ix,jy,ks,kp,nunc,nage)= &
+ wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)
+ else
+ if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
- w=wx*wy
+ w=wx*wy
wetgridunc(ix,jy,ks,kp,nunc,nage)= &
wetgridunc(ix,jy,ks,kp,nunc,nage)+deposit(ks)*w
- endif
+ endif
- if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
+ if ((ixp.ge.0).and.(jyp.ge.0).and.(ixp.le.numxgrid-1).and. &
(jyp.le.numygrid-1)) then
- w=(1.-wx)*(1.-wy)
+ w=(1.-wx)*(1.-wy)
wetgridunc(ixp,jyp,ks,kp,nunc,nage)= &
wetgridunc(ixp,jyp,ks,kp,nunc,nage)+deposit(ks)*w
- endif
+ endif
- if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
+ if ((ixp.ge.0).and.(jy.ge.0).and.(ixp.le.numxgrid-1).and. &
(jy.le.numygrid-1)) then
- w=(1.-wx)*wy
+ w=(1.-wx)*wy
wetgridunc(ixp,jy,ks,kp,nunc,nage)= &
wetgridunc(ixp,jy,ks,kp,nunc,nage)+deposit(ks)*w
- endif
+ endif
- if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
+ if ((ix.ge.0).and.(jyp.ge.0).and.(ix.le.numxgrid-1).and. &
(jyp.le.numygrid-1)) then
- w=wx*(1.-wy)
+ w=wx*(1.-wy)
wetgridunc(ix,jyp,ks,kp,nunc,nage)= &
wetgridunc(ix,jyp,ks,kp,nunc,nage)+deposit(ks)*w
- endif
+ endif
+ endif
end do
end if
diff --git a/src/writeprecip.f90 b/src/writeprecip.f90
new file mode 100644
index 0000000000000000000000000000000000000000..13bd1f83b4713e6794d9df0d307ba013972a8a81
--- /dev/null
+++ b/src/writeprecip.f90
@@ -0,0 +1,82 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
+! *
+! This file is part of FLEXPART. *
+! *
+! FLEXPART is free software: you can redistribute it and/or modify *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or *
+! (at your option) any later version. *
+! *
+! FLEXPART is distributed in the hope that it will be useful, *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+! GNU General Public License for more details. *
+! *
+! You should have received a copy of the GNU General Public License *
+! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
+!**********************************************************************
+
+subroutine writeprecip(itime,imem)
+
+ !*****************************************************************************
+ ! *
+ ! This routine produces a file containing total precipitation for each *
+ ! releases point. *
+ ! *
+ ! Author: S. Eckhardt *
+ ! 7 Mai 2017 *
+ !*****************************************************************************
+
+ use point_mod
+ use par_mod
+ use com_mod
+
+ implicit none
+
+ integer :: jjjjmmdd,ihmmss,itime,i
+ real(kind=dp) :: jul
+ character :: adate*8,atime*6
+
+ integer :: ix,jy,imem
+ real :: xp1,yp1
+
+
+ if (itime.eq.0) then
+ open(unitprecip,file=path(2)(1:length(2))//'wetscav_precip.txt', &
+ form='formatted',err=998)
+ else
+ open(unitprecip,file=path(2)(1:length(2))//'wetscav_precip.txt', &
+ ACCESS='APPEND',form='formatted',err=998)
+ endif
+
+ jul=bdate+real(itime,kind=dp)/86400._dp
+ call caldate(jul,jjjjmmdd,ihmmss)
+ write(adate,'(i8.8)') jjjjmmdd
+ write(atime,'(i6.6)') ihmmss
+
+ do i=1,numpoint
+ xp1=xpoint1(i)*dx+xlon0 !lat, long (real) coord
+ yp1=ypoint1(i)*dy+ylat0 !lat, long (real) coord
+ ix=int((xpoint1(i)+xpoint2(i))/2.)
+ jy=int((ypoint1(i)+ypoint2(i))/2.)
+ write(unitprecip,*) jjjjmmdd, ihmmss, &
+ xp1,yp1,lsprec(ix,jy,1,imem),convprec(ix,jy,1,imem) !time is the same as in the ECMWF windfield
+! units mm/h, valid for the time given in the windfield
+ end do
+
+ close(unitprecip)
+
+ return
+
+
+998 write(*,*) ' #### FLEXPART MODEL ERROR! THE FILE #### '
+ write(*,*) ' #### '//path(2)(1:length(2))//'header_txt'//' #### '
+ write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS #### '
+ write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### '
+ write(*,*) ' #### THE PROGRAM AGAIN. #### '
+ stop
+
+end subroutine writeprecip