From 694a1809708acd5ab93a8cf4b287f771bc7c4f0c Mon Sep 17 00:00:00 2001 From: Shixuan Zhang Date: Tue, 20 Feb 2018 20:48:45 -0800 Subject: [PATCH] add "limiter #5 for simple condensation model" let qme=0 when f = 0 and df/dt = 0 To be tested. --- .../cam/src/physics/cam/phys_control.F90 | 9 ++- components/cam/src/physics/cam/physpkg.F90 | 6 +- .../physics/cam/simple_condensation_model.F90 | 67 ++++++++++++++++--- 3 files changed, 69 insertions(+), 13 deletions(-) diff --git a/components/cam/src/physics/cam/phys_control.F90 b/components/cam/src/physics/cam/phys_control.F90 index 192f14b16e75..aa4c3750d687 100644 --- a/components/cam/src/physics/cam/phys_control.F90 +++ b/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 diff --git a/components/cam/src/physics/cam/physpkg.F90 b/components/cam/src/physics/cam/physpkg.F90 index 587c01f392c2..78d660a96e29 100644 --- a/components/cam/src/physics/cam/physpkg.F90 +++ b/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) diff --git a/components/cam/src/physics/cam/simple_condensation_model.F90 b/components/cam/src/physics/cam/simple_condensation_model.F90 index 929de032e1f3..5511ed0a174b 100644 --- a/components/cam/src/physics/cam/simple_condensation_model.F90 +++ b/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,16 +146,16 @@ 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 @@ -158,6 +163,9 @@ subroutine simple_RKZ_tend(state, ptend, tcwat, qcwat, lcwat, ast, qmeold, & 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) @@ -165,8 +173,8 @@ subroutine simple_RKZ_tend(state, ptend, tcwat, qcwat, lcwat, ast, qmeold, & ! 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