first version of the backward scavenging

src/com_mod.f90

@@ -575,7 +575,7 @@ 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
+  logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,OHREA,ASSSPEC,SCAVDEP
 
   ! 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,6 +594,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
   !                    (i.e. transfer of mass between these two occurs
 
 
@@ -664,6 +665,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
@@ -684,7 +686,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
+ 
 
 
   !*******************************************************

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       *
@@ -60,6 +60,12 @@ subroutine conccalc(itime,weight)
   real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
 
 
+  integer :: usekernel
+
+  usekernel=0
+  if (usekernel.ne.1) then
+     write (*,*) 'NOT USING THE KERNEL!'
+  endif
   ! For forward simulations, make a loop over the number of species;
   ! for backward simulations, make an additional loop over the
   ! releasepoints
@@ -180,16 +186,24 @@ subroutine conccalc(itime,weight)
   ! If a particle is close to the domain boundary, do not use the kernel either.
   !*****************************************************************************
 
-      if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
+      if (((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)).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
-          do ks=1,nspec
-            gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+          if (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*weight*xscav_frac1(i,ks)
+             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
@@ -219,46 +233,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 (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*w*weight*xscav_frac1(i,ks)
+               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 (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks)
+               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 (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*w*weight*xscav_frac1(i,ks)
+               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 (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks)
+               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
@@ -284,11 +328,19 @@ 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
-            do ks=1,nspec
-              griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*weight*xscav_frac1(i,ks)
+               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 +370,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 (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks)
+                 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 (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks)
+                 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,24 +407,39 @@ 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 (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks)
+                 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 (SCAVDEP) 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)/(1-xscav_frac1(i,ks))/rhoi*weight*w*xscav_frac1(i,ks)
+                 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

src/concoutput.f90

@@ -263,6 +263,12 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
       write(unitoutgridppt) itime
     endif
 
+    if (iout.eq.6) then !scavdep output
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_scav_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+      write(unitoutgrid) itime
+     endif
+
     do kp=1,maxpointspec_act
       do nage=1,nageclass
 
@@ -341,7 +347,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

src/readcommand.f90

@@ -338,6 +338,21 @@ 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
   !******************
 
@@ -383,9 +398,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, 5 OR 6 FOR     #### '
     write(*,*) ' #### STANDARD FLEXPART OUTPUT OR  9 - 13     #### '
     write(*,*) ' #### FOR NETCDF OUTPUT                       #### '
     stop

src/releaseparticles.f90

@@ -56,6 +56,7 @@ 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
@@ -85,6 +86,9 @@ 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
 
@@ -175,11 +179,24 @@ 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 (SCAVDEP) 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
   !**************************************
             end do
@@ -371,6 +388,12 @@ 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
 
@@ -378,12 +401,22 @@ subroutine releaseparticles(itime)
             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
+      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

@@ -375,7 +375,7 @@ subroutine timemanager
   !*****************************************************
 
       if ((itime.eq.loutend).and.(outnum.gt.0.)) then
-        if ((iout.le.3.).or.(iout.eq.5)) then
+        if ((iout.le.3.).or.(iout.eq.5).or.(iout.eq.6)) then
           if (surf_only.ne.1) then 
             if (lnetcdfout.eq.1) then 
               call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)

src/wetdepo.f90

@@ -173,7 +173,20 @@ subroutine wetdepo(itime,ltsample,loutnext)
     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
+!  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
@@ -413,6 +426,22 @@ subroutine wetdepo(itime,ltsample,loutnext)
         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