add "limiter #5 for simple condensation model"

let qme=0 when f = 0 and df/dt = 0

To be tested.

components/cam/src/physics/cam/phys_control.F90

@@ -165,6 +165,7 @@ module phys_control
 logical :: l_rkz_lmt_2       = .false.
 logical :: l_rkz_lmt_3       = .false.
 logical :: l_rkz_lmt_4       = .true.
+logical :: l_rkz_lmt_5       = .false.
 
 !======================================================================= 
 contains
@@ -202,7 +203,7 @@ subroutine phys_ctl_readnl(nlfile)
       l_tracer_aero, l_vdiff, l_rayleigh, l_gw_drag, l_ac_energy_chk, &
       l_bc_energy_fix, l_dry_adj, l_st_mac, l_st_mic, l_rad, &
       simple_macrop_opt, rkz_cldfrc_opt, rkz_term_A_opt, rkz_term_B_opt, rkz_term_C_opt, &
-      rkz_term_C_ql_opt, rkz_term_C_fmin, l_rkz_lmt_2, l_rkz_lmt_3, l_rkz_lmt_4, &
+      rkz_term_C_ql_opt, rkz_term_C_fmin, l_rkz_lmt_2, l_rkz_lmt_3, l_rkz_lmt_4, l_rkz_lmt_5, &
       prc_coef1,prc_exp,prc_exp1,cld_sed,mg_prc_coeff_fix, &
       rrtmg_temp_fix
    !-----------------------------------------------------------------------------
@@ -294,6 +295,7 @@ subroutine phys_ctl_readnl(nlfile)
    call mpibcast(l_rkz_lmt_2,                     1 , mpilog,  0, mpicom)
    call mpibcast(l_rkz_lmt_3,                     1 , mpilog,  0, mpicom)
    call mpibcast(l_rkz_lmt_4,                     1 , mpilog,  0, mpicom)
+   call mpibcast(l_rkz_lmt_5,                     1 , mpilog,  0, mpicom)
 
    call mpibcast(cld_macmic_num_steps,            1 , mpiint,  0, mpicom)
    call mpibcast(prc_coef1,                       1 , mpir8,   0, mpicom)
@@ -447,7 +449,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, mi
                        ,l_bc_energy_fix_out, l_dry_adj_out, l_st_mac_out, l_st_mic_out, l_rad_out  &
                        ,simple_macrop_opt_out, rkz_cldfrc_opt_out, rkz_term_A_opt_out, rkz_term_B_opt_out &
                        ,rkz_term_C_opt_out, rkz_term_C_ql_opt_out, rkz_term_C_fmin_out &
-                       ,l_rkz_lmt_2_out, l_rkz_lmt_3_out, l_rkz_lmt_4_out &
+                       ,l_rkz_lmt_2_out, l_rkz_lmt_3_out, l_rkz_lmt_4_out, l_rkz_lmt_5_out &
                        ,prc_coef1_out,prc_exp_out,prc_exp1_out, cld_sed_out,mg_prc_coeff_fix_out,rrtmg_temp_fix_out)
 
 !-----------------------------------------------------------------------
@@ -525,6 +527,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, mi
    logical,           intent(out), optional :: l_rkz_lmt_2_out
    logical,           intent(out), optional :: l_rkz_lmt_3_out
    logical,           intent(out), optional :: l_rkz_lmt_4_out
+   logical,           intent(out), optional :: l_rkz_lmt_5_out
 
    logical,           intent(out), optional :: mg_prc_coeff_fix_out
    logical,           intent(out), optional :: rrtmg_temp_fix_out
@@ -598,6 +601,8 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, mi
    if ( present(l_rkz_lmt_2_out         ) ) l_rkz_lmt_2_out       = l_rkz_lmt_2
    if ( present(l_rkz_lmt_3_out         ) ) l_rkz_lmt_3_out       = l_rkz_lmt_3
    if ( present(l_rkz_lmt_4_out         ) ) l_rkz_lmt_4_out       = l_rkz_lmt_4
+   if ( present(l_rkz_lmt_5_out         ) ) l_rkz_lmt_5_out       = l_rkz_lmt_5
+
    if ( present(cld_macmic_num_steps_out) ) cld_macmic_num_steps_out = cld_macmic_num_steps
    if ( present(prc_coef1_out           ) ) prc_coef1_out            = prc_coef1
    if ( present(prc_exp_out             ) ) prc_exp_out              = prc_exp

components/cam/src/physics/cam/physpkg.F90

@@ -2033,7 +2033,7 @@ subroutine tphysbc (ztodt,               &
     logical :: l_rkz_lmt_2
     logical :: l_rkz_lmt_3
     logical :: l_rkz_lmt_4
-
+    logical :: l_rkz_lmt_5
 
     call phys_getopts( microp_scheme_out      = microp_scheme, &
                        macrop_scheme_out      = macrop_scheme, &
@@ -2056,6 +2056,7 @@ subroutine tphysbc (ztodt,               &
                       ,l_rkz_lmt_2_out        = l_rkz_lmt_2 &
                       ,l_rkz_lmt_3_out        = l_rkz_lmt_3 &
                       ,l_rkz_lmt_4_out        = l_rkz_lmt_4 &
+                      ,l_rkz_lmt_5_out        = l_rkz_lmt_5 &
                       )
 
     !-----------------------------------------------------------------------
@@ -2488,7 +2489,8 @@ subroutine tphysbc (ztodt,               &
                                      rkz_term_C_fmin, &
                                      l_rkz_lmt_2, &
                                      l_rkz_lmt_3, &
-                                     l_rkz_lmt_4  &
+                                     l_rkz_lmt_4, &
+                                     l_rkz_lmt_5  &
                                      )
 
                call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol)          

components/cam/src/physics/cam/simple_condensation_model.F90

@@ -74,6 +74,9 @@ subroutine simple_RKZ_init()
     call addfld ('RKZ_zqme', (/'lev'/), 'I', 'kg/kg/s', 'condensation rate before limiters in the simple RKZ scheme')
     call addfld ('RKZ_qme',  (/'lev'/), 'I', 'kg/kg/s', 'condensation rate after limiters in the simple RKZ scheme')
 
+    call addfld ('RKZ_qme_lm4',  (/'lev'/), 'I', 'kg/kg/s', 'condensation rate difference before and after limiter 4 in the simple RKZ scheme')
+    call addfld ('RKZ_qme_lm5',  (/'lev'/), 'I', 'kg/kg/s', 'condensation rate difference before and after limiter 5 in the simple RKZ scheme')
+
   end subroutine simple_RKZ_init
 
   !------------------------------------------------------------------------
@@ -90,7 +93,8 @@ subroutine simple_RKZ_tend(state, ptend, tcwat, qcwat, lcwat, ast, qmeold, &
                              rkz_term_C_fmin,   & 
                              l_rkz_lmt_2, &
                              l_rkz_lmt_3, &
-                             l_rkz_lmt_4  &
+                             l_rkz_lmt_4, &
+                             l_rkz_lmt_5  &
                              )
 
   use shr_kind_mod, only: r8=>shr_kind_r8
@@ -129,6 +133,7 @@ subroutine simple_RKZ_tend(state, ptend, tcwat, qcwat, lcwat, ast, qmeold, &
   logical,  intent(in) :: l_rkz_lmt_2
   logical,  intent(in) :: l_rkz_lmt_3
   logical,  intent(in) :: l_rkz_lmt_4
+  logical,  intent(in) :: l_rkz_lmt_5
 
   ! tmp work arrays
 
@@ -141,32 +146,35 @@ subroutine simple_RKZ_tend(state, ptend, tcwat, qcwat, lcwat, ast, qmeold, &
   logical  :: lq(pcnst)              ! logical array used when calling subroutine physics_ptend_init indicating 
                                      ! which tracers are affected by this parameterization
 
-  real(r8) ::     ast_old(pcols,pver)! cloud fraction of previous time step
+  real(r8) :: ast_old(pcols,pver)! cloud fraction of previous time step
 
-  real(r8) ::        qsat(pcols,pver)! saturation specific humidity
-  real(r8) ::         esl(pcols,pver)! saturation vapor pressure (output from subroutine qsat_water, not used)
-  real(r8) ::     dqsatdT(pcols,pver)! dqsat/dT
-  real(r8) ::         gam(pcols,pver)! L/cpair * dqsat/dT
+  real(r8) :: qsat(pcols,pver)! saturation specific humidity
+  real(r8) :: esl(pcols,pver)! saturation vapor pressure (output from subroutine qsat_water, not used)
+  real(r8) :: dqsatdT(pcols,pver)! dqsat/dT
+  real(r8) :: gam(pcols,pver)! L/cpair * dqsat/dT
 
   real(r8) :: rhu00                  ! threshold grid-box-mean RH used in the diagnostic cldfrc scheme
   real(r8) :: rhgbm(pcols,pver)      ! grid box mean relative humidity
-  real(r8) ::     dastdRH(pcols,pver)! df/dRH where f is the cloud fraction and RH the relative humidity
+  real(r8) :: dastdRH(pcols,pver)    ! df/dRH where f is the cloud fraction and RH the relative humidity
 
   real(r8) :: qtend(pcols,pver)      ! Moisture tendency caused by other processes
   real(r8) :: ltend(pcols,pver)      ! liquid condensate tendency caused by other processes
   real(r8) :: ttend(pcols,pver)      ! Temperature tendency caused by other processes
 
   real(r8) :: qme   (pcols,pver)     ! total condensation rate
 
+  real(r8) :: qmebf   (pcols,pver)   ! total condensation rate before appling limiter 
+  real(r8) :: qmedf   (pcols,pver)   ! difference of total condensation rate before and after limiter application
+
   real(r8) :: term_A (pcols,pver)
   real(r8) :: term_B (pcols,pver)
   real(r8) :: term_C (pcols,pver)
   real(r8) :: rdenom (pcols,pver)    ! 1/denominator in the final expression of total grid-box mean condensation rate 
                                      ! when option 2 is used for term C. Note that in this case, 1/denominator will be 
                                      ! multipled to all three terms (A, B, and C).
 
-  real(r8) ::  ql_incld(pcols,pver)  ! in-cloud liquid concentration
-  real(r8) ::  dfdt    (pcols,pver)  ! df/dt where f is the cloud fraction
+  real(r8) :: ql_incld(pcols,pver)  ! in-cloud liquid concentration
+  real(r8) :: dfdt    (pcols,pver)  ! df/dt where f is the cloud fraction
   real(r8) :: zforcing (pcols,pver)
   real(r8) :: zc3      (pcols,pver)
 
@@ -175,6 +183,7 @@ subroutine simple_RKZ_tend(state, ptend, tcwat, qcwat, lcwat, ast, qmeold, &
   real(r8) :: zlim (pcols,pver)
   real(r8),parameter :: zsmall = 1e-12_r8
 
+
   !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
   rdtime = 1._r8/dtime
   ncol  = state%ncol
@@ -443,6 +452,8 @@ subroutine simple_RKZ_tend(state, ptend, tcwat, qcwat, lcwat, ast, qmeold, &
      !---------------------------------------------------------------------------
      if (l_rkz_lmt_4) then
 
+        qmebf(:ncol,:pver)  = qme(:ncol,:pver)
+
         ! Avoid negative qv
 
         zqvnew(:ncol,:pver) = state%q(:ncol,:pver,1) - dtime*qme(:ncol,:pver)
@@ -457,6 +468,44 @@ subroutine simple_RKZ_tend(state, ptend, tcwat, qcwat, lcwat, ast, qmeold, &
            qme(:ncol,:pver) = ( zsmall - state%q(:ncol,:pver,ixcldliq) )*rdtime
         end where
 
+        qmedf(:ncol,:pver)  = qme(:ncol,:pver) - qmebf(:ncol,:pver)
+        call outfld('RKZ_qme_lm4',   qmedf,    pcols, lchnk)
+
+     end if
+
+     !---------------------------------------------------------------------------
+     ! Limit condensation/evaporation rate to avoid negative water
+     ! concentrations. Also, let qme=0 when f=0 and df/dt=0.
+     ! (say, without any changes in cloud fraction and its tendency, the
+     ! condensation does not play any roles, we think qme should be zero, thus,
+     ! we add this limiter to avoid nonzero qme with f=0 and df/dt=0) . 
+     !---------------------------------------------------------------------------
+     if (l_rkz_lmt_5) then
+
+        qmebf(:ncol,:pver)  = qme(:ncol,:pver)
+
+        ! Avoid negative qv
+
+        zqvnew(:ncol,:pver) = state%q(:ncol,:pver,1) - dtime*qme(:ncol,:pver)
+        where( zqvnew(:ncol,:pver).lt.zsmall )
+           qme(:ncol,:pver) = ( state%q(:ncol,:pver,1) - zsmall )*rdtime
+        end where
+
+        ! Avoid negative ql (note that qme could be negative, which would mean
+        ! evaporation)
+
+        zqlnew(:ncol,:pver) = state%q(:ncol,:pver,ixcldliq) + dtime*qme(:ncol,:pver)
+        where( zqlnew(:ncol,:pver).lt.zsmall )
+           qme(:ncol,:pver) = ( zsmall - state%q(:ncol,:pver,ixcldliq) )*rdtime
+        end where
+
+        where ((ast(:ncol,:pver).eq.0._r8).and.(dfdt(:ncol,:pver).eq.0._r8))
+          qme(:ncol,:pver) = 0._r8
+        end where
+
+        qmedf(:ncol,:pver)  = qme(:ncol,:pver) - qmebf(:ncol,:pver)
+        call outfld('RKZ_qme_lm5',   qmedf,    pcols, lchnk)
+
      end if
 
      ! Send limited condensation rate to output