Code for cloud water should be correct if the total cw is stored in field clwc...

Code for cloud water should be correct if the total cw is stored in field clwc (old scheme) or in field qc (new scheme). Minor edits in some files.

src/advance.f90

@@ -549,19 +549,6 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
   !if (mod(itime,10800).ne.0) dump=.true.
   !!! CHANGE
       
-  !!!----- TEST OF THE WELL-MIXED CRITERION: modified by mc,  not to be included in final version mc
-      ! if (zt.lt.h) then
-      !     i_well=int(zt/h*25.)+1
-      !     well_mixed_vector(i_well)=well_mixed_vector(i_well)+dt
-      !     well_mixed_norm=well_mixed_norm+dt
-      !     avg_air_dens(i_well)=avg_air_dens(i_well)+rhoa*dt
-      !     avg_wst=avg_wst+wst*dt
-      !     avg_ol=avg_ol+ol*dt
-      !     avg_h=avg_h+h*dt
-      ! end if
-      ! h_well=h
-  !------- END TEST
-      
   ! Determine probability of deposition
   !************************************
 

src/conccalc.f90

@@ -20,7 +20,7 @@
 !**********************************************************************
 
 subroutine conccalc(itime,weight)
-  !                      i     i
+  !                 i     i
   !*****************************************************************************
   !                                                                            *
   !     Calculation of the concentrations on a regular grid using volume       *

src/concoutput_mpi.f90

@@ -320,7 +320,6 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
 ! CONCENTRATION OR MIXING RATIO
             do kz=1,numzgrid
               do l=1,nclassunc
-!                auxgrid(l)=gridunc0(ix,jy,kz,ks,kp,l,nage)
                 auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
               end do
               call mean(auxgrid,grid(ix,jy,kz), &

src/gridcheck.f90

@@ -83,7 +83,7 @@ subroutine gridcheck
   integer :: gotGrid
   real(kind=4) :: xaux1,xaux2,yaux1,yaux2
   real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
-  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parID
+  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
   !HSO  end
   integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
   real :: sizesouth,sizenorth,xauxa,pint,conversion_factor
@@ -197,9 +197,16 @@ subroutine gridcheck
 !ZHG FOR CLOUDS FROM GRIB
   elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
     isec1(6)=246         ! indicatorOfParameter
+    ! readclouds=.true.
+    ! sumclouds=.false.
   elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
     isec1(6)=247         ! indicatorOfParameter
 !ZHG end
+! ESO qc(=clwc+ciwc)
+  elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
+    isec1(6)=201031      ! indicatorOfParameter
+    ! readclouds=.true.
+    ! sumclouds=.true.
   elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
     isec1(6)=134         ! indicatorOfParameter
   elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot

src/makefile

@@ -417,6 +417,7 @@ timemanager.o: com_mod.o flux_mod.o netcdf_output_mod.o oh_mod.o outg_mod.o \
 	par_mod.o point_mod.o unc_mod.o xmass_mod.o
 timemanager_mpi.o: com_mod.o flux_mod.o mpi_mod.o oh_mod.o outg_mod.o \
 	par_mod.o point_mod.o unc_mod.o xmass_mod.o
+unc_mod.o: par_mod.o
 verttransform.o: cmapf_mod.o com_mod.o par_mod.o 
 verttransform_gfs.o: cmapf_mod.o com_mod.o par_mod.o
 verttransform_nests.o: com_mod.o par_mod.o

src/par_mod.f90

@@ -41,11 +41,14 @@ module par_mod
 
   !****************************************************************
   ! Parameter defining KIND parameter for double/single precision
+  ! 'cp' sets the precision for some internal calculations (sp or dp),
+  ! sp is default but dp can be used for increased precision.
+  ! WARNING: for now, do not use 'dp'
   !****************************************************************
 
   integer,parameter :: dp=selected_real_kind(P=15)
   integer,parameter :: sp=selected_real_kind(6)
-
+  integer,parameter :: cp=sp ! sp is default
 
   !***********************************************************
   ! Number of directories/files used for FLEXPART input/output

src/readwind_mpi.f90

@@ -457,7 +457,7 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
     call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
 !hg
     call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
-    call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
+    if (.not.sumclouds) call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
 !hg end
 
   endif

src/timemanager.f90

@@ -112,6 +112,7 @@ subroutine timemanager
   ! integer(kind=2) :: cbt(maxpart)
   real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc
   real :: drygridtotalunc,xold,yold,zold,xmassfract
+  real, parameter :: e_inv = 1.0/exp(1.0)
   !double precision xm(maxspec,maxpointspec_act),
   !    +                 xm_depw(maxspec,maxpointspec_act),
   !    +                 xm_depd(maxspec,maxpointspec_act)
@@ -593,7 +594,7 @@ subroutine timemanager
                    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
+                   ! if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.inv_e.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))

src/timemanager_mpi.f90

@@ -115,6 +115,7 @@ subroutine timemanager
 
   real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc
   real :: drygridtotalunc,xold,yold,zold,xmassfract
+  real, parameter :: e_inv = 1.0/exp(1.0)
 !double precision xm(maxspec,maxpointspec_act),
 !    +                 xm_depw(maxspec,maxpointspec_act),
 !    +                 xm_depd(maxspec,maxpointspec_act)
@@ -754,7 +755,7 @@ subroutine timemanager
                    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
+                   ! if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.inv_e.and.checklifetime(j,ks).eq.0.)then
                        !Mass below 1% of initial >register lifetime
                    !     checklifetime(j,ks)=abs(itra1(j)-itramem(j))
 

src/unc_mod.f90

@@ -21,20 +21,22 @@
 
 module unc_mod
 
+  use par_mod, only:cp
+
   implicit none
 
   real,allocatable ,dimension (:,:,:,:,:,:,:) :: gridunc
   real,allocatable, dimension (:,:,:,:,:,:,:) :: griduncn
-  real,allocatable, dimension (:,:,:,:,:,:) :: drygridunc
-  real,allocatable, dimension (:,:,:,:,:,:) :: drygriduncn
-  real,allocatable, dimension (:,:,:,:,:,:) :: wetgridunc
-  real,allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn
+  real(cp),allocatable, dimension (:,:,:,:,:,:) :: drygridunc
+  real(cp),allocatable, dimension (:,:,:,:,:,:) :: drygriduncn
+  real(cp),allocatable, dimension (:,:,:,:,:,:) :: wetgridunc
+  real(cp),allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn
 
 ! For sum of individual contributions, used for the MPI version
-  real,allocatable, dimension (:,:,:,:,:,:) :: drygridunc0
-  real,allocatable, dimension (:,:,:,:,:,:) :: drygriduncn0
-  real,allocatable, dimension (:,:,:,:,:,:) :: wetgridunc0
-  real,allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn0
+  real(cp),allocatable, dimension (:,:,:,:,:,:) :: drygridunc0
+  real(cp),allocatable, dimension (:,:,:,:,:,:) :: drygriduncn0
+  real(cp),allocatable, dimension (:,:,:,:,:,:) :: wetgridunc0
+  real(cp),allocatable, dimension (:,:,:,:,:,:) :: wetgriduncn0
 
   real,allocatable, dimension (:,:,:,:,:) :: init_cond