compiles after merge scavenging into test dev

src/advance.f90

@@ -116,6 +116,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
   real :: dcwsave
   real :: usigold,vsigold,wsigold,r,rs
   real :: uold,vold,wold,vdepo(maxspec)
+  real :: h1(2)
   !real uprof(nzmax),vprof(nzmax),wprof(nzmax)
   !real usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
   !real rhoprof(nzmax),rhogradprof(nzmax)
@@ -222,18 +223,44 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
   jyp=jy+1
 
 
+ ! Determine the lower left corner and its distance to the current position
+  !*************************************************************************
+
+  ddx=xt-real(ix)
+  ddy=yt-real(jy)
+  rddx=1.-ddx
+  rddy=1.-ddy
+  p1=rddx*rddy
+  p2=ddx*rddy
+  p3=rddx*ddy
+  p4=ddx*ddy
+
+ ! Calculate variables for time interpolation
+  !*******************************************
+
+  dt1=real(itime-memtime(1))
+  dt2=real(memtime(2)-itime)
+  dtt=1./(dt1+dt2)
+
   ! Compute maximum mixing height around particle position
   !*******************************************************
 
   h=0.
   if (ngrid.le.0) then
     do k=1,2
-      mind=memind(k) ! eso: compatibility with 3-field version
-      do j=jy,jyp
-        do i=ix,ixp
-          if (hmix(i,j,1,mind).gt.h) h=hmix(i,j,1,mind)
-        end do
-      end do
+       mind=memind(k) ! eso: compatibility with 3-field version
+       if (interpolhmix) then
+             h1(k)=p1*hmix(ix ,jy ,1,mind) &
+                 + p2*hmix(ixp,jy ,1,mind) &
+                 + p3*hmix(ix ,jyp,1,mind) &
+                 + p4*hmix(ixp,jyp,1,mind)
+        else
+          do j=jy,jyp
+            do i=ix,ixp
+               if (hmix(i,j,1,mind).gt.h) h=hmix(i,j,1,mind)
+            end do
+          end do
+        endif
     end do
     tropop=tropopause(nix,njy,1,1)
   else
@@ -248,6 +275,7 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
     tropop=tropopausen(nix,njy,1,1,ngrid)
   endif
 
+  if (interpolhmix) h=(h1(1)*dt2+h1(2)*dt1)*dtt 
   zeta=zt/h
 
 
@@ -445,6 +473,14 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
           delz=wp*dtf
         endif
 
+        if (turboff) then
+!sec switch off turbulence
+          up=0.0
+          vp=0.0
+          wp=0.0
+          delz=0.
+        endif
+
   !****************************************************
   ! Compute turbulent vertical displacement of particle
   !****************************************************
@@ -646,6 +682,12 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
     nrand=nrand+1
   endif
 
+  if (turboff) then
+!sec switch off turbulence
+    ux=0.0
+    vy=0.0
+    wp=0.0
+  endif
 
   ! If particle represents only a single species, add gravitational settling
   ! velocity. The settling velocity is zero for gases

src/com_mod.f90

@@ -576,8 +576,10 @@ module com_mod
   integer :: numxgridn,numygridn
   real :: dxoutn,dyoutn,outlon0n,outlat0n,xoutshiftn,youtshiftn
   !real outheight(maxzgrid),outheighthalf(maxzgrid)
+
   logical :: DEP,DRYDEP,DRYDEPSPEC(maxspec),WETDEP,WETDEPSPEC(maxspec),&
        & OHREA,ASSSPEC
+  logical :: DRYBKDEP,WETBKDEP
 
   ! numxgrid,numygrid       number of grid points in x,y-direction
   ! numxgridn,numygridn     number of grid points in x,y-direction for nested output grid
@@ -597,6 +599,7 @@ module com_mod
   ! WETDEPSPEC              .true., if wet deposition is switched on for that species
   ! OHREA                   .true., if OH reaction is switched on
   ! ASSSPEC                 .true., if there are two species asscoiated
+  ! DRYBKDEP,WETBKDEP        .true., for bkwd runs, where mass deposited and source regions is calculated - either for dry or for wet deposition
   !                    (i.e. transfer of mass between these two occurs
 
 
@@ -667,6 +670,7 @@ module com_mod
   real(kind=dp), allocatable, dimension(:) :: xtra1, ytra1
   real, allocatable, dimension(:) :: ztra1 
   real, allocatable, dimension(:,:) :: xmass1
+  real, allocatable, dimension(:,:) :: xscav_frac1
 
   ! eso: Moved from timemanager
   real, allocatable, dimension(:) :: uap,ucp,uzp,us,vs,ws
@@ -687,7 +691,8 @@ module com_mod
   ! numparticlecount        counts the total number of particles that have been released
   ! xtra1,ytra1,ztra1       spatial positions of the particles
   ! xmass1 [kg]             particle masses
-  
+  ! xscav_frac1             fraction of particle masse which has been scavenged at receptor
+ 
 
 
   !*******************************************************
@@ -750,6 +755,11 @@ module com_mod
   integer :: mpi_mode=0 ! .gt. 0 if running MPI version
   logical :: lroot=.true. ! true if serial version, or if MPI .and. root process
   
+  logical :: usekernel=.false.    ! true if the output kernel shall be switched on
+  logical :: interpolhmix=.false. ! true if the hmix shall be interpolated
+  logical :: turboff=.false.       ! true if the turbulence shall be switched off
+  
+  
 contains
   subroutine com_mod_allocate_part(nmpart)
   !*******************************************************************************    

src/conccalc.f90

@@ -20,7 +20,7 @@
 !**********************************************************************
 
 subroutine conccalc(itime,weight)
-  !                 i     i
+  !                      i     i
   !*****************************************************************************
   !                                                                            *
   !     Calculation of the concentrations on a regular grid using volume       *
@@ -58,13 +58,13 @@ subroutine conccalc(itime,weight)
   real :: rhoprof(2),rhoi
   real :: xl,yl,wx,wy,w
   real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
-
+!  integer xscav_count
 
   ! For forward simulations, make a loop over the number of species;
   ! for backward simulations, make an additional loop over the
   ! releasepoints
   !***************************************************************************
-
+!  xscav_count=0
   do i=1,numpart
     if (itra1(i).ne.itime) goto 20
 
@@ -75,7 +75,8 @@ subroutine conccalc(itime,weight)
     end do
 33   continue
 
-
+!  if (xscav_frac1(i,1).lt.0) xscav_count=xscav_count+1
+           
   ! For special runs, interpolate the air density to the particle position
   !************************************************************************
   !***********************************************************************
@@ -171,7 +172,6 @@ subroutine conccalc(itime,weight)
       jy=int(yl)
       if (yl.lt.0.) jy=jy-1
 
-  ! if (i.eq.10000) write(*,*) itime,xtra1(i),ytra1(i),ztra1(i),xl,yl
 
 
   ! For particles aged less than 3 hours, attribute particle mass to grid cell
@@ -182,17 +182,25 @@ subroutine conccalc(itime,weight)
 
       if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
            (xl.gt.real(numxgrid-1)-0.5).or. &
-           (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
+           (yl.gt.real(numygrid-1)-0.5))) then             ! no kernel, direct attribution to grid cell
         if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
              (jy.le.numygrid-1)) then
-          do ks=1,nspec
-            gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+          if (DRYBKDEP.or.WETBKDEP) then
+             do ks=1,nspec
+               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                 xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
+             end do
+          else
+             do ks=1,nspec
+               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                  gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                  xmass1(i,ks)/rhoi*weight
-          end do
+             end do
+          endif
         endif
 
-      else                                 ! attribution via uniform kernel
+      else                                 ! attribution via uniform kernel 
 
         ddx=xl-real(ix)                   ! distance to left cell border
         ddy=yl-real(jy)                   ! distance to lower cell border
@@ -219,46 +227,76 @@ subroutine conccalc(itime,weight)
         if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
           if ((jy.ge.0).and.(jy.le.numygrid-1)) then
             w=wx*wy
-            do ks=1,nspec
-              gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+               do ks=1,nspec
+                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
+               end do
+            else
+               do ks=1,nspec
+                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                    gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                    xmass1(i,ks)/rhoi*weight*w
-            end do
+               end do
+            endif
           endif
 
           if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
             w=wx*(1.-wy)
-            do ks=1,nspec
-              gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+              do ks=1,nspec
+                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+               end do
+             else
+              do ks=1,nspec
+                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                    gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                    xmass1(i,ks)/rhoi*weight*w
-            end do
+               end do
+             endif
           endif
-        endif
+        endif !ix ge 0
 
 
         if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
           if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
             w=(1.-wx)*(1.-wy)
-            do ks=1,nspec
-              gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+               do ks=1,nspec
+                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
+               end do
+            else
+               do ks=1,nspec
+                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                    gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                    xmass1(i,ks)/rhoi*weight*w
-            end do
+               end do
+            endif
           endif
 
           if ((jy.ge.0).and.(jy.le.numygrid-1)) then
             w=(1.-wx)*wy
-            do ks=1,nspec
-              gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+               do ks=1,nspec
+                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+               end do
+            else
+               do ks=1,nspec
+                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                    gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                    xmass1(i,ks)/rhoi*weight*w
-            end do
+               end do
+            endif
           endif
-        endif
-      endif
-
-
+        endif !ixp ge 0
+     endif
 
   !************************************
   ! Do everything for the nested domain
@@ -281,14 +319,22 @@ subroutine conccalc(itime,weight)
 
         if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
              (xl.gt.real(numxgridn-1)-0.5).or. &
-             (yl.gt.real(numygridn-1)-0.5)) then             ! no kernel, direct attribution to grid cell
+             (yl.gt.real(numygridn-1)-0.5).or.(.not.usekernel)) then             ! no kernel, direct attribution to grid cell
           if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
                (jy.le.numygridn-1)) then
-            do ks=1,nspec
-              griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+            if (DRYBKDEP.or.WETBKDEP) then
+               do ks=1,nspec
+                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
+               end do
+            else
+               do ks=1,nspec
+                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                    griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                    xmass1(i,ks)/rhoi*weight
-            end do
+               end do
+            endif
           endif
 
         else                                 ! attribution via uniform kernel
@@ -318,20 +364,36 @@ subroutine conccalc(itime,weight)
           if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
             if ((jy.ge.0).and.(jy.le.numygridn-1)) then
               w=wx*wy
-              do ks=1,nspec
-                griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+              if (DRYBKDEP.or.WETBKDEP) then
+                 do ks=1,nspec
+                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+                 end do
+              else
+                do ks=1,nspec
+                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                      griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                      xmass1(i,ks)/rhoi*weight*w
-              end do
+                 end do
+              endif
             endif
 
             if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
               w=wx*(1.-wy)
-              do ks=1,nspec
-                griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+              if (DRYBKDEP.or.WETBKDEP) then
+                 do ks=1,nspec
+                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+                 end do
+              else
+                 do ks=1,nspec
+                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                      griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                      xmass1(i,ks)/rhoi*weight*w
-              end do
+                 end do
+              endif
             endif
           endif
 
@@ -339,28 +401,44 @@ subroutine conccalc(itime,weight)
           if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
             if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
               w=(1.-wx)*(1.-wy)
-              do ks=1,nspec
-                griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+              if (DRYBKDEP.or.WETBKDEP) then
+                 do ks=1,nspec
+                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
+                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+                 end do
+              else
+                 do ks=1,nspec
+                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                      griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                      xmass1(i,ks)/rhoi*weight*w
-              end do
+                 end do
+              endif
             endif
 
             if ((jy.ge.0).and.(jy.le.numygridn-1)) then
               w=(1.-wx)*wy
-              do ks=1,nspec
-                griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+              if (DRYBKDEP.or.WETBKDEP) then
+                 do ks=1,nspec
+                   griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
+                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
+                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
+                 end do
+              else
+                 do ks=1,nspec
+                    griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                      griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                      xmass1(i,ks)/rhoi*weight*w
-              end do
+                 end do
+              endif
             endif
           endif
         endif
-
       endif
     endif
 20  continue
   end do
+!  write(*,*) 'xscav count:',xscav_count
 
   !***********************************************************************
   ! 2. Evaluate concentrations at receptor points, using the kernel method

src/concoutput.f90

@@ -245,23 +245,32 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
   do ks=1,nspec
 
     write(anspec,'(i3.3)') ks
-    if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
-      if (ldirect.eq.1) then
-        open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
+
+    if (DRYBKDEP.or.WETBKDEP) then !scavdep output
+      if (DRYBKDEP) & 
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+      if (WETBKDEP) & 
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
+           atime//'_'//anspec,form='unformatted')