Index: src/advance.f90
===================================================================
--- src/advance.f90	(revision a652cd5abe62c42ac456c3c62c1950b31af9ab94)
+++ src/advance.f90	(revision 94735e2096ac99bd12a8f8117f7920fd93e872e2)
@@ -117,4 +117,5 @@
   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)
@@ -223,4 +224,23 @@
 
 
+ ! 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
   !*******************************************************
@@ -229,10 +249,17 @@
   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)
@@ -249,4 +276,5 @@
   endif
 
+  if (interpolhmix) h=(h1(1)*dt2+h1(2)*dt1)*dtt 
   zeta=zt/h
 
@@ -446,4 +474,12 @@
         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
@@ -647,4 +683,10 @@
   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
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision c8fc72499e835239b052197045c24d51657bb79f)
+++ src/com_mod.f90	(revision 58882983dbcc65e87e1aa78634bff47e668268e9)
@@ -577,6 +577,8 @@
   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
@@ -598,4 +600,5 @@
   ! 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
 
@@ -668,4 +671,5 @@
   real, allocatable, dimension(:) :: ztra1 
   real, allocatable, dimension(:,:) :: xmass1
+  real, allocatable, dimension(:,:) :: xscav_frac1
 
   ! eso: Moved from timemanager
@@ -688,5 +692,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
+ 
 
 
@@ -750,4 +755,9 @@
   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
Index: src/conccalc.f90
===================================================================
--- src/conccalc.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ src/conccalc.f90	(revision 9287c012848e4e4bb0f37eb6dad76b516ed063d4)
@@ -21,5 +21,5 @@
 
 subroutine conccalc(itime,weight)
-  !                 i     i
+  !                      i     i
   !*****************************************************************************
   !                                                                            *
@@ -59,5 +59,5 @@
   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;
@@ -65,5 +65,5 @@
   ! releasepoints
   !***************************************************************************
-
+!  xscav_count=0
   do i=1,numpart
     if (itra1(i).ne.itime) goto 20
@@ -76,5 +76,6 @@
 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
   !************************************************************************
@@ -172,5 +173,4 @@
       if (yl.lt.0.) jy=jy-1
 
-  ! if (i.eq.10000) write(*,*) itime,xtra1(i),ytra1(i),ztra1(i),xl,yl
 
 
@@ -183,15 +183,23 @@
       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
@@ -220,20 +228,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 (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
-          endif
-        endif
+               end do
+             endif
+          endif
+        endif !ix ge 0
 
 
@@ -241,23 +265,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 (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
-          endif
-        endif
-      endif
-
-
+               end do
+            endif
+          endif
+        endif !ixp ge 0
+     endif
 
   !************************************
@@ -282,12 +320,20 @@
         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
 
@@ -319,18 +365,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 (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
@@ -340,26 +402,42 @@
             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
 
   !***********************************************************************
Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision 16b61a5e2d3b1020be1cff080e479ba26f3f7063)
+++ src/concoutput.f90	(revision 16b61a5e2d3b1020be1cff080e479ba26f3f7063)
@@ -246,21 +246,30 @@
 
     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
@@ -342,5 +351,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/drydepokernel.f90
===================================================================
--- src/drydepokernel.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ src/drydepokernel.f90	(revision 1c0d5e63f47af4e13fdf119fed31e0558305d2ac)
@@ -102,34 +102,38 @@
    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
Index: src/ecmwf_mod.f90
===================================================================
--- src/ecmwf_mod.f90	(revision 62e65c73ae7e3d8c4b8b884065ede76e07211ec8)
+++ src/ecmwf_mod.f90	(revision 842074e7416f59098921221d89faf6edc3639279)
@@ -43,6 +43,8 @@
 
 !  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
 
Index: src/get_vdep_prob.f90
===================================================================
--- src/get_vdep_prob.f90	(revision 92a74b26a4c380bde7d83f918a7d164605d5b263)
+++ src/get_vdep_prob.f90	(revision 92a74b26a4c380bde7d83f918a7d164605d5b263)
@@ -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
+
Index: src/get_wetscav.f90
===================================================================
--- src/get_wetscav.f90	(revision c7e771d7de021b45aeef132bb2857bf7e012e216)
+++ src/get_wetscav.f90	(revision c7e771d7de021b45aeef132bb2857bf7e012e216)
@@ -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
Index: src/getfields.f90
===================================================================
--- src/getfields.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/getfields.f90	(revision d1a8707797c5c3fabdbd530210b84b85fb18b258)
@@ -138,4 +138,8 @@
 40  indmin=indj
 
+   if (WETBKDEP) then
+        call writeprecip(itime,memind(1))
+   endif
+
   else
 
@@ -169,4 +173,8 @@
 60  indmin=indj
 
+   if (WETBKDEP) then
+        call writeprecip(itime,memind(1))
+   endif
+
   endif
 
Index: src/interpol_rain.f90
===================================================================
--- src/interpol_rain.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/interpol_rain.f90	(revision c7e771d7de021b45aeef132bb2857bf7e012e216)
@@ -21,5 +21,5 @@
 
 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
@@ -77,4 +77,5 @@
   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
 
 
@@ -113,35 +114,37 @@
   !***********************
 
-  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
Index: src/interpol_vdep.f90
===================================================================
--- src/interpol_vdep.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/interpol_vdep.f90	(revision 5844866080c2cda8fb3b9c8aa6db8739cd125deb)
@@ -54,5 +54,5 @@
 
   ! a) Bilinear horizontal interpolation
-
+! write(*,*) 'interpol: ',dt1,dt2,dtt,lsynctime,ix,jy
   do m=1,2
     indexh=memind(m)
Index: src/makefile
===================================================================
--- src/makefile	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ src/makefile	(revision d1a8707797c5c3fabdbd530210b84b85fb18b258)
@@ -68,5 +68,5 @@
 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
@@ -144,4 +144,5 @@
 advance.o		initialize.o		\
 writeheader.o		writeheader_txt.o  	\
+writeprecip.o \
 writeheader_surf.o     	assignland.o\
 part0.o 		gethourlyOH.o\
@@ -152,5 +153,6 @@
 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   	\
@@ -271,4 +273,6 @@
 
 ## 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
@@ -424,8 +428,10 @@
 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
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ src/par_mod.f90	(revision d1a8707797c5c3fabdbd530210b84b85fb18b258)
@@ -194,4 +194,5 @@
   integer,parameter :: maxpart=10000000
   integer,parameter :: maxspec=4
+
   real,parameter :: minmass=0.0001
 
@@ -258,5 +259,5 @@
   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
Index: src/readcommand.f90
===================================================================
--- src/readcommand.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ src/readcommand.f90	(revision 8ee24a5772873f9cec6d380cad07cb952da5d0f2)
@@ -298,4 +298,6 @@
   !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
@@ -320,9 +322,26 @@
      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
 
@@ -386,7 +405,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 OR 5 FOR        #### '
     write(*,*) ' #### STANDARD FLEXPART OUTPUT OR  9 - 13     #### '
     write(*,*) ' #### FOR NETCDF OUTPUT                       #### '
Index: src/readreleases.f90
===================================================================
--- src/readreleases.f90	(revision c8fc72499e835239b052197045c24d51657bb79f)
+++ src/readreleases.f90	(revision c8fc72499e835239b052197045c24d51657bb79f)
@@ -521,4 +521,18 @@
   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
   !*********************************************************************
Index: src/releaseparticles.f90
===================================================================
--- src/releaseparticles.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/releaseparticles.f90	(revision 75a4dedbe36ebfba4e9c67fce8af8641dafcd8eb)
@@ -176,9 +176,16 @@
   ! 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
   !**************************************
@@ -202,5 +209,4 @@
 
             ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux
-
   ! Interpolation of topography and density
   !****************************************
@@ -330,5 +336,5 @@
   !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
@@ -375,14 +381,13 @@
             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
 
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision 18adf6049d9d47356934a37eb7784506f795e220)
+++ src/timemanager.f90	(revision d1a8707797c5c3fabdbd530210b84b85fb18b258)
@@ -105,7 +105,7 @@
 ! 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)
@@ -114,5 +114,7 @@
   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),
@@ -145,4 +147,7 @@
 !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
 
@@ -541,10 +546,53 @@
         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
Index: src/verttransform.f90
===================================================================
--- src/verttransform.f90	(revision 341f4b7eb1c45f4749b96fc8dbea2dfff30140f2)
+++ src/verttransform.f90	(revision 6481010989731330559aee3129d788f71fc7948f)
@@ -91,4 +91,7 @@
 
   logical :: init = .true.
+  logical :: init_w = .false.
+  logical :: init_r = .false.
+
 
   !ZHG SEP 2014 tests  
@@ -103,6 +106,6 @@
   ! CHARACTER(LEN=3)  :: aspec
   ! integer :: virr=0
-  real :: tot_cloud_h
-  real :: dbg_height(nzmax) 
+  !real :: tot_cloud_h
+  !real :: dbg_height(nzmax) 
 !ZHG
 
@@ -123,4 +126,17 @@
   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
@@ -161,4 +177,14 @@
     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
@@ -178,5 +204,5 @@
     init=.false.
 
-    dbg_height = height
+!    dbg_height = height
 
   endif
@@ -598,5 +624,5 @@
         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
@@ -604,5 +630,5 @@
 ! 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]
Index: src/wetdepo.f90
===================================================================
--- src/wetdepo.f90	(revision c8fc72499e835239b052197045c24d51657bb79f)
+++ src/wetdepo.f90	(revision c8fc72499e835239b052197045c24d51657bb79f)
@@ -42,18 +42,11 @@
 !                                                                            *
 ! 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               *
@@ -71,27 +64,10 @@
 
   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
@@ -119,282 +95,37 @@
     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
-!**********************************************************
+! Determine age class of the particle - nage is used for the kernel
+!******************************************************************
+     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
-      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
-!******************************************************************
-        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
-
-
-        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
@@ -417,6 +148,6 @@
       endif
 
-
-    end do !all species
+!    endif ! no deposition
+    end do ! loop over species
 
 ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
Index: src/wetdepokernel.f90
===================================================================
--- src/wetdepokernel.f90	(revision 4c64400c622119505b289ef83be130b0b81b9590)
+++ src/wetdepokernel.f90	(revision 1c0d5e63f47af4e13fdf119fed31e0558305d2ac)
@@ -95,31 +95,36 @@
   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
Index: src/writeprecip.f90
===================================================================
--- src/writeprecip.f90	(revision None)
+++ src/writeprecip.f90	(revision None)
@@ -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
