bugfix: MPI version gave wrong wet deposition when using ECMWF cloud water...

bugfix: MPI version gave wrong wet deposition when using ECMWF cloud water fields. Cloud water in ECMWF fields now uses parameter qc, or reads clwc+ciwc. Added minmass variable as limit for terminating particles.

src/com_mod.f90

@@ -350,17 +350,17 @@ module com_mod
   real :: tt(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: qv(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 :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !liquid   [kg/kg]
+  real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !ice      [kg/kg]
+  real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0  !combined [m3/m3]
 
   real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: drhodz(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
   real :: tth(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
   real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
-  real :: clwch(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
-  real :: ciwch(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
+  real :: clwch(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)=0.0
+  real :: ciwch(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)=0.0
 
   real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
   !scavenging NIK, PS
@@ -370,6 +370,7 @@ module com_mod
 !ZHG Sep 2015  
    real :: icloud_stats(0:nxmax-1,0:nymax-1,5,numwfmem)
    real :: clh(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
+   real :: clw4(0:nxmax-1,0:nymax-1,numwfmem) ! eso: =icloud_stats(:,:,4,:)
 
 
   ! uu,vv,ww [m/2]       wind components in x,y and z direction

src/mpi_mod.f90

@@ -118,7 +118,7 @@ module mpi_mod
   logical, parameter :: mp_dev_mode = .false.
   logical, parameter :: mp_dbg_out = .false.
   logical, parameter :: mp_time_barrier=.true.
-  logical, parameter :: mp_measure_time=.false.
+  logical, parameter :: mp_measure_time=.true.
   logical, parameter :: mp_exact_numpart=.true.
 
 ! for measuring CPU/Wall time
@@ -906,7 +906,7 @@ contains
 ! Send variables from getfield process (id_read) to other processes
 !**********************************************************************
 
-! The non-reader processes need to know if clouds were read.
+! The non-reader processes need to know if cloud water 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 
@@ -962,9 +962,10 @@ contains
 
 ! cloud water/ice:
     if (readclouds) then
-      call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2_size1*5,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
+      ! call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2s1*5,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
+      ! if (mp_ierr /= 0) goto 600
+      call MPI_Bcast(clw4(:,:,li:ui),d2s1,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)
@@ -1335,8 +1336,11 @@ contains
 ! Send cloud water if it exists. Increment counter always (as on receiving end)
       if (readclouds) then
         i=i+1
-        call MPI_Isend(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,dest,tm1,&
+        ! call MPI_Isend(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,dest,tm1,&
+        !      &MPI_COMM_WORLD,reqs(i),mp_ierr)
+        call MPI_Isend(clw4(:,:,mind),d2s1,mp_pp,dest,tm1,&
              &MPI_COMM_WORLD,reqs(i),mp_ierr)
+
         if (mp_ierr /= 0) goto 600
 
         ! call MPI_Isend(clwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
@@ -1542,7 +1546,9 @@ contains
     if (readclouds) then
       j=j+1
 
-      call MPI_Irecv(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,id_read,MPI_ANY_TAG,&
+      ! call MPI_Irecv(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,id_read,MPI_ANY_TAG,&
+      !      &MPI_COMM_WORLD,reqs(j),mp_ierr)
+      call MPI_Irecv(clw4(:,:,mind),d2s1*5,mp_pp,id_read,MPI_ANY_TAG,&
            &MPI_COMM_WORLD,reqs(j),mp_ierr)
       if (mp_ierr /= 0) goto 600
 

src/par_mod.f90

@@ -210,13 +210,12 @@ module par_mod
   !**************************************************
 
   integer,parameter :: maxpart=40000000
-!  integer,parameter :: maxpart=60000000
-!  integer,parameter :: maxpart=120000000
   integer,parameter :: maxspec=6
+  integer,parameter :: minmass=0.0 !0.0001
 
   ! maxpart                 Maximum number of particles
   ! maxspec                 Maximum number of chemical species per release
-
+  ! minmass                 Terminate particles carrying less mass
 
   ! maxpoint is also set dynamically during runtime
   ! maxpoint                Maximum number of release locations

src/readwind.f90

@@ -189,13 +189,14 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
       isec1(6)=132         ! indicatorOfParameter
     elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
       isec1(6)=133         ! indicatorOfParameter
-!ZHG ! READ CLOUD FIELD : assume these to be added together to one variable
+! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field qc 
     elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
       isec1(6)=246         ! indicatorOfParameter
-! 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.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
+      isec1(6)=247         ! indicatorOfParameter
+! ESO qc(=clwc+ciwc):
+    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
+      isec1(6)=201031         ! indicatorOfParameter
     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
@@ -359,16 +360,24 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
       if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
       if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
 !ZHG READING CLOUD FIELDS ASWELL
+! ESO TODO: add check for if one of clwc/ciwc missing (error),
+! also if all 3 cw fields present, use qc and disregard the others
       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.
-!if (clwch(i,j,nlev_ec-k+2,n) .gt. 0)        write(*,*) 'readwind: found water!', clwch(i,j,nlev_ec-k+2,n)
+        readclouds=.true.
+        sumclouds=.false.
+      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)
       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)
-        !write(*,*) 'found ice!'
-!      endif
 !ZHG end
+!ESO read qc (=clwc+ciwc)
+      if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
+        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        readclouds=.true.
+        sumclouds=.true.
+!if (clwch(i,j,nlev_ec-k+2,n) .gt. 0)        write(*,*) 'readwind: found water!', clwch(i,j,nlev_ec-k+2,n)
+      endif
 
     end do
   end do
@@ -426,7 +435,7 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
     call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
 !ZHG
     call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
-    call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
+    if (.not.sumclouds) call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
 !ZHG end
 
   endif

src/readwind_mpi.f90

@@ -82,9 +82,9 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
   integer :: gotGrid
 !HSO  end
 
-  real(kind=4), intent(inout) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real(kind=4), intent(inout) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
-  real(kind=4), intent(inout) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+  real(sp), intent(inout) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real(sp), intent(inout) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real(sp), intent(inout) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
   integer, intent(in) :: indj,n
   integer :: i,j,k,levdiff2,ifield,iumax,iwmax
 
@@ -115,7 +115,7 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
 
   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
@@ -199,13 +199,14 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
       isec1(6)=132         ! indicatorOfParameter
     elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
       isec1(6)=133         ! indicatorOfParameter
-!ZHG ! READ CLOUD FIELD : assume these to be added together to one variable
+! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field qc 
     elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
       isec1(6)=246         ! indicatorOfParameter
-! 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.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
+      isec1(6)=247         ! indicatorOfParameter
+! ESO qc(=clwc+ciwc):
+    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
+      isec1(6)=201031         ! indicatorOfParameter
     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
@@ -372,13 +373,18 @@ subroutine readwind(indj,n,uuh,vvh,wwh)
       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.
-!if (clwch(i,j,nlev_ec-k+2,n) .gt. 0)        write(*,*) 'readwind: found water!', clwch(i,j,nlev_ec-k+2,n)
+        sumclouds=.false.
+      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)
       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)
-        !write(*,*) 'found ice!'
-!      endif
 !ZHG end
+!ESO read qc (=clwc+ciwc)
+      if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
+        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        readclouds=.true.
+        sumclouds=.true.
+      endif
 
     end do
   end do

src/timemanager.f90

@@ -606,7 +606,7 @@ subroutine timemanager
             endif
           end do
 
-          if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
+          if (xmassfract.lt.minmass) then   ! terminate all particles carrying less mass
             itra1(j)=-999999999
             if (verbosity.gt.0) then
               print*,'terminated particle ',j,' for small mass'
@@ -654,21 +654,6 @@ subroutine timemanager
   ! they are compromising the final result (or not)
     if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl  
           
-!!------------------------------------------------------------------------------
-! these lines below to test the well-mixed condition, modified by  mc, not to be
-! included in final release:
-    ! if (itime.eq.0) then
-    ! open(551,file='/home/mc/test_cbl/out/WELLMIXEDTEST_CBL_lonlat_9_33_100_3hours_3htp_cd.DAT')
-    ! open(552,file='/home/mc/test_cbl/out/avg_ol_h_wst_lonlat_9_33_100_3hours_3htp_cd.DAT')
-    ! end if
-    ! write(552,'(5F16.7)')itime*1./3600.,avg_wst/well_mixed_norm,avg_ol/well_mixed_norm,avg_h/well_mixed_norm
-    ! do j=1,25
-    !      !write(551,*))itime*1.,h_well/50.*j,well_mixed_vector(j)/well_mixed_norm*50.
-    !     avg_air_dens(j)=avg_air_dens(j)/well_mixed_vector(j)
-    !     write(551,'(5F16.7)')itime*1.,h_well/25.*j,well_mixed_vector(j)/well_mixed_norm*25.,avg_air_dens(j),0.04*j
-    ! end do 
-!!------------------------------------------------------------------------------
-
   end do
 
 

src/timemanager_mpi.f90

@@ -102,7 +102,6 @@ subroutine timemanager
 
   implicit none
 
-  logical :: fc_mp=.true., fc_sync=.true.
   logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete
   integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0,mind
 ! integer :: ksp
@@ -351,7 +350,7 @@ subroutine timemanager
     if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
       do ks=1,nspec
         do kp=1,maxpointspec_act
-          if (decay(ks).gt.0.) then
+          if (decay(ks).gt.0.) then ! TODO move this statement up 2 levels
             do nage=1,nageclass
               do l=1,nclassunc
 ! Mother output grid
@@ -484,7 +483,6 @@ subroutine timemanager
           else 
             if (lroot) then
               if (lnetcdfout.eq.1) then
-
                 call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
                      &drygridtotalunc)
               else
@@ -771,7 +769,7 @@ subroutine timemanager
             endif
           end do
 
-          if (xmassfract.lt.0.00005 .and. sum(real(npart(npoint(j)))*xmass1(j,:)).lt.1.0) then   ! terminate all particles carrying less mass
+          if (xmassfract.lt.minmass) then ! .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

src/verttransform.f90

@@ -250,7 +250,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
   qv(:,:,1,n)=qvh(:,:,1,n)
 !hg adding the cloud water 
   clwc(:,:,1,n)=clwch(:,:,1,n)
-  ciwc(:,:,1,n)=ciwch(:,:,1,n)   
+  if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)   
 !hg 
   pv(:,:,1,n)=pvh(:,:,1)
   rho(:,:,1,n)=rhoh(:,:,1)
@@ -261,7 +261,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
 
 !hg adding the cloud water
   clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
-  ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
+  if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
 !hg
   pv(:,:,nz,n)=pvh(:,:,nuvz)
   rho(:,:,nz,n)=rhoh(:,:,nuvz)
@@ -279,7 +279,7 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
           qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
 !hg adding the cloud water
           clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
-          ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
+          if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
 !hg
           pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
           rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
@@ -309,7 +309,8 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
                +qvh(ix,jy,kz,n)*dz1)/dz
 !hg adding the cloud water
           clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
-          ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
+          if (.not.sumclouds) &
+               &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
 !hg
           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
           rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
@@ -577,9 +578,13 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
 ! to include future cloud processing by non-precipitating-clouds. 
 !***********************************************************************************
     write(*,*) 'using cloud water from ECMWF'
-    clw(:,:,:,n)=0
-    icloud_stats(:,:,:,n)=0
+    clw(:,:,:,n)=0.0
+    icloud_stats(:,:,:,n)=0.0
     clouds(:,:,:,n)=0
+! If water/ice are read separately into clwc and ciwc, store sum in clwc
+    if (.not.sumclouds) then 
+      clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
+    end if
     do jy=0,nymin1
       do ix=0,nxmin1
         lsp=lsprec(ix,jy,1,n)
@@ -634,7 +639,6 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
 !**************************************************************************
 !   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'
     write(*,*) 'using cloud water from Parameterization'
     do jy=0,nymin1
       do ix=0,nxmin1
@@ -676,6 +680,9 @@ subroutine verttransform(n,uuh,vvh,wwh,pvh)
     end do
   endif !readclouds
 
+! eso: copy the relevant data to clw4 to reduce amount of communicated data for MPI
+  clw4(:,:,n) = icloud_stats(:,:,4,n)
+
      !********* TEST ***************
      ! WRITE OUT SOME TEST VARIABLES
      !********* TEST ************'**

src/wetdepo.f90

@@ -316,7 +316,9 @@ subroutine wetdepo(itime,ltsample,loutnext)
 
           !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)
+!            cl =icloud_stats(ix,jy,4,n)*(grfraction(1)/cc)
+! ESO: stop using icloud_stats
+            cl =clw4(ix,jy,n)*(grfraction(1)/cc)
           else                                  !parameterize cloudwater m2/m3
             !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
             cl=1.6E-6*prec(1)**0.36