Index: options/COMMAND
===================================================================
--- options/COMMAND	(revision 60403cdfd1dfc5def4c83aff4c9c55a236804c93)
+++ options/COMMAND	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -1,8 +1,8 @@
 &COMMAND
  LDIRECT=          1,
- IBDATE=   20110310,
- IBTIME=     120000,
- IEDATE=   20110311,
- IETIME=     040000,
+ IBDATE=   20140911,
+ IBTIME=     050000,
+ IEDATE=   20140911,
+ IETIME=     230000,
  LOUTSTEP=       3600,
  LOUTAVER=       3600,
Index: options/OUTGRID_NEST
===================================================================
--- options/OUTGRID_NEST	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ options/OUTGRID_NEST	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -7,9 +7,9 @@
 
 1.  ------.----       4X,F11.4
-       -15.0000       GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
+       -85.0000       GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
     OUTLONLEFT        (left boundary of the first grid cell - not its centre)
 
 2.  ------.----       4X,F11.4
-        35.0000       GEOGRAFICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
+        30.0000       GEOGRAFICAL LATITUDE OF LOWER LEFT CORNER OF OUTPUT GRID
     OUTLATLOWER       (lower boundary of the first grid cell - not its centre)
 
Index: options/RELEASES
===================================================================
--- options/RELEASES	(revision 60403cdfd1dfc5def4c83aff4c9c55a236804c93)
+++ options/RELEASES	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -4,28 +4,12 @@
  /
 &RELEASE
- IDATE1=   20110310,
- ITIME1=     120000,
- IDATE2=   20110310,
- ITIME2=     130000,
- LON1=  15.004    ,
- LON2=  15.054    ,
- LAT1=  44.734    ,
- LAT2=  44.784    ,
- Z1=  1500.000000    ,
- Z2=  3000.000000    ,
- ZKIND=     1,
- MASS=  2.1000E8    ,
- PARTS=   2000000
- COMMENT="TEST1                          ",
- /
-&RELEASE
- IDATE1=   20110310,
- ITIME1=     120000,
- IDATE2=   20110310,
- ITIME2=     130000,
- LON1=  20.004    ,
- LON2=  20.054    ,
- LAT1=  44.734    ,
- LAT2=  44.784    ,
+ IDATE1=   20140911,
+ ITIME1=     070000,
+ IDATE2=   20140911,
+ ITIME2=     080000,
+ LON1=  -75.250    ,
+ LON2=  -75.250    ,
+ LAT1=  38.000    ,
+ LAT2=  38.000    ,
  Z1=  1500.000000    ,
  Z2=  3000.000000    ,
@@ -33,4 +17,20 @@
  MASS=  2.0000E8    ,
  PARTS=   2000000
+ COMMENT="TEST1                          ",
+ /
+&RELEASE
+ IDATE1=   20140911,
+ ITIME1=     070000,
+ IDATE2=   20140911,
+ ITIME2=     080000,
+ LON1=  -79.250   ,
+ LON2=  -79.250   ,
+ LAT1=  38.100    ,
+ LAT2=  38.100    ,
+ Z1=  100.000000    ,
+ Z2=  1000.000000    ,
+ ZKIND=     1,
+ MASS=  2.0000E8    ,
+ PARTS=   2000000
  COMMENT="TEST2                          ",
  /
Index: src/FLEXPART_MPI.f90
===================================================================
--- src/FLEXPART_MPI.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/FLEXPART_MPI.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -424,10 +424,23 @@
 
 ! NIK 16.02.2005 
-  write(*,*) '**********************************************'
-  write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count
-  write(*,*) 'Total number of occurences of in-cloud    scavenging', tot_inc_count
-  write(*,*) '**********************************************'
-
   if (lroot) then
+! eso TODO: do MPI_Reduce (sum) of total occurences across processes
+    call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, 1, mp_pp, MPI_SUM, id_root, &
+         & mp_comm_used, mp_ierr)
+    call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, 1, mp_pp, MPI_SUM, id_root, &
+         & mp_comm_used, mp_ierr)
+  else
+    call MPI_Reduce(tot_blc_count, tot_blc_count, 1, mp_pp, MPI_SUM, id_root, &
+         & mp_comm_used, mp_ierr)
+    call MPI_Reduce(tot_inc_count, tot_inc_count, 1, mp_pp, MPI_SUM, id_root, &
+         & mp_comm_used, mp_ierr)
+  end if
+
+  if (lroot) then
+    write(*,*) '**********************************************'
+    write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count
+    write(*,*) 'Total number of occurences of in-cloud    scavenging', tot_inc_count
+    write(*,*) '**********************************************'
+
     write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
          &XPART MODEL RUN!'
Index: src/advance.f90
===================================================================
--- src/advance.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/advance.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -229,6 +229,5 @@
   if (ngrid.le.0) then
     do k=1,2
-! eso: compatibility with 3-field version
-      mind=memind(k)
+      mind=memind(k) ! eso: compatibility with 3-field version
       do j=jy,jyp
         do i=ix,ixp
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/com_mod.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -269,6 +269,5 @@
   !                         disc or on tape
 
-  integer :: memtime(numwfmem),memind(numwfmem) ! eso: or memind(3) and change
-                                                ! interpol_rain 
+  integer :: memtime(numwfmem),memind(3) ! eso: or memind(numwfmem) 
 
   ! memtime [s]             validation times of wind fields in memory
Index: src/conccalc_mpi.f90
===================================================================
--- src/conccalc_mpi.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/conccalc_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -62,5 +62,6 @@
   real :: xl,yl,wx,wy,w
   real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
-
+  integer :: mind2
+  ! mind2        eso: pointer to 2nd windfield in memory
 
   ! For forward simulations, make a loop over the number of species;
@@ -70,4 +71,6 @@
 
   if (mp_measure_time) call mpif_mtime('conccalc',0)
+
+  mind2=memind(2)
 
   do i=1,numpart
@@ -131,8 +134,13 @@
   !*****************************************************************************
       do ind=indz,indzp
-        rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
-             +p2*rho(ixp,jy ,ind,2) &
-             +p3*rho(ix ,jyp,ind,2) &
-             +p4*rho(ixp,jyp,ind,2)
+        ! rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
+        !      +p2*rho(ixp,jy ,ind,2) &
+        !      +p3*rho(ix ,jyp,ind,2) &
+        !      +p4*rho(ixp,jyp,ind,2)
+        rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,mind2) &
+             +p2*rho(ixp,jy ,ind,mind2) &
+             +p3*rho(ix ,jyp,ind,mind2) &
+             +p4*rho(ixp,jyp,ind,mind2)
+
       end do
       rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
Index: src/concoutput.f90
===================================================================
--- src/concoutput.f90	(revision f13406cf18bf27ac76cdb58b9abebb11329d6757)
+++ src/concoutput.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -99,4 +99,6 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  integer :: mind
+! mind        eso:added to ensure identical results between 2&3-fields versions
 
 
@@ -140,4 +142,5 @@
   !*******************************************************************
 
+  mind=memind(2)
   do kz=1,numzgrid
     if (kz.eq.1) then
@@ -162,17 +165,20 @@
         iix=max(min(nint(xl),nxmin1),0)
         jjy=max(min(nint(yl),nymin1),0)
-        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
-             rho(iix,jjy,kzz-1,2)*dz2)/dz
+        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
+        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
+        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
+             rho(iix,jjy,kzz-1,mind)*dz2)/dz
       end do
     end do
   end do
 
-    do i=1,numreceptor
-      xl=xreceptor(i)
-      yl=yreceptor(i)
-      iix=max(min(nint(xl),nxmin1),0)
-      jjy=max(min(nint(yl),nymin1),0)
-      densityoutrecept(i)=rho(iix,jjy,1,2)
-    end do
+  do i=1,numreceptor
+    xl=xreceptor(i)
+    yl=yreceptor(i)
+    iix=max(min(nint(xl),nxmin1),0)
+    jjy=max(min(nint(yl),nymin1),0)
+    !densityoutrecept(i)=rho(iix,jjy,1,2)
+    densityoutrecept(i)=rho(iix,jjy,1,mind)
+  end do
 
 
Index: src/concoutput_mpi.f90
===================================================================
--- src/concoutput_mpi.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/concoutput_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -104,4 +104,6 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  integer :: mind
+! mind        eso:added to ensure identical results between 2&3-fields versions
 
 ! Measure execution time
@@ -148,4 +150,5 @@
 !*******************************************************************
 
+  mind=memind(2)
   do kz=1,numzgrid
     if (kz.eq.1) then
@@ -170,6 +173,8 @@
         iix=max(min(nint(xl),nxmin1),0)
         jjy=max(min(nint(yl),nymin1),0)
-        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
-             rho(iix,jjy,kzz-1,2)*dz2)/dz
+        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
+        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
+        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
+             rho(iix,jjy,kzz-1,mind)*dz2)/dz
       end do
     end do
@@ -181,5 +186,6 @@
     iix=max(min(nint(xl),nxmin1),0)
     jjy=max(min(nint(yl),nymin1),0)
-    densityoutrecept(i)=rho(iix,jjy,1,2)
+    !densityoutrecept(i)=rho(iix,jjy,1,2)
+    densityoutrecept(i)=rho(iix,jjy,1,mind)
   end do
 
Index: src/concoutput_nest.f90
===================================================================
--- src/concoutput_nest.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/concoutput_nest.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -95,4 +95,6 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  integer :: mind 
+! mind        eso:added to ensure identical results between 2&3-fields versions
 
 
@@ -134,4 +136,5 @@
   !*******************************************************************
 
+  mind=memind(2)
   do kz=1,numzgrid
     if (kz.eq.1) then
@@ -156,17 +159,20 @@
         iix=max(min(nint(xl),nxmin1),0)
         jjy=max(min(nint(yl),nymin1),0)
-        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
-             rho(iix,jjy,kzz-1,2)*dz2)/dz
+        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
+        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
+        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
+             rho(iix,jjy,kzz-1,mind)*dz2)/dz
       end do
     end do
   end do
 
-    do i=1,numreceptor
-      xl=xreceptor(i)
-      yl=yreceptor(i)
-      iix=max(min(nint(xl),nxmin1),0)
-      jjy=max(min(nint(yl),nymin1),0)
-      densityoutrecept(i)=rho(iix,jjy,1,2)
-    end do
+  do i=1,numreceptor
+    xl=xreceptor(i)
+    yl=yreceptor(i)
+    iix=max(min(nint(xl),nxmin1),0)
+    jjy=max(min(nint(yl),nymin1),0)
+    !densityoutrecept(i)=rho(iix,jjy,1,2)
+    densityoutrecept(i)=rho(iix,jjy,1,mind)
+  end do
 
 
Index: src/concoutput_nest_mpi.f90
===================================================================
--- src/concoutput_nest_mpi.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/concoutput_nest_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -99,4 +99,6 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  integer :: mind 
+! mind        eso:added to ensure identical results between 2&3-fields versions
 
 ! Measure execution time
@@ -149,4 +151,5 @@
   !*******************************************************************
 
+  mind=memind(2)
   do kz=1,numzgrid
     if (kz.eq.1) then
@@ -171,6 +174,8 @@
         iix=max(min(nint(xl),nxmin1),0)
         jjy=max(min(nint(yl),nymin1),0)
-        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
-             rho(iix,jjy,kzz-1,2)*dz2)/dz
+        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
+        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
+        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
+             rho(iix,jjy,kzz-1,mind)*dz2)/dz
       end do
     end do
@@ -182,5 +187,6 @@
     iix=max(min(nint(xl),nxmin1),0)
     jjy=max(min(nint(yl),nymin1),0)
-    densityoutrecept(i)=rho(iix,jjy,1,2)
+    !densityoutrecept(i)=rho(iix,jjy,1,2)
+    densityoutrecept(i)=rho(iix,jjy,1,mind)
   end do
 
Index: src/concoutput_surf_mpi.f90
===================================================================
--- src/concoutput_surf_mpi.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/concoutput_surf_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -105,16 +105,10 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  integer :: mind
+! mind        eso:added to get consistent results between 2&3-fields versions
 
 ! Measure execution time
-  if (mp_measure_time) then
-    call cpu_time(mp_root_time_beg)
-    mp_root_wtime_beg = mpi_wtime()
-  end if
-
-  if (verbosity.eq.1) then
-     print*,'inside concoutput_surf '
-     CALL SYSTEM_CLOCK(count_clock)
-     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-  endif
+  if (mp_measure_time) call mpif_mtime('rootonly',0)
+
 
   ! Determine current calendar date, needed for the file name
@@ -165,4 +159,5 @@
   !*******************************************************************
 
+  mind=memind(2)
   do kz=1,numzgrid
     if (kz.eq.1) then
@@ -187,17 +182,20 @@
         iix=max(min(nint(xl),nxmin1),0)
         jjy=max(min(nint(yl),nymin1),0)
-        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
-             rho(iix,jjy,kzz-1,2)*dz2)/dz
+        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
+        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
+        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
+             rho(iix,jjy,kzz-1,mind)*dz2)/dz
       end do
     end do
   end do
 
-    do i=1,numreceptor
-      xl=xreceptor(i)
-      yl=yreceptor(i)
-      iix=max(min(nint(xl),nxmin1),0)
-      jjy=max(min(nint(yl),nymin1),0)
-      densityoutrecept(i)=rho(iix,jjy,1,2)
-    end do
+  do i=1,numreceptor
+    xl=xreceptor(i)
+    yl=yreceptor(i)
+    iix=max(min(nint(xl),nxmin1),0)
+    jjy=max(min(nint(yl),nymin1),0)
+    !densityoutrecept(i)=rho(iix,jjy,1,2)
+    densityoutrecept(i)=rho(iix,jjy,1,mind)
+  end do
 
 
@@ -655,10 +653,5 @@
   end do
 
-  if (mp_measure_time) then
-    call cpu_time(mp_root_time_end)
-    mp_root_wtime_end = mpi_wtime()
-    mp_root_time_total = mp_root_time_total + (mp_root_time_end - mp_root_time_beg)
-    mp_root_wtime_total = mp_root_wtime_total + (mp_root_wtime_end - mp_root_wtime_beg)
-  end if
-
+  if (mp_measure_time) call mpif_mtime('rootonly',1)
+  
 end subroutine concoutput_surf
Index: src/concoutput_surf_nest_mpi.f90
===================================================================
--- src/concoutput_surf_nest_mpi.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/concoutput_surf_nest_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -98,4 +98,6 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  integer :: mind 
+! mind        eso:added to get consistent results between 2&3-fields versions
 
 ! Measure execution time
@@ -148,4 +150,5 @@
   !*******************************************************************
 
+  mind=memind(2)
   do kz=1,numzgrid
     if (kz.eq.1) then
@@ -170,17 +173,20 @@
         iix=max(min(nint(xl),nxmin1),0)
         jjy=max(min(nint(yl),nymin1),0)
-        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
-             rho(iix,jjy,kzz-1,2)*dz2)/dz
+        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
+        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
+        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
+             rho(iix,jjy,kzz-1,mind)*dz2)/dz
       end do
     end do
   end do
 
-    do i=1,numreceptor
-      xl=xreceptor(i)
-      yl=yreceptor(i)
-      iix=max(min(nint(xl),nxmin1),0)
-      jjy=max(min(nint(yl),nymin1),0)
-      densityoutrecept(i)=rho(iix,jjy,1,2)
-    end do
+  do i=1,numreceptor
+    xl=xreceptor(i)
+    yl=yreceptor(i)
+    iix=max(min(nint(xl),nxmin1),0)
+    jjy=max(min(nint(yl),nymin1),0)
+    !densityoutrecept(i)=rho(iix,jjy,1,2)
+    densityoutrecept(i)=rho(iix,jjy,1,mind)
+  end do
 
 
Index: src/getfields_mpi.f90
===================================================================
--- src/getfields_mpi.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/getfields_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -149,5 +149,5 @@
 ! 3 fields in memory
 ! Note that the reader process never needs to keep 3 fields in memory,
-! so it is possible to save some memory here
+! (2 is enough) so it is possible to save some memory
 !*********************************************************************
       if (mind3.eq.32.or.mind3.eq.19) then
@@ -245,5 +245,5 @@
         end if
         memtime(2)=wftime(indj+1)
-! :DEV: not used
+! DEV: not used
         if (.not.lmp_sync) memind(3)=3 ! indicate position of next read
         nstop = 1
Index: src/initial_cond_calc.f90
===================================================================
--- src/initial_cond_calc.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/initial_cond_calc.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -44,4 +44,6 @@
   real :: ddx,ddy
   real :: rhoprof(2),rhoi,xl,yl,wx,wy,w
+  integer :: mind2
+  ! mind2        eso: pointer to 2nd windfield in memory
 
 
@@ -88,9 +90,11 @@
   ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
   !*****************************************************************************
+    mind2=memind(2)
+
     do ind=indz,indzp
-      rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
-           +p2*rho(ixp,jy ,ind,2) &
-           +p3*rho(ix ,jyp,ind,2) &
-           +p4*rho(ixp,jyp,ind,2)
+      rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,mind2) &
+           +p2*rho(ixp,jy ,ind,mind2) &
+           +p3*rho(ix ,jyp,ind,mind2) &
+           +p4*rho(ixp,jyp,ind,mind2)
     end do
     rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
Index: src/makefile
===================================================================
--- src/makefile	(revision cd8513828ea365795e9076588d88b13779b34083)
+++ src/makefile	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -62,5 +62,5 @@
 
 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) -mtune=native $(FUSER) # -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 -fdump-core -Warray-bounds -ffpe-trap=invalid,overflow,denormal -Wall -fcheck=all $(FUSER)  # ,underflow,zero
+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 -Warray-bounds  -Wall -fcheck=all $(FUSER)  # -ffpe-trap=invalid,overflow,denormal,underflow,zero -fdump-core
 
 LDFLAGS  = $(FFLAGS) -L$(LIBPATH1) -L$(LIBPATH2) $(LIBS)
Index: src/mpi_mod.f90
===================================================================
--- src/mpi_mod.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/mpi_mod.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -44,4 +44,6 @@
 ! mp_pid                  Process ID of each MPI process                     *
 ! mp_seed                 Parameter for random number seed                   *
+! read_grp_min            Minimum number of processes at which one will be   *
+!                         used as reader                                     *
 ! numpart_mpi,            Number of particles per node                       *
 ! maxpart_mpi                                                                *
@@ -50,5 +52,5 @@
 !                         loops over particles. Will be all processes        *
 !                         unless a dedicated process runs getfields/readwind *
-!                                                                            *
+! lmp_sync                If .false., use asynchronous MPI                   *
 !                                                                            *
 !                                                                            *
@@ -69,5 +71,5 @@
 ! Set aside a process for reading windfields if using at least these many processes
 !==================================================
-  integer, parameter, private :: read_grp_min=99
+  integer, parameter, private :: read_grp_min=4
 !==================================================
 
@@ -86,6 +88,6 @@
 
 ! MPI tags/requests for send/receive operation
-  integer :: tm1!=100
-  integer, parameter :: nvar_async=27
+  integer :: tm1
+  integer, parameter :: nvar_async=27 !29 :DBG:
   !integer, dimension(:), allocatable :: tags
   integer, dimension(:), allocatable :: reqs
@@ -106,5 +108,5 @@
 !===============================================================================
 
-! mp_dbg_mode       MPI related output for debugging etc.
+! mp_dbg_mode       Used for debugging MPI.
 ! mp_dev_mode       various checks related to debugging the parallel code
 ! mp_dbg_out        write some arrays to data file for debugging
@@ -197,5 +199,15 @@
            & 'numwfmem should be set to 2 for syncronous reading'
       write(*,FMT='(80("#"))')
-    end if
+! Force "syncronized" version if all processes will call getfields
+    else if (.not.lmp_sync.and.mp_np.lt.read_grp_min) then
+      if (lroot) then
+        write(*,FMT='(80("#"))')
+        write(*,*) '#### mpi_mod::mpif_init> WARNING: ', &
+             & 'all procs call getfields. Setting lmp_sync=.true.'
+        write(*,FMT='(80("#"))')
+      end if
+      lmp_sync=.true. ! :DBG: eso fix this...
+    end if
+! TODO: Add warnings for unimplemented flexpart features
 
 ! Set ID of process that calls getfield/readwind. 
@@ -211,5 +223,5 @@
     call MPI_Comm_group (MPI_COMM_WORLD, world_group_id, mp_ierr)
 
-! Create a MPI group/communiactor that will calculate trajectories. 
+! Create a MPI group/communicator that will calculate trajectories. 
 ! Skip this step if program is run with only a few processes
 !************************************************************************
@@ -228,5 +240,5 @@
       lmp_use_reader = .true.
 
-! Map the subgroup IDs= 0,1,2,...,mp_np-2, skipping 'readwind' process
+! Map the subgroup IDs= 0,1,2,...,mp_np-2, skipping reader process
       j=-1 
       do i=0, mp_np-2 !loop over all (subgroup) IDs
@@ -289,12 +301,17 @@
     end if
     if (mp_dev_mode) write(*,*) 'PID, mp_seed: ',mp_pid, mp_seed
+    if (mp_dbg_mode) then
+! :DBG: For debugging, set all seed to 0 and maxrand to e.g. 20
+      mp_seed=0
+      if (lroot) write(*,*) 'MPI: setting seed=0'
+    end if
 
 ! Allocate request array for asynchronous MPI
     if (.not.lmp_sync) then
       allocate(reqs(0:nvar_async*mp_np-1))
-      reqs=MPI_REQUEST_NULL
+      reqs(:)=MPI_REQUEST_NULL
     else
       allocate(reqs(0:1))
-      reqs=MPI_REQUEST_NULL
+      reqs(:)=MPI_REQUEST_NULL
     end if
 
@@ -389,5 +406,5 @@
     integer :: add_part=0
 
-    call MPI_BCAST(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
+    call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
 
 ! MPI subgroup does the particle-loop 
@@ -815,5 +832,5 @@
 !*******************************************************************************
 ! DESCRIPTION
-!   Distribute 'getfield' variables from reader process to all processes.
+!   Distribute 'getfield' variables from reader process
 ! 
 !   Called from timemanager
@@ -887,4 +904,9 @@
 !**********************************************************************
 
+! The non-reader processes need to know if clouds were read.
+! TODO: only at first step or always?
+    call MPI_Bcast(readclouds,1,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+
 ! Static fields/variables sent only at startup
     if (first_call) then
@@ -904,5 +926,5 @@
       if (mp_ierr /= 0) goto 600 
       call MPI_Bcast(height,nzmax,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
-    if (mp_ierr /= 0) goto 600
+      if (mp_ierr /= 0) goto 600
 
       first_call=.false.
@@ -934,5 +956,13 @@
     if (mp_ierr /= 0) goto 600 
     call MPI_Bcast(clouds(:,:,:,li:ui),d3s1,MPI_INTEGER1,id_read,MPI_COMM_WORLD,mp_ierr)
-    if (mp_ierr /= 0) goto 600 
+    if (mp_ierr /= 0) goto 600
+
+! cloud water/ice:
+    if (readclouds) then
+      call MPI_Bcast(clwc(:,:,:,li:ui),d3s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
+      if (mp_ierr /= 0) goto 600
+      call MPI_Bcast(ciwc(:,:,:,li:ui),d3s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
+      if (mp_ierr /= 0) goto 600
+    end if
 
 ! 2D fields
@@ -969,5 +999,4 @@
     if (mp_ierr /= 0) goto 600 
 
-
     if (mp_measure_time) call mpif_mtime('commtime',1)
 
@@ -995,4 +1024,6 @@
 !   step
 !
+! TODO
+!   Transfer cloud water/ice if and when available for nested
 !
 !***********************************************************************
@@ -1140,6 +1171,5 @@
 ! DESCRIPTION
 !   Distribute 'getfield' variables from reader process to all processes.
-! 
-!   Called from timemanager by root process only
+!   Called from timemanager by reader process only
 !
 ! NOTE
@@ -1154,6 +1184,7 @@
 !   memstat -- input, for resolving pointer to windfield index being read
 !   mind    -- index where to place new fields
-!   !!! Under development, don't use yet !!!
-!
+!
+! TODO
+!   Transfer cloud water/ice
 !
 !*******************************************************************************
@@ -1173,13 +1204,5 @@
     integer :: d1s1 = maxwf
 
-    !integer :: d3s1,d3s2,d2s1,d2s2
-
 !*******************************************************************************
-
-! :TODO: don't need these
-    ! d3s1=d3_size1
-    ! d3s2=d3_size2
-    ! d2s1=d2_size1 
-    ! d2s2=d2_size2
 
 ! At the time the send is posted, the reader process is one step further
@@ -1220,5 +1243,5 @@
 
     do dest=0,mp_np-1 ! mp_np-2 will also work if last proc reserved for reading
-                      ! :TODO: use mp_partgroup_np here
+                      ! TODO: use mp_partgroup_np here
       if (dest.eq.id_read) cycle
       i=dest*nvar_async
@@ -1305,4 +1328,19 @@
       if (mp_ierr /= 0) goto 600 
 
+! Send cloud water if it exists. Increment counter always (as on receiving end)
+      if (readclouds) then
+        i=i+1
+        call MPI_Isend(clwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
+             &MPI_COMM_WORLD,reqs(i),mp_ierr)
+        if (mp_ierr /= 0) goto 600
+        i=i+1
+
+        call MPI_Isend(ciwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
+             &MPI_COMM_WORLD,reqs(i),mp_ierr)
+        if (mp_ierr /= 0) goto 600
+      ! else
+      !   i=i+2
+      end if
+
     end do
 
@@ -1311,5 +1349,5 @@
     goto 601
 
-600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
+600 write(*,*) "#### mpi_mod::mpif_gf_send_vars_async> mp_ierr \= 0", mp_ierr
     stop
 
@@ -1320,7 +1358,6 @@
 !*******************************************************************************
 ! DESCRIPTION
-!   Receive 'getfield' variables from reader process to all processes.
-! 
-!   Called from timemanager by all processes except root 
+!   Receive 'getfield' variables from reader process.
+!   Called from timemanager by all processes except reader 
 !
 ! NOTE
@@ -1330,4 +1367,6 @@
 !   memstat -- input, used to resolve windfield index being received
 !
+! TODO
+!   Transfer cloud water/ice
 !
 !*******************************************************************************
@@ -1388,87 +1427,127 @@
 ! Get MPI tags/requests for communications
     j=mp_pid*nvar_async
-    call MPI_Irecv(uu(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(vv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(uupol(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(vvpol(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(ww(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
+    call MPI_Irecv(uu(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(vv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(uupol(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(vvpol(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(ww(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     if (mp_ierr /= 0) goto 600
     j=j+1
-    call MPI_Irecv(tt(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(rho(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(drhodz(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(tth(:,:,:,mind),d3s2,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(qvh(:,:,:,mind),d3s2,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-
-    call MPI_Irecv(qv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
+    call MPI_Irecv(tt(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(rho(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(drhodz(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(tth(:,:,:,mind),d3s2,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(qvh(:,:,:,mind),d3s2,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+
+    call MPI_Irecv(qv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     if (mp_ierr /= 0) goto 600
-j=j+1
-    call MPI_Irecv(pv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(clouds(:,:,:,mind),d3s1,MPI_INTEGER1,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)    
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-
-    call MPI_Irecv(cloudsh(:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(vdep(:,:,:,mind),d2s2,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-
-    call MPI_Irecv(ps(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
+    j=j+1
+    call MPI_Irecv(pv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(clouds(:,:,:,mind),d3s1,MPI_INTEGER1,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)    
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+
+    call MPI_Irecv(cloudsh(:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(vdep(:,:,:,mind),d2s2,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(ps(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     if (mp_ierr /= 0) goto 600
     j=j+1
-    call MPI_Irecv(sd(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
+    call MPI_Irecv(sd(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     if (mp_ierr /= 0) goto 600
     j=j+1
-    call MPI_Irecv(tcc(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
+    call MPI_Irecv(tcc(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     if (mp_ierr /= 0) goto 600
     j=j+1
-    call MPI_Irecv(tt2(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(td2(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(lsprec(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(convprec(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-
-    call MPI_Irecv(ustar(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(wstar(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(hmix(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-    j=j+1
-    call MPI_Irecv(tropopause(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
+    call MPI_Irecv(tt2(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(td2(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(lsprec(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(convprec(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+
+    call MPI_Irecv(ustar(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(wstar(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(hmix(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+    j=j+1
+    call MPI_Irecv(tropopause(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     if (mp_ierr /= 0) goto 600
     j=j+1
-    call MPI_Irecv(oli(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
-    if (mp_ierr /= 0) goto 600 
+    call MPI_Irecv(oli(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
+         &MPI_COMM_WORLD,reqs(j),mp_ierr)
+    if (mp_ierr /= 0) goto 600 
+
+
+! Post request for clwc. These data are possibly not sent, request must then be cancelled
+! For now assume that data at all steps either have or do not have water 
+    if (readclouds) then
+      j=j+1
+      call MPI_Irecv(clwc(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+           &MPI_COMM_WORLD,reqs(j),mp_ierr)    
+      if (mp_ierr /= 0) goto 600 
+      j=j+1
+      call MPI_Irecv(ciwc(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
+           &MPI_COMM_WORLD,reqs(j),mp_ierr)    
+      if (mp_ierr /= 0) goto 600 
+    end if
 
 
@@ -1477,5 +1556,5 @@
     goto 601
 
-600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
+600 write(*,*) "#### mpi_mod::mpif_gf_recv_vars_async> MPI ERROR ", mp_ierr
     stop
 
@@ -1494,17 +1573,41 @@
 !
 !*******************************************************************************
+    use com_mod, only: readclouds
+
     implicit none
 
-!    if (mp_measure_time) call mpif_mtime('commtime',0)
+
+    integer :: n_req
+    integer :: i
+
+!***********************************************************************
+
+    n_req=nvar_async*mp_np
+
+    if (mp_measure_time) call mpif_mtime('commtime',0)
 
 !    call MPI_Wait(rm1,MPI_STATUS_IGNORE,mp_ierr)
-    call MPI_Waitall(nvar_async*mp_np,reqs,MPI_STATUSES_IGNORE,mp_ierr)
-    if (mp_ierr /= 0) goto 600 
-
-!    if (mp_measure_time) call mpif_mtime('commtime',1)
+
+! TODO: cancel recv request if readclouds=.false.
+!    if (readclouds) then
+    call MPI_Waitall(n_req,reqs,MPI_STATUSES_IGNORE,mp_ierr)
+!    endif
+    ! else
+    !   do i = 0, nvar_async*mp_np-1
+    !     if (mod(i,27).eq.0 .or. mod(i,28).eq.0) then
+    !       call MPI_Cancel(reqs(i),mp_ierr)
+    !       cycle
+    !     end if
+    !     call MPI_Wait(reqs(i),MPI_STATUS_IGNORE,mp_ierr)
+    !   end do
+    ! end if
+
+    if (mp_ierr /= 0) goto 600 
+
+    if (mp_measure_time) call mpif_mtime('commtime',1)
 
     goto 601
 
-600 write(*,*) "#### mpi_mod::mpif_gf_request> ERROR, mp_ierr \= 0 ", mp_ierr
+600 write(*,*) "#### mpi_mod::mpif_gf_request> MPI ERROR ", mp_ierr
     stop
 
@@ -1606,4 +1709,5 @@
 
 !**********************************************************************
+
     grid_size3d=numxgridn*numygridn*numzgrid*maxspec* &
          & maxpointspec_act*nclassunc*maxageclass
@@ -1666,4 +1770,6 @@
     character(LEN=*), intent(in) :: ident
     integer, intent(in) :: imode
+
+!***********************************************************************
 
     select case(ident)
@@ -1800,4 +1906,6 @@
 
     integer :: ip,j,r
+
+!***********************************************************************
 
     if (mp_measure_time) then
@@ -1857,5 +1965,5 @@
 ! In the implementation with 3 fields, the processes may have posted
 ! MPI_Irecv requests that should be cancelled here
-!! :TODO:
+!! TODO:
     ! if (.not.lmp_sync) then
     !   r=mp_pid*nvar_async
@@ -1868,5 +1976,5 @@
     call MPI_FINALIZE(mp_ierr)
     if (mp_ierr /= 0) then
-      write(*,*) '#### mpif_finalize::MPI_FINALIZE> ERROR ####'
+      write(*,*) '#### mpif_finalize::MPI_FINALIZE> MPI ERROR ', mp_ierr, ' ####'
       stop
     end if
@@ -1887,4 +1995,6 @@
     integer, save :: free_lun=100
     logical :: exists, iopen
+
+!***********************************************************************
 
     loop1: do
@@ -1912,4 +2022,6 @@
     character(LEN=40) :: fn_1, fn_2
 
+!***********************************************************************
+
     write(c_ts, FMT='(I8.8,BZ)') tstep
     fn_1='-'//trim(adjustl(c_ts))//'-'//trim(ident)
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/par_mod.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -216,6 +216,6 @@
   ! ---------
   integer,parameter :: maxwf=50000, maxtable=1000, numclass=13, ni=11
-  !integer,parameter :: numwfmem=2 ! Serial version/MPI with 2 fields
-  integer,parameter :: numwfmem=3 ! MPI with 3 fields
+  integer,parameter :: numwfmem=2 ! Serial version/MPI with 2 fields
+  !integer,parameter :: numwfmem=3 ! MPI with 3 fields
 
   ! maxwf                   maximum number of wind fields to be used for simulation
Index: src/pathnames
===================================================================
--- src/pathnames	(revision 0aded104fb66fe1f03a999e11f39e51a2c3a0a6d)
+++ src/pathnames	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -2,5 +2,5 @@
 ../output/
 /
-/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global
+/xnilu_wrk/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global
 ============================================
 
Index: src/readspecies.f90
===================================================================
--- src/readspecies.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/readspecies.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -216,50 +216,66 @@
 ! Check scavenging parameters given in SPECIES file
   if (dquer(pos_spec).gt.0) then !is particle
-    write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter A &
-         &(Rain collection efficiency)  ', weta(pos_spec)
-    write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter B &
-         &(Snow collection efficiency)  ', wetb(pos_spec)
+    if (lroot) then
+      write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter A &
+           &(Rain collection efficiency)  ', weta(pos_spec)
+      write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter B &
+           &(Snow collection efficiency)  ', wetb(pos_spec)
+    end if
     if (weta(pos_spec).gt.1.0 .or. wetb(pos_spec).gt.1.0) then
-      write(*,*) '*******************************************'
-      write(*,*) ' WARNING: Particle below-cloud scavenging parameter A or B &
-           &is out of likely range'
-      write(*,*) '          Likely range is 0.0-1.0'
-      write(*,*) '*******************************************'
+      if (lroot) then
+        write(*,*) '*******************************************'
+        write(*,*) ' WARNING: Particle below-cloud scavenging parameter A or B &
+             &is out of likely range'
+        write(*,*) '          Likely range is 0.0-1.0'
+        write(*,*) '*******************************************'
+      end if
     endif
-    write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Ai &
-         &(CCN efficiency)  ', weta_in(pos_spec)
-    write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Bi &
-         &(IN efficiency)  ', wetb_in(pos_spec)
-    if (weta_in(pos_spec).gt.1.0 .or. weta_in(pos_spec).lt.0.7 )then
-      write(*,*) '*******************************************'
-      write(*,*) ' WARNING: Particle in-cloud scavenging parameter A is out of likely range'
-      write(*,*) '          The default value is 0.9 for CCN '
-      write(*,*) '*******************************************'
+    if (lroot) then
+      write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Ai &
+           &(CCN efficiency)  ', weta_in(pos_spec)
+      write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Bi &
+           &(IN efficiency)  ', wetb_in(pos_spec)
+    end if
+    if (weta_in(pos_spec).gt.1.0 .or. weta_in(pos_spec).lt.0.7) then
+      if (lroot) then
+        write(*,*) '*******************************************'
+        write(*,*) ' WARNING: Particle in-cloud scavenging parameter A is out of likely range'
+        write(*,*) '          The default value is 0.9 for CCN '
+        write(*,*) '*******************************************'
+      end if
     endif
     if (wetb_in(pos_spec).gt.0.2 .or. wetb_in(pos_spec).lt.0.01) then
-      write(*,*) '*******************************************'
-      write(*,*) ' WARNING: Particle in-cloud scavenging parameter B is out of likely range'
-      write(*,*) '          The default value is 0.1 for IN '
-      write(*,*) '*******************************************'
+      if (lroot) then
+        write(*,*) '*******************************************'
+        write(*,*) ' WARNING: Particle in-cloud scavenging parameter B is out of likely range'
+        write(*,*) '          The default value is 0.1 for IN '
+        write(*,*) '*******************************************'
+      end if
     endif
 
   else !is gas
-    write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter A  ', weta(pos_spec)
-    write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter B  ', wetb(pos_spec)
-    write(*,*) ' Gas in-cloud scavenging uses default values as in Hertel et al 1995'
+    if (lroot) then
+      write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter A  ', weta(pos_spec)
+      write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter B  ', wetb(pos_spec)
+      write(*,*) ' Gas in-cloud scavenging uses default values as in Hertel et al 1995'
+    end if
+    if (weta(pos_spec).gt.0.) then !if wet deposition is turned on
+      if (weta(pos_spec).gt.1E-04 .or. weta(pos_spec).lt.1E-09) then
+        if (lroot) then
+          write(*,*) '*******************************************'
+          write(*,*) ' WARNING: Gas below-cloud scavenging parameter A is out of likely range'
+          write(*,*) '          Likely range is 1E-04 to 1E-08 (see Hertel et al 1995)'
+          write(*,*) '*******************************************'
+        end if
+      endif
+    end if
     if (wetb(pos_spec).gt.0.) then !if wet deposition is turned on
       if (wetb(pos_spec).gt.0.8 .or. wetb(pos_spec).lt.0.6) then
-        write(*,*) '*******************************************'
-        write(*,*) ' WARNING: Gas below-cloud scavenging parameter B is out of likely range'
-        write(*,*) '          Likely range is 0.6 to 0.8 (see Hertel et al 1995)'
-        write(*,*) '*******************************************'
-      endif
-    end if
-    if (weta(pos_spec).gt.0.) then !if wet deposition is turned on
-      if (weta(pos_spec).gt.1E-04 .or. weta(pos_spec).lt.1E-09) then
-        write(*,*) '*******************************************'
-        write(*,*) ' WARNING: Gas below-cloud scavenging parameter A is out of likely range'
-        write(*,*) '          Likely range is 1E-04 to 1E-08 (see Hertel et al 1995)'
-        write(*,*) '*******************************************'
+        if (lroot) then
+          write(*,*) '*******************************************'
+          write(*,*) ' WARNING: Gas below-cloud scavenging parameter B is out of likely range'
+          write(*,*) '          Likely range is 0.6 to 0.8 (see Hertel et al 1995)'
+          write(*,*) '*******************************************'
+        end if
       endif
     end if
Index: src/readwind.f90
===================================================================
--- src/readwind.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/readwind.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -77,7 +77,7 @@
 !HSO  end
 
-  real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+  real(sp) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real(sp) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real(sp) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
   integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
 
@@ -91,12 +91,12 @@
 
   integer :: isec1(56),isec2(22+nxmax+nymax)
-  real(kind=4) :: zsec4(jpunp)
-  real(kind=4) :: xaux,yaux,xaux0,yaux0
-  real(kind=8) :: xauxin,yauxin
-  real,parameter :: eps=1.e-4
-  real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
-  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
-
-  logical :: hflswitch,strswitch, readcloud
+  real(sp) :: zsec4(jpunp)
+  real(sp) :: xaux,yaux
+  real(dp) :: xauxin,yauxin
+  real(sp),parameter :: eps=1.e-4
+  real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
+  real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
+
+  logical :: hflswitch,strswitch!,readcloud
 
 !HSO  grib api error messages
@@ -107,5 +107,5 @@
   strswitch=.false.
 !hg
-  readcloud=.false.
+!  readcloud=.false.
 !hg end
   levdiff2=nlev_ec-nwz+1
@@ -362,5 +362,4 @@
         !write(*,*) 'found water!'
       endif
-
       if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
         ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
@@ -379,5 +378,4 @@
 
 50 call grib_close_file(ifile)
-
 
 !error message if no fields found with correct first longitude in it
Index: src/readwind_mpi.f90
===================================================================
--- src/readwind_mpi.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/readwind_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -41,5 +41,7 @@
 !  eso 2014:
 !    This version for reading windfields in advance by a separate
-!    MPI process. 
+!    MPI process.
+!    TODO: can be merged with readwind.f90 if using numwfmem in
+!          shift_field
 !**********************************************************************
 !                                                                     *
@@ -97,5 +99,5 @@
   integer :: isec1(56),isec2(22+nxmax+nymax)
   real(kind=4) :: zsec4(jpunp)
-  real(kind=4) :: xaux,yaux,xaux0,yaux0
+  real(kind=4) :: xaux,yaux
   real(kind=8) :: xauxin,yauxin
   real,parameter :: eps=1.e-4
@@ -103,5 +105,5 @@
   real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
 
-  logical :: hflswitch,strswitch, readcloud
+  logical :: hflswitch,strswitch !,readcloud
 
 !HSO  grib api error messages
@@ -115,5 +117,5 @@
   strswitch=.false.
 !hg
-  readcloud=.false.
+!  readcloud=.false.
 !hg end
   levdiff2=nlev_ec-nwz+1
@@ -234,7 +236,9 @@
     elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
       isec1(6)=176         ! indicatorOfParameter
-    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
+!    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS --wrong
+    elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS --correct
       isec1(6)=180         ! indicatorOfParameter
-    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
+!    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS --wrong
+    elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS --correct
       isec1(6)=181         ! indicatorOfParameter
     elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
Index: src/releaseparticles_mpi.f90
===================================================================
--- src/releaseparticles_mpi.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/releaseparticles_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -69,5 +69,6 @@
   real(kind=dp) :: juldate,julmonday,jul,jullocal,juldiff
   real,parameter :: eps=nxmax/3.e5,eps2=1.e-6
-  
+  integer :: mind2
+! mind2        eso: pointer to 2nd windfield in memory
 
   integer :: idummy = -7
@@ -291,21 +292,21 @@
               do kz=1,nz
                 if (ngrid.gt.0) then
-                  r=p1*rhon(ix ,jy ,kz,2,ngrid) &
-                       +p2*rhon(ixp,jy ,kz,2,ngrid) &
-                       +p3*rhon(ix ,jyp,kz,2,ngrid) &
-                       +p4*rhon(ixp,jyp,kz,2,ngrid)
-                  t=p1*ttn(ix ,jy ,kz,2,ngrid) &
-                       +p2*ttn(ixp,jy ,kz,2,ngrid) &
-                       +p3*ttn(ix ,jyp,kz,2,ngrid) &
-                       +p4*ttn(ixp,jyp,kz,2,ngrid)
+                  r=p1*rhon(ix ,jy ,kz,mind2,ngrid) &
+                       +p2*rhon(ixp,jy ,kz,mind2,ngrid) &
+                       +p3*rhon(ix ,jyp,kz,mind2,ngrid) &
+                       +p4*rhon(ixp,jyp,kz,mind2,ngrid)
+                  t=p1*ttn(ix ,jy ,kz,mind2,ngrid) &
+                       +p2*ttn(ixp,jy ,kz,mind2,ngrid) &
+                       +p3*ttn(ix ,jyp,kz,mind2,ngrid) &
+                       +p4*ttn(ixp,jyp,kz,mind2,ngrid)
                 else
-                  r=p1*rho(ix ,jy ,kz,2) &
-                       +p2*rho(ixp,jy ,kz,2) &
-                       +p3*rho(ix ,jyp,kz,2) &
-                       +p4*rho(ixp,jyp,kz,2)
-                  t=p1*tt(ix ,jy ,kz,2) &
-                       +p2*tt(ixp,jy ,kz,2) &
-                       +p3*tt(ix ,jyp,kz,2) &
-                       +p4*tt(ixp,jyp,kz,2)
+                  r=p1*rho(ix ,jy ,kz,mind2) &
+                       +p2*rho(ixp,jy ,kz,mind2) &
+                       +p3*rho(ix ,jyp,kz,mind2) &
+                       +p4*rho(ixp,jyp,kz,mind2)
+                  t=p1*tt(ix ,jy ,kz,mind2) &
+                       +p2*tt(ixp,jy ,kz,mind2) &
+                       +p3*tt(ix ,jyp,kz,mind2) &
+                       +p4*tt(ixp,jyp,kz,mind2)
                 endif
                 press=r*r_air*t/100.
@@ -377,15 +378,15 @@
               if (ngrid.gt.0) then
                 do n=1,2
-                  rhoaux(n)=p1*rhon(ix ,jy ,indz+n-1,2,ngrid) &
-                       +p2*rhon(ixp,jy ,indz+n-1,2,ngrid) &
-                       +p3*rhon(ix ,jyp,indz+n-1,2,ngrid) &
-                       +p4*rhon(ixp,jyp,indz+n-1,2,ngrid)
+                  rhoaux(n)=p1*rhon(ix ,jy ,indz+n-1,mind2,ngrid) &
+                       +p2*rhon(ixp,jy ,indz+n-1,mind2,ngrid) &
+                       +p3*rhon(ix ,jyp,indz+n-1,mind2,ngrid) &
+                       +p4*rhon(ixp,jyp,indz+n-1,mind2,ngrid)
                 end do
               else
                 do n=1,2
-                  rhoaux(n)=p1*rho(ix ,jy ,indz+n-1,2) &
-                       +p2*rho(ixp,jy ,indz+n-1,2) &
-                       +p3*rho(ix ,jyp,indz+n-1,2) &
-                       +p4*rho(ixp,jyp,indz+n-1,2)
+                  rhoaux(n)=p1*rho(ix ,jy ,indz+n-1,mind2) &
+                       +p2*rho(ixp,jy ,indz+n-1,mind2) &
+                       +p3*rho(ix ,jyp,indz+n-1,mind2) &
+                       +p4*rho(ixp,jyp,indz+n-1,mind2)
                 end do
               endif
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/timemanager_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -180,5 +180,5 @@
     if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) then
 ! readwind process skips this step
-      if (.not.(lmpreader.and.lmp_use_reader)) then !.or..not.lmp_use_reader)
+      if (.not.(lmpreader.and.lmp_use_reader)) then
         call ohreaction(itime,lsynctime,loutnext)
       endif
@@ -223,6 +223,6 @@
       call mpif_gf_send_vars(memstat)
       call mpif_gf_send_vars_nest(memstat)
-! This version is also used whenever 2 fields are needed (in this case,
-! async send/recv is impossible)
+! Version 2  (lmp_sync=.false.) is also used whenever 2 new fields are read (in which
+! case async send/recv is impossible.
     else if (.not.lmp_sync.and.lmp_use_reader.and.memstat.ge.32) then
       call mpif_gf_send_vars(memstat)
@@ -231,5 +231,5 @@
 
 ! Version 2 (lmp_sync=.false.) is for holding three fields in memory. Uses a
-! read-ahead process where sending/receiving of the 3rd fields are done in
+! read-ahead process where sending/receiving of the 3rd fields is done in
 ! the background in parallel with performing computations with fields 1&2
 !********************************************************************************
@@ -241,16 +241,16 @@
       end if
 
-
 ! COMPLETION CHECK:
 ! Issued at start of each new field period. 
       if (memstat.ne.0.and.memstat.lt.32.and.lmp_use_reader) then
-! :DEV: z0(7) changes with time, so should be dimension (numclass,2)
-! to allow transfer of the future value in the background
+! TODO: z0(7) changes with time, so should be dimension (numclass,2) to
+! allow transfer of the future value in the background
+        call MPI_Bcast(z0,numclass,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
         call mpif_gf_request
-        call MPI_Bcast(z0,numclass,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
       end if
 
-
 ! RECVEIVING PROCESS(ES):
+      ! eso TODO: at this point we do not know if clwc/ciwc will be available
+      ! at next time step. Issue receive request anyway, cancel at mpif_gf_request
       if (memstat.gt.0.and.lmp_use_reader.and..not.lmpreader) then
         call mpif_gf_recv_vars_async(memstat)
@@ -755,4 +755,5 @@
     if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl  
 
+! TODO: delete for release version?
 !!-------------------------------------------------------------------------------
 ! These lines below to test the well-mixed condition, modified by  mc, not to
@@ -778,5 +779,5 @@
 
 ! eso :TODO: this not implemented yet (transfer particles to PID 0 or rewrite)
-! the tools to do this is in mpi_mod.f90
+! the tools to do this are already in mpi_mod.f90
   if (lroot) then 
     do j=1,numpart
@@ -787,8 +788,8 @@
 
   if (ipout.eq.2) then
-! MPI: process 0 creates the file, the other processes append to it
+! MPI process 0 creates the file, the other processes append to it
     do ip=0, mp_partgroup_np-1
       if (ip.eq.mp_partid) then 
-        if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
+        !if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
         call partoutput(itime)    ! dump particle positions
       end if
Index: src/wetdepo.f90
===================================================================
--- src/wetdepo.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
+++ src/wetdepo.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
@@ -79,12 +79,12 @@
   real :: wetdeposit(maxspec),restmass
   real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
-  save lfr,cfr
-
-  real :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
-  real :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
+  !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 :: bclr(6) = (/274.35758, 332839.59273, 226656.57259, 58005.91340, 6588.38582, 0.244984/) !rain (Laakso et al 2003)
-  real :: bcls(6) = (/ 22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
+  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
 
@@ -159,5 +159,5 @@
   !********************************************************************
     interp_time=nint(itime-0.5*ltsample)
-
+    
     if (ngrid.eq.0) then
       call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
@@ -171,5 +171,5 @@
 
 !  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
+    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
 
   ! get the level were the actual particle is in
@@ -231,5 +231,5 @@
 
 
-!ZHG calculated for 1) both 2) lsp 3) convp 
+  !ZHG calculated for 1) both 2) lsp 3) convp 
     grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
     grfraction(2)=max(0.05,cc*(lfr(i)))
@@ -252,9 +252,9 @@
       wetscav=0.   
 
-     if (ngrid.gt.0) then
-       act_temp=ttn(ix,jy,hz,n,ngrid)
-     else
-       act_temp=tt(ix,jy,hz,n)
-     endif
+      if (ngrid.gt.0) then
+        act_temp=ttn(ix,jy,hz,n,ngrid)
+      else
+        act_temp=tt(ix,jy,hz,n)
+      endif
      
 
@@ -264,37 +264,37 @@
       if (clouds_v.ge.4) then !below cloud
 
-         if (weta(ks).gt.0. .or. wetb(ks).gt.0) then !if positive below-cloud parameters given in SPECIES file (either A or B)
-            blc_count=blc_count+1
+        if (weta(ks).gt.0. .or. wetb(ks).gt.0) then !if positive below-cloud parameters given in SPECIES file (either A or B)
+          blc_count=blc_count+1
 
 !GAS
-           if (dquer(ks) .le. 0.) then  !is gas
-            !Gas scavenging coefficient based on Hertel et al 1995 using the below-cloud scavenging parameters A (=weta) and B (=wetb) from SPECIES file
+          if (dquer(ks) .le. 0.) then  !is gas
+  ! Gas scavenging coefficient based on Hertel et al 1995 using the below-cloud scavenging parameters A (=weta) and B (=wetb) from SPECIES file
             wetscav=weta(ks)*prec(1)**wetb(ks) 
 
-!AEROSOL
-           else !is particle
-             !NIK 17.02.2015
-             ! For the calculation here particle size needs to be in meter and not um as dquer is changed to in readreleases
-             dquer_m=dquer(ks)/1000000. !conversion from um to m
-
-             !ZHG snow or rain removal is applied based on the temperature.
-             if (act_temp .ge. 273)  then !Rain 
-  
-             !Particle RAIN scavenging coefficient based on Laakso et al 2003, the below-cloud scavenging (rain efficienty) parameter A (=weta) from SPECIES file
-                 wetscav= weta(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))
-
-             elseif (act_temp .lt. 273)  then !snow
-  
-             !Particle SNOW scavenging coefficient based on Kyro et al 2009, the below-cloud scavenging (Snow efficiency) parameter B (=wetb) from SPECIES file
-                wetscav= wetb(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         
+  !AEROSOL
+          else !is particle
+  !NIK 17.02.2015
+  ! For the calculation here particle size needs to be in meter and not um as dquer is changed to in readreleases
+            dquer_m=dquer(ks)/1000000. !conversion from um to m
+           
+  !ZHG snow or rain removal is applied based on the temperature.
+            if (act_temp .ge. 273)  then !Rain 
+
+  !Particle RAIN scavenging coefficient based on Laakso et al 2003, the below-cloud scavenging (rain efficienty) parameter A (=weta) from SPECIES file
+              wetscav= weta(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))
+
+            elseif (act_temp .lt. 273)  then !snow
+
+  !Particle SNOW scavenging coefficient based on Kyro et al 2009, the below-cloud scavenging (Snow efficiency) parameter B (=wetb) from SPECIES file
+              wetscav= wetb(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 !gas or particle
+        endif ! positive below-cloud scavenging parameters given in Species file
       endif !end below-cloud
 
@@ -305,40 +305,40 @@
       if (clouds_v.lt.4) then !in-cloud
 
-          inc_count=inc_count+1
-
-          !ZHG liquid water parameterization (CLWC+CIWC) 
-          if (readclouds) then !get cloud water clwc & ciwc units Kg/Kg
-             cl=clwc(ix,jy,hz,n)+ciwc(ix,jy,hz,n)
-          else !parameterize cloudwater
-             cl=2E-7*prec(1)**0.36
+        inc_count=inc_count+1
+
+  !ZHG liquid water parameterization (CLWC+CIWC) 
+        if (readclouds) then !get cloud water clwc & ciwc units Kg/Kg
+          cl=clwc(ix,jy,hz,n)+ciwc(ix,jy,hz,n)
+        else !parameterize cloudwater
+          cl=2E-7*prec(1)**0.36
+        endif
+
+  ! AEROSOL
+        if (dquer(ks).gt. 0.) then ! is particle
+          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
           endif
 
-! AEROSOL
-          if (dquer(ks).gt. 0.) then ! is particle
-            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
-            endif
-
-            !ZHG  calculate the activated fraction based on the In-cloud scavenging parameters Ai (=weta_in) and Bi (=wetb_in) from SPECIES file
-            ! frac_act is the combined IN and CCN efficiency
-            ! The default values are 0.9 for CCN and 0.1 IN
-            ! This parameterization is based on Verheggen et al. (2007) & Cozich et al. (2006)
-            frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
+! ZHG  calculate the activated fraction based on the In-cloud scavenging parameters Ai (=weta_in) and Bi (=wetb_in) from SPECIES file
+! frac_act is the combined IN and CCN efficiency
+! The default values are 0.9 for CCN and 0.1 IN
+! This parameterization is based on Verheggen et al. (2007) & Cozich et al. (2006)
+          frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
  
-            !ZHG Use the activated fraction and the liqid water to calculate the washout
-            S_i= frac_act/cl
-
-! GAS
-          else ! is gas
-            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
-            S_i=1/cle
-          endif
-
-          ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
-          wetscav=S_i*(prec(1)/3.6E6)/clouds_h
+  !ZHG Use the activated fraction and the liqid water to calculate the washout
+          S_i= frac_act/cl
+
+  ! GAS
+        else ! is gas
+          cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+          S_i=1/cle
+        endif
+
+  ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
+        wetscav=S_i*(prec(1)/3.6E6)/clouds_h
 
 !          write(*,*) 'in-cloud, act_temp=',act_temp,',prec=',prec(1),',wetscav=',wetscav,',jpart=',jpart,',clouds_h=,', &
@@ -346,5 +346,5 @@
 
       endif !incloud
-
+     
   !**************************************************
   ! CALCULATE DEPOSITION 
@@ -355,5 +355,5 @@
       if (wetscav.gt.0.) then
         wetdeposit(ks)=xmass1(jpart,ks)* &
-        (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
+             (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
       else ! if no scavenging
         wetdeposit(ks)=0.
@@ -389,7 +389,7 @@
     if (ldirect.eq.1) then
       call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
-      real(ytra1(jpart)),nage,kp)
+           real(ytra1(jpart)),nage,kp)
       if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
-        wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
+           wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
     endif
 
@@ -397,7 +397,7 @@
   end do ! all particles
 
-! count the total number of below-cloud and in-cloud occurences:
-   tot_blc_count=tot_blc_count+blc_count
-   tot_inc_count=tot_inc_count+inc_count
+  ! count the total number of below-cloud and in-cloud occurences:
+  tot_blc_count=tot_blc_count+blc_count
+  tot_inc_count=tot_inc_count+inc_count
 
 end subroutine wetdepo
