diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_common.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_common.F90 index 1ae654edf..d50713fa1 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_common.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_common.F90 @@ -14,7 +14,7 @@ module GFS_PBL_generic_common subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & imp_physics_thompson, ltaerosol,mraerosol, & imp_physics_mg, ntgl, imp_physics_gfdl, & - imp_physics_zhao_carr, imp_physics_nssl,& + imp_physics_nssl, & nssl_hail_on, nssl_ccn_on, nssl_3moment, kk, & errmsg, errflg) implicit none @@ -22,7 +22,7 @@ subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & integer, intent(in ) :: imp_physics, imp_physics_wsm6, & imp_physics_thompson, & imp_physics_mg, ntgl, imp_physics_gfdl, & - imp_physics_zhao_carr,imp_physics_nssl + imp_physics_nssl logical, intent(in ) :: ltaerosol, mraerosol, nssl_hail_on, nssl_ccn_on, nssl_3moment integer, intent(out) :: kk character(len=*), intent(out) :: errmsg @@ -53,9 +53,6 @@ subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & elseif (imp_physics == imp_physics_gfdl) then ! GFDL MP kk = 7 - elseif (imp_physics == imp_physics_zhao_carr) then -! Zhao/Carr/Sundqvist - kk = 3 elseif (imp_physics == imp_physics_nssl) then IF ( nssl_hail_on ) THEN kk = 16 diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_post.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_post.F90 index 01033f4d6..1e0c66ff4 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_post.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_post.F90 @@ -11,7 +11,7 @@ module GFS_PBL_generic_post subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, & ntqv, ntcw, ntiw, ntrw, ntsw, ntlnc, ntinc, ntrnc, ntsnc, ntgnc, ntwa, ntia, ntgl, ntoz, ntke, ntkev,nqrimef, & trans_aero, ntchs, ntchm, ntccn, nthl, nthnc, ntgv, nthv, ntrz, ntgz, nthz, & - imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, imp_physics_mg, & + imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, imp_physics_mg, & imp_physics_fer_hires, imp_physics_nssl, nssl_ccn_on, ltaerosol, mraerosol, nssl_hail_on, nssl_3moment, & cplflx, cplaqm, cplchm, lssav, flag_for_pbl_generic_tend, ldiag3d, lsidea, hybedmf, do_shoc, satmedmf, & shinhong, do_ysu, dvdftra, dusfc1, dvsfc1, dtsfc1, dqsfc1, dtf, dudt, dvdt, dtdt, htrsw, htrlw, xmu, & @@ -33,7 +33,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, integer, intent(in) :: ntccn, nthl, nthnc, ntgv, nthv, ntrz, ntgz, nthz logical, intent(in) :: trans_aero integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6 - integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires + integer, intent(in) :: imp_physics_mg, imp_physics_fer_hires integer, intent(in) :: imp_physics_nssl logical, intent(in) :: nssl_ccn_on, nssl_hail_on, nssl_3moment logical, intent(in) :: ltaerosol, cplflx, cplaqm, cplchm, lssav, ldiag3d, lsidea, use_med_flux, mraerosol @@ -109,7 +109,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, call set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & imp_physics_thompson, ltaerosol,mraerosol, & imp_physics_mg, ntgl, imp_physics_gfdl, & - imp_physics_zhao_carr, imp_physics_nssl,& + imp_physics_nssl,& nssl_hail_on, nssl_ccn_on, nssl_3moment, kk, & errmsg, errflg) if (errflg /= 0) return @@ -244,14 +244,6 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, dqdt(i,k,ntoz) = dvdftra(i,k,7) enddo enddo - elseif (imp_physics == imp_physics_zhao_carr) then - do k=1,levs - do i=1,im - dqdt(i,k,1) = dvdftra(i,k,1) - dqdt(i,k,ntcw) = dvdftra(i,k,2) - dqdt(i,k,ntoz) = dvdftra(i,k,3) - enddo - enddo elseif (imp_physics == imp_physics_nssl ) then ! nssl IF ( nssl_hail_on ) THEN diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_post.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_post.meta index 057d061a4..ea1a1b1d8 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_post.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_post.meta @@ -260,13 +260,6 @@ dimensions = () type = integer intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in [imp_physics_mg] standard_name = identifier_for_morrison_gettelman_microphysics_scheme long_name = choice of Morrison-Gettelman microphysics scheme diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.F90 index eab767147..027e0dd4c 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.F90 @@ -14,7 +14,7 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, ntwa, ntia, ntgl, ntoz, ntke, ntkev, nqrimef, trans_aero, ntchs, ntchm, & ntccn, nthl, nthnc, ntgv, nthv, ntrz, ntgz, nthz, & imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, & - imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires, imp_physics_nssl, & + imp_physics_mg, imp_physics_fer_hires, imp_physics_nssl, & ltaerosol, mraerosol, nssl_ccn_on, nssl_hail_on, nssl_3moment, & hybedmf, do_shoc, satmedmf, qgrs, vdftra, save_u, save_v, save_t, save_q, & flag_for_pbl_generic_tend, ldiag3d, qdiag3d, lssav, ugrs, vgrs, tgrs, errmsg, errflg) @@ -32,7 +32,7 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, integer, intent(in) :: ntccn, nthl, nthnc, ntgv, nthv, ntrz, ntgz, nthz logical, intent(in) :: trans_aero, ldiag3d, qdiag3d, lssav integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6 - integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires + integer, intent(in) :: imp_physics_mg, imp_physics_fer_hires logical, intent(in) :: ltaerosol, hybedmf, do_shoc, satmedmf, flag_for_pbl_generic_tend, mraerosol integer, intent(in) :: imp_physics_nssl logical, intent(in) :: nssl_hail_on, nssl_ccn_on, nssl_3moment @@ -191,16 +191,6 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, enddo enddo rtg_ozone_index = 7 - elseif (imp_physics == imp_physics_zhao_carr) then -! Zhao/Carr/Sundqvist - do k=1,levs - do i=1,im - vdftra(i,k,1) = qgrs(i,k,ntqv) - vdftra(i,k,2) = qgrs(i,k,ntcw) - vdftra(i,k,3) = qgrs(i,k,ntoz) - enddo - enddo - rtg_ozone_index = 3 elseif (imp_physics == imp_physics_nssl ) then ! nssl IF ( nssl_hail_on ) THEN @@ -275,7 +265,7 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, call set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & imp_physics_thompson, ltaerosol,mraerosol, & imp_physics_mg, ntgl, imp_physics_gfdl, & - imp_physics_zhao_carr, imp_physics_nssl,& + imp_physics_nssl, & nssl_hail_on, nssl_ccn_on, nssl_3moment, kk, & errmsg, errflg) if (errflg /= 0) return diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.meta index 7a8e72bba..b605e9754 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.meta @@ -266,13 +266,6 @@ dimensions = () type = integer intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in [imp_physics_mg] standard_name = identifier_for_morrison_gettelman_microphysics_scheme long_name = choice of Morrison-Gettelman microphysics scheme diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.fv3.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.fv3.F90 index cbd660414..851d0180d 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.fv3.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.fv3.F90 @@ -17,7 +17,7 @@ module GFS_rad_time_vary !! subroutine GFS_rad_time_vary_timestep_init (lrseeds, rseeds, & lslwr, lsswr, isubc_lw, isubc_sw, icsdsw, icsdlw, cnx, cny, isc, jsc, & - imap, jmap, sec, kdt, imp_physics, imp_physics_zhao_carr, ipsd0, ipsdlim,& + imap, jmap, sec, kdt, imp_physics, ipsd0, ipsdlim, & ps_2delt, ps_1delt, t_2delt, t_1delt, qv_2delt, qv_1delt, t, qv, ps, & errmsg, errflg) @@ -31,7 +31,7 @@ subroutine GFS_rad_time_vary_timestep_init (lrseeds, rseeds, logical, intent(in) :: lrseeds integer, intent(in), optional :: rseeds(:,:) integer, intent(in) :: isubc_lw, isubc_sw, cnx, cny, isc, jsc, kdt - integer, intent(in) :: imp_physics, imp_physics_zhao_carr, ipsd0, ipsdlim + integer, intent(in) :: imp_physics, ipsd0, ipsdlim logical, intent(in) :: lslwr, lsswr integer, intent(inout), optional :: icsdsw(:), icsdlw(:) integer, intent(in) :: imap(:), jmap(:) @@ -82,17 +82,6 @@ subroutine GFS_rad_time_vary_timestep_init (lrseeds, rseeds, end if !lrseeds endif ! isubc_lw and isubc_sw - if (imp_physics == imp_physics_zhao_carr) then - if (kdt == 1) then - t_2delt = t - t_1delt = t - qv_2delt = qv - qv_1delt = qv - ps_2delt = ps - ps_1delt = ps - endif - endif - endif end subroutine GFS_rad_time_vary_timestep_init diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.fv3.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.fv3.meta index a7ac1381c..0a03e1592 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.fv3.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.fv3.meta @@ -131,13 +131,6 @@ dimensions = () type = integer intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in [ipsd0] standard_name = initial_seed_for_mcica long_name = initial permutation seed for mcica radiation diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.scm.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.scm.F90 index 46585c9da..9c84bf994 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.scm.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.scm.F90 @@ -17,7 +17,7 @@ module GFS_rad_time_vary !! subroutine GFS_rad_time_vary_timestep_init (lrseeds, rseeds, & lslwr, lsswr, isubc_lw, isubc_sw, icsdsw, icsdlw, cnx, cny, isc, jsc, & - imap, jmap, sec, kdt, imp_physics, imp_physics_zhao_carr, ipsd0, ipsdlim,& + imap, jmap, sec, kdt, imp_physics, ipsd0, ipsdlim, & ps_2delt, ps_1delt, t_2delt, t_1delt, qv_2delt, qv_1delt, t, qv, ps, & errmsg, errflg) @@ -31,7 +31,7 @@ subroutine GFS_rad_time_vary_timestep_init (lrseeds, rseeds, logical, intent(in) :: lrseeds integer, intent(in) :: rseeds(:,:) integer, intent(in) :: isubc_lw, isubc_sw, cnx, cny, isc, jsc, kdt - integer, intent(in) :: imp_physics, imp_physics_zhao_carr, ipsd0, ipsdlim + integer, intent(in) :: imp_physics, ipsd0, ipsdlim logical, intent(in) :: lslwr, lsswr integer, intent(inout) :: icsdsw(:), icsdlw(:) integer, intent(in) :: imap(:), jmap(:) @@ -82,17 +82,6 @@ subroutine GFS_rad_time_vary_timestep_init (lrseeds, rseeds, end if !lrseeds endif ! isubc_lw and isubc_sw - if (imp_physics == imp_physics_zhao_carr) then - if (kdt == 1) then - t_2delt = t - t_1delt = t - qv_2delt = qv - qv_1delt = qv - ps_2delt = ps - ps_1delt = ps - endif - endif - endif end subroutine GFS_rad_time_vary_timestep_init diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.scm.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.scm.meta index a7ac1381c..0a03e1592 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.scm.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rad_time_vary.scm.meta @@ -131,13 +131,6 @@ dimensions = () type = integer intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in [ipsd0] standard_name = initial_seed_for_mcica long_name = initial permutation seed for mcica radiation diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.F90 index 754fe12bb..426c5cf13 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.F90 @@ -27,8 +27,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& ntss3, ntss4, ntss5, ntsu, ntbcb, ntbcl, ntocb, ntocl, ntchm, & imp_physics,imp_physics_nssl, nssl_ccn_on, nssl_invertccn, & imp_physics_thompson, imp_physics_tempo, imp_physics_gfdl, & - imp_physics_zhao_carr, & - imp_physics_zhao_carr_pdf, imp_physics_mg, imp_physics_wsm6, & + imp_physics_mg, imp_physics_wsm6, & imp_physics_fer_hires, iovr, iovr_rand, iovr_maxrand, iovr_max, & iovr_dcorr, iovr_exp, iovr_exprand, idcor, idcor_con, idcor_hogan, & idcor_oreopoulos, dcorr_con, julian, yearlen, lndp_var_list, lsswr, & @@ -124,8 +123,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& imp_physics_thompson, & imp_physics_tempo, & imp_physics_gfdl, & - imp_physics_zhao_carr, & - imp_physics_zhao_carr_pdf, & imp_physics_mg, imp_physics_wsm6, & imp_physics_nssl, & imp_physics_fer_hires, & @@ -725,7 +722,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& ! if (ntcw > 0) then ! prognostic cloud schemes ccnd = 0.0_kind_phys - if (ncnd == 1) then ! Zhao_Carr_Sundqvist + if (ncnd == 1) then do k=1,LMK do i=1,IM ccnd(i,k,1) = tracer1(i,k,ntcw) ! liquid water/ice @@ -994,15 +991,13 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& ! --- add suspended convective cloud water to grid-scale cloud water ! only for cloud fraction & radiation computation ! it is to enhance cloudiness due to suspended convec cloud water -! for zhao/moorthi's (imp_phys=99) & -! ferrier's (imp_phys=5) microphysics schemes +! for ferrier's (imp_phys=5) microphysics schemes - if ((num_p3d == 4) .and. (npdf3d == 3)) then ! same as imp_physics = imp_physics_zhao_carr_pdf + if ((num_p3d == 4) .and. (npdf3d == 3)) then do k=1,lm k1 = k + kd do i=1,im !GJF: this is not consistent with GFS_typedefs, - ! but it looks like the Zhao-Carr-PDF scheme is not in the CCPP deltaq(i,k1) = 0.0!Tbd%phy_f3d(i,k,5) !GJF: this variable is not in phy_f3d anymore cnvw (i,k1) = cnvw_in(i,k) cnvc (i,k1) = cnvc_in(i,k) @@ -1027,10 +1022,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& enddo endif - if (imp_physics == imp_physics_zhao_carr) then - ccnd(1:IM,1:LMK,1) = ccnd(1:IM,1:LMK,1) + cnvw(1:IM,1:LMK) - endif - !> - Call radiation_clouds_prop() to calculate cloud properties. call radiation_clouds_prop & & ( plyr, plvl, tlyr, tvly, qlyr, qstl, rhly, & ! --- inputs: @@ -1041,7 +1032,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,& & imp_physics, imp_physics_nssl, imp_physics_fer_hires, & & imp_physics_gfdl, imp_physics_thompson, & & imp_physics_wsm6, imp_physics_tempo, & - & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, & & imp_physics_mg, iovr, iovr_rand, iovr_maxrand, iovr_max, & & iovr_dcorr, iovr_exp, iovr_exprand, idcor, idcor_con, & & idcor_hogan, idcor_oreopoulos, lcrick, lcnorm, & diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.meta index 697aadfb5..aa19a4c06 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmg_pre.meta @@ -491,20 +491,6 @@ dimensions = () type = integer intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in -[imp_physics_zhao_carr_pdf] - standard_name = identifier_for_zhao_carr_pdf_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme with PDF clouds - units = flag - dimensions = () - type = integer - intent = in [imp_physics_mg] standard_name = identifier_for_morrison_gettelman_microphysics_scheme long_name = choice of Morrison-Gettelman microphysics scheme diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 index a335f56a4..ae035b5b9 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 @@ -5,7 +5,7 @@ module GFS_rrtmgp_cloud_mp use machine, only: kind_phys use radiation_tools, only: check_error_msg - use module_radiation_clouds, only: progcld_thompson + use module_radiation_clouds, only: progcld_thompson, cld_frac_XuRandall use rrtmgp_lw_cloud_optics, only: & radliq_lwr => radliq_lwrLW, radliq_upr => radliq_uprLW,& radice_lwr => radice_lwrLW, radice_upr => radice_uprLW @@ -371,9 +371,9 @@ subroutine cloud_mp_GF(nCol, nLev, lsmask, t_lay, p_lev, p_lay, qs_lay, relhum, !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b) if(qi > 1.E-8) cld_cnv_reice(iCol,iLay) = max(173.45 + 2.14*(t_lay(iCol,iLay)-273.15), 20.) - ! Xu-Randall (1996) cloud-fraction. - cld_cnv_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay), & - qs_lay(iCol,iLay), relhum(iCol,iLay), qc+qi, alpha0) + ! Xu-Randall (1996) cloud-fraction, lambda = 0.5 + cld_cnv_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay)*0.01, & + qs_lay(iCol,iLay), relhum(iCol,iLay), qc+qi, alpha0, 0.5, 1.0) endif enddo enddo @@ -503,9 +503,10 @@ subroutine cloud_mp_SAMF(nCol, nLev, t_lay, p_lev, p_lay, qs_lay, relhum, cld_cnv_reliq(iCol,iLay) = reliq_def cld_cnv_reice(iCol,iLay) = reice_def - ! Xu-Randall (1996) cloud-fraction. - cld_cnv_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay), & - qs_lay(iCol,iLay), relhum(iCol,iLay), cnv_mixratio(iCol,iLay), alpha0) + ! Xu-Randall (1996) cloud-fraction, lambda = 0.5 + cld_cnv_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay)*0.01, & + qs_lay(iCol,iLay), relhum(iCol,iLay), cnv_mixratio(iCol,iLay), & + alpha0, 0.5, 1.0) endif enddo enddo @@ -738,10 +739,12 @@ subroutine cloud_mp_thompson(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldice, i_c cld_swp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,4) * tem1 * deltaP) ! Xu-Randall (1996) cloud-fraction. **Additionally, Conditioned on relative-humidity** + ! lambda = 0.5 cld_mr = cld_condensate(iCol,iLay,1) + cld_condensate(iCol,iLay,2) + & cld_condensate(iCol,iLay,3) + cld_condensate(iCol,iLay,4) - cld_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay), & - qs_lay(iCol,iLay), relhum(iCol,iLay), cld_mr, alpha0, cond_cfrac_onRH) + cld_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay)*0.01, & + qs_lay(iCol,iLay), relhum(iCol,iLay), cld_mr, alpha0, & + 0.5, 1.0, cond_cfrac_onRH) enddo enddo @@ -769,52 +772,6 @@ subroutine cloud_mp_thompson(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldice, i_c end subroutine cloud_mp_thompson -!> This function computes the cloud-fraction following -!! Xu-Randall(1996) \cite xu_and_randall_1996 -!! - function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha, cond_cfrac_onRH) - implicit none - ! Inputs - logical, intent(in), optional :: & - cond_cfrac_onRH ! If true, cloud-fracion set to unity when rh>99% - - real(kind_phys), intent(in) :: & - p_lay, & !< Pressure (Pa) - qs_lay, & !< Saturation vapor-pressure (Pa) - relhum, & !< Relative humidity - cld_mr, & !< Total cloud mixing ratio - alpha !< Scheme parameter (default=100) - - ! Outputs - real(kind_phys) :: cld_frac_XuRandall - - ! Locals - real(kind_phys) :: clwt, clwm, onemrh, tem1, tem2, tem3 - - ! Parameters - real(kind_phys) :: & - lambda = 0.50, & ! - P = 0.25 - - clwt = 1.0e-8 * (p_lay*0.001) - if (cld_mr > clwt) then - if(present(cond_cfrac_onRH) .and. relhum > 0.99) then - cld_frac_XuRandall = 1. - else - onemrh = max(1.e-10, 1.0 - relhum) - tem1 = alpha / min(max((onemrh*qs_lay)**lambda,0.0001),1.0) - tem2 = max(min(tem1*(cld_mr - clwt), 50.0 ), 0.0 ) - tem3 = sqrt(sqrt(relhum)) ! This assumes "p" = 0.25. Identical, but cheaper than relhum**p - ! - cld_frac_XuRandall = max( tem3*(1.0-exp(-tem2)), 0.0 ) - endif - else - cld_frac_XuRandall = 0.0 - endif - - return - end function - !> This routine is a wrapper to update the Thompson effective particle sizes used by the !! RRTMGP radiation scheme. subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc, & diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.F90 index df3d69bdd..2ba1029f3 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.F90 @@ -35,9 +35,9 @@ module GFS_rrtmgp_setup !! \htmlinclude GFS_rrtmgp_setup_init.html !! subroutine GFS_rrtmgp_setup_init(do_RRTMGP, imp_physics, imp_physics_fer_hires, & - imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, & - imp_physics_zhao_carr_pdf, imp_physics_mg, si, levr, ictm, isol, ico2, iaer, & - ntcw, ntoz, iovr, isubc_sw, isubc_lw, lalw1bd, idate, & + imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, & + imp_physics_mg, si, levr, ictm, isol, ico2, iaer, & + ntcw, ntoz, iovr, isubc_sw, isubc_lw, lalw1bd, idate, & me, aeros_file, iaermdl, iaerflg, con_pi, con_t0c, con_c, con_boltz, con_plnk, & solar_file, con_solr_2008, con_solr_2002, co2usr_file, co2cyc_file, ipsd0, & errmsg, errflg) @@ -50,8 +50,6 @@ subroutine GFS_rrtmgp_setup_init(do_RRTMGP, imp_physics, imp_physics_fer_hires, imp_physics_gfdl, & !< Flag for gfdl scheme imp_physics_thompson, & !< Flag for thompsonscheme imp_physics_wsm6, & !< Flag for wsm6 scheme - imp_physics_zhao_carr, & !< Flag for zhao-carr scheme - imp_physics_zhao_carr_pdf, & !< Flag for zhao-carr+PDF scheme imp_physics_mg !< Flag for MG scheme real(kind_phys), intent(in) :: & con_pi, con_t0c, con_c, con_boltz, con_plnk, con_solr_2008, con_solr_2002 diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.meta index 763790cb9..d9fe8359e 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.meta @@ -53,20 +53,6 @@ dimensions = () type = integer intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in -[imp_physics_zhao_carr_pdf] - standard_name = identifier_for_zhao_carr_pdf_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme with PDF clouds - units = flag - dimensions = () - type = integer - intent = in [imp_physics_mg] standard_name = identifier_for_morrison_gettelman_microphysics_scheme long_name = choice of Morrison-Gettelman microphysics scheme diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.F90 index b3d59c095..e7749e8c4 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.F90 @@ -17,7 +17,6 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, & xlon, xlat, gt0, gq0, sigmain,sigmaout,qmicro, & omegain,omegaout,imp_physics, imp_physics_mg, & - imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, & imp_physics_gfdl, imp_physics_thompson, dtidx, ntlnc, & imp_physics_wsm6, imp_physics_fer_hires, prsi, ntinc, & imp_physics_nssl, imp_physics_tempo, & @@ -33,7 +32,7 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & ! interface variables logical, intent(in) :: otsptflag(:)! on/off switch for tracer transport (size ntrac) integer, intent(in ) :: im, levs, nn, ntrac, ntcw, ntiw, ntclamt, ntrw, ntsw,& - ntrnc, ntsnc, ntgl, ntgnc, imp_physics, imp_physics_mg, imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, & + ntrnc, ntsnc, ntgl, ntgnc, imp_physics, imp_physics_mg, & imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6,imp_physics_fer_hires, & imp_physics_nssl, imp_physics_tempo, me, index_of_process_conv_trans integer, intent(in ), dimension(:) :: islmsk, kpbl, kinver @@ -57,7 +56,6 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & real(kind=kind_phys), intent(inout ), dimension(:,:), optional :: sigmain, omegain real(kind=kind_phys), intent(inout ), dimension(:,:), optional :: sigmaout, qmicro, omegaout real(kind=kind_phys), intent(inout), dimension(:,:) :: rhc, save_qc - ! save_qi is not allocated for Zhao-Carr MP real(kind=kind_phys), intent(inout), dimension(:,:) :: save_qi real(kind=kind_phys), intent(inout), dimension(:,:) :: save_tcp real(kind=kind_phys), intent(inout), dimension(:,:,:) :: clw @@ -192,19 +190,7 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & rhc(:,:) = 1.0 endif - if (imp_physics == imp_physics_zhao_carr .or. imp_physics == imp_physics_zhao_carr_pdf) then ! zhao-carr microphysics - !GF* move to GFS_MP_generic_pre (from gscond/precpd) - ! do i=1,im - ! psautco_l(i) = Model%psautco(1)*work1(i) + Model%psautco(2)*work2(i) - ! prautco_l(i) = Model%prautco(1)*work1(i) + Model%prautco(2)*work2(i) - ! enddo - !*GF - do k=1,levs - do i=1,im - clw(i,k,1) = gq0(i,k,ntcw) - enddo - enddo - elseif (imp_physics == imp_physics_gfdl) then + if (imp_physics == imp_physics_gfdl) then clw(1:im,:,1) = gq0(1:im,:,ntcw) elseif (imp_physics == imp_physics_thompson .or. & imp_physics == imp_physics_tempo) then diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta index 4cf339e4d..72c4daf4c 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_3.meta @@ -302,20 +302,6 @@ dimensions = () type = integer intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in -[imp_physics_zhao_carr_pdf] - standard_name = identifier_for_zhao_carr_pdf_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme with PDF clouds - units = flag - dimensions = () - type = integer - intent = in [imp_physics_gfdl] standard_name = identifier_for_gfdl_microphysics_scheme long_name = choice of GFDL microphysics scheme diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.F90 index f9a2b76ea..7ca706a81 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.F90 @@ -11,7 +11,7 @@ module GFS_suite_interstitial_4 subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntrac, ntcw, ntiw, ntclamt, & ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, ntccn, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & imp_physics_nssl, imp_physics_tempo, nssl_invertccn, nssl_ccn_on, & - imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, convert_dry_rho, dtf, save_qc, save_qi, con_pi, dtidx, dtend,& + convert_dry_rho, dtf, save_qc, save_qi, con_pi, dtidx, dtend, & index_of_process_conv_trans, gq0, clw, prsl, save_tcp, con_rd, con_eps, nssl_cccn, nwfa, spechum, ldiag3d, & qdiag3d, save_lnc, save_inc, ntk, ntke, otsptflag, errmsg, errflg) @@ -31,14 +31,13 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr logical, intent(in) :: otsptflag(:)! on/off switch for tracer transport by updraft and integer, intent(in ) :: im, levs, tracers_total, ntrac, ntcw, ntiw, ntclamt, ntrw, & ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, ntccn, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & - imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, imp_physics_nssl, imp_physics_tempo + imp_physics_nssl, imp_physics_tempo logical, intent(in) :: ltaerosol, convert_dry_rho logical, intent(in) :: nssl_ccn_on, nssl_invertccn real(kind=kind_phys), intent(in ) :: con_pi, dtf real(kind=kind_phys), intent(in ), dimension(:,:) :: save_qc - ! save_qi is not allocated for Zhao-Carr MP real(kind=kind_phys), intent(in ), dimension(:,:) :: save_qi, save_lnc, save_inc ! dtend and dtidx are only allocated if ldiag3d @@ -75,12 +74,10 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr errflg = 0 ! This code was previously in GFS_SCNV_generic_post, but it really belongs - ! here, because it fixes the convective transportable_tracers mess for Zhao-Carr - ! and GFDL MP from GFS_suite_interstitial_3. This whole code around clw(:,:,2) - ! being set to -999 for Zhao-Carr MP (which doesn't have cloud ice) and GFDL-MP - ! (which does have cloud ice, but for some reason it was decided to code it up - ! in the same way as for Zhao-Carr, nowadays unnecessary and confusing) needs - ! to be cleaned up. The convection schemes doing something different internally + ! here, because it fixes the convective transportable_tracers mess for + ! GFDL MP from GFS_suite_interstitial_3. This whole code around clw(:,:,2) + ! being set to -999 for GFDL-MP + ! The convection schemes doing something different internally ! based on clw(i,k,2) being -999.0 or not is not a good idea. do k=1,levs do i=1,im @@ -96,9 +93,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr endif endif if(ntcw>0) then - if (imp_physics == imp_physics_zhao_carr .or. & - imp_physics == imp_physics_zhao_carr_pdf .or. & - imp_physics == imp_physics_gfdl) then + if (imp_physics == imp_physics_gfdl) then idtend=dtidx(100+ntcw,index_of_process_conv_trans) if(idtend>=1) then dtend(:,:,idtend) = dtend(:,:,idtend) + clw(:,:,1)+clw(:,:,2) - gq0(:,:,ntcw) @@ -155,9 +150,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr if (ntcw > 0) then ! for microphysics - if (imp_physics == imp_physics_zhao_carr .or. & - imp_physics == imp_physics_zhao_carr_pdf .or. & - imp_physics == imp_physics_gfdl) then + if (imp_physics == imp_physics_gfdl) then gq0(1:im,:,ntcw) = clw(1:im,:,1) + clw(1:im,:,2) elseif (ntiw > 0) then diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.meta index 718b6ab95..681e550d1 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_suite_interstitial_4.meta @@ -165,20 +165,6 @@ dimensions = () type = integer intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in -[imp_physics_zhao_carr_pdf] - standard_name = identifier_for_zhao_carr_pdf_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme with PDF clouds - units = flag - dimensions = () - type = integer - intent = in [convert_dry_rho] standard_name = flag_for_converting_hydrometeors_from_moist_to_dry_air long_name = flag for converting hydrometeors from moist to dry air @@ -400,4 +386,4 @@ units = 1 dimensions = () type = integer - intent = out \ No newline at end of file + intent = out diff --git a/physics/MP/Zhao_Carr/zhaocarr_gscond.f b/physics/MP/Zhao_Carr/zhaocarr_gscond.f deleted file mode 100644 index 1d22c09ac..000000000 --- a/physics/MP/Zhao_Carr/zhaocarr_gscond.f +++ /dev/null @@ -1,530 +0,0 @@ -!> \file zhaocarr_gscond.f -!! This file contains the subroutine that calculates grid-scale -!! condensation and evaporation for use in Zhao and Carr (1997) -!! \cite zhao_and_carr_1997 scheme. - -!> This module contains the CCPP-compliant zhao_carr_gscond scheme. - module zhaocarr_gscond - - implicit none - public :: zhaocarr_gscond_init, zhaocarr_gscond_run - private - logical :: is_initialized = .False. - contains - - -! \brief Brief description of the subroutine -! -!> \section arg_table_zhaocarr_gscond_init Argument Table -!! - subroutine zhaocarr_gscond_init (imp_physics, & - & imp_physics_zhao_carr, & - & imp_physics_zhao_carr_pdf, & - & errmsg, errflg) - implicit none - - ! Interface variables - integer, intent(in ) :: imp_physics - integer, intent(in ) :: imp_physics_zhao_carr, & - & imp_physics_zhao_carr_pdf - - ! CCPP error handling - character(len=*), intent( out) :: errmsg - integer, intent( out) :: errflg - - ! Initialize the CCPP error handling variables - errmsg = '' - errflg = 0 - - if (is_initialized) return - - ! Consistency checks - if (imp_physics/=imp_physics_zhao_carr .and. & - & imp_physics/=imp_physics_zhao_carr_pdf) then - write(errmsg,'(*(a))') "Logic error: namelist choice of & - & microphysics is different from Zhao-Carr MP" - errflg = 1 - return - end if - - is_initialized = .true. - end subroutine zhaocarr_gscond_init - - -!> \defgroup condense GFS gscond Main -!> @{ -!! This subroutine computes grid-scale condensation and evaporation of -!! cloud condensate. -!! -#if 0 -!> \section arg_table_zhaocarr_gscond_run Argument Table -!! \htmlinclude zhaocarr_gscond_run.html -!! -#endif -!> \section general_gscond GFS gscond Scheme General Algorithm -!! -# Calculate ice-water identification number \f$IW\f$ in order to make a distinction between -!! cloud water and cloud ice (table2 of Zhao and Carr (1997) \cite zhao_and_carr_1997). -!! -# Calculate the changes in \f$t\f$, \f$q\f$ and \f$p\f$ due to all the processes except microphysics. -!! -# Calculate cloud evaporation rate (\f$E_c\f$, eq. 19 of Zhao and Carr (1997)\cite zhao_and_carr_1997). -!! -# Calculate cloud condensation rate (\f$C_g\f$, eq.8 of Zhao and Carr (1997)\cite zhao_and_carr_1997). -!! -# Update \f$t\f$, \f$q\f$, \f$cwm\f$ due to cloud evaporation and condensation processes. -!> \section Zhao-Carr_cond_detailed GFS gscond Scheme Detailed Algorithm -!> @{ - subroutine zhaocarr_gscond_run (im,km,dt,dtf,prsl,ps,q,clw1 & - &, clw2, cwm, t, tp, qp, psp & - &, psat,hvap,grav,hfus,ttp,rd,cp,eps,epsm1,rv & - &, tp1, qp1, psp1, u, lprnt, ipr, errmsg, errflg) - -! -! ****************************************************************** -! * * -! * subroutine for grid-scale condensation & evaporation * -! * for the mrf model at ncep. * -! * * -! ****************************************************************** -! * * -! * created by: q. zhao jan. 1995 * -! * modified by: h.-l. pan sep. 1998 * -! * modified by: s. moorthi aug. 1998, 1999, 2000 * -! * * -! * references: * -! * * -! ****************************************************************** -! - use machine , only : kind_phys - use funcphys , only : fpvs -! use namelist_def, only: nsdfi,fhdfi - implicit none -! -! Interface variables - integer, intent(in) :: im, km, ipr - real(kind=kind_phys), intent(in) :: dt, dtf - real(kind=kind_phys), intent(in) :: prsl(:,:), ps(:) - real(kind=kind_phys), intent(inout) :: q(:,:), t(:,:) - real(kind=kind_phys), intent(in) :: clw1(:,:), clw2(:,:) - real(kind=kind_phys), intent(out) :: cwm(:,:) - real(kind=kind_phys), intent(inout) :: & - &, tp(:,:), qp(:,:), psp(:) & - &, tp1(:,:), qp1(:,:), psp1(:) - real(kind=kind_phys), intent(in) :: u(:,:) - logical, intent(in) :: lprnt - real(kind=kind_phys), intent(in) :: psat, hvap, grav, hfus & - &, ttp, rd, cp, eps, epsm1, rv -! - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg -! -! Local variables - real (kind=kind_phys) h1, d00, elwv, eliv - &, epsq, r, cpr, rcp -! - parameter (h1=1.e0, d00=0.e0, epsq=2.e-12) -! - real(kind=kind_phys), parameter :: cons_0=0.0, cons_m15=-15.0 -! - real (kind=kind_phys) qi(im), qint(im), ccrik, e0 - &, cond, rdt, us, cclimit, climit - &, tmt0, tmt15, qik, cwmik - &, ai, qw, u00ik, tik, pres, pp0, fi - &, at, aq, ap, fiw, elv, qc, rqik - &, rqikk, tx1, tx2, tx3, es, qs - &, tsq, delq, condi, cone0, us00, ccrik1 - &, aa, ab, ac, ad, ae, af, ag - &, el2orc, albycp -! real (kind=kind_phys) vprs(im) - integer iw(im,km), i, k, iwik -! - ! Initialize CCPP error handling variables - errmsg = '' - errflg = 0 -! -!-----------------GFS interstitial in driver ---------------------------- - do i = 1,im - do k= 1,km - cwm(i,k) = clw1(i,k)+clw2(i,k) - enddo - enddo -!-----------------prepare constants for later uses----------------- -! - elwv = hvap - eliv = hvap+hfus - r = rd - cpr = cp*r - rcp = h1/cp - el2orc = hvap*hvap / (rv*cp) - albycp = hvap / cp -! - rdt = h1/dt - us = h1 - cclimit = 1.0e-3 - climit = 1.0e-20 -! - do i = 1, im - iw(i,km) = d00 - enddo -! -! check for first time step -! -! if (tp(1,1) < 1.) then -! do k = 1, km -! do i = 1, im -! tp(i,k) = t(i,k) -! qp(i,k) = max(q(i,k),epsq) -! tp1(i,k) = t(i,k) -! qp1(i,k) = max(q(i,k),epsq) -! enddo -! enddo -! do i = 1, im -! psp(i) = ps(i) -! psp1(i) = ps(i) -! enddo -! endif -! -!************************************************************* -!> -# Begining of grid-scale condensation/evaporation loop (start of -!! k-loop, i-loop) -!************************************************************* -! -! do k = km-1,2,-1 - do k = km,1,-1 -! vprs(:) = 0.001 * fpvs(t(:,k)) ! fpvs in pa -!----------------------------------------------------------------------- -!------------------qw, qi and qint-------------------------------------- - do i = 1, im - tmt0 = t(i,k)-273.16 - tmt15 = min(tmt0,cons_m15) - qik = max(q(i,k),epsq) - cwmik = max(cwm(i,k),climit) -! -! ai = 0.008855 -! bi = 1.0 -! if (tmt0 .lt. -20.0) then -! ai = 0.007225 -! bi = 0.9674 -! end if -! -! the global qsat computation is done in pa - pres = prsl(i,k) -! -! qw = vprs(i) - qw = min(pres, fpvs(t(i,k))) -! - qw = eps * qw / (pres + epsm1 * qw) - qw = max(qw,epsq) -! qi(i) = qw *(bi+ai*min(tmt0,cons_0)) -! qint(i) = qw *(1.-0.00032*tmt15*(tmt15+15.)) - qi(i) = qw - qint(i) = qw -! if (tmt0 .le. -40.) qint(i) = qi(i) - -!> -# Compute ice-water identification number IW. -!!\n The distinction between cloud water and cloud ice is made by the -!! cloud identification number IW, which is zero for cloud water and -!! unity for cloud ice (Table 2 in Zhao and Carr (1997) -!! \cite zhao_and_carr_1997): -!! - All clouds are defined to consist of liquid water below the -!! freezing level (\f$T\geq 0^oC\f$) and of ice particles above the -!! \f$T=-15^oC\f$ level. -!! - In the temperature region between \f$-15^oC\f$ and \f$0^oC\f$, -!! clouds may be composed of liquid water or ice. If there are cloud -!! ice particles above this point at the previous or current time step, -!! or if the cloud at this point at the previous time step consists of -!! ice particles, then the cloud substance at this point is considered -!! to be ice particles because of the cloud seeding effect and the -!! memory of its content. Otherwise, all clouds in this region are -!! considered to contain supercooled cloud water. - -!-------------------ice-water id number iw------------------------------ - if(tmt0.lt.-15.0) then - u00ik = u(i,k) - fi = qik - u00ik*qi(i) - if(fi > d00.or.cwmik > climit) then - iw(i,k) = 1 - else - iw(i,k) = 0 - end if - end if -! - if(tmt0.ge.0.0) then - iw(i,k) = 0 - end if -! - if (tmt0 < 0.0 .and. tmt0 >= -15.0) then - iw(i,k) = 0 - if (k < km) then - if (iw(i,k+1) == 1 .and. cwmik > climit) iw(i,k) = 1 - endif - end if - enddo -!> -# Condensation and evaporation of cloud -!--------------condensation and evaporation of cloud-------------------- - do i = 1, im -!> - Compute the changes in t, q and p (\f$A_{t}\f$,\f$A_{q}\f$ and -!! \f$A_{p}\f$) caused by all the processes except grid-scale -!! condensation and evaporation. -!!\f[ -!! A_{t}=(t-tp)/dt -!!\f] -!!\f[ -!! A_{q}=(q-qp)/dt -!!\f] -!!\f[ -!! A_{p}=(prsl-\frac{prsl}{ps} \times psp)/dt -!!\f] -!------------------------at, aq and dp/dt------------------------------- - qik = max(q(i,k),epsq) - cwmik = max(cwm(i,k),climit) - iwik = iw(i,k) - u00ik = u(i,k) - tik = t(i,k) - pres = prsl(i,k) - pp0 = (pres / ps(i)) * psp(i) - at = (tik-tp(i,k)) * rdt - aq = (qik-qp(i,k)) * rdt - ap = (pres-pp0) * rdt -!> - Calculate the saturation specific humidity \f$q_{s}\f$ and the -!! relative humidity \f$f\f$ using IW. -!----------------the satuation specific humidity------------------------ - fiw = float(iwik) - elv = (h1-fiw)*elwv + fiw*eliv - qc = (h1-fiw)*qint(i) + fiw*qi(i) -! if (lprnt) print *,' qc=',qc,' qint=',qint(i),' qi=',qi(i) -!----------------the relative humidity---------------------------------- - if(qc.le.1.0e-10) then - rqik=d00 - else - rqik = qik/qc - endif - -!> - According to Sundqvist et al. (1989) \cite sundqvist_et_al_1989, -!! estimate cloud fraction \f$b\f$ at a grid point from relative -!! humidity \f$f\f$ using the equation -!!\f[ -!! b=1-\left ( \frac{f_{s}-f}{f_{s}-u} \right )^{1/2} -!!\f] -!! for \f$f>u\f$; and \f$b=0\f$ for \f$f1.0\times10^{-3}\f$, condense water vapor -!! into cloud condensate (\f$C_{g}\f$). -!!\n Using \f$q=fq_{s}\f$, \f$q_{s}=\epsilon e_{s}/p\f$, and the -!! Clausius-Clapeyron equation \f$de_{s}/dT=\epsilon Le_{s}/RT^{2}\f$, -!! where \f$q_{s}\f$ is the saturation specific humidity,\f$e_{s}\f$ -!! is the saturation vapor pressure, \f$R\f$ is the specific gas -!! constant for dry air, \f$f\f$ is the relative humidity, and -!! \f$\epsilon=0.622\f$, the expression for \f$C_{g}\f$ has the form -!!\f[ -!! C_{g}=\frac{M-q_{s}f_{t}}{1+(f\epsilon L^{2}q_{s}/RC_{p}T^{2})}+E_{c} -!!\f] -!! where -!!\f[ -!! M=A_{q}-\frac{f\epsilon Lq_{s}}{RT^{2}}A_{t}+\frac{fq_{s}}{p}A_{p} -!!\f] -!! To close the system, an equation for the relative humidity tendency -!! \f$f_{t}\f$ was derived by Sundqvist et al.(1989) -!! \cite sundqvist_et_al_1989 using the hypothesis that the quantity -!! \f$M+E_{c}\f$ is divided into one part,\f$bM\f$,which condenses -!! in the already cloudy portion of a grid square, and another part, -!! \f$(1-b)M+E_{c}\f$,which is used to increase the relative humidity -!! of the cloud-free portion and the cloudiness in the square. The -!! equation is written as -!!\f[ -!! f_{t}=\frac{2(1-b)(f_{s}-u)[(1-b)M+E_{c}]}{2q_{s}(1-b)(f_{s}-u)+cwm/b} -!!\f] -!! - Check and correct if over condensation occurs. -!! - Update t, q and cwm (according to Eqs(6) and (7) in Zhao and Carr (1997) -!! \cite zhao_and_carr_1997) -!!\f[ -!! cwm=cwm+(C_{g}-E_{c})\times dt -!!\f] -!!\f[ -!! q=q-(C_{g}-E_{c})\times dt -!!\f] -!!\f[ -!! t=t+\frac{L}{C_{p}}(C_{g}-E_{c})\times dt -!!\f] -!!\n where \f$L\f$ is the latent heat of condensation/deposition, and -!! \f$C_{p}\f$ is the specific heat of air at constant pressure. - -!----------------cloud cover ratio ccrik-------------------------------- - if (rqik .lt. u00ik) then - ccrik = d00 - elseif(rqik.ge.us) then - ccrik = us - else - rqikk = min(us,rqik) - ccrik = h1-sqrt((us-rqikk)/(us-u00ik)) - endif -!-----------correct ccr if it is too small in large cwm regions-------- -! if(ccrik.ge.0.01.and.ccrik.le.0.2.and -! & .cwmik.ge.0.2e-3) then -! ccrik=min(1.0,cwmik*1.0e3) -! end if -!---------------------------------------------------------------------- -! if no cloud exists then evaporate any existing cloud condensate -!----------------evaporation of cloud water----------------------------- - e0 = d00 - if (ccrik <= cclimit.and. cwmik > climit) then -! -! first iteration - increment halved -! - tx1 = tik - tx3 = qik -! - es = min(pres, fpvs(tx1)) - qs = u00ik * eps * es / (pres + epsm1*es) - tsq = tx1 * tx1 - delq = 0.5 * (qs - tx3) * tsq / (tsq + el2orc * qs) -! - tx2 = delq - tx1 = tx1 - delq * albycp - tx3 = tx3 + delq -! -! second iteration -! - es = min(pres, fpvs(tx1)) - qs = u00ik * eps * es / (pres + epsm1*es) - tsq = tx1 * tx1 - delq = (qs - tx3) * tsq / (tsq + el2orc * qs) -! - tx2 = tx2 + delq - tx1 = tx1 - delq * albycp - tx3 = tx3 + delq -! -! third iteration -! - es = min(pres, fpvs(tx1)) - qs = u00ik * eps * es / (pres + epsm1*es) - tsq = tx1 * tx1 - delq = (qs - tx3) * tsq / (tsq + el2orc * qs) - tx2 = tx2 + delq -! - e0 = max(tx2*rdt, cons_0) -! if (lprnt .and. i .eq. ipr .and. k .eq. 34) -! & print *,' tx2=',tx2,' qc=',qc,' u00ik=',u00ik,' rqik=',rqik -! &,' cwmik=',cwmik,' e0',e0 - -! e0 = max(qc*(u00ik-rqik)*rdt, cons_0) - e0 = min(cwmik*rdt, e0) - e0 = max(cons_0,e0) - end if -! if cloud cover > 0.2 condense water vapor in to cloud condensate -!-----------the eqs. for cond. has been reorganized to reduce cpu------ - cond = d00 -! if (ccrik .gt. 0.20 .and. qc .gt. epsq) then - if (ccrik .gt. cclimit .and. qc .gt. epsq) then - us00 = us - u00ik - ccrik1 = 1.0 - ccrik - aa = eps*elv*pres*qik - ab = ccrik*ccrik1*qc*us00 - ac = ab + 0.5*cwmik - ad = ab * ccrik1 - ae = cpr*tik*tik - af = ae * pres - ag = aa * elv - ai = cp * aa - cond = (ac-ad)*(af*aq-ai*at+ae*qik*ap)/(ac*(af+ag)) -!-----------check & correct if over condensation occurs----------------- - condi = (qik -u00ik *qc*1.0)*rdt - cond = min(cond, condi) -!----------check & correct if supersatuation is too high---------------- -! qtemp=qik-max(0.,(cond-e0))*dt -! if(qc.le.1.0e-10) then -! rqtmp=0.0 -! else -! rqtmp=qtemp/qc -! end if -! if(rqtmp.ge.1.10) then -! cond=(qik-1.10*qc)*rdt -! end if -!----------------------------------------------------------------------- - cond = max(cond, d00) -!-------------------update of t, q and cwm------------------------------ - end if - cone0 = (cond-e0) * dt - cwm(i,k) = cwm(i,k) + cone0 -! if (lprnt .and. i .eq. ipr) print *,' t=',t(i,k),' cone0',cone0 -! &,' cond=',cond,' e0=',e0,' elv=',elv,' rcp=',rcp,' k=',k -! &,' cwm=',cwm(i,k) - t(i,k) = t(i,k) + elv*rcp*cone0 - q(i,k) = q(i,k) - cone0 - enddo ! end of i-loop! - enddo ! end of k-loop! -! -!********************************************************************* -!> -# End of the condensation/evaporation loop (end of i-loop,k-loop). -!********************************************************************* -! -!> -# Store \f$t\f$, \f$q\f$, \f$ps\f$ for next time step. - - if (dt > dtf+0.001) then ! three time level - do k = 1, km - do i = 1, im - tp(i,k) = tp1(i,k) - qp(i,k) = qp1(i,k) -! - tp1(i,k) = t(i,k) - qp1(i,k) = max(q(i,k),epsq) - enddo - enddo - do i = 1, im - psp(i) = psp1(i) - psp1(i) = ps(i) - enddo - else ! two time level scheme - tp1, qp1, psp1 not used - do k = 1, km -! write(0,*)' in gscond k=',k,' im=',im,' km=',km - do i = 1, im -! write(0,*)' in gscond i=',i - tp(i,k) = t(i,k) - qp(i,k) = max(q(i,k),epsq) -! qp(i,k) = q(i,k) - tp1(i,k) = tp(i,k) - qp1(i,k) = qp(i,k) - enddo - enddo - do i = 1, im - psp(i) = ps(i) - psp1(i) = ps(i) - enddo - endif -!----------------------------------------------------------------------- - return - end subroutine zhaocarr_gscond_run -!> @} -!> @} - end module zhaocarr_gscond diff --git a/physics/MP/Zhao_Carr/zhaocarr_gscond.meta b/physics/MP/Zhao_Carr/zhaocarr_gscond.meta deleted file mode 100644 index ed57ca909..000000000 --- a/physics/MP/Zhao_Carr/zhaocarr_gscond.meta +++ /dev/null @@ -1,300 +0,0 @@ -[ccpp-table-properties] - name = zhaocarr_gscond - type = scheme - dependencies = ../../tools/funcphys.f90,../../hooks/machine.F,../../hooks/physcons.F90 - -######################################################################## -[ccpp-arg-table] - name = zhaocarr_gscond_init - type = scheme -[imp_physics] - standard_name = control_for_microphysics_scheme - long_name = choice of microphysics scheme - units = flag - dimensions = () - type = integer - intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in -[imp_physics_zhao_carr_pdf] - standard_name = identifier_for_zhao_carr_pdf_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme with PDF clouds - units = flag - dimensions = () - type = integer - intent = in -[errmsg] - standard_name = ccpp_error_message - long_name = error message for error handling in CCPP - units = none - dimensions = () - type = character - kind = len=* - intent = out -[errflg] - standard_name = ccpp_error_code - long_name = error code for error handling in CCPP - units = 1 - dimensions = () - type = integer - intent = out -######################################################################## -[ccpp-arg-table] - name = zhaocarr_gscond_run - type = scheme -[im] - standard_name = horizontal_loop_extent - long_name = horizontal loop extent - units = count - dimensions = () - type = integer - intent = in -[km] - standard_name = vertical_layer_dimension - long_name = vertical layer dimension - units = count - dimensions = () - type = integer - intent = in -[dt] - standard_name = timestep_for_physics - long_name = physics time step - units = s - dimensions = () - type = real - kind = kind_phys - intent = in -[dtf] - standard_name = timestep_for_dynamics - long_name = dynamics time step - units = s - dimensions = () - type = real - kind = kind_phys - intent = in -[prsl] - standard_name = air_pressure - long_name = layer mean air pressure - units = Pa - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = in -[ps] - standard_name = surface_air_pressure - long_name = surface pressure - units = Pa - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in -[q] - standard_name = specific_humidity_of_new_state - long_name = water vapor specific humidity - units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout -[clw1] - standard_name = ice_water_mixing_ratio_convective_transport_tracer - long_name = ratio of mass of ice water to mass of dry air plus vapor (without condensates) in the convectively transported tracer array - units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = in -[clw2] - standard_name = cloud_condensed_water_mixing_ratio_convective_transport_tracer - long_name = ratio of mass of cloud water to mass of dry air plus vapor (without condensates) in the convectively transported tracer array - units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = in -[cwm] - standard_name = cloud_liquid_water_mixing_ratio_of_new_state - long_name = moist cloud condensed water mixing ratio - units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = out -[t] - standard_name = air_temperature_of_new_state - long_name = layer mean air temperature - units = K - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout -[tp] - standard_name = air_temperature_two_timesteps_back - long_name = air temperature two timesteps back - units = K - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout -[qp] - standard_name = specific_humidity_two_timesteps_back - long_name = water vapor specific humidity two timesteps back - units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout -[psp] - standard_name = surface_air_pressure_two_timesteps_back - long_name = surface air pressure two timesteps back - units = Pa - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = inout -[psat] - standard_name = saturation_pressure_at_triple_point_of_water - long_name = saturation pressure at triple point of water - units = Pa - dimensions = () - type = real - kind = kind_phys - intent = in -[hvap] - standard_name = latent_heat_of_vaporization_of_water_at_0C - long_name = latent heat of evaporation/sublimation - units = J kg-1 - dimensions = () - type = real - kind = kind_phys - intent = in -[grav] - standard_name = gravitational_acceleration - long_name = gravitational acceleration - units = m s-2 - dimensions = () - type = real - kind = kind_phys - intent = in -[hfus] - standard_name = latent_heat_of_fusion_of_water_at_0C - long_name = latent heat of fusion - units = J kg-1 - dimensions = () - type = real - kind = kind_phys - intent = in -[ttp] - standard_name = triple_point_temperature_of_water - long_name = triple point temperature of water - units = K - dimensions = () - type = real - kind = kind_phys - intent = in -[rd] - standard_name = gas_constant_of_dry_air - long_name = ideal gas constant for dry air - units = J kg-1 K-1 - dimensions = () - type = real - kind = kind_phys - intent = in -[cp] - standard_name = specific_heat_of_dry_air_at_constant_pressure - long_name = specific heat of dry air at constant pressure - units = J kg-1 K-1 - dimensions = () - type = real - kind = kind_phys - intent = in -[eps] - standard_name = ratio_of_dry_air_to_water_vapor_gas_constants - long_name = rd/rv - units = none - dimensions = () - type = real - kind = kind_phys - intent = in -[epsm1] - standard_name = ratio_of_dry_air_to_water_vapor_gas_constants_minus_one - long_name = (rd/rv) - 1 - units = none - dimensions = () - type = real - kind = kind_phys - intent = in -[rv] - standard_name = gas_constant_water_vapor - long_name = ideal gas constant for water vapor - units = J kg-1 K-1 - dimensions = () - type = real - kind = kind_phys - intent = in -[tp1] - standard_name = air_temperature_on_previous_timestep_in_xyz_dimensioned_restart_array - long_name = air temperature at previous timestep - units = K - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout -[qp1] - standard_name = specific_humidity_on_previous_timestep_in_xyz_dimensioned_restart_array - long_name = water vapor specific humidity at previous timestep - units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout -[psp1] - standard_name = surface_air_pressure_on_previous_timestep - long_name = surface air surface pressure at previous timestep - units = Pa - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = inout -[u] - standard_name = critical_relative_humidity - long_name = critical relative humidity - units = frac - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = in -[lprnt] - standard_name = flag_print - long_name = flag for printing diagnostics to output - units = flag - dimensions = () - type = logical - intent = in -[ipr] - standard_name = horizontal_index_of_printed_column - long_name = horizontal index of printed column - units = index - dimensions = () - type = integer - intent = in -[errmsg] - standard_name = ccpp_error_message - long_name = error message for error handling in CCPP - units = none - dimensions = () - type = character - kind = len=* - intent = out -[errflg] - standard_name = ccpp_error_code - long_name = error code for error handling in CCPP - units = 1 - dimensions = () - type = integer - intent = out diff --git a/physics/MP/Zhao_Carr/zhaocarr_precpd.f b/physics/MP/Zhao_Carr/zhaocarr_precpd.f deleted file mode 100644 index 922fedfcc..000000000 --- a/physics/MP/Zhao_Carr/zhaocarr_precpd.f +++ /dev/null @@ -1,741 +0,0 @@ -!> \file zhaocarr_precpd.f -!! This file contains the subroutine that calculates precipitation -!! processes from suspended cloud water/ice. - -!> This module contains the CCPP-compliant zhao_carr_precpd scheme. - module zhaocarr_precpd - - implicit none - public :: zhaocarr_precpd_init, zhaocarr_precpd_run - private - logical :: is_initialized = .False. - contains - - subroutine zhaocarr_precpd_init (imp_physics, & - & imp_physics_zhao_carr, & - & imp_physics_zhao_carr_pdf, & - & errmsg, errflg) - implicit none - - ! Interface variables - integer, intent(in ) :: imp_physics - integer, intent(in ) :: imp_physics_zhao_carr, & - & imp_physics_zhao_carr_pdf - ! CCPP error handling - character(len=*), intent( out) :: errmsg - integer, intent( out) :: errflg - - ! Initialize the CCPP error handling variables - errmsg = '' - errflg = 0 - - if (is_initialized) return - - ! Consistency checks - if (imp_physics/=imp_physics_zhao_carr .and. & - & imp_physics/=imp_physics_zhao_carr_pdf) then - write(errmsg,'(*(a))') "Logic error: namelist choice of & - & microphysics is different from Zhao-Carr MP" - errflg = 1 - return - end if - - is_initialized = .true. - end subroutine zhaocarr_precpd_init - -!> \defgroup precip GFS precpd Main -!! \brief This subroutine computes the conversion from condensation to -!! precipitation (snow or rain) or evaporation of rain. -!! -!! \section arg_table_zhaocarr_precpd_run Argument Table -!! \htmlinclude zhaocarr_precpd_run.html -!! -!> \section general_precpd GFS precpd Scheme General Algorithm -!! The following two equations can be used to calculate the -!! precipitation rates of rain and snow at each model level: -!!\f[ -!! P_{r}(\eta)=\frac{p_{s}-p_{t}}{g\eta_{s}}\int_{\eta}^{\eta_{t}}(P_{raut}+P_{racw}+P_{sacw}+P_{sm1}+P_{sm2}-E_{rr})d\eta -!! \f] -!! and -!! \f[ -!! P_{s}(\eta)=\frac{p_{s}-p_{t}}{g\eta_{s}}\int_{\eta}^{\eta_{t}}(P_{saut}+P_{saci}-P_{sm1}-P_{sm2}-E_{rs})d\eta -!! \f] -!! where \f$p_{s}\f$ and\f$p_{t}\f$ are the surface pressure and the -!! pressure at the top of model domain, respectively, and \f$g\f$ is -!! gravity. The implementation of the precipitation scheme also -!! includes a simplified procedure of computing \f$P_{r}\f$ -!! and \f$P_{s}\f$ (Zhao and Carr (1997) \cite zhao_and_carr_1997). -!! -!! The calculation is as follows: -!! -# Calculate precipitation production by auto conversion and accretion (\f$P_{saut}\f$, \f$P_{saci}\f$, \f$P_{raut}\f$). -!! - The accretion of cloud water by rain, \f$P_{racw}\f$, is not included in the current operational scheme. -!! -# Calculate evaporation of precipitation (\f$E_{rr}\f$ and \f$E_{rs}\f$). -!! -# Calculate melting of snow (\f$P_{sm1}\f$ and \f$P_{sm2}\f$, \f$P_{sacw}\f$). -!! -# Update t and q due to precipitation (snow or rain) production. -!! -# Calculate precipitation at surface (\f$rn\f$) and fraction of frozen precipitation (\f$sr\f$). -!! \section Zhao-Carr_precip_detailed GFS precpd Scheme Detailed Algorithm -!> @{ - subroutine zhaocarr_precpd_run (im,km,dt,del,prsl,q,cwm,t,rn & - &, grav, hvap, hfus, ttp, cp, eps, epsm1 & - &, sr,rainp,u00k,psautco,prautco,evpco,wminco & - &, wk1,lprnt,jpr,errmsg,errflg) - -! -! ****************************************************************** -! * * -! * subroutine for precipitation processes * -! * from suspended cloud water/ice * -! * * -! ****************************************************************** -! * * -! * originally created by q. zhao jan. 1995 * -! * ------- * -! * modified and rewritten by shrinivas moorthi oct. 1998 * -! * ----------------- * -! * and hua-lu pan * -! * ---------- * -! * * -! * references: * -! * * -! * zhao and carr (1997), monthly weather review (august) * -! * sundqvist et al., (1989) monthly weather review. (august) * -! * chuang 2013, modify sr to define frozen precipitation fraction* -! ****************************************************************** -! -! in this code vertical indexing runs from surface to top of the -! model -! -! argument list: -! -------------- -! im : inner dimension over which calculation is made -! km : number of vertical levels -! dt : time step in seconds -! del(km) : pressure layer thickness (bottom to top) -! prsl(km) : pressure values for model layers (bottom to top) -! q(im,km) : specific humidity (updated in the code) -! cwm(im,km) : condensate mixing ratio (updated in the code) -! t(im,km) : temperature (updated in the code) -! rn(im) : precipitation over one time-step dt (m/dt) -!old sr(im) : index (=-1 snow, =0 rain/snow, =1 rain) -!new sr(im) : "snow ratio", ratio of snow to total precipitation -! cll(im,km) : cloud cover -!hchuang rn(im) unit in m per time step -! precipitation rate conversion 1 mm/s = 1 kg/m2/s -! - use machine , only : kind_phys - use funcphys , only : fpvs - - implicit none -! include 'constant.h' -! -! Interface variables - integer, intent(in) :: im, km, jpr - real (kind=kind_phys), intent(in) :: grav, hvap, hfus, ttp, cp, & - & eps, epsm1 - real (kind=kind_phys), intent(in) :: dt - real (kind=kind_phys), intent(in) :: del(:,:), prsl(:,:) - real (kind=kind_phys), intent(inout) :: q(:,:), t(:,:), & - & cwm(:,:) - real (kind=kind_phys), intent(out) :: rn(:), sr(:), rainp(:,:) - real (kind=kind_phys), intent(in) :: u00k(:,:) - real (kind=kind_phys), intent(in) :: psautco(:), prautco(:), & - & evpco, wminco(:), wk1(:) - logical, intent(in) :: lprnt - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg -! -! Local variables - real (kind=kind_phys) g, h1, h1000 - &, d00 - &, elwv, eliv, row - &, epsq, eliw - &, rcp, rrow - parameter ( h1=1.e0, h1000=1000.0 & - &, d00=0.e0, row=1.e3 & - &, epsq=2.e-12) -! - - real(kind=kind_phys), parameter :: cons_0=0.0, cons_p01=0.01 & - &, cons_20=20.0 & - &, cons_m30=-30.0, cons_50=50.0 -! - real (kind=kind_phys) rnp(im), psautco_l(im), prautco_l(im) & - &, wk2(im) -! - real (kind=kind_phys) err(im), ers(im), precrl(im) & - &, precsl(im), precrl1(im), precsl1(im) & - &, rq(im), condt(im) & - &, conde(im), rconde(im), tmt0(im) & - &, wmin(im,km), wmink(im), pres(im) & - &, wmini(im,km), ccr(im) & - &, tt(im), qq(im), ww(im) & - &, zaodt - real (kind=kind_phys) cclim(km) -! - integer iw(im,km), ipr(im), iwl(im), iwl1(im) -! - logical comput(im) -! - real (kind=kind_phys) ke, rdt, us, climit, cws, csm1 - &, crs1, crs2, cr, aa2, dtcp, c00, cmr - &, tem, c1, c2, wwn -! &, tem, c1, c2, u00b, u00t, wwn - &, precrk, precsk, pres1, qk, qw, qi - &, qint, fiw, wws, cwmk, expf - &, psaut, psaci, amaxcm, tem1, tem2 - &, tmt0k, psm1, psm2, ppr - &, rprs, erk, pps, sid, rid, amaxps - &, praut, fi, qc, amaxrq, rqkll - integer i, k, ihpr, n -! - ! Initialize CCPP error handling variables - errmsg = '' - errflg = 0 -!-------------- GFS psautco/prautco interstitial ---------------- - do i=1, im - wk2(i) = 1.0-wk1(i) - psautco_l(i) = psautco(1)*wk1(i) + psautco(2)*wk2(i) - prautco_l(i) = prautco(1)*wk1(i) + prautco(2)*wk2(i) - enddo -!-----------------------preliminaries --------------------------------- -! -! do k=1,km -! do i=1,im -! cll(i,k) = 0.0 -! enddo -! enddo -! - g=grav - elwv=hvap - eliv=hvap+hfus - eliw=eliv-elwv - rcp=h1/cp - rrow=h1/row - rdt = h1 / dt -! ke = 2.0e-5 ! commented on 09/10/99 -- opr value -! ke = 2.0e-6 -! ke = 1.0e-5 -!!! ke = 5.0e-5 -!! ke = 7.0e-5 - ke = evpco -! ke = 7.0e-5 - us = h1 - climit = 1.0e-20 - cws = 0.025 -! - zaodt = 800.0 * rdt -! - csm1 = 5.0000e-8 * zaodt - crs1 = 5.00000e-6 * zaodt - crs2 = 6.66600e-10 * zaodt - cr = 5.0e-4 * zaodt - aa2 = 1.25e-3 * zaodt -! - ke = ke * sqrt(rdt) -! ke = ke * sqrt(zaodt) -! - dtcp = dt * rcp -! -! c00 = 1.5e-1 * dt -! c00 = 10.0e-1 * dt -! c00 = 3.0e-1 * dt !05/09/2000 -! c00 = 1.0e-4 * dt !05/09/2000 -! c00 = prautco * dt !05/09/2000 - cmr = 1.0 / 3.0e-4 -! cmr = 1.0 / 5.0e-4 -! c1 = 100.0 - c1 = 300.0 - c2 = 0.5 -! -! -!--------calculate c0 and cmr using lc at previous step----------------- -! - do k=1,km - do i=1,im - tem = (prsl(i,k)*0.00001) -! tem = sqrt(tem) - iw(i,k) = 0.0 -! wmin(i,k) = 1.0e-5 * tem -! wmini(i,k) = 1.0e-5 * tem ! testing for ras -! - - wmin(i,k) = wminco(1) * tem - wmini(i,k) = wminco(2) * tem - - - rainp(i,k) = 0.0 - - enddo - enddo - do i=1,im -! c0(i) = 1.5e-1 -! cmr(i) = 3.0e-4 -! - iwl1(i) = 0 - precrl1(i) = d00 - precsl1(i) = d00 - comput(i) = .false. - rn(i) = d00 - sr(i) = d00 - ccr(i) = d00 -! - rnp(i) = d00 - enddo -!> -# Select columns where rain can be produced, where -!!\f[ -!! cwm > \min (wmin, wmini) -!!\f] -!! where the cloud water and ice conversion threshold: -!! \f[ -!! wmin=wminco(1)\times prsl\times 10^{-5} -!! \f] -!! \f[ -!! wmini=wminco(2)\times prsl\times 10^{-5} -!! \f] - -!------------select columns where rain can be produced-------------- - do k=1, km-1 - do i=1,im - tem = min(wmin(i,k), wmini(i,k)) - if (cwm(i,k) > tem) comput(i) = .true. - enddo - enddo - ihpr = 0 - do i=1,im - if (comput(i)) then - ihpr = ihpr + 1 - ipr(ihpr) = i - endif - enddo -!*********************************************************************** -!-----------------begining of precipitation calculation----------------- -!*********************************************************************** -! do k=km-1,2,-1 - do k=km,1,-1 - do n=1,ihpr - precrl(n) = precrl1(n) - precsl(n) = precsl1(n) - err (n) = d00 - ers (n) = d00 - iwl (n) = 0 -! - i = ipr(n) - tt(n) = t(i,k) - qq(n) = q(i,k) - ww(n) = cwm(i,k) - wmink(n) = wmin(i,k) - pres(n) = prsl(i,k) -! - precrk = max(cons_0, precrl1(n)) - precsk = max(cons_0, precsl1(n)) - wwn = max(ww(n), climit) -! if (wwn .gt. wmink(n) .or. (precrk+precsk) .gt. d00) then - if (wwn > climit .or. (precrk+precsk) > d00) then - comput(n) = .true. - else - comput(n) = .false. - endif - enddo -! -! es(1:ihpr) = fpvs(tt(1:ihpr)) - do n=1,ihpr - if (comput(n)) then - i = ipr(n) - conde(n) = (dt/g) * del(i,k) - condt(n) = conde(n) * rdt - rconde(n) = h1 / conde(n) - qk = max(epsq, qq(n)) - tmt0(n) = tt(n) - 273.16 - wwn = max(ww(n), climit) -! -! pl = pres(n) * 0.01 -! call qsatd(tt(n), pl, qc) -! rq(n) = max(qq(n), epsq) / max(qc, 1.0e-10) -! rq(n) = max(1.0e-10, rq(n)) ! -- relative humidity--- -! -! the global qsat computation is done in pa - pres1 = pres(n) -! qw = es(n) - qw = min(pres1, fpvs(tt(n))) - qw = eps * qw / (pres1 + epsm1 * qw) - qw = max(qw,epsq) -! -! tmt15 = min(tmt0(n), cons_m15) -! ai = 0.008855 -! bi = 1.0 -! if (tmt0(n) .lt. -20.0) then -! ai = 0.007225 -! bi = 0.9674 -! endif -! qi = qw * (bi + ai*min(tmt0(n),cons_0)) -! qint = qw * (1.-0.00032*tmt15*(tmt15+15.)) -! - qi = qw - qint = qw -! if (tmt0(n).le.-40.) qint = qi -! -!-------------------ice-water id number iw------------------------------ -!> -# Calculate ice-water identification number IW (see algorithm in -!! \ref condense). - if(tmt0(n) < -15.) then - fi = qk - u00k(i,k)*qi - if(fi > d00 .or. wwn > climit) then - iwl(n) = 1 - else - iwl(n) = 0 - endif -! endif - elseif (tmt0(n) >= 0.) then - iwl(n) = 0 -! -! if(tmt0(n).lt.0.0.and.tmt0(n).ge.-15.0) then - else - iwl(n) = 0 - if(iwl1(n) == 1 .and. wwn > climit) iwl(n) = 1 - endif -! -! if(tmt0(n).ge.0.) then -! iwl(n) = 0 -! endif -!----------------the satuation specific humidity------------------------ - fiw = float(iwl(n)) - qc = (h1-fiw)*qint + fiw*qi -!----------------the relative humidity---------------------------------- - if(qc <= 1.0e-10) then - rq(n) = d00 - else - rq(n) = qk / qc - endif -!----------------cloud cover ratio ccr---------------------------------- -!> -# Calculate cloud fraction \f$b\f$ (see algorithm in \ref condense) - if(rq(n) < u00k(i,k)) then - ccr(n) = d00 - elseif(rq(n) >= us) then - ccr(n) = us - else - rqkll = min(us,rq(n)) - ccr(n) = h1-sqrt((us-rqkll)/(us-u00k(i,k))) - endif -! - endif - enddo -!-------------------ice-water id number iwl------------------------------ -! do n=1,ihpr -! if (comput(n) .and. (ww(n) .gt. climit)) then -! if (tmt0(n) .lt. -15.0 -! * .or. (tmt0(n) .lt. 0.0 .and. iwl1(n) .eq. 1)) -! * iwl(n) = 1 -! cll(ipr(n),k) = 1.0 ! cloud cover! -! cll(ipr(n),k) = min(1.0, ww(n)*cclim(k)) ! cloud cover! -! endif -! enddo -! -!> -# Precipitation production by auto conversion and accretion -!! - The autoconversion of cloud ice to snow (\f$P_{saut}\f$) is simulated -!! using the equation from Lin et al.(1983)\cite lin_et_al_1983 -!!\f[ -!! P_{saut}=a_{1}(cwm-wmini) -!!\f] -!! Since snow production in this process is caused by the increase in -!! size of cloud ice particles due to depositional growth and -!! aggregation of small ice particles, \f$P_{saut}\f$ is a function of -!! temperature as determined by coefficient \f$a_{1}\f$, given by -!! \f[ -!! a_{1}=psautco \times dt \times exp\left[ 0.025\left(T-273.15\right)\right] -!! \f] -!! -!! - The accretion of cloud ice by snow (\f$P_{saci}\f$) in the -!! regions where cloud ice exists is simulated by -!!\f[ -!! P_{saci}=C_{s}cwm P_{s} -!!\f] -!! where \f$P_{s}\f$ is the precipitation rate of snow. The collection -!! coefficient \f$C_{s}\f$ is a function of temperature since the open -!! structures of ice crystals at relative warm temperatures are more -!! likely to stick, given a collision, than crystals of other shapes -!! (Rogers (1979) \cite rogers_1979). Above the freezing level, -!! \f$C_{s}\f$ is expressed by -!!\f[ -!! C_{s}=c_{1}exp\left[ 0.025\left(T-273.15\right)\right] -!!\f] -!! where \f$c_{1}=1.25\times 10^{-3} m^{2}kg^{-1}s^{-1}\f$ are used. -!! \f$C_{s}\f$ is set to zero below the freezing level. -!! -!--- precipitation production -- auto conversion and accretion -! - do n=1,ihpr - if (comput(n) .and. ccr(n) > 0.0) then - wws = ww(n) - cwmk = max(cons_0, wws) - i = ipr(n) -! amaxcm = max(cons_0, cwmk - wmink(n)) - if (iwl(n) == 1) then ! ice phase - amaxcm = max(cons_0, cwmk - wmini(i,k)) - expf = dt * exp(0.025*tmt0(n)) - psaut = min(cwmk, psautco_l(i)*expf*amaxcm) - ww(n) = ww(n) - psaut - cwmk = max(cons_0, ww(n)) -! cwmk = max(cons_0, ww(n)-wmini(i,k)) - psaci = min(cwmk, aa2*expf*precsl1(n)*cwmk) - - ww(n) = ww(n) - psaci - precsl(n) = precsl(n) + (wws - ww(n)) * condt(n) - else ! liquid water -! -!> - Following Sundqvist et al. (1989)\cite sundqvist_et_al_1989, -!! the autoconversion of cloud water to rain (\f$P_{raut}\f$) can be -!! parameterized from the cloud water mixing ratio \f$m\f$ and cloud -!! coverage \f$b\f$, that is, -!!\f[ -!! P_{raut}=(prautco \times dt )\times (cwm-wmin)\left\{1-exp[-(\frac{cwm-wmin}{m_{r}b})^{2}]\right\} -!!\f] -!! where \f$m_{r}\f$ is \f$3.0\times 10^{-4}\f$. -! for using sundqvist precip formulation of rain -! - amaxcm = max(cons_0, cwmk - wmink(n)) -!! amaxcm = cwmk - tem1 = precsl1(n) + precrl1(n) - tem2 = min(max(cons_0, 268.0-tt(n)), cons_20) - tem = (1.0+c1*sqrt(tem1*rdt)) * (1+c2*sqrt(tem2)) -! - tem2 = amaxcm * cmr * tem / max(ccr(n),cons_p01) - tem2 = min(cons_50, tem2*tem2) -! praut = c00 * tem * amaxcm * (1.0-exp(-tem2)) - praut = (prautco_l(i)*dt) * tem * amaxcm - & * (1.0-exp(-tem2)) - praut = min(praut, cwmk) - ww(n) = ww(n) - praut -! -! - Calculate the accretion of cloud water by rain \f$P_{racw}\f$, -! can be expressed using the cloud mixing ratio \f$cwm\f$ and rainfall -! rate \f$P_{r}\f$: -!\f[ -! P_{racw}=C_{r}cwmP_{r} -!\f] -! where \f$C_{r}=5.0\times10^{-4}m^{2}kg^{-1}s^{-1}\f$ is the -! collection coeffiecient. Note that this process is not included in -! current operational physcics. -! below is for zhao's precip formulation (water) -! -! amaxcm = max(cons_0, cwmk - wmink(n)) -! praut = min(cwmk, c00*amaxcm*amaxcm) -! ww(n) = ww(n) - praut -! -! cwmk = max(cons_0, ww(n)) -! tem1 = precsl1(n) + precrl1(n) -! pracw = min(cwmk, cr*dt*tem1*cwmk) -! ww(n) = ww(n) - pracw -! - precrl(n) = precrl(n) + (wws - ww(n)) * condt(n) -! -!hchuang code change [+1l] : add record to record information in vertical -! turn rnp in unit of ww (cwm and q, kg/kg ???) - rnp(n) = rnp(n) + (wws - ww(n)) - endif - endif - enddo -!> -# Evaporation of precipitation (\f$E_{rr}\f$ and \f$E_{rs}\f$) -!!\n Evaporation of precipitation is an important process that moistens -!! the layers below cloud base. Through this process, some of the -!! precipitating water is evaporated back to the atmosphere and the -!! precipitation efficiency is reduced. -!! - Evaporation of rain is calculated using the equation (Sundqvist(1988)\cite sundqvist_1988): -!!\f[ -!! E_{rr}= evpco \times (u-f)(P_{r})^{\beta} -!!\f] -!! where \f$u\f$ is u00k, \f$f\f$ is the relative humidity. -!! \f$\beta = 0.5\f$ are empirical parameter. -!! - Evaporation of snow is calculated using the equation: -!!\f[ -!! E_{rs}=[C_{rs1}+C_{rs2}(T-273.15)](\frac{u-f}{u})P_{s} -!!\f] -!! where \f$C_{rs1}=5\times 10^{-6}m^{2}kg^{-1}s^{-1}\f$ and -!! \f$C_{rs2}=6.67\times 10^{-10}m^{2}kg^{-1}K^{-1}s^{-1}\f$. The -!! evaporation of melting snow below the freezing level is ignored in -!! this scheme because of the difficulty in the latent heat treatment -!! since the surface of a melting snowflake is usually covered by a -!! thin layer of liquid water. -! -!-----evaporation of precipitation------------------------- -!**** err & ers positive--->evaporation-- negtive--->condensation -! - do n=1,ihpr - if (comput(n)) then - i = ipr(n) - qk = max(epsq, qq(n)) - tmt0k = max(cons_m30, tmt0(n)) - precrk = max(cons_0, precrl(n)) - precsk = max(cons_0, precsl(n)) - amaxrq = max(cons_0, u00k(i,k)-rq(n)) * conde(n) -!---------------------------------------------------------------------- -! increase the evaporation for strong/light prec -!---------------------------------------------------------------------- - ppr = ke * amaxrq * sqrt(precrk) -! ppr = ke * amaxrq * sqrt(precrk*rdt) - if (tmt0(n) .ge. 0.) then - pps = 0. - else - pps = (crs1+crs2*tmt0k) * amaxrq * precsk / u00k(i,k) - end if -!---------------correct if over-evapo./cond. occurs-------------------- - erk=precrk+precsk - if(rq(n).ge.1.0e-10) erk = amaxrq * qk * rdt / rq(n) - if (ppr+pps .gt. abs(erk)) then - rprs = erk / (precrk+precsk) - ppr = precrk * rprs - pps = precsk * rprs - endif - ppr = min(ppr, precrk) - pps = min(pps, precsk) - err(n) = ppr * rconde(n) - ers(n) = pps * rconde(n) - precrl(n) = precrl(n) - ppr -!hchuang code change [+1l] : add record to record information in vertical -! use err for kg/kg/dt not the ppr (mm/dt=kg/m2/dt) -! - rnp(n) = rnp(n) - err(n) -! - precsl(n) = precsl(n) - pps - endif - enddo -!> -# Melting of snow (\f$P_{sm1}\f$ and \f$P_{sm2}\f$) -!!\n In this scheme, we allow snow melting to take place in certain -!! temperature regions below the freezing level in two ways. In both -!! cases, the melted snow is assumed to become raindrops. -!! - One is the continuous melting of snow due to the increase in -!! temperature as it falls down through the freezing level. This -!! process is parameterized as a function of temperature and snow -!! precipitation rate, that is, -!!\f[ -!! P_{sm1}=C_{sm}(T-273.15)^{2}P_{s} -!!\f] -!! where \f$C_{sm}=5\times 10^{-8}m^{2}kg^{-1}K^{-2}s^{-1}\f$ -!! cause the falling snow to melt almost completely before it reaches -!! the \f$T=278.15 K\f$ level. -!! - Another is the immediate melting of melting snow by collection of -!! the cloud water below the freezing level. In order to calculate the -!! melting rate, the collection rate of cloud water by melting snow is -!! computed first. Similar to the collection of cloud water by rain, -!! the collection of cloud water by melting snow can be parameterized -!! to be proportional to the cloud water mixing ratio \f$m\f$ and the -!! precipitation rate of snow \f$P_{s}\f$: -!!\f[ -!! P_{sacw}=C_{r}cwmP_{s} -!!\f] -!! where \f$C_{r}\f$ is the collection coefficient, -!! \f$C_{r}=5.0\times 10^{-4}m^{2}kg^{-1}s^{-1}\f$ . The melting rate -!! of snow then can be computed from -!!\f[ -!! P_{sm2}=C_{ws}P_{sacw} -!!\f] -!! where \f$C_{ws}=0.025\f$. -!--------------------melting of the snow-------------------------------- - do n=1,ihpr - if (comput(n)) then - if (tmt0(n) .gt. 0.) then - amaxps = max(cons_0, precsl(n)) - psm1 = csm1 * tmt0(n) * tmt0(n) * amaxps - psm2 = cws * cr * max(cons_0, ww(n)) * amaxps - ppr = (psm1 + psm2) * conde(n) - if (ppr .gt. amaxps) then - ppr = amaxps - psm1 = amaxps * rconde(n) - endif - precrl(n) = precrl(n) + ppr -! -!hchuang code change [+1l] : add record to record information in vertical -! turn ppr (mm/dt=kg/m2/dt) to kg/kg/dt -> ppr/air density (kg/m3) - rnp(n) = rnp(n) + ppr * rconde(n) -! - precsl(n) = precsl(n) - ppr - else - psm1 = d00 - endif -! -!---------------update t and q------------------------------------------ -!> - Update t and q. -!!\f[ -!! t=t-\frac{L}{C_{p}}(E_{rr}+E_{rs}+P_{sm1})\times dt -!!\f] -!!\f[ -!! q=q+(E_{rr}+E_{rs})\times dt -!!\f] - - tt(n) = tt(n) - dtcp * (elwv*err(n)+eliv*ers(n)+eliw*psm1) - qq(n) = qq(n) + dt * (err(n)+ers(n)) - endif - enddo -! - do n=1,ihpr - iwl1(n) = iwl(n) - precrl1(n) = max(cons_0, precrl(n)) - precsl1(n) = max(cons_0, precsl(n)) - i = ipr(n) - t(i,k) = tt(n) - q(i,k) = qq(n) - cwm(i,k) = ww(n) - iw(i,k) = iwl(n) -!hchuang code change [+1l] : add record to record information in vertical -! rnp = precrl1*rconde(n) unit in kg/kg/dt -! - rainp(i,k) = rnp(n) - enddo -! -! move water from vapor to liquid should the liquid amount be negative -! - do i = 1, im - if (cwm(i,k) < 0.) then - tem = q(i,k) + cwm(i,k) - if (tem >= 0.0) then - q(i,k) = tem - t(i,k) = t(i,k) - elwv * rcp * cwm(i,k) - cwm(i,k) = 0. - elseif (q(i,k) > 0.0) then - cwm(i,k) = tem - t(i,k) = t(i,k) + elwv * rcp * q(i,k) - q(i,k) = 0.0 - endif - endif - enddo -! - enddo ! k loop ends here! -!********************************************************************** -!-----------------------end of precipitation processes----------------- -!********************************************************************** -! -!> -# Calculate precipitation at surface (\f$rn\f$)and determine -!! fraction of frozen precipitation (\f$sr\f$). -!!\f[ -!! rn= (P_{r}(\eta_{sfc})+P_{s}(\eta_{sfc}))/10^3 -!!\f] -!!\f[ -!! sr=\frac{P_{s}(\eta_{sfc})}{P_{s}(\eta_{sfc})+P_{r}(\eta_{sfc})} -!!\f] - do n=1,ihpr - i = ipr(n) - rn(i) = (precrl1(n) + precsl1(n)) * rrow ! precip at surface -! -!----sr=1 if sfc prec is rain ; ----sr=-1 if sfc prec is snow -!----sr=0 for both of them or no sfc prec -! -! rid = 0. -! sid = 0. -! if (precrl1(n) .ge. 1.e-13) rid = 1. -! if (precsl1(n) .ge. 1.e-13) sid = -1. -! sr(i) = rid + sid ! sr=1 --> rain, sr=-1 -->snow, sr=0 -->both -! chuang, june 2013: change sr to define fraction of frozen precipitation instead -! because wpc uses it in their winter experiment - - rid = precrl1(n) + precsl1(n) - if (rid < 1.e-13) then - sr(i) = 0. - else - sr(i) = precsl1(n)/rid - endif - enddo -! - return - end subroutine zhaocarr_precpd_run -!> @} - - end module zhaocarr_precpd diff --git a/physics/MP/Zhao_Carr/zhaocarr_precpd.meta b/physics/MP/Zhao_Carr/zhaocarr_precpd.meta deleted file mode 100644 index 86e6c7d67..000000000 --- a/physics/MP/Zhao_Carr/zhaocarr_precpd.meta +++ /dev/null @@ -1,270 +0,0 @@ -[ccpp-table-properties] - name = zhaocarr_precpd - type = scheme - dependencies = ../../tools/funcphys.f90,../../hooks/machine.F,../../hooks/physcons.F90 - -######################################################################## -[ccpp-arg-table] - name = zhaocarr_precpd_init - type = scheme -[imp_physics] - standard_name = control_for_microphysics_scheme - long_name = choice of microphysics scheme - units = flag - dimensions = () - type = integer - intent = in -[imp_physics_zhao_carr] - standard_name = identifier_for_zhao_carr_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme - units = flag - dimensions = () - type = integer - intent = in -[imp_physics_zhao_carr_pdf] - standard_name = identifier_for_zhao_carr_pdf_microphysics_scheme - long_name = choice of Zhao-Carr microphysics scheme with PDF clouds - units = flag - dimensions = () - type = integer - intent = in -[errmsg] - standard_name = ccpp_error_message - long_name = error message for error handling in CCPP - units = none - dimensions = () - type = character - kind = len=* - intent = out -[errflg] - standard_name = ccpp_error_code - long_name = error code for error handling in CCPP - units = 1 - dimensions = () - type = integer - intent = out - -######################################################################## -[ccpp-arg-table] - name = zhaocarr_precpd_run - type = scheme -[im] - standard_name = horizontal_loop_extent - long_name = horizontal loop extent - units = count - dimensions = () - type = integer - intent = in -[km] - standard_name = vertical_layer_dimension - long_name = vertical layer dimension - units = count - dimensions = () - type = integer - intent = in -[dt] - standard_name = timestep_for_physics - long_name = physics time step - units = s - dimensions = () - type = real - kind = kind_phys - intent = in -[del] - standard_name = air_pressure_difference_between_midlayers - long_name = pressure level thickness - units = Pa - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = in -[prsl] - standard_name = air_pressure - long_name = layer mean pressure - units = Pa - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = in -[q] - standard_name = specific_humidity_of_new_state - long_name = water vapor specific humidity - units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout -[cwm] - standard_name = cloud_liquid_water_mixing_ratio_of_new_state - long_name = moist cloud condensed water mixing ratio - units = kg kg-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout -[t] - standard_name = air_temperature_of_new_state - long_name = layer mean air temperature - units = K - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = inout -[rn] - standard_name = lwe_thickness_of_explicit_precipitation_amount - long_name = explicit precipitation amount on physics timestep - units = m - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out -[grav] - standard_name = gravitational_acceleration - long_name = gravitational acceleration - units = m s-2 - dimensions = () - type = real - kind = kind_phys - intent = in -[hvap] - standard_name = latent_heat_of_vaporization_of_water_at_0C - long_name = latent heat of evaporation/sublimation - units = J kg-1 - dimensions = () - type = real - kind = kind_phys - intent = in -[hfus] - standard_name = latent_heat_of_fusion_of_water_at_0C - long_name = latent heat of fusion - units = J kg-1 - dimensions = () - type = real - kind = kind_phys - intent = in -[ttp] - standard_name = triple_point_temperature_of_water - long_name = triple point temperature of water - units = K - dimensions = () - type = real - kind = kind_phys - intent = in -[cp] - standard_name = specific_heat_of_dry_air_at_constant_pressure - long_name = specific heat of dry air at constant pressure - units = J kg-1 K-1 - dimensions = () - type = real - kind = kind_phys - intent = in -[eps] - standard_name = ratio_of_dry_air_to_water_vapor_gas_constants - long_name = rd/rv - units = none - dimensions = () - type = real - kind = kind_phys - intent = in -[epsm1] - standard_name = ratio_of_dry_air_to_water_vapor_gas_constants_minus_one - long_name = (rd/rv) - 1 - units = none - dimensions = () - type = real - kind = kind_phys - intent = in -[sr] - standard_name = ratio_of_snowfall_to_rainfall - long_name = ratio of snowfall to large-scale rainfall - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out -[rainp] - standard_name = tendency_of_rain_water_mixing_ratio_due_to_microphysics - long_name = tendency of rain water mixing ratio due to microphysics - units = kg kg-1 s-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = out -[u00k] - standard_name = critical_relative_humidity - long_name = critical relative humidity - units = frac - dimensions = (horizontal_loop_extent,vertical_layer_dimension) - type = real - kind = kind_phys - intent = in -[psautco] - standard_name = autoconversion_to_snow_coefficient - long_name = conversion coefficient from cloud ice to snow - units = none - dimensions = (2) - type = real - kind = kind_phys - intent = in -[prautco] - standard_name = autoconversion_to_rain_coefficient - long_name = conversion coefficient from cloud water to rain - units = none - dimensions = (2) - type = real - kind = kind_phys - intent = in -[evpco] - standard_name = precipitation_evaporation_coefficient - long_name = coefficient for evaporation of rainfall - units = none - dimensions = () - type = real - kind = kind_phys - intent = in -[wminco] - standard_name = cloud_condensate_autoconversion_threshold_coefficient - long_name = conversion coefficient from cloud liquid and ice to precipitation - units = none - dimensions = (2) - type = real - kind = kind_phys - intent = in -[wk1] - standard_name = grid_size_related_coefficient_used_in_scale_sensitive_schemes - long_name = grid size related coefficient used in scale-sensitive schemes - units = none - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in -[lprnt] - standard_name = flag_print - long_name = flag for printing diagnostics to output - units = flag - dimensions = () - type = logical - intent = in -[jpr] - standard_name = horizontal_index_of_printed_column - long_name = horizontal index of printed column - units = index - dimensions = () - type = integer - intent = in -[errmsg] - standard_name = ccpp_error_message - long_name = error message for error handling in CCPP - units = none - dimensions = () - type = character - kind = len=* - intent = out -[errflg] - standard_name = ccpp_error_code - long_name = error code for error handling in CCPP - units = 1 - dimensions = () - type = integer - intent = out - diff --git a/physics/Radiation/radiation_clouds.f b/physics/Radiation/radiation_clouds.f index d779d56c2..3582dd8ba 100644 --- a/physics/Radiation/radiation_clouds.f +++ b/physics/Radiation/radiation_clouds.f @@ -28,7 +28,6 @@ ! ntrac, ntcw, ntiw, ntrw, ntsw, ntgl, ntclamt, ! ! imp_physics, imp_physics_nssl, imp_physics_fer_hires, ! ! imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, ! -! imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, ! ! imp_physics_mg, iovr, iovr_rand, iovr_maxrand, iovr_max, ! ! iovr_dcorr, iovr_exp, iovr_exprand, idcor, idcor_con, ! ! idcor_hogan, idcor_oreopoulos, ! @@ -44,8 +43,6 @@ ! clds,mtop,mbot,de_lgth,alpha) ! ! ! ! internal/external accessable subroutines: ! -! 'progcld_zhao_carr' --- zhao/moorthi prognostic cloud scheme ! -! 'progcld_zhao_carr_pdf' --- zhao/moorthi prognostic cloud + pdfcld ! ! 'progcld_gfdl_lin' --- GFDL-Lin cloud microphysics ! ! 'progcld_fer_hires' --- Ferrier-Aligo cloud microphysics ! ! 'progcld_thompson_wsm6' --- Thompson/wsm6 cloud microphysics (EMC) ! @@ -105,8 +102,7 @@ ! boundary (include bl-cld). h,m,l clds domain boundaries are ! ! adjusted for better agreement with observations. ! ! jan 2011, yu-tai hou - changed virtual temperature ! -! as input variable instead of originally computed inside the ! -! two prognostic cld schemes 'progcld_zhao_carr' ! +! as input variable ! ! aug 2012, yu-tai hou - modified subroutine cld_init ! ! to pass all fixed control variables at the start. and set ! ! their correponding internal module variables to be used by ! @@ -226,13 +222,13 @@ module module_radiation_clouds & 161.503, 168.262, 175.248, 182.473, 189.952, 197.699, & & 205.728, 214.055, 222.694, 231.661, 240.971, 250.639/) - public progcld_zhao_carr, progcld_zhao_carr_pdf, & - & progcld_gfdl_lin, progclduni, progcld_fer_hires, & + public progcld_gfdl_lin, progclduni, progcld_fer_hires, & & cld_init, radiation_clouds_prop, & & progcld_thompson_wsm6, progcld_thompson, cal_cldfra3, & & find_cloudLayers, adjust_cloudIce, adjust_cloudH2O, & & adjust_cloudFinal, gethml + public cld_frac_XuRandall ! ================= contains ! ================= @@ -242,8 +238,6 @@ module module_radiation_clouds !!\param si model vertical sigma layer interface !!\param NLAY vertical layer number !!\param imp_physics cloud microphysics scheme control flag -!!\n =99: Zhao-Carr/Sundqvist microphysics cloud -!!\n =98: Zhao-Carr/Sundqvist microphysics cloud = pdfcld !!\n =11: GFDL microphysics cloud !!\n =8: Thompson microphysics !!\n =6: WSM6 microphysics @@ -302,11 +296,7 @@ subroutine cld_init & if (me == 0) then print *, VTAGCLD !print out version tag print *,' - Using Prognostic Cloud Method' - if (imp_physics == 99) then - print *,' --- Zhao/Carr/Sundqvist microphysics' - elseif (imp_physics == 98) then - print *,' --- zhao/carr/sundqvist + pdf cloud' - elseif (imp_physics == 11) then + if (imp_physics == 11) then print *,' --- GFDL Lin cloud microphysics' elseif (imp_physics == 8) then print *,' --- Thompson cloud microphysics' @@ -346,7 +336,6 @@ subroutine radiation_clouds_prop & & imp_physics, imp_physics_nssl, imp_physics_fer_hires, & & imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, & & imp_physics_tempo, & - & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, & & imp_physics_mg, iovr, iovr_rand, iovr_maxrand, iovr_max, & & iovr_dcorr, iovr_exp, iovr_exprand, idcor, idcor_con, & & idcor_hogan, idcor_oreopoulos, lcrick, lcnorm, & @@ -365,8 +354,7 @@ subroutine radiation_clouds_prop & ! ================= subprogram documentation block ================ ! ! ! -! subprogram: radiation_clouds_prop computes cloud related quantities using ! -! zhao/moorthi's prognostic cloud microphysics scheme. ! +! subprogram: radiation_clouds_prop computes cloud related quantities ! ! ! ! abstract: this program computes cloud fractions from cloud ! ! condensates, calculates liquid/ice cloud droplet effective radius, ! @@ -379,8 +367,6 @@ subroutine radiation_clouds_prop & ! ! ! subprograms called: ! ! ! -! 'progcld_zhao_carr' --- zhao/moorthi prognostic cloud scheme ! -! 'progcld_zhao_carr_pdf' --- zhao/moorthi prognostic cloud + pdfcld ! ! 'progcld_gfdl_lin' --- GFDL-Lin cloud microphysics ! ! 'progcld_fer_hires' --- Ferrier-Aligo cloud microphysics ! ! 'progcld_thompson_wsm6' --- Thompson/wsm6 cloud microphysics (EMC) ! @@ -437,8 +423,6 @@ subroutine radiation_clouds_prop & ! imp_physics_gfdl : GFDL microphysics scheme ! ! imp_physics_thompson : Thompson microphysics scheme ! ! imp_physics_wsm6 : WSMG microphysics scheme ! -! imp_physics_zhao_carr : Zhao-Carr microphysics scheme ! -! imp_physics_zhao_carr_pdf : Zhao-Carr microphysics scheme with PDF clouds ! imp_physics_mg : Morrison-Gettelman microphysics scheme ! ! iovr : choice of cloud-overlap ! ! iovr_rand : flag of cloud-overlap: random (=0) ! @@ -523,8 +507,6 @@ subroutine radiation_clouds_prop & & imp_physics_thompson, ! Flag for thompsonscheme & imp_physics_tempo, ! Flag for TEMPO scheme & imp_physics_wsm6, ! Flag for wsm6 scheme - & imp_physics_zhao_carr, ! Flag for zhao-carr scheme - & imp_physics_zhao_carr_pdf, ! Flag for zhao-carr+PDF scheme & imp_physics_mg ! Flag for MG scheme integer, intent(in) :: & @@ -618,9 +600,8 @@ subroutine radiation_clouds_prop & end do end do - if (imp_physics == imp_physics_zhao_carr .or. & - & imp_physics == imp_physics_mg) then ! zhao/moorthi's p - ! or unified cloud and/or with MG microphysics + if (imp_physics == imp_physics_mg) then ! + ! unified cloud and/or with MG microphysics if (uni_cld .and. ncndl >= 2) then call progclduni (plyr, plvl, tlyr, tvly, ccnd, ncndl, & ! --- inputs @@ -632,30 +613,8 @@ subroutine radiation_clouds_prop & & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs & cld_reice,cld_rwp, cld_rerain,cld_swp, & & cld_resnow) - else - call progcld_zhao_carr (plyr ,plvl, tlyr, tvly, qlyr, & ! --- inputs - & qstl, rhly, ccnd(1:IX,1:NLAY,1), xlat, xlon, & - & slmsk, dz, delp, IX, NLAY, NLP1, uni_cld, & - & lmfshal, lmfdeep2, xr_con, xr_exp, & - & cldcov, effrl, effri, effrr, effrs, effr_in, & - & dzlay, & - & cldtot, cldcnv, lcrick, lcnorm, con_ttp, & ! inout - & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs - & cld_reice,cld_rwp, cld_rerain,cld_swp, & - & cld_resnow) endif - elseif(imp_physics == imp_physics_zhao_carr_pdf) then ! zhao/moorthi's prognostic cloud+pdfcld - - call progcld_zhao_carr_pdf (plyr, plvl, tlyr, tvly, qlyr, & ! --- inputs - & qstl, rhly, ccnd(1:IX,1:NLAY,1), cnvw, cnvc, & - & xlat, xlon, slmsk, dz, delp, IX, NLAY, NLP1, & - & deltaq, sup, kdt, me, dzlay, & - & cldtot, cldcnv, lcrick, lcnorm, con_thgni, & ! inout - & con_ttp, cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs - & cld_reice,cld_rwp, cld_rerain,cld_swp, & - & cld_resnow) - elseif (imp_physics == imp_physics_gfdl) then ! GFDL cloud scheme if (.not. lgfdlmprad) then @@ -756,715 +715,128 @@ subroutine radiation_clouds_prop & & cldcov(:,1:LM), effrl, effri, effrs, & & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & & dzlay, gridkm, top_at_1, & - & cldtot, cldcnv, & ! inout - & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs - & cld_reice,cld_rwp, cld_rerain,cld_swp, & - & cld_resnow) - else - - !-- MYNN PBL or convective GF - !-- use cloud fractions with SGS clouds - do k=1,NLAY - do i=1,IX - cld_frac(i,k) = clouds1(i,k) - enddo - enddo - - ! --- use clduni as with the GFDL microphysics. - ! --- make sure that effr_in=.true. in the input.nml! - call progclduni (plyr, plvl, tlyr, tvly, ccnd, ncndl, & ! --- inputs - & xlat, xlon, slmsk, dz, delp, IX, NLAY, NLP1, & - & cld_frac, & - & effrl, effri, effrr, effrs, effr_in , & - & dzlay, & - & cldtot, cldcnv, lcrick, lcnorm, con_ttp, & ! inout - & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs - & cld_reice,cld_rwp, cld_rerain,cld_swp, & - & cld_resnow) - endif - - else - ! MYNN PBL or GF convective are not used - - if (icloud == 3) then - call progcld_thompson (plyr,plvl,tlyr,qlyr,qstl,rhly, & ! --- inputs - & tracer1,xlat,xlon,slmsk,dz,delp, & - & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & - & ntsw-1,ntgl-1, & - & IX, LM, NLP1, uni_cld, lmfshal, lmfdeep2, & - & cldcov(:,1:LM), effrl, effri, effrs, & - & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & - & dzlay, gridkm, top_at_1, & - & cldtot, cldcnv, & ! inout - & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs - & cld_reice,cld_rwp, cld_rerain,cld_swp, & - & cld_resnow) - - else - call progcld_thompson_wsm6 (plyr,plvl,tlyr,qlyr,qstl, & ! --- inputs - & rhly,tracer1,xlat,xlon,slmsk,dz,delp, & - & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & - & ntsw-1,ntgl-1,con_ttp,xr_cnvcld, & - & IX, NLAY, NLP1, xr_con, xr_exp, uni_cld, & - & lmfshal, lmfdeep2, & - & cldcov(:,1:NLAY), cnvw, effrl, effri, effrs, & - & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & - & dzlay, & - & cldtot, cldcnv, lcnorm, & ! inout - & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs - & cld_reice,cld_rwp, cld_rerain,cld_swp, & - & cld_resnow) - endif - endif ! MYNN PBL or GF - - endif ! end if_imp_physics - -!> - Compute SFC/low/middle/high cloud top pressure for each cloud -!! domain for given latitude. -! ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,L,m,h; -! --- i=1,2 are low-lat (<45 degree) and pole regions) - - do i =1, IX - rxlat(i) = abs( xlat(i) / con_pi ) ! if xlat in pi/2 -> -pi/2 range -! rxlat(i) = abs(0.5 - xlat(i)/con_pi) ! if xlat in 0 -> pi range - enddo - - do id = 1, 4 - tem1 = ptopc(id,2) - ptopc(id,1) - - do i =1, IX - ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*rxlat(i)-1.0 ) - enddo - enddo - - ! Compute cloud decorrelation length - if (idcor == idcor_hogan) then - call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) - endif - if (idcor == idcor_oreopoulos) then - call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) - endif - if (idcor == idcor_con) then - de_lgth(:) = dcorr_con - endif - - ! Call subroutine get_alpha_exper to define alpha parameter for exponential cloud overlap options - if ( iovr == iovr_dcorr .or. iovr == iovr_exp & - & .or. iovr == iovr_exprand) then - call get_alpha_exper(ix, nLay, iovr, iovr_exprand, dzlay, & - & de_lgth, cld_frac, alpha) - else - de_lgth(:) = 0. - alpha(:,:) = 0. - endif - -!> - Call gethml() to compute low,mid,high,total, and boundary layer -!! cloud fractions and clouds top/bottom layer indices for low, mid, -!! and high clouds. -! --- compute low, mid, high, total, and boundary layer cloud fractions -! and clouds top/bottom layer indices for low, mid, and high clouds. -! The three cloud domain boundaries are defined by ptopc. The cloud -! overlapping method is defined by control flag 'iovr', which may -! be different for lw and sw radiation programs. - - call gethml & -! --- inputs: - & ( plyr, ptop1, cldtot, cldcnv, dz, de_lgth, alpha, & - & IX, NLAY, iovr, iovr_rand, iovr_maxrand, iovr_max, & - & iovr_dcorr, iovr_exp, iovr_exprand, top_at_1, si, & -! --- outputs: - & clds, mtop, mbot & - & ) - -!................................... - end subroutine radiation_clouds_prop - -!> This subroutine computes cloud related quantities using -!! zhao/moorthi's prognostic cloud microphysics scheme. -!>\section progcld_zhao_carr General Algorithm - subroutine progcld_zhao_carr & - & ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, & ! --- inputs: - & xlat,xlon,slmsk,dz,delp, IX, NLAY, NLP1, & - & uni_cld, lmfshal, lmfdeep2, xr_con, xr_exp, cldcov, & - & effrl,effri,effrr,effrs,effr_in, & - & dzlay, cldtot, cldcnv, lcrick, lcnorm, con_ttp, & - & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs - & cld_reice,cld_rwp, cld_rerain,cld_swp, cld_resnow & - & ) - -! ================= subprogram documentation block ================ ! -! ! -! subprogram: progcld_zhao_carr computes cloud related quantities using ! -! zhao/moorthi's prognostic cloud microphysics scheme. ! -! ! -! abstract: this program computes cloud fractions from cloud ! -! condensates, calculates liquid/ice cloud droplet effective radius, ! -! and computes the low, mid, high, total and boundary layer cloud ! -! fractions and the vertical indices of low, mid, and high cloud ! -! top and base. the three vertical cloud domains are set up in the ! -! initial subroutine "cld_init". ! -! ! -! usage: call progcld_zhao_carr ! -! ! -! attributes: ! -! language: fortran 90 ! -! machine: ibm-sp, sgi ! -! ! -! ! -! ==================== definition of variables ==================== ! -! ! -! input variables: ! -! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! -! plvl (IX,NLP1) : model level pressure in mb (100Pa) ! -! tlyr (IX,NLAY) : model layer mean temperature in k ! -! tvly (IX,NLAY) : model layer virtual temperature in k ! -! qlyr (IX,NLAY) : layer specific humidity in gm/gm ! -! qstl (IX,NLAY) : layer saturate humidity in gm/gm ! -! rhly (IX,NLAY) : layer relative humidity (=qlyr/qstl) ! -! clw (IX,NLAY) : layer cloud condensate amount ! -! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2! -! range, otherwise see in-line comment ! -! xlon (IX) : grid longitude in radians (not used) ! -! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) ! -! dz (ix,nlay) : layer thickness (km) ! -! delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! -! IX : horizontal dimention ! -! NLAY,NLP1 : vertical layer/level dimensions ! -! uni_cld : logical - true for cloud fraction from shoc ! -! lmfshal : logical - true for mass flux shallow convection ! -! lmfdeep2 : logical - true for mass flux deep convection ! -! cldcov : layer cloud fraction (used when uni_cld=.true. ! -! effrl : effective radius for liquid water -! effri : effective radius for ice water -! effrr : effective radius for rain water -! effrs : effective radius for snow water -! effr_in : logical, if .true. use input effective radii -! dzlay(ix,nlay) : thickness between model layer centers (km) ! -! lmfshal : mass-flux shallow conv scheme flag ! -! lmfdeep2 : scale-aware mass-flux deep conv scheme flag ! -! lcrick : control flag for eliminating CRICK ! -! =t: apply layer smoothing to eliminate CRICK ! -! =f: do not apply layer smoothing ! -! lcnorm : control flag for in-cld condensate ! -! =t: normalize cloud condensate ! -! =f: not normalize cloud condensate ! -! ! -! output variables: ! -! cloud profiles: ! -! cld_frac (:,:) - layer total cloud fraction ! -! cld_lwp (:,:) - layer cloud liq water path (g/m**2) ! -! cld_reliq (:,:) - mean eff radius for liq cloud (micron) ! -! cld_iwp (:,:) - layer cloud ice water path (g/m**2) ! -! cld_reice (:,:) - mean eff radius for ice cloud (micron) ! -! cld_rwp (:,:) - layer rain drop water path not assigned ! -! cld_rerain(:,:) - mean eff radius for rain drop (micron) ! -! *** cld_swp (:,:) - layer snow flake water path not assigned ! -! cld_resnow(:,:) - mean eff radius for snow flake (micron) ! -! ! -! ==================== end of description ===================== ! -! - implicit none - -! --- inputs - integer, intent(in) :: IX, NLAY, NLP1 - - logical, intent(in) :: uni_cld, lmfshal, lmfdeep2, effr_in, & - & lcrick, lcnorm - - real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & - & tlyr, tvly, qlyr, qstl, rhly, clw, cldcov, delp, dz, & - & effrl, effri, effrr, effrs, dzlay - - real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & - & slmsk - real (kind=kind_phys), intent(in) :: con_ttp, xr_con, xr_exp - -! --- inputs/outputs - - real (kind=kind_phys), dimension(:,:), intent(inout) :: & - & cld_frac, cld_lwp, cld_reliq, cld_iwp, cld_reice, & - & cld_rwp, cld_rerain, cld_swp, cld_resnow - -! --- local variables: - real (kind=kind_phys), dimension(IX,NLAY) :: cldtot, cldcnv, & - & cwp, cip, crp, csp, rew, rei, res, rer, tem2d, clwf - - real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & - & tem1, tem2, tem3 - - integer :: i, k, id, nf - -! -!===> ... begin here -! -!> - Assgin liquid/ice/rain/snow cloud effective radius from input or predefined values. - if(effr_in) then - do k = 1, NLAY - do i = 1, IX - cldtot(i,k) = 0.0 - cldcnv(i,k) = 0.0 - cwp (i,k) = 0.0 - cip (i,k) = 0.0 - crp (i,k) = 0.0 - csp (i,k) = 0.0 - rew (i,k) = effrl (i,k) - rei (i,k) = effri (i,k) - rer (i,k) = effrr (i,k) - res (i,k) = effrs (i,k) - tem2d (i,k) = min(1.0, max(0.0,(con_ttp-tlyr(i,k))*0.05)) - clwf(i,k) = 0.0 - enddo - enddo - else - do k = 1, NLAY - do i = 1, IX - cldtot(i,k) = 0.0 - cldcnv(i,k) = 0.0 - cwp (i,k) = 0.0 - cip (i,k) = 0.0 - crp (i,k) = 0.0 - csp (i,k) = 0.0 - rew (i,k) = reliq_def ! default liq radius to 10 micron - rei (i,k) = reice_def ! default ice radius to 50 micron - rer (i,k) = rrain_def ! default rain radius to 1000 micron - res (i,k) = rsnow_def ! default snow radius to 250 micron - tem2d (i,k) = min(1.0, max(0.0, (con_ttp-tlyr(i,k))*0.05)) - clwf(i,k) = 0.0 - enddo - enddo - endif -! - if ( lcrick ) then - do i = 1, IX - clwf(i,1) = 0.75*clw(i,1) + 0.25*clw(i,2) - clwf(i,nlay) = 0.75*clw(i,nlay) + 0.25*clw(i,nlay-1) - enddo - do k = 2, NLAY-1 - do i = 1, IX - clwf(i,K) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1) - enddo - enddo - else - do k = 1, NLAY - do i = 1, IX - clwf(i,k) = clw(i,k) - enddo - enddo - endif - -!> - Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . - - do k = 1, NLAY - do i = 1, IX - clwt = max(0.0, clwf(i,k)) * gfac * delp(i,k) - cip(i,k) = clwt * tem2d(i,k) - cwp(i,k) = clwt - cip(i,k) - enddo - enddo - -!> - Compute effective liquid cloud droplet radius over land. - - if(.not. effr_in) then - do i = 1, IX - if (nint(slmsk(i)) == 1) then - do k = 1, NLAY - rew(i,k) = 5.0 + 5.0 * tem2d(i,k) - enddo - endif - enddo - endif - - if (uni_cld) then ! use unified sgs clouds generated outside - do k = 1, NLAY - do i = 1, IX - cldtot(i,k) = cldcov(i,k) - enddo - enddo - - else - -!> - Compute layer cloud fraction. - - - if (.not. lmfshal) then - call cloud_fraction_XuRandall & - & ( IX, NLAY, xr_con, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) & ! --- outputs - else - call cloud_fraction_mass_flx_1 & - & ( IX, NLAY, lmfdeep2, xr_con, xr_exp, plyr, clwf, rhly, & - & qstl, & ! --- inputs - & cldtot ) - endif - - endif ! if (uni_cld) then - - do k = 1, NLAY - do i = 1, IX - if (cldtot(i,k) < climit) then - cldtot(i,k) = 0.0 - cwp(i,k) = 0.0 - cip(i,k) = 0.0 - crp(i,k) = 0.0 - csp(i,k) = 0.0 - endif - enddo - enddo - - if ( lcnorm ) then - do k = 1, NLAY - do i = 1, IX - if (cldtot(i,k) >= climit) then - tem1 = 1.0 / max(climit2, cldtot(i,k)) - cwp(i,k) = cwp(i,k) * tem1 - cip(i,k) = cip(i,k) * tem1 - crp(i,k) = crp(i,k) * tem1 - csp(i,k) = csp(i,k) * tem1 - endif - enddo - enddo - endif - -!> - Compute effective ice cloud droplet radius following Heymsfield -!! and McFarquhar (1996) \cite heymsfield_and_mcfarquhar_1996. - - if(.not.effr_in) then - do k = 1, NLAY - do i = 1, IX - tem2 = tlyr(i,k) - con_ttp - - if (cip(i,k) > 0.0) then - tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k)) - - if (tem2 < -50.0) then - rei(i,k) = (1250.0/9.917) * tem3 ** 0.109 - elseif (tem2 < -40.0) then - rei(i,k) = (1250.0/9.337) * tem3 ** 0.08 - elseif (tem2 < -30.0) then - rei(i,k) = (1250.0/9.208) * tem3 ** 0.055 - else - rei(i,k) = (1250.0/9.387) * tem3 ** 0.031 - endif -! rei(i,k) = max(20.0, min(rei(i,k), 300.0)) -! rei(i,k) = max(10.0, min(rei(i,k), 100.0)) - rei(i,k) = max(10.0, min(rei(i,k), 150.0)) -! rei(i,k) = max(5.0, min(rei(i,k), 130.0)) - endif - enddo - enddo - endif - -! - do k = 1, NLAY - do i = 1, IX - cld_frac(i,k) = cldtot(i,k) - cld_lwp(i,k) = cwp(i,k) - cld_reliq(i,k) = rew(i,k) - cld_iwp(i,k) = cip(i,k) - cld_reice(i,k) = rei(i,k) -! cld_rwp(i,k) = 0.0 - cld_rerain(i,k) = rer(i,k) -! cld_swp(i,k) = 0.0 - cld_resnow(i,k) = res(i,k) - enddo - enddo -! -!................................... - end subroutine progcld_zhao_carr -!----------------------------------- -!----------------------------------- - -!> This subroutine computes cloud related quantities using -!! zhao/moorthi's prognostic cloud microphysics scheme + pdfcld. -!>\section progcld_zhao_carr_pdf General Algorithm - subroutine progcld_zhao_carr_pdf & - & ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw,cnvw,cnvc, & ! --- inputs: - & xlat,xlon,slmsk, dz, delp, & - & ix, nlay, nlp1, & - & deltaq,sup,kdt,me, & - & dzlay, cldtot, cldcnv, lcrick, lcnorm, con_thgni, con_ttp, & - & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs - & cld_reice,cld_rwp, cld_rerain,cld_swp, cld_resnow & - & ) - -! ================= subprogram documentation block ================ ! -! ! -! subprogram: progcld_zhao_carr_pdf computes cloud related quantities using ! -! zhao/moorthi's prognostic cloud microphysics scheme. ! -! ! -! abstract: this program computes cloud fractions from cloud ! -! condensates, calculates liquid/ice cloud droplet effective radius, ! -! and computes the low, mid, high, total and boundary layer cloud ! -! fractions and the vertical indices of low, mid, and high cloud ! -! top and base. the three vertical cloud domains are set up in the ! -! initial subroutine "cld_init". ! -! ! -! usage: call progcld_zhao_carr_pdf ! -! ! -! attributes: ! -! language: fortran 90 ! -! machine: ibm-sp, sgi ! -! ! -! ! -! ==================== defination of variables ==================== ! -! ! -! input variables: ! -! plyr (ix,nlay) : model layer mean pressure in mb (100pa) ! -! plvl (ix,nlp1) : model level pressure in mb (100pa) ! -! tlyr (ix,nlay) : model layer mean temperature in k ! -! tvly (ix,nlay) : model layer virtual temperature in k ! -! qlyr (ix,nlay) : layer specific humidity in gm/gm ! -! qstl (ix,nlay) : layer saturate humidity in gm/gm ! -! rhly (ix,nlay) : layer relative humidity (=qlyr/qstl) ! -! clw (ix,nlay) : layer cloud condensate amount ! -! xlat (ix) : grid latitude in radians, default to pi/2 -> -pi/2! -! range, otherwise see in-line comment ! -! xlon (ix) : grid longitude in radians (not used) ! -! slmsk (ix) : sea/land mask array (sea:0,land:1,sea-ice:2) ! -! dz (ix,nlay) : layer thickness (km) ! -! delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! -! ix : horizontal dimention ! -! nlay,nlp1 : vertical layer/level dimensions ! -! cnvw (ix,nlay) : layer convective cloud condensate ! -! cnvc (ix,nlay) : layer convective cloud cover ! -! deltaq(ix,nlay) : half total water distribution width ! -! sup : supersaturation ! -! dzlay(ix,nlay) : thickness between model layer centers (km) ! -! lcrick : control flag for eliminating crick ! -! =t: apply layer smoothing to eliminate crick ! -! =f: do not apply layer smoothing ! -! lcnorm : control flag for in-cld condensate ! -! =t: normalize cloud condensate ! -! =f: not normalize cloud condensate ! -! ! -! output variables: ! -! cloud profiles: ! -! cld_frac (:,:) - layer total cloud fraction ! -! cld_lwp (:,:) - layer cloud liq water path (g/m**2) ! -! cld_reliq (:,:) - mean eff radius for liq cloud (micron) ! -! cld_iwp (:,:) - layer cloud ice water path (g/m**2) ! -! cld_reice (:,:) - mean eff radius for ice cloud (micron) ! -! cld_rwp (:,:) - layer rain drop water path not assigned ! -! cld_rerain(:,:) - mean eff radius for rain drop (micron) ! -! *** cld_swp (:,:) - layer snow flake water path not assigned ! -! cld_resnow(:,:) - mean eff radius for snow flake (micron) ! -! ! -! ==================== end of description ===================== ! -! - implicit none - -! --- inputs - integer, intent(in) :: ix, nlay, nlp1,kdt - logical, intent(in) :: lcrick, lcnorm - real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & - & tlyr, tvly, qlyr, qstl, rhly, clw, dz, delp, dzlay -! & tlyr, tvly, qlyr, qstl, rhly, clw, cnvw, cnvc -! real (kind=kind_phys), dimension(:,:), intent(in) :: deltaq - real (kind=kind_phys), intent(in) :: con_thgni, con_ttp - real (kind=kind_phys), dimension(:,:) :: deltaq, cnvw, cnvc - real (kind=kind_phys) qtmp,qsc,rhs - real (kind=kind_phys), intent(in) :: sup - real (kind=kind_phys), parameter :: epsq = 1.0e-12 - - real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & - & slmsk - integer :: me - -! --- inputs/outputs - - real (kind=kind_phys), dimension(:,:), intent(inout) :: & - & cld_frac, cld_lwp, cld_reliq, cld_iwp, cld_reice, & - & cld_rwp, cld_rerain, cld_swp, cld_resnow - -! --- local variables: - real (kind=kind_phys), dimension(ix,nlay) :: cldtot, cldcnv, & - & cwp, cip, crp, csp, rew, rei, res, rer, tem2d, clwf - - real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & - & tem1, tem2, tem3 - - integer :: i, k, id, nf - -! -!===> ... begin here -! - do k = 1, nlay - do i = 1, ix - cldtot(i,k) = 0.0 - cldcnv(i,k) = 0.0 - cwp (i,k) = 0.0 - cip (i,k) = 0.0 - crp (i,k) = 0.0 - csp (i,k) = 0.0 - rew (i,k) = reliq_def ! default liq radius to 10 micron - rei (i,k) = reice_def ! default ice radius to 50 micron - rer (i,k) = rrain_def ! default rain radius to 1000 micron - res (i,k) = rsnow_def ! default snow radius to 250 micron - tem2d (i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) ) - clwf(i,k) = 0.0 - enddo - enddo -! - if ( lcrick ) then - do i = 1, ix - clwf(i,1) = 0.75*clw(i,1) + 0.25*clw(i,2) - clwf(i,nlay) = 0.75*clw(i,nlay) + 0.25*clw(i,nlay-1) - enddo - do k = 2, nlay-1 - do i = 1, ix - clwf(i,k) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1) - enddo - enddo - else - do k = 1, nlay - do i = 1, ix - clwf(i,k) = clw(i,k) - enddo - enddo - endif - - if(kdt==1) then - do k = 1, nlay - do i = 1, ix - deltaq(i,k) = (1.-0.95)*qstl(i,k) - enddo - enddo - endif - -!> -# Calculate liquid/ice condensate path in \f$ g/m^2 \f$ + & cldtot, cldcnv, & ! inout + & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs + & cld_reice,cld_rwp, cld_rerain,cld_swp, & + & cld_resnow) + else - do k = 1, nlay - do i = 1, ix - clwt = max(0.0,(clwf(i,k)+cnvw(i,k))) * gfac * delp(i,k) - cip(i,k) = clwt * tem2d(i,k) - cwp(i,k) = clwt - cip(i,k) - enddo - enddo + !-- MYNN PBL or convective GF + !-- use cloud fractions with SGS clouds + do k=1,NLAY + do i=1,IX + cld_frac(i,k) = clouds1(i,k) + enddo + enddo -!> -# Calculate effective liquid cloud droplet radius over land. + ! --- use clduni as with the GFDL microphysics. + ! --- make sure that effr_in=.true. in the input.nml! + call progclduni (plyr, plvl, tlyr, tvly, ccnd, ncndl, & ! --- inputs + & xlat, xlon, slmsk, dz, delp, IX, NLAY, NLP1, & + & cld_frac, & + & effrl, effri, effrr, effrs, effr_in , & + & dzlay, & + & cldtot, cldcnv, lcrick, lcnorm, con_ttp, & ! inout + & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs + & cld_reice,cld_rwp, cld_rerain,cld_swp, & + & cld_resnow) + endif - do i = 1, ix - if (nint(slmsk(i)) == 1) then - do k = 1, nlay - rew(i,k) = 5.0 + 5.0 * tem2d(i,k) - enddo - endif - enddo + else + ! MYNN PBL or GF convective are not used -!> -# Calculate layer cloud fraction. + if (icloud == 3) then + call progcld_thompson (plyr,plvl,tlyr,qlyr,qstl,rhly, & ! --- inputs + & tracer1,xlat,xlon,slmsk,dz,delp, & + & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & + & ntsw-1,ntgl-1, & + & IX, LM, NLP1, uni_cld, lmfshal, lmfdeep2, & + & cldcov(:,1:LM), effrl, effri, effrs, & + & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & + & dzlay, gridkm, top_at_1, & + & cldtot, cldcnv, & ! inout + & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs + & cld_reice,cld_rwp, cld_rerain,cld_swp, & + & cld_resnow) - do k = 1, nlay - do i = 1, ix - tem1 = tlyr(i,k) - 273.16 - if(tem1 < (con_thgni - 273.16)) then ! for pure ice, has to be consistent with gscond - qsc = sup * qstl(i,k) - rhs = sup else - qsc = qstl(i,k) - rhs = 1.0 + call progcld_thompson_wsm6 (plyr,plvl,tlyr,qlyr,qstl, & ! --- inputs + & rhly,tracer1,xlat,xlon,slmsk,dz,delp, & + & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & + & ntsw-1,ntgl-1,con_ttp,xr_cnvcld, & + & IX, NLAY, NLP1, xr_con, xr_exp, uni_cld, & + & lmfshal, lmfdeep2, & + & cldcov(:,1:NLAY), cnvw, effrl, effri, effrs, & + & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & + & dzlay, & + & cldtot, cldcnv, lcnorm, & ! inout + & cld_frac, cld_lwp, cld_reliq, cld_iwp, & ! --- outputs + & cld_reice,cld_rwp, cld_rerain,cld_swp, & + & cld_resnow) endif - if(rhly(i,k) >= rhs) then - cldtot(i,k) = 1.0 - else - qtmp = qlyr(i,k) + clwf(i,k) - qsc - if(deltaq(i,k) > epsq) then -! if(qtmp <= -deltaq(i,k) .or. cwmik < epsq) then - if(qtmp <= -deltaq(i,k)) then - cldtot(i,k) = 0.0 - elseif(qtmp >= deltaq(i,k)) then - cldtot(i,k) = 1.0 - else - cldtot(i,k) = 0.5*qtmp/deltaq(i,k) + 0.5 - cldtot(i,k) = max(cldtot(i,k),0.0) - cldtot(i,k) = min(cldtot(i,k),1.0) - endif - else - if(qtmp > 0.) then - cldtot(i,k) = 1.0 - else - cldtot(i,k) = 0.0 - endif - endif - endif - cldtot(i,k) = cnvc(i,k) + (1-cnvc(i,k))*cldtot(i,k) - cldtot(i,k) = max(cldtot(i,k),0.) - cldtot(i,k) = min(cldtot(i,k),1.) + endif ! MYNN PBL or GF - enddo - enddo + endif ! end if_imp_physics - do k = 1, nlay - do i = 1, ix - if (cldtot(i,k) < climit) then - cldtot(i,k) = 0.0 - cwp(i,k) = 0.0 - cip(i,k) = 0.0 - crp(i,k) = 0.0 - csp(i,k) = 0.0 - endif - enddo +!> - Compute SFC/low/middle/high cloud top pressure for each cloud +!! domain for given latitude. +! ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,L,m,h; +! --- i=1,2 are low-lat (<45 degree) and pole regions) + + do i =1, IX + rxlat(i) = abs( xlat(i) / con_pi ) ! if xlat in pi/2 -> -pi/2 range +! rxlat(i) = abs(0.5 - xlat(i)/con_pi) ! if xlat in 0 -> pi range enddo - if ( lcnorm ) then - do k = 1, nlay - do i = 1, ix - if (cldtot(i,k) >= climit) then - tem1 = 1.0 / max(climit2, cldtot(i,k)) - cwp(i,k) = cwp(i,k) * tem1 - cip(i,k) = cip(i,k) * tem1 - crp(i,k) = crp(i,k) * tem1 - csp(i,k) = csp(i,k) * tem1 - endif - enddo + do id = 1, 4 + tem1 = ptopc(id,2) - ptopc(id,1) + + do i =1, IX + ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*rxlat(i)-1.0 ) enddo + enddo + + ! Compute cloud decorrelation length + if (idcor == idcor_hogan) then + call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) + endif + if (idcor == idcor_oreopoulos) then + call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) + endif + if (idcor == idcor_con) then + de_lgth(:) = dcorr_con endif -!> -# Calculate effective ice cloud droplet radius following Heymsfield -!! and McFarquhar (1996) \cite heymsfield_and_mcfarquhar_1996. - - do k = 1, nlay - do i = 1, ix - tem2 = tlyr(i,k) - con_ttp + ! Call subroutine get_alpha_exper to define alpha parameter for exponential cloud overlap options + if ( iovr == iovr_dcorr .or. iovr == iovr_exp & + & .or. iovr == iovr_exprand) then + call get_alpha_exper(ix, nLay, iovr, iovr_exprand, dzlay, & + & de_lgth, cld_frac, alpha) + else + de_lgth(:) = 0. + alpha(:,:) = 0. + endif - if (cip(i,k) > 0.0) then -! tem3 = gord * cip(i,k) * (plyr(i,k)/delp(i,k)) / tvly(i,k) - tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k)) +!> - Call gethml() to compute low,mid,high,total, and boundary layer +!! cloud fractions and clouds top/bottom layer indices for low, mid, +!! and high clouds. +! --- compute low, mid, high, total, and boundary layer cloud fractions +! and clouds top/bottom layer indices for low, mid, and high clouds. +! The three cloud domain boundaries are defined by ptopc. The cloud +! overlapping method is defined by control flag 'iovr', which may +! be different for lw and sw radiation programs. - if (tem2 < -50.0) then - rei(i,k) = (1250.0/9.917) * tem3 ** 0.109 - elseif (tem2 < -40.0) then - rei(i,k) = (1250.0/9.337) * tem3 ** 0.08 - elseif (tem2 < -30.0) then - rei(i,k) = (1250.0/9.208) * tem3 ** 0.055 - else - rei(i,k) = (1250.0/9.387) * tem3 ** 0.031 - endif -! rei(i,k) = max(20.0, min(rei(i,k), 300.0)) -! rei(i,k) = max(10.0, min(rei(i,k), 100.0)) - rei(i,k) = max(10.0, min(rei(i,k), 150.0)) -! rei(i,k) = max(5.0, min(rei(i,k), 130.0)) - endif - enddo - enddo + call gethml & +! --- inputs: + & ( plyr, ptop1, cldtot, cldcnv, dz, de_lgth, alpha, & + & IX, NLAY, iovr, iovr_rand, iovr_maxrand, iovr_max, & + & iovr_dcorr, iovr_exp, iovr_exprand, top_at_1, si, & +! --- outputs: + & clds, mtop, mbot & + & ) -! - do k = 1, NLAY - do i = 1, IX - cld_frac(i,k) = cldtot(i,k) - cld_lwp(i,k) = cwp(i,k) - cld_reliq(i,k) = rew(i,k) - cld_iwp(i,k) = cip(i,k) - cld_reice(i,k) = rei(i,k) -! cld_rwp(i,k) = 0.0 - cld_rerain(i,k) = rer(i,k) -! cld_swp(i,k) = 0.0 - cld_resnow(i,k) = res(i,k) - enddo - enddo -! !................................... - end subroutine progcld_zhao_carr_pdf -!----------------------------------- + end subroutine radiation_clouds_prop !----------------------------------- @@ -1822,6 +1194,8 @@ subroutine progcld_fer_hires & real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & & tem1, tem2, tem3 + real (kind=kind_phys) :: xrc3 + integer :: i, k, id, nf ! @@ -1893,14 +1267,22 @@ subroutine progcld_fer_hires & !> - Calculate layer cloud fraction. if (.not. lmfshal) then - call cloud_fraction_XuRandall & - & ( IX, NLAY, xr_con, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) & ! --- outputs + xrc3 = xr_con + do k = 1, NLAY-1 + do i = 1, IX + cldtot(i,k) = cld_frac_XuRandall(plyr(i,k), qstl(i,k), & + & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0.) + end do + end do else - call cloud_fraction_mass_flx_1 & - & ( IX, NLAY, lmfdeep2, xr_con, xr_exp, plyr, clwf, rhly, & - & qstl, & ! --- inputs - & cldtot ) & ! --- outputs + xrc3 = 100. + if (lmfdeep2) xrc3 = xr_con + do k = 1, NLAY-1 + do i = 1, IX + cldtot(i,k) = cld_frac_XuRandall(plyr(i,k), qstl(i,k), & + & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0.) + end do + end do endif endif ! if (uni_cld) then @@ -2072,8 +1454,12 @@ subroutine progcld_thompson_wsm6 & real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & & tem1, tem2, tem3 + real (kind=kind_phys) :: xrc3 + integer :: i, k, id, nf + logical :: cond_cfrac_onRH + ! --- constant values real (kind=kind_phys), parameter :: snow2ice = 0.25 real (kind=kind_phys), parameter :: coef_t = 0.025 @@ -2196,14 +1582,24 @@ subroutine progcld_thompson_wsm6 & !> - Calculate layer cloud fraction. if (.not. lmfshal) then - call cloud_fraction_XuRandall & - & ( IX, NLAY, xr_con, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) & ! --- outputs + xrc3 = xr_con + do k = 1, NLAY-1 + do i = 1, IX + cldtot(i,k) = cld_frac_XuRandall(plyr(i,k), qstl(i,k), & + & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0.) + end do + end do else - call cloud_fraction_mass_flx_2 & - & ( IX, NLAY, lmfdeep2, xr_con, xr_exp, plyr, clwf, rhly, & - & qstl, & ! --- inputs - & cldtot ) + xrc3 = 100. + if (lmfdeep2) xrc3 = xr_con + cond_cfrac_onRH = .true. + do k = 1, NLAY-1 + do i = 1, IX + cldtot(i,k) = cld_frac_XuRandall(plyr(i,k), qstl(i,k), & + & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0., & + & cond_cfrac_onRH) + end do + end do endif endif ! if (uni_cld) then @@ -2261,6 +1657,7 @@ subroutine progcld_thompson_wsm6 & enddo enddo + !............................................ end subroutine progcld_thompson_wsm6 !............................................ @@ -2545,6 +1942,7 @@ subroutine progcld_thompson & iwp_ex(i) = iwp_ex(i)*1.E-3 enddo ! + !............................................ end subroutine progcld_thompson !............................................ @@ -3730,155 +3128,56 @@ SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte) END SUBROUTINE adjust_cloudFinal -!> This subroutine computes the Xu-Randall cloud fraction scheme. - subroutine cloud_fraction_XuRandall & - & ( IX, NLAY, xr_con, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) & ! --- outputs - -! --- inputs: - integer, intent(in) :: IX, NLAY - real (kind=kind_phys), intent(in) :: xr_con, xr_exp - real (kind=kind_phys), dimension(:,:), intent(in) :: plyr, clwf, & - & rhly, qstl - -! --- outputs - real (kind=kind_phys), dimension(:,:), intent(inout) :: cldtot - -! --- local variables: - - real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & - & tem1, tem2 - integer :: i, k - -!> - Compute layer cloud fraction. - - clwmin = 0.0 - do k = 1, NLAY - do i = 1, IX - clwt = 1.0e-6 * (plyr(i,k)*0.001) - - if (clwf(i,k) > clwt) then - - onemrh= max( 1.e-10, 1.0-rhly(i,k) ) - clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) - - tem1 = min(max((onemrh*qstl(i,k))**xr_exp,0.0001),1.0) - tem1 = xr_con / tem1 - - value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) - tem2 = sqrt(sqrt(rhly(i,k))) - - cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) - endif - enddo - enddo - - end subroutine cloud_fraction_XuRandall - -!> - subroutine cloud_fraction_mass_flx_1 & - & ( IX, NLAY, lmfdeep2, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) & ! --- outputs - -! --- inputs: - integer, intent(in) :: IX, NLAY - real (kind=kind_phys), intent(in) :: xrc3, xr_exp - real (kind=kind_phys), dimension(:,:), intent(in) :: plyr, clwf, & - & rhly, qstl - logical, intent(in) :: lmfdeep2 - -! --- outputs - real (kind=kind_phys), dimension(:,:), intent(inout) :: cldtot - -! --- local variables: - - real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & - & tem1, tem2 - integer :: i, k - -!> - Compute layer cloud fraction. - - clwmin = 0.0 - do k = 1, NLAY - do i = 1, IX - clwt = 1.0e-6 * (plyr(i,k)*0.001) -! clwt = 2.0e-6 * (plyr(i,k)*0.001) - - if (clwf(i,k) > clwt) then - onemrh= max( 1.e-10, 1.0-rhly(i,k) ) - clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) -! - tem1 = min(max((onemrh*qstl(i,k))**xr_exp,0.0001),1.0) !jhan - if (lmfdeep2) then - tem1 = xrc3 / tem1 - else - tem1 = 100.0 / tem1 - endif -! - value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) - tem2 = sqrt( sqrt(rhly(i,k)) ) - - cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) - endif - enddo - enddo - - end subroutine cloud_fraction_mass_flx_1 - -!> - subroutine cloud_fraction_mass_flx_2 & - & ( IX, NLAY, lmfdeep2, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) & ! --- outputs - -! --- inputs: - integer, intent(in) :: IX, NLAY - real (kind=kind_phys), intent(in) :: xrc3, xr_exp - real (kind=kind_phys), dimension(:,:), intent(in) :: plyr, clwf, & - & rhly, qstl - logical, intent(in) :: lmfdeep2 - -! --- outputs - real (kind=kind_phys), dimension(:,:), intent(inout) :: cldtot - -! --- local variables: - - real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & - & tem1, tem2 - integer :: i, k - -!> - Compute layer cloud fraction. - - clwmin = 0.0 - do k = 1, NLAY-1 - do i = 1, IX - clwt = 1.0e-6 * (plyr(i,k)*0.001) - - if (clwf(i,k) > clwt) then - if(rhly(i,k) > 0.99) then - cldtot(i,k) = 1. - else - onemrh= max( 1.e-10, 1.0-rhly(i,k) ) - clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) - - tem1 = min(max((onemrh*qstl(i,k))**xr_exp,0.0001),1.0) - if (lmfdeep2) then - tem1 = xrc3 / tem1 - else - tem1 = 100.0 / tem1 - endif - - value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) - tem2 = sqrt( sqrt(rhly(i,k)) ) - - cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) - endif +!> This function computes the cloud-fraction following +!! Xu-Randall(1996) \cite xu_and_randall_1996 +!! + function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha, & ! --- inputs + & lambda, factor, cond_cfrac_onRH) + implicit none + ! Inputs + logical, intent(in), optional :: & + & cond_cfrac_onRH ! If true, cloud-fracion set to unity when rh>99% + + real(kind_phys), intent(in) :: & + & p_lay, & !< Pressure (100Pa) + & qs_lay, & !< Saturation vapor-pressure (Pa) + & relhum, & !< Relative humidity + & cld_mr, & !< Total cloud mixing ratio + & alpha, & !< Scheme parameter (default=100) + & lambda, & + & factor ! factor=1.0 for RRTMGP, factor=0 for RRTMG + + ! Outputs + real(kind_phys) :: cld_frac_XuRandall + + ! Locals + real(kind_phys) :: clwt, clwm, onemrh, tem1, tem2, tem3 + +! ! Parameters +! real(kind_phys) :: & +! lambda = 0.50 ! , & +! P = 0.25 + + clwt = 1.0e-6 * (p_lay*0.001) + clwm = clwt * factor + if (cld_mr > clwt) then + if(present(cond_cfrac_onRH) .and. relhum > 0.99) then + cld_frac_XuRandall = 1. else - cldtot(i,k) = 0.0 + onemrh = max(1.e-10, 1.0 - relhum) + tem1 = alpha/min(max((onemrh*qs_lay)**lambda,0.0001),1.0) + tem2 = max(min(tem1*(cld_mr - clwm), 50.0 ), 0.0 ) + tem3 = sqrt(sqrt(relhum)) ! This assumes "p" = 0.25. Identical, but cheaper than relhum**p + ! + cld_frac_XuRandall = max( tem3*(1.0-exp(-tem2)), 0.0 ) endif - enddo - enddo + else + cld_frac_XuRandall = 0.0 + endif + + return + end function - end subroutine cloud_fraction_mass_flx_2 !........................................! end module module_radiation_clouds !>@}