Index: src/FLEXPART_MPI.f90
===================================================================
--- src/FLEXPART_MPI.f90	(revision df967a995b6f35cdba22c7f6992c9bc67442d6b7)
+++ src/FLEXPART_MPI.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -165,8 +165,9 @@
     write(*,FMT='(80("#"))')
     write(*,*) '#### FLEXPART_MPI> ERROR: ', &
-         & 'MPI version not (yet) working with backward runs. ',&
-         & 'Use the serial version instead.'
+         & 'MPI version not (yet) working with backward runs. '
+    write(*,*) '#### Use the serial version instead.'
     write(*,FMT='(80("#"))')
-    stop
+    ! call mpif_finalize
+    ! stop
   end if
 
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision 6b22af900a57ac3fc1a984db76f4ef2cb0436af0)
+++ src/com_mod.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -130,5 +130,5 @@
   ! gdomainfill             .T., if domain-filling is global, .F. if not
 
-!hg
+!ZHG SEP 2015 wheather or not to read clouds from GRIB
   logical :: readclouds
 
@@ -347,7 +347,9 @@
   real :: tt(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: qv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
-!hg adding cloud water 
-  real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
-  real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
+!ZHG adding cloud water 
+  real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem) !liquid   [kg/kg]
+  real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem) !ice      [kg/kg]
+  real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)  !combined [m3/m3]
+
   real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
@@ -362,6 +364,8 @@
   integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   integer :: cloudsh(0:nxmax-1,0:nymax-1,numwfmem)
-  ! integer :: icloudbot(0:nxmax-1,0:nymax-1,numwfmem)
-  ! integer :: icloudthck(0:nxmax-1,0:nymax-1,numwfmem)
+
+!ZHG Sep 2015  
+   real :: icloud_stats(0:nxmax-1,0:nymax-1,5,numwfmem)
+   real :: clh(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
 
 
@@ -405,6 +409,6 @@
   real :: tropopause(0:nxmax-1,0:nymax-1,1,numwfmem)
   real :: oli(0:nxmax-1,0:nymax-1,1,numwfmem)
-! real :: diffk(0:nxmax-1,0:nymax-1,1,numwfmem) this is not in use?
-
+! real :: diffk(0:nxmax-1,0:nymax-1,1,numwfmem) ESO: this is not in use?
+! logical :: beneath_cloud=.true.
   ! ps                   surface pressure
   ! sd                   snow depth
@@ -667,4 +671,8 @@
 
 
+  !CGZ-lifetime
+  real, allocatable, dimension(:,:) ::checklifetime, species_lifetime
+  !CGZ-lifetime
+
   ! numpart                 actual number of particles in memory
   ! itra1 (maxpart) [s]     temporal positions of the particles
@@ -759,5 +767,7 @@
          & idt(nmpart),itramem(nmpart),itrasplit(nmpart),&
          & xtra1(nmpart),ytra1(nmpart),ztra1(nmpart),&
-         & xmass1(nmpart, maxspec))
+         & xmass1(nmpart, maxspec),&
+         & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
+
 
     allocate(uap(nmpart),ucp(nmpart),uzp(nmpart),us(nmpart),&
Index: src/gridcheck.f90
===================================================================
--- src/gridcheck.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ src/gridcheck.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -195,10 +195,10 @@
   elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
     isec1(6)=133         ! indicatorOfParameter
-!hg
+!ZHG FOR CLOUDS FROM GRIB
   elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
     isec1(6)=246         ! indicatorOfParameter
   elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
     isec1(6)=247         ! indicatorOfParameter
-!hg end
+!ZHG end
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     isec1(6)=134         ! indicatorOfParameter
Index: src/mpi_mod.f90
===================================================================
--- src/mpi_mod.f90	(revision a1f4dd6e486a6db079a9512ecc10c867a4fc09d4)
+++ src/mpi_mod.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -89,5 +89,5 @@
 ! MPI tags/requests for send/receive operation
   integer :: tm1
-  integer, parameter :: nvar_async=27 !29 :DBG:
+  integer, parameter :: nvar_async=26 !27 !29 :DBG:
 !integer, dimension(:), allocatable :: tags
   integer, dimension(:), allocatable :: reqs
@@ -119,5 +119,5 @@
   logical, parameter :: mp_dbg_out = .false.
   logical, parameter :: mp_time_barrier=.true.
-  logical, parameter :: mp_measure_time=.true.
+  logical, parameter :: mp_measure_time=.false.
   logical, parameter :: mp_exact_numpart=.true.
 
@@ -207,5 +207,5 @@
       write(*,FMT='(80("#"))')
 ! Force "syncronized" version if all processes will call getfields
-    else if (.not.lmp_sync.and.mp_np.lt.read_grp_min) then
+    else if ((.not.lmp_sync.and.mp_np.lt.read_grp_min).or.(mp_np.eq.1)) then
       if (lroot) then
         write(*,FMT='(80("#"))')
@@ -963,8 +963,11 @@
 ! 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
+      call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2_size1*5,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
+      if (mp_ierr /= 0) goto 600
+
+      ! 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
 
@@ -1189,5 +1192,4 @@
 !
 ! TODO
-!   Transfer cloud water/ice
 !
 !*******************************************************************************
@@ -1334,16 +1336,18 @@
       if (readclouds) then
         i=i+1
-        call MPI_Isend(clwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
+        call MPI_Isend(icloud_stats(:,:,:,mind),d2s1*5,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
+
+        ! 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
+
       end if
-
     end do
 
@@ -1371,5 +1375,4 @@
 !
 ! TODO
-!   Transfer cloud water/ice
 !
 !*******************************************************************************
@@ -1390,10 +1393,4 @@
 !    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 this immediate receive is posted, memstat is the state of
@@ -1545,11 +1542,17 @@
     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 
+
+      call MPI_Irecv(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,id_read,MPI_ANY_TAG,&
+           &MPI_COMM_WORLD,reqs(j),mp_ierr)
+      if (mp_ierr /= 0) goto 600
+
+      ! 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
 
@@ -2089,5 +2092,7 @@
     implicit none
 
-    integer, parameter :: li=1, ui=3 ! wfmem indices (i.e, operate on entire field)
+    integer :: li=1, ui=2 ! wfmem indices (i.e, operate on entire field)
+
+    if (.not.lmp_sync) ui=3
 
 
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision 4c4261c79d2d3b541ca53a77d926381a82b1814b)
+++ src/par_mod.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -76,5 +76,5 @@
   real,parameter :: hmixmin=100., hmixmax=4500., turbmesoscale=0.16
   real,parameter :: d_trop=50., d_strat=0.1
-
+  real,parameter :: rho_water=1000 !ZHG 2015 [kg/m3]
   ! karman                  Karman's constant
   ! href [m]                Reference height for dry deposition
@@ -226,6 +226,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 de7afbd99ab0adaeab4e79d9be82610ba6277d59)
+++ src/pathnames	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -2,6 +2,8 @@
 ../output/
 /
-/xnilu_wrk/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global/
+/xnilu_wrk/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global
 ============================================
+
+/xnilu_wrk/flex_wrk/zhg/ECMWF_CLWC/Availables
 
 
Index: src/readwind.f90
===================================================================
--- src/readwind.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ src/readwind.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -106,7 +106,7 @@
   hflswitch=.false.
   strswitch=.false.
-!hg
+!ZHG test the grib fields that have lcwc without using them
 !  readcloud=.false.
-!hg end
+
   levdiff2=nlev_ec-nwz+1
   iumax=0
@@ -190,10 +190,11 @@
     elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
       isec1(6)=133         ! indicatorOfParameter
-!hg
+!ZHG ! READ CLOUD FIELD : assume these to be added together to one variable
     elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
       isec1(6)=246         ! indicatorOfParameter
-    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
-      isec1(6)=247         ! indicatorOfParameter
- !hg end
+! ICE AND WATER IS ADDED TOGETHER IN NEW WINDFIELDS
+!    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
+!      isec1(6)=247         ! indicatorOfParameter
+!ZHG end
     elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
       isec1(6)=134         ! indicatorOfParameter
@@ -226,7 +227,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
@@ -356,15 +359,15 @@
       if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
       if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
-!hg READING CLOUD FIELDS ASWELL
-      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content
+!ZHG READING CLOUD FIELDS ASWELL
+      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
         clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
         readclouds = .true.
-        !write(*,*) 'found water!'
+!if (clwch(i,j,nlev_ec-k+2,n) .gt. 0)        write(*,*) 'readwind: found water!', clwch(i,j,nlev_ec-k+2,n)
       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)
+!      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)
         !write(*,*) 'found ice!'
-      endif
-!hg end
+!      endif
+!ZHG end
 
     end do
@@ -422,8 +425,8 @@
     call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
     call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
-!hg
+!ZHG
     call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
     call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
-!hg end
+!ZHG end
 
   endif
Index: src/readwind_mpi.f90
===================================================================
--- src/readwind_mpi.f90	(revision 5f9d14aa526df192efa539f87fc409fdeeeee06c)
+++ src/readwind_mpi.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -200,10 +200,11 @@
     elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
       isec1(6)=133         ! indicatorOfParameter
-!hg
+!ZHG ! READ CLOUD FIELD : assume these to be added together to one variable
     elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
       isec1(6)=246         ! indicatorOfParameter
-    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
-      isec1(6)=247         ! indicatorOfParameter
- !hg end
+! ICE AND WATER IS ADDED TOGETHER IN NEW WINDFIELDS
+!    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
+!      isec1(6)=247         ! indicatorOfParameter
+!ZHG end
     elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
       isec1(6)=134         ! indicatorOfParameter
@@ -368,16 +369,15 @@
       if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
       if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
-!hg READING CLOUD FIELDS ASWELL
-      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content
+!ZHG READING CLOUD FIELDS ASWELL
+      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
         clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
         readclouds = .true.
-        !write(*,*) 'found water!'
+!if (clwch(i,j,nlev_ec-k+2,n) .gt. 0)        write(*,*) 'readwind: found water!', clwch(i,j,nlev_ec-k+2,n)
       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)
+!      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)
         !write(*,*) 'found ice!'
-      endif
-!hg end
+!      endif
+!ZHG end
 
     end do
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision 5f42c27a3ad632d81e1e6ecd2f81d821b0c11373)
+++ src/timemanager.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -136,4 +136,12 @@
   !**********************************************************************
 
+!ZHG 2015
+!CGZ-lifetime: set lifetime to 0
+  checklifetime(:,:)=0
+  species_lifetime(:,:)=0
+  print*, 'Initialized lifetime'
+!CGZ-lifetime: set lifetime to 0
+  
+
 
   if (verbosity.gt.0) then
@@ -411,4 +419,13 @@
         if (iflux.eq.1) call fluxoutput(itime)
         write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
+ 
+        !CGZ-lifetime: output species lifetime
+!ZHG
+        write(*,*) 'Overview species lifetime in days', &
+             real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
+        write(*,*) 'all info:',species_lifetime
+!ZHG
+        !CGZ-lifetime: output species lifetime
+
         !write(*,46) float(itime)/3600,itime,numpart
 45      format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES:    Uncertainty: ',3f7.3)
@@ -575,4 +592,14 @@
                    xmassfract=max(xmassfract,real(npart(npoint(j)))* &
                    xmass1(j,ks)/xmass(npoint(j),ks))
+!ZHG 2015
+                  !CGZ-lifetime: Check mass fraction left/save lifetime
+                   if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
+                       !Mass below 1% of initial >register lifetime
+                       checklifetime(j,ks)=abs(itra1(j)-itramem(j))
+                       species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
+                       species_lifetime(ks,2)= species_lifetime(ks,2)+1
+                   endif
+                   !CGZ-lifetime: Check mass fraction left/save lifetime
+!ZHG 2015
             else
               xmassfract=1.
Index: src/timemanager_mpi.f90
===================================================================
--- src/timemanager_mpi.f90	(revision a1f4dd6e486a6db079a9512ecc10c867a4fc09d4)
+++ src/timemanager_mpi.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -139,4 +139,23 @@
 !**********************************************************************
 
+
+!  itime=0
+  if (lroot) then
+  !  write(*,45) itime,numpart*mp_partgroup_np,gridtotalunc,wetgridtotalunc,drygridtotalunc
+    write(*,46) float(itime)/3600,itime,numpart*mp_partgroup_np
+    
+    if (verbosity.gt.0) then
+      write (*,*) 'timemanager> starting simulation'
+    end if
+  end if ! (lroot)
+
+!CGZ-lifetime: set lifetime to 0
+  checklifetime(:,:)=0
+  species_lifetime(:,:)=0
+  print*, 'Initialized lifetime'
+!CGZ-lifetime: set lifetime to 0
+
+
+
   do itime=0,ideltas,lsynctime
 
@@ -544,11 +563,15 @@
         endif
         
-        if (lroot) then
+        !CGZ-lifetime: output species lifetime
+        ! if (lroot) then
+        !   write(*,*) 'Overview species lifetime in days', &
+        !        real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
+        !   write(*,*) 'all info:',species_lifetime
           write(*,45) itime,numpart_tot_mpi,gridtotalunc,&
                &wetgridtotalunc,drygridtotalunc
-          if (verbosity.gt.0) then
-            write (*,*) 'timemanager> starting simulation'
-          end if
-        end if
+        !   if (verbosity.gt.0) then
+        !     write (*,*) 'timemanager> starting simulation'
+        !   end if
+        ! end if
 
 45      format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES:    Uncertainty: ',3f7.3)
@@ -729,7 +752,19 @@
 
             if (mdomainfill.eq.0) then
-              if (xmass(npoint(j),ks).gt.0.) &
+              if (xmass(npoint(j),ks).gt.0.)then 
                    xmassfract=max(xmassfract,real(npart(npoint(j)))* &
                    xmass1(j,ks)/xmass(npoint(j),ks))
+                   
+                   !CGZ-lifetime: Check mass fraction left/save lifetime
+                   if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.0.01.and.checklifetime(j,ks).eq.0.)then
+                       !Mass below 1% of initial >register lifetime
+                       checklifetime(j,ks)=abs(itra1(j)-itramem(j))
+
+                       species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
+                       species_lifetime(ks,2)= species_lifetime(ks,2)+1
+                   endif
+                   !CGZ-lifetime: Check mass fraction left/save lifetime
+                   
+              endif
             else
               xmassfract=1.
@@ -737,5 +772,7 @@
           end do
 
-          if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
+          if (xmassfract.lt.0.00005 .and. sum(real(npart(npoint(j)))*xmass1(j,:)).lt.1.0) then   ! terminate all particles carrying less mass
+          !            print*,'terminated particle ',j,' for small mass (', sum(real(npart(npoint(j)))* &
+          !         xmass1(j,:)), ' of ', sum(xmass(npoint(j),:)),')'
             itra1(j)=-999999999
           endif
@@ -758,4 +795,5 @@
                  call initial_cond_calc(itime+lsynctime,j)
             itra1(j)=-999999999
+            !print*, 'terminated particle ',j,'for age'
           endif
         endif
@@ -828,5 +866,5 @@
   endif
   deallocate(gridunc)
-  deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
+  deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass, checklifetime)
   deallocate(ireleasestart,ireleaseend,npart,kindz)
   deallocate(xmasssave)
Index: src/verttransform.f90
===================================================================
--- src/verttransform.f90	(revision 4d456391fa86b9dc93fe32bac700d3c926f265cb)
+++ src/verttransform.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -97,11 +97,17 @@
   logical :: init = .true.
 
-!hg 
+  !ZHG SEP 2014 tests  
   integer :: cloud_ver,cloud_min, cloud_max 
+  integer ::teller(5), convpteller=0, lspteller=0
   real :: cloud_col_wat, cloud_water
-!hg temporary variables for testing
+  !ZHG 2015 temporary variables for testing
   real :: rcw(0:nxmax-1,0:nymax-1)
   real :: rpc(0:nxmax-1,0:nymax-1)
-!hg 
+  character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/'
+  character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
+  CHARACTER(LEN=3)  :: aspec
+  integer :: virr=0
+  real :: tot_cloud_h
+!ZHG
 
 !*************************************************************************
@@ -561,10 +567,73 @@
 
 
+!***********************************************************************************  
+if (readclouds) then !HG METHOD
+! The method is loops all grids vertically and constructs the 3D matrix for clouds
+! Cloud top and cloud bottom gid cells are assigned as well as the total column 
+! cloud water. For precipitating grids, the type and whether it is in or below 
+! cloud scavenging are assigned with numbers 2-5 (following the old metod).
+! Distinction is done for lsp and convp though they are treated the same in regards
+! to scavenging. Also clouds that are not precipitating are defined which may be 
+! to include future cloud processing by non-precipitating-clouds. 
+!***********************************************************************************
   if (readclouds) write(*,*) 'using cloud water from ECMWF'
-!  if (.not.readclouds) write(*,*) 'using cloud water from parameterization'
-
-  rcw(:,:)=0
-  rpc(:,:)=0
-
+  clw(:,:,:,n)=0
+  icloud_stats(:,:,:,n)=0
+  clouds(:,:,:,n)=0
+  do jy=0,nymin1
+   do ix=0,nxmin1
+    lsp=lsprec(ix,jy,1,n)
+    convp=convprec(ix,jy,1,n)
+    prec=lsp+convp
+    tot_cloud_h=0
+    ! Find clouds in the vertical
+    do kz=1, nz-1 !go from top to bottom
+     if (clwc(ix,jy,kz,n).gt.0) then      
+      ! 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)) 
+      icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
+      icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
+!ZHG 2015 extra for testing
+      clh(ix,jy,kz,n)=height(kz+1)-height(kz)
+      icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
+      icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
+!ZHG
+    endif
+    end do
+
+   ! If Precipitation. Define removal type in the vertical
+  if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+
+    do kz=nz,1,-1 !go Bottom up!
+     if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
+        cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) 
+        clouds(ix,jy,kz,n)=1                               ! is a cloud
+        if (lsp.ge.convp) then
+           clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
+        else
+          clouds(ix,jy,kz,n)=2                             ! convp in-cloud
+        endif                                              ! convective or large scale
+     elseif((clw(ix,jy,kz,n).le.0) .and. (icloud_stats(ix,jy,3,n).ge.height(kz)) ) then ! is below cloud
+        if (lsp.ge.convp) then
+          clouds(ix,jy,kz,n)=5                             ! lsp dominated washout
+        else
+          clouds(ix,jy,kz,n)=4                             ! convp dominated washout
+        endif                                              ! convective or large scale 
+     endif
+
+     if (height(kz).ge. 19000) then                        ! set a max height for removal
+        clouds(ix,jy,kz,n)=0
+     endif !clw>0
+    end do !nz
+   endif ! precipitation
+  end do
+ end do
+!**************************************************************************
+ else       ! use old definitions
+!**************************************************************************
+!   create a cloud and rainout/washout field, clouds occur where rh>80%
+!   total cloudheight is stored at level 0
+ if (.not.readclouds) write(*,*) 'using cloud water from Parameterization'
   do jy=0,nymin1
     do ix=0,nxmin1
@@ -603,4 +672,89 @@
       end do
 !END OLD METHOD
+      end do
+     end do
+endif !readclouds
+
+     !********* TEST ***************
+     ! WRITE OUT SOME TEST VARIABLES
+     !********* TEST ************'**
+!teller(:)=0
+!virr=virr+1
+!WRITE(aspec, '(i3.3)'), virr
+
+!if (readclouds) then
+!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement.txt'
+!else
+!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement_old.txt'
+!endif
+!
+!OPEN(UNIT=118, FILE=fnameH,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!do kz_inv=1,nz-1
+!  kz=nz-kz_inv+1
+!  !kz=91
+!  do jy=0,nymin1
+!     do ix=0,nxmin1
+!          if (clouds(ix,jy,kz,n).eq.1) teller(1)=teller(1)+1 ! no precipitation cloud
+!          if (clouds(ix,jy,kz,n).eq.2) teller(2)=teller(2)+1 ! convp dominated rainout
+!          if (clouds(ix,jy,kz,n).eq.3) teller(3)=teller(3)+1 ! lsp dominated rainout
+!          if (clouds(ix,jy,kz,n).eq.4) teller(4)=teller(4)+1 ! convp dominated washout
+!          if (clouds(ix,jy,kz,n).eq.5) teller(5)=teller(5)+1 ! lsp dominated washout
+!          
+!        !  write(*,*) height(kz),teller
+!     end do
+!  end do
+!  write(118,*) height(kz),teller
+!  teller(:)=0
+!end do
+!teller(:)=0
+!write(*,*) teller 
+!write(*,*) aspec
+!
+!fnameA=trim(zhgpath)//trim(aspec)//'cloudV.txt'
+!fnameB=trim(zhgpath)//trim(aspec)//'cloudT.txt'
+!fnameC=trim(zhgpath)//trim(aspec)//'cloudB.txt'
+!fnameD=trim(zhgpath)//trim(aspec)//'cloudW.txt'
+!fnameE=trim(zhgpath)//trim(aspec)//'old_cloudV.txt'
+!fnameF=trim(zhgpath)//trim(aspec)//'lsp.txt'
+!fnameG=trim(zhgpath)//trim(aspec)//'convp.txt'
+!if (readclouds) then
+!OPEN(UNIT=111, FILE=fnameA,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=112, FILE=fnameB,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=113, FILE=fnameC,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=114, FILE=fnameD,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!else
+!OPEN(UNIT=115, FILE=fnameE,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=116, FILE=fnameF,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=117, FILE=fnameG,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!endif
+!
+!do ix=0,nxmin1
+!if (readclouds) then
+!write(111,*) (icloud_stats(ix,jy,1,n),jy=0,nymin1)
+!write(112,*) (icloud_stats(ix,jy,2,n),jy=0,nymin1)
+!write(113,*) (icloud_stats(ix,jy,3,n),jy=0,nymin1)
+!write(114,*) (icloud_stats(ix,jy,4,n),jy=0,nymin1)
+!else
+!write(115,*) (cloudsh(ix,jy,n),jy=0,nymin1)    !integer
+!write(116,*) (lsprec(ix,jy,1,n),jy=0,nymin1)   !7.83691406E-02 
+!write(117,*) (convprec(ix,jy,1,n),jy=0,nymin1) !5.38330078E-02
+!endif
+!end do
+!
+!if (readclouds) then
+!CLOSE(111)
+!CLOSE(112)
+!CLOSE(113)
+!CLOSE(114)
+!else
+!CLOSE(115)
+!CLOSE(116)
+!CLOSE(117)
+!endif
+!
+!END ********* TEST *************** END
+! WRITE OUT SOME TEST VARIABLES
+!END ********* TEST *************** END
+
 
 ! PS 2012
@@ -613,32 +767,4 @@
 !            lconvectprec = .true.
 !           endif
-!HG METHOD
-!readclouds =.true.
-!      if (readclouds) then 
-!hg added APR 2014  Cloud Water=clwc(ix,jy,kz,n)  Cloud Ice=ciwc(ix,jy,kz,n)
-!hg Use the cloud water variables to determine existence of clouds. This makes the PS code obsolete 
-!        cloud_min=99999
-!        cloud_max=-1
-!        cloud_col_wat=0
-
-!        do kz=1, nz
-!         !clw & ciw are given in kg/kg  
-          ! cloud_water=clwc(ix,jy,kz,n)+ciwc(ix,jy,kz,n)
-          ! if (cloud_water .gt. 0) then
-          !   cloud_min=min(nint(height(kz)),cloud_min) !hg needs reset each grid
-          !   cloud_max=max(nint(height(kz)),cloud_max) !hg needs reset each grid
-          !   cloud_col_wat=cloud_col_wat+cloud_water !hg needs reset each grid 
-          ! endif
-          ! cloud_ver=max(0,cloud_max-cloud_min)
-
-          ! icloudbot(ix,jy,n)=cloud_min
-          ! icloudthck(ix,jy,n)=cloud_ver
-          ! rcw(ix,jy)=cloud_col_wat
-          ! rpc(ix,jy)=prec
-!write(*,*) 'Using clouds from ECMWF' !hg END Henrik Code
-!END HG METHOD
-
-
-
 !      else ! windfields does not contain cloud data 
 !          rhmin = 0.90 ! standard condition for presence of clouds
@@ -698,30 +824,6 @@
 !      endif !hg read clouds
 
-    end do
-  end do
-
-!do 102 kz=1,nuvz
-!write(an,'(i02)') kz+10
-!write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
-!open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
-!do 101 jy=0,nymin1
-!    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
-!101   continue
-! close(4)
-!102   continue
-
-! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
-! do 103 jy=0,nymin1
-!     write (4,*)
-!+       (height(kz),kz=1,nuvz)
-!103   continue
-! close(4)
-
-!open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
-! do 104 jy=0,nymin1
-!     write (4,*)
-!+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
-!104   continue
-! close(4)
+
+
 
 !eso measure CPU time
Index: src/wetdepo.f90
===================================================================
--- src/wetdepo.f90	(revision 0f20c31f55c80f501e55c69aad5a2dd4081628f1)
+++ src/wetdepo.f90	(revision d6a09777de02ba48864df5a70d9160b51767c086)
@@ -90,4 +90,5 @@
 
   integer :: blc_count, inc_count
+  real    :: Si_dummy, wetscav_dummy
 
 
@@ -190,4 +191,5 @@
       clouds_h=cloudsh(ix,jy,n)
     else
+      ! new removal not implemented for nests yet 
       clouds_v=cloudsn(ix,jy,hz,n,ngrid)
       clouds_h=cloudsnh(ix,jy,n,ngrid)
@@ -231,5 +233,7 @@
 
 
-  !ZHG calculated for 1) both 2) lsp 3) convp 
+  !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)))
@@ -252,4 +256,5 @@
       wetscav=0.   
 
+      !ZHG test if it nested?
       if (ngrid.gt.0) then
         act_temp=ttn(ix,jy,hz,n,ngrid)
@@ -259,5 +264,5 @@
      
 
-  !****i*******************
+  !***********************
   ! BELOW CLOUD SCAVENGING
   !***********************  
@@ -267,26 +272,25 @@
           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
-            wetscav=weta(ks)*prec(1)**wetb(ks) 
-
-  !AEROSOL
+            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.
+!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
+! 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
             if (act_temp .ge. 273 .and. weta(ks).gt.0.)  then !Rain 
-
-  !Particle RAIN scavenging coefficient based on Laakso et al 2003, the below-cloud scavenging (rain efficienty) parameter A (=weta) from SPECIES file
+              ! ZHG 2014 : 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 .and. wetb(ks).gt.0.)  then !snow
-
-  !Particle SNOW scavenging coefficient based on Kyro et al 2009, the below-cloud scavenging (Snow efficiency) parameter B (=wetb) from SPECIES file
+            elseif (act_temp .lt. 273 .and. wetb(ks).gt.0.)  then ! Snow
+              ! ZHG 2014 : 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))
@@ -298,28 +302,26 @@
           endif !gas or particle
         endif ! positive below-cloud scavenging parameters given in Species file
-      endif !end below-cloud
+      endif !end BELOW
 
 
   !********************
   ! IN CLOUD SCAVENGING
-  !********************
-      if (clouds_v.lt.4) then !in-cloud
-
+      !******************************************************
+      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
-        if (weta_in(ks).gt.0. .or. wetb_in(ks).gt.0.) then !if positive in-cloud parameters given in SPECIES file (either Ai or Bi)
-
+        if (weta_in(ks).gt.0. .or. wetb_in(ks).gt.0.) then 
 ! if negative coefficients (turned off) set to zero for use in equation
           if (weta_in(ks).lt.0.) weta_in(ks)=0.
           if (wetb_in(ks).lt.0.) wetb_in(ks)=0.
 
-          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
+          !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
+          if (readclouds) then                  !icloud_stats(ix,jy,4,n) has units kg/m2
+            cl =icloud_stats(ix,jy,4,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
@@ -329,9 +331,5 @@
               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)
+           ! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
             frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
  
@@ -339,23 +337,60 @@
 
   ! AEROSOL
+          !**************************************************
           if (dquer(ks).gt. 0.) then ! is particle
 
             S_i= frac_act/cl
 
+          !*********************
   ! GAS
           else ! is gas
                
             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
-            S_i=frac_act/cle
-
+            !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) then
+           wetscav=S_i*(prec(1)/3.6E6)
+          else
           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=,', &
-!          clouds_h,',cl=',cl, 'diff to old scheme=', cl-2E-7*prec(1)**0.36
+          endif
+
+!ZHG 2015 TEST
+!          Si_dummy=frac_act/2E-7*prec(1)**0.36
+!           wetscav_dummy=Si_dummy*(prec(1)/3.6E6)/clouds_h
+!           if (clouds_v.lt.4) then
+!           talltest=talltest+1
+!if(talltest .eq. 1) OPEN(UNIT=199, FILE=utfil,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!if(talltest .lt. 100001)  write(199,*) prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90
+!if(talltest .lt. 100001)  write(199,*) wetscav, wetscav_dummy,prec(1),ytra1(jpart)-90,clouds_v,cl
+!if(talltest .eq. 100001) CLOSE(199)
+!if(talltest .eq. 100001) STOP
+!
+!write(*,*)  'PREC kg/m2s CLOUD kg/m2', (prec(1)/3.6E6), cl !, '2E-7*prec(1)**0.36',  2E-7*prec(1)**0.36,'2E-7*prec(1)**0.36*clouds_h',2E-7*prec(1)**0.36*clouds_h
+!write(*,*)  'PREC kg/m2s LSP+convp kg/m2', prec(1), convp+lsp
+!write(*,*)  wetscav, wetscav_dummy
+!write(*,*) cc, grfraction(1), cc/grfraction(1)
+
+!write(*,*)  'Lmbda_old', (prec(1)/3.6E6)/(clouds_h*2E-7*prec(1)**0.36)
+
+
+!write(*,*) '**************************************************'
+!write(*,*)  'clouds_h', clouds_h, 'clouds_v',clouds_v,'abs(ltsample)', abs(ltsample)
+!write(*,*)  'readclouds', readclouds, 'wetscav',wetscav, 'wetscav_dummy', wetscav_dummy
+!write(*,*)  'S_i', S_i , 'Si_dummy', Si_dummy, 'prec(1)', prec(1) 
+
+
+!           write(*,*) 'PRECIPITATION ,cl  ECMWF , cl PARAMETIZED, clouds_v, lat' &
+!                      ,prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90
+
+!endif 
+
         endif ! positive in-cloud scavenging parameters given in Species file
       endif !incloud
+!END ZHG TEST
      
   !**************************************************
@@ -368,4 +403,8 @@
         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.
