Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision 62e65c73ae7e3d8c4b8b884065ede76e07211ec8)
+++ src/com_mod.f90	(revision 462f74bead132c09c0766e7de3ced7c1ef50dc58)
@@ -576,5 +576,5 @@
   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
@@ -595,4 +595,5 @@
   ! 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
 
@@ -665,4 +666,5 @@
   real, allocatable, dimension(:) :: ztra1 
   real, allocatable, dimension(:,:) :: xmass1
+  real, allocatable, dimension(:,:) :: xscav_frac1
 
   ! eso: Moved from timemanager
@@ -685,5 +687,6 @@
   ! xtra1,ytra1,ztra1       spatial positions of the particles
   ! xmass1 [kg]             particle masses
-  
+  ! xscav_frac1             fraction of particle masse which has been scavenged at receptor
+ 
 
 
Index: src/conccalc.f90
===================================================================
--- src/conccalc.f90	(revision fdc0f03d0fd969385acd98d3a3206e9d802172bb)
+++ src/conccalc.f90	(revision 462f74bead132c09c0766e7de3ced7c1ef50dc58)
@@ -21,5 +21,5 @@
 
 subroutine conccalc(itime,weight)
-  !                 i     i
+  !                      i     i
   !*****************************************************************************
   !                                                                            *
@@ -61,4 +61,10 @@
 
 
+  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
@@ -181,14 +187,22 @@
   !*****************************************************************************
 
-      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
 
@@ -220,20 +234,36 @@
           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
-          endif
-        endif
+               end do
+             endif
+          endif
+        endif !ix ge 0
 
 
@@ -241,23 +271,37 @@
           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
-          endif
-        endif
-      endif
-
-
+               end do
+            endif
+          endif
+        endif !ixp ge 0
+     endif
 
   !************************************
@@ -285,9 +329,17 @@
           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
 
@@ -319,18 +371,34 @@
             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
@@ -340,22 +408,37 @@
             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
Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision db712a8ffdd9c478928d0173162055a7c109b178)
+++ src/concoutput.f90	(revision 462f74bead132c09c0766e7de3ced7c1ef50dc58)
@@ -264,4 +264,10 @@
     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
@@ -342,5 +348,5 @@
 ! 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
Index: src/readcommand.f90
===================================================================
--- src/readcommand.f90	(revision c04b739b567dd89ec6a4491a728fbe13d3bd0b68)
+++ src/readcommand.f90	(revision 462f74bead132c09c0766e7de3ced7c1ef50dc58)
@@ -339,4 +339,19 @@
   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
   !******************
@@ -384,7 +399,7 @@
   !**********************************************************************
 
-  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                       #### '
Index: src/releaseparticles.f90
===================================================================
--- src/releaseparticles.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/releaseparticles.f90	(revision 462f74bead132c09c0766e7de3ced7c1ef50dc58)
@@ -57,4 +57,5 @@
   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
@@ -86,4 +87,7 @@
   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
@@ -176,9 +180,22 @@
   ! 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
   !**************************************
@@ -372,4 +389,10 @@
               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
@@ -379,10 +402,20 @@
             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
 
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision 18adf6049d9d47356934a37eb7784506f795e220)
+++ src/timemanager.f90	(revision 462f74bead132c09c0766e7de3ced7c1ef50dc58)
@@ -376,5 +376,5 @@
 
       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 
Index: src/wetdepo.f90
===================================================================
--- src/wetdepo.f90	(revision 05cf28de9ff073ff7e3cdc064f0ed59c61a6dd09)
+++ src/wetdepo.f90	(revision 462f74bead132c09c0766e7de3ced7c1ef50dc58)
@@ -174,5 +174,18 @@
 
 !  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
@@ -414,4 +427,20 @@
       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
