Fixes for CALIPSO simulator and MG precipitation fluxes

Merge branch 'wlin/atm/calipsofix-mgflx' into next (PR #1333)

This fixes CALIPSO simulator that shows abnormal spatial distribution.
The root cause is that the pressure cutoff for MG calculation leads to
precipitation flux arrays containing undefined values above the cut off
level. It was not an issue with original L30 configuration.

Fixes for precipitation flux calculations in MG microphysics are also included.

modified: components/cam/src/physics/cam/cospsimulator_intr.F90
modified: components/cam/src/physics/cam/micro_mg1_0.F90
modified: components/cam/src/physics/cam/micro_mg2_0.F90
modified: components/cam/src/physics/cam/micro_mg_cam.F90

[BFB]

components/cam/src/physics/cam/cospsimulator_intr.F90

@@ -2456,6 +2456,32 @@ subroutine cospsimulator_intr_run(state,pbuf, cam_in,emis,coszrs,cliqwp_in,cicew
           reff_cosp(1:ncol,1:pver,8) = 0._r8
       end where
 
+      ! reset negative values to avoid COSP warning
+      where (rel(1:ncol,1:pver) .lt. 0.0_r8)
+          reff_cosp(1:ncol,1:pver,1) = 0.0_r8
+      end where
+      where (rei(1:ncol,1:pver) .lt. 0.0_r8)
+          reff_cosp(1:ncol,1:pver,2) = 0.0_r8
+      end where
+      where (ls_reffrain(1:ncol,1:pver) .lt. 0.0_r8)
+          reff_cosp(1:ncol,1:pver,3) = 0.0_r8
+      end where
+      where (ls_reffsnow(1:ncol,1:pver) .lt. 0.0_r8)
+          reff_cosp(1:ncol,1:pver,4) = 0.0_r8
+      end where
+      where (cv_reffliq(1:ncol,1:pver) .lt. 0.0_r8)
+          reff_cosp(1:ncol,1:pver,5) = 0.0_r8
+      end where
+      where (cv_reffice(1:ncol,1:pver) .lt. 0.0_r8)
+          reff_cosp(1:ncol,1:pver,6) = 0.0_r8
+      end where
+      where (ls_reffrain(1:ncol,1:pver) .lt. 0.0_r8)
+          reff_cosp(1:ncol,1:pver,7) = 0.0_r8
+      end where
+      where (ls_reffsnow(1:ncol,1:pver) .lt. 0.0_r8)
+          reff_cosp(1:ncol,1:pver,8) = 0.0_r8
+      end where
+
    end if  !!if cam5
 
    !! Make sure interpolated values are not less than 0 - COSP was complaining and resetting small negative values to zero.

components/cam/src/physics/cam/micro_mg1_0.F90

@@ -1449,6 +1449,8 @@ subroutine micro_mg_tend ( &
       nstend(i,1:pver)=0._r8
       prect(i)=0._r8
       preci(i)=0._r8
+      rflx(i,1:pver+1)=0._r8
+      sflx(i,1:pver+1)=0._r8
       qniic(i,1:pver)=0._r8
       qric(i,1:pver)=0._r8
       nsic(i,1:pver)=0._r8
@@ -1502,6 +1504,10 @@ subroutine micro_mg_tend ( &
       prect(i)=0._r8
       preci(i)=0._r8
 
+      ! initialize fluxes
+      rflx(i,1:pver+1)=0._r8
+      sflx(i,1:pver+1)=0._r8
+
       do k=top_lev,pver
 
          qcvar=relvar(i,k)
@@ -2822,6 +2828,10 @@ subroutine micro_mg_tend ( &
             end if
          end if
 
+         ! Precip Flux Calculation (Diagnostic)
+         rflx(i,k+1)=(prect(i)-preci(i)) * rhow
+         sflx(i,k+1)=preci(i) * rhow
+
          ! if rain/snow mix ratio is zero so should number concentration
 
          if (qniic(i,k).lt.qsmall) then
@@ -2963,15 +2973,6 @@ subroutine micro_mg_tend ( &
             rercld(i,k)=rercld(i,k)/arcld(i,k)
          end if
 
-         !calculate precip fluxes and adding them to summing sub-stepping variables
-         !! flux is zero at top interface
-         rflx(i,1)=0.0_r8
-         sflx(i,1)=0.0_r8
-
-         !! calculating the precip flux (kg/m2/s) as mixingratio(kg/kg)*airdensity(kg/m3)*massweightedfallspeed(m/s)
-         rflx(i,k+1)=qrout(i,k)*rho(i,k)*umr(k)
-         sflx(i,k+1)=qsout(i,k)*rho(i,k)*ums(k)
-
          !! add to summing sub-stepping variable
          rflx1(i,k+1)=rflx1(i,k+1)+rflx(i,k+1)
          sflx1(i,k+1)=sflx1(i,k+1)+sflx(i,k+1)
@@ -3311,6 +3312,12 @@ subroutine micro_mg_tend ( &
       prect(i) = prect(i)+(faloutc(pver)+falouti(pver))/g/nstep/1000._r8  
       preci(i) = preci(i)+(falouti(pver))/g/nstep/1000._r8
 
+      ! Add fallout to Precip Flux: note unit change m/s *kg/m3 = kg/m2
+      do k = top_lev,pver
+         rflx(i,k+1)=rflx(i,k+1)+(faloutc(k))/g/nstep/1000._r8 * rhow
+         sflx(i,k+1)=sflx(i,k+1)+(falouti(k))/g/nstep/1000._r8 * rhow
+      end do
+
    end do   !! nstep loop
 
    ! end sedimentation

components/cam/src/physics/cam/micro_mg2_0.F90

@@ -366,6 +366,7 @@ subroutine micro_mg_tend ( &
      cmeout,                       deffi,                        &
      pgamrad,                      lamcrad,                      &
      qsout,                        dsout,                        &
+     lflx,               iflx,                                   &
      rflx,               sflx,               qrout,              &
      reff_rain,                    reff_snow,                    &
      qcsevap,            qisevap,            qvres,              &
@@ -501,6 +502,8 @@ subroutine micro_mg_tend ( &
   real(r8), intent(out) :: lamcrad(:,:)      ! slope of droplet distribution for optics (radiation) (1/m)
   real(r8), intent(out) :: qsout(:,:)        ! snow mixing ratio (kg/kg)
   real(r8), intent(out) :: dsout(:,:)        ! snow diameter (m)
+  real(r8), intent(out) :: lflx(:,:)         ! grid-box average liquid condensate flux (kg m^-2 s^-1)
+  real(r8), intent(out) :: iflx(:,:)         ! grid-box average ice condensate flux (kg m^-2 s^-1)
   real(r8), intent(out) :: rflx(:,:)         ! grid-box average rain flux (kg m^-2 s^-1)
   real(r8), intent(out) :: sflx(:,:)         ! grid-box average snow flux (kg m^-2 s^-1)
   real(r8), intent(out) :: qrout(:,:)        ! grid-box average rain mixing ratio (kg/kg)
@@ -951,6 +954,8 @@ subroutine micro_mg_tend ( &
 
   rflx=0._r8
   sflx=0._r8
+  lflx=0._r8
+  iflx=0._r8
 
   ! initialize precip output
 
@@ -1950,13 +1955,6 @@ subroutine micro_mg_tend ( &
   qsout = qs
   nsout = ns * rho
 
-  ! calculate precip fluxes
-  ! calculate the precip flux (kg/m2/s) as mixingratio(kg/kg)*airdensity(kg/m3)*massweightedfallspeed(m/s)
-  ! ---------------------------------------------------------------------
-
-  rflx(:,2:) = rflx(:,2:) + (qric*rho*umr*precip_frac)
-  sflx(:,2:) = sflx(:,2:) + (qsic*rho*ums*precip_frac)
-
   ! calculate n0r and lamr from rain mass and number
   ! divide by precip fraction to get in-precip (local) values of
   ! rain mass and number, divide by rhow to get rain number in kg^-1
@@ -2217,6 +2215,11 @@ subroutine micro_mg_tend ( &
 
         end do
 
+        ! Ice flux
+        do k = 1,nlev
+          iflx(i,k+1) = iflx(i,k+1) + falouti(k) / g / real(nstep)
+        end do
+
         ! units below are m/s
         ! sedimentation flux at surface is added to precip flux at surface
         ! to get total precip (cloud + precip water) rate
@@ -2283,6 +2286,11 @@ subroutine micro_mg_tend ( &
 
         end do
 
+        !Liquid condensate flux here
+        do k = 1,nlev
+           lflx(i,k+1) = lflx(i,k+1) + faloutc(k) / g / real(nstep)
+        end do
+
         prect(i) = prect(i)+faloutc(nlev)/g/real(nstep)/1000._r8
 
      end do
@@ -2334,6 +2342,11 @@ subroutine micro_mg_tend ( &
 
         end do
 
+        ! Rain Flux
+        do k = 1,nlev
+           rflx(i,k+1) = rflx(i,k+1) + faloutr(k) / g / real(nstep)
+        end do
+
         prect(i) = prect(i)+faloutr(nlev)/g/real(nstep)/1000._r8
 
      end do
@@ -2385,6 +2398,11 @@ subroutine micro_mg_tend ( &
 
         end do   !! k loop
 
+        ! Snow Flux
+        do k = 1,nlev
+           sflx(i,k+1) = sflx(i,k+1) + falouts(k) / g / real(nstep)
+        end do
+
         prect(i) = prect(i)+falouts(nlev)/g/real(nstep)/1000._r8
         preci(i) = preci(i)+falouts(nlev)/g/real(nstep)/1000._r8
 

components/cam/src/physics/cam/micro_mg_cam.F90

@@ -1120,6 +1120,8 @@ subroutine micro_mg_cam_tend(state, ptend, dtime, pbuf)
    real(r8), target :: prodsnow(state%psetcols,pver)   ! Local production of snow
    real(r8), target :: cmeice(state%psetcols,pver)     ! Rate of cond-evap of ice within the cloud
    real(r8), target :: qsout(state%psetcols,pver)      ! Snow mixing ratio
+   real(r8), target :: cflx(state%psetcols,pverp)      ! grid-box avg liq condensate flux (kg m^-2 s^-1)
+   real(r8), target :: iflx(state%psetcols,pverp)      ! grid-box avg ice condensate flux (kg m^-2 s^-1)
    real(r8), target :: rflx(state%psetcols,pverp)      ! grid-box average rain flux (kg m^-2 s^-1)
    real(r8), target :: sflx(state%psetcols,pverp)      ! grid-box average snow flux (kg m^-2 s^-1)
    real(r8), target :: qrout(state%psetcols,pver)      ! Rain mixing ratio
@@ -1246,6 +1248,8 @@ subroutine micro_mg_cam_tend(state, ptend, dtime, pbuf)
    real(r8), allocatable, target :: packed_prodsnow(:,:)
    real(r8), allocatable, target :: packed_cmeout(:,:)
    real(r8), allocatable, target :: packed_qsout(:,:)
+   real(r8), allocatable, target :: packed_cflx(:,:)
+   real(r8), allocatable, target :: packed_iflx(:,:)
    real(r8), allocatable, target :: packed_rflx(:,:)
    real(r8), allocatable, target :: packed_sflx(:,:)
    real(r8), allocatable, target :: packed_qrout(:,:)
@@ -1817,6 +1821,10 @@ subroutine micro_mg_cam_tend(state, ptend, dtime, pbuf)
    call post_proc%add_field(p(cmeice), p(packed_cmeout))
    allocate(packed_qsout(mgncol,nlev))
    call post_proc%add_field(p(qsout), p(packed_qsout))
+   allocate(packed_cflx(mgncol,nlev+1))
+   call post_proc%add_field(p(cflx), p(packed_cflx))
+   allocate(packed_iflx(mgncol,nlev+1))
+   call post_proc%add_field(p(iflx), p(packed_iflx))
    allocate(packed_rflx(mgncol,nlev+1))
    call post_proc%add_field(p(rflx), p(packed_rflx))
    allocate(packed_sflx(mgncol,nlev+1))
@@ -2144,6 +2152,7 @@ subroutine micro_mg_cam_tend(state, ptend, dtime, pbuf)
                  packed_cmeout,          packed_dei,             &
                  packed_mu,              packed_lambdac,         &
                  packed_qsout,           packed_des,             &
+                 packed_cflx,    packed_iflx,                    &
                  packed_rflx,    packed_sflx,    packed_qrout,   &
                  reff_rain_dum,          reff_snow_dum,          &
                  packed_qcsevap, packed_qisevap, packed_qvres,   &
@@ -2264,14 +2273,26 @@ subroutine micro_mg_cam_tend(state, ptend, dtime, pbuf)
       end do
    end do
 
+   ! array must be zeroed beyond trop_cloud_top_pre otherwise undefined values will be used in cosp.
+   mgflxprc(:ncol,1:top_lev) = 0.0_r8
+   mgflxsnw(:ncol,1:top_lev) = 0.0_r8
+
    mgflxprc(:ncol,top_lev:pverp) = rflx(:ncol,top_lev:pverp) + sflx(:ncol,top_lev:pverp)
    mgflxsnw(:ncol,top_lev:pverp) = sflx(:ncol,top_lev:pverp)
 
+!ADD CONDENSATE FLUXES FOR MG2 (ice and snow already added for MG1)
+   if (micro_mg_version .ge. 2) then
+      mgflxprc(:ncol,top_lev:pverp) = mgflxprc(:ncol,top_lev:pverp) + iflx(:ncol,top_lev:pverp) + cflx(:ncol,top_lev:pverp)
+      mgflxsnw(:ncol,top_lev:pverp) = mgflxsnw(:ncol,top_lev:pverp) + iflx(:ncol,top_lev:pverp)
+   end if
+
    mgmrprc(:ncol,top_lev:pver) = qrout(:ncol,top_lev:pver) + qsout(:ncol,top_lev:pver)
    mgmrsnw(:ncol,top_lev:pver) = qsout(:ncol,top_lev:pver)
 
    !! calculate effective radius of convective liquid and ice using dcon and deicon (not used by code, not useful for COSP)
    !! hard-coded as average of hard-coded values used for deep/shallow convective detrainment (near line 1502/1505)
+   ! this needs to be replaced by clubb_liq_deep and clubb_ice+deep accordingly
+
    cvreffliq(:ncol,top_lev:pver) = 9.0_r8
    cvreffice(:ncol,top_lev:pver) = 37.0_r8
 
@@ -2358,6 +2379,7 @@ subroutine micro_mg_cam_tend(state, ptend, dtime, pbuf)
                cldfsnow(i,k) = 0.25_r8
             endif
          endif
+
          ! Calculate in-cloud snow water path
          icswp(i,k) = qsout(i,k) / max( mincld, cldfsnow(i,k) ) * state_loc%pdel(i,k) / gravit
       end do