From 6cd1e1a5bcb8e2efbde475a38e10c2a92dac4596 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Mon, 5 Jan 2026 16:18:02 +0000 Subject: [PATCH 01/19] remove Zhao-Carr cloud schemes and combine Xu-Randall cloud fraction calculations --- physics/Radiation/radiation_clouds.f | 980 ++++----------------------- 1 file changed, 134 insertions(+), 846 deletions(-) diff --git a/physics/Radiation/radiation_clouds.f b/physics/Radiation/radiation_clouds.f index d779d56c2..8cce6550b 100644 --- a/physics/Radiation/radiation_clouds.f +++ b/physics/Radiation/radiation_clouds.f @@ -44,8 +44,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) ! @@ -226,8 +224,7 @@ 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, & @@ -242,8 +239,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 +297,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' @@ -379,8 +370,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) ! @@ -618,9 +607,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 +620,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 @@ -759,712 +725,125 @@ subroutine radiation_clouds_prop & & 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$ + & 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 +1201,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 +1274,16 @@ subroutine progcld_fer_hires & !> - Calculate layer cloud fraction. if (.not. lmfshal) then + xrc3 = xr_con call cloud_fraction_XuRandall & - & ( IX, NLAY, xr_con, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) & ! --- outputs + & ( IX, NLAY, xrc3, 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 ) & ! --- outputs + xrc3 = 100. + if (lmfdeep2) xrc3 = xr_con + call cloud_fraction_XuRandall & + & ( IX, NLAY, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs + & cldtot ) endif endif ! if (uni_cld) then @@ -2072,6 +1455,8 @@ 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 ! --- constant values @@ -2196,14 +1581,16 @@ subroutine progcld_thompson_wsm6 & !> - Calculate layer cloud fraction. if (.not. lmfshal) then + xrc3 = xr_con call cloud_fraction_XuRandall & - & ( IX, NLAY, xr_con, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) & ! --- outputs + & ( IX, NLAY, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs + & cldtot ) ! --- outputs 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 + call cloud_fraction_XuRandall & + & ( IX, NLAY, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs + & cldtot, cond_cfrac_onRH = .true.) endif endif ! if (uni_cld) then @@ -2261,6 +1648,7 @@ subroutine progcld_thompson_wsm6 & enddo enddo + !............................................ end subroutine progcld_thompson_wsm6 !............................................ @@ -2545,6 +1933,7 @@ subroutine progcld_thompson & iwp_ex(i) = iwp_ex(i)*1.E-3 enddo ! + !............................................ end subroutine progcld_thompson !............................................ @@ -3732,110 +3121,15 @@ 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 + & ( IX, NLAY, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs + & cldtot, cond_cfrac_onRH) ! --- outputs ! --- inputs: integer, intent(in) :: IX, NLAY - real (kind=kind_phys), intent(in) :: xrc3, xr_exp + real (kind=kind_phys), intent(in) :: xrc3, xr_exp real (kind=kind_phys), dimension(:,:), intent(in) :: plyr, clwf, & & rhly, qstl - logical, intent(in) :: lmfdeep2 + logical, intent(in), optional :: cond_cfrac_onRH ! --- outputs real (kind=kind_phys), dimension(:,:), intent(inout) :: cldtot @@ -3854,31 +3148,25 @@ subroutine cloud_fraction_mass_flx_2 & clwt = 1.0e-6 * (plyr(i,k)*0.001) if (clwf(i,k) > clwt) then - if(rhly(i,k) > 0.99) then + if(present(cond_cfrac_onRH) .and. 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 + tem1 = xrc3 / 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 - else - cldtot(i,k) = 0.0 endif enddo enddo - end subroutine cloud_fraction_mass_flx_2 + end subroutine cloud_fraction_XuRandall !........................................! end module module_radiation_clouds !>@} From 83191f3056b9357df7d28107f2f1a7e059ed9f41 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Wed, 7 Jan 2026 15:49:00 +0000 Subject: [PATCH 02/19] change Xu-Randall cloud fraction scheme to be elemental --- physics/Radiation/radiation_clouds.f | 134 +++++++++++++++------------ 1 file changed, 76 insertions(+), 58 deletions(-) diff --git a/physics/Radiation/radiation_clouds.f b/physics/Radiation/radiation_clouds.f index 8cce6550b..54c673e91 100644 --- a/physics/Radiation/radiation_clouds.f +++ b/physics/Radiation/radiation_clouds.f @@ -1275,15 +1275,21 @@ subroutine progcld_fer_hires & if (.not. lmfshal) then xrc3 = xr_con - call cloud_fraction_XuRandall & - & ( IX, NLAY, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) ! --- outputs + 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 xrc3 = 100. if (lmfdeep2) xrc3 = xr_con - call cloud_fraction_XuRandall & - & ( IX, NLAY, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) + 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 @@ -1459,6 +1465,8 @@ subroutine progcld_thompson_wsm6 & 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 @@ -1582,15 +1590,23 @@ subroutine progcld_thompson_wsm6 & if (.not. lmfshal) then xrc3 = xr_con - call cloud_fraction_XuRandall & - & ( IX, NLAY, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot ) ! --- outputs + 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 xrc3 = 100. if (lmfdeep2) xrc3 = xr_con - call cloud_fraction_XuRandall & - & ( IX, NLAY, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot, cond_cfrac_onRH = .true.) + 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 @@ -3119,54 +3135,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, xrc3, xr_exp, plyr, clwf, rhly, qstl, & ! --- inputs - & cldtot, cond_cfrac_onRH) ! --- 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), optional :: cond_cfrac_onRH - -! --- 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(present(cond_cfrac_onRH) .and. 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) - tem1 = xrc3 / 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 +!> 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 + 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_XuRandall !........................................! end module module_radiation_clouds !>@} From 3a66ae6c2d500ab50908f4b76d8207250bfd1af1 Mon Sep 17 00:00:00 2001 From: "Clara.Draoper-NOAA" Date: Wed, 7 Jan 2026 16:37:41 +0000 Subject: [PATCH 03/19] Fixed too-high lower bound for SMC Deleted code associated with already removed print statements Deleted excess comments. --- physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 | 81 ++++---------------- 1 file changed, 17 insertions(+), 64 deletions(-) diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index a5f855f11..7070efb87 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -187,9 +187,9 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & integer :: lsoil_incr integer, allocatable :: mask_tile(:) - integer,allocatable :: stc_updated(:), slc_updated(:) + integer,allocatable :: stc_updated(:) logical :: soil_freeze, soil_ice - integer :: soiltype, n_stc, n_slc + integer :: soiltype real(kind=kind_phys) :: slc_new integer :: i, j, ij, l, k, ib @@ -198,10 +198,6 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & real(kind=kind_phys) :: smp !< for computing supercooled water real(kind=kind_phys) :: hc_incr - integer :: nother, nsnowupd - integer :: nstcupd, nslcupd, nfrozen, nfrozen_upd - logical :: print_update_stats = .False. - ! --- Initialize CCPP error handling variables errmsg = '' errflg = 0 @@ -223,7 +219,7 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & endif if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) then - print*, "adding land iau increments " + print*, "adding land iau increments CSD-edited" endif if (Land_IAU_Control%lsoil .ne. km) then @@ -236,11 +232,9 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & allocate(stc_inc_flat(Land_IAU_Control%nx * Land_IAU_Control%ny, km)) !GFS_Control%ncols allocate(slc_inc_flat(Land_IAU_Control%nx * Land_IAU_Control%ny, km)) !GFS_Control%ncols allocate(stc_updated(Land_IAU_Control%nx * Land_IAU_Control%ny)) - allocate(slc_updated(Land_IAU_Control%nx * Land_IAU_Control%ny)) !copy background stc stc_updated = 0 - slc_updated = 0 ib = 1 do j = 1, Land_IAU_Control%ny do k = 1, km @@ -260,53 +254,39 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & lensfc = Land_IAU_Control%nx * Land_IAU_Control%ny if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) print*,' adjusting first ', lsoil_incr, ' surface layers only, delt ', delt - ! initialize variables for counts statitics to be zeros - nother = 0 ! grid cells not land - nsnowupd = 0 ! grid cells with snow (temperature not yet updated) - nstcupd = 0 ! grid cells that are updated stc - nslcupd = 0 ! grid cells that are updated slc - nfrozen = 0 ! not update as frozen soil - nfrozen_upd = 0 ! not update as frozen soil - -!TODO---if only fv3 increment files are used, this can be read from file + allocate(mask_tile(lensfc)) call calculate_landinc_mask(weasd, vegtype, soiltyp, lensfc, isice_table, mask_tile) - !IAU increments are in units of 1/sec !Land_IAU_Control%dtp - !* only updating soil temp for now + dz(1) = -zsoil(1) + do k = 2, km + dz(k) = -zsoil(k) + zsoil(k-1) + enddo + !IAU increments are in units of 1/sec ij_loop : do ij = 1, lensfc ! mask: 1 - soil, 2 - snow, 0 - land-ice, -1 - not land if (mask_tile(ij) == 1) then soil_freeze=.false. soil_ice=.false. - do k = 1, lsoil_incr ! k = 1, km + do k = 1, lsoil_incr if ( stc(ij,k) < con_t0c) soil_freeze=.true. if ( smc(ij,k) - slc(ij,k) > 0.001 ) soil_ice=.true. if (Land_IAU_Control%upd_stc) then - stc(ij,k) = stc(ij,k) + stc_inc_flat(ij,k)*delt !Land_IAU_Control%dtp if (k==1) then stc_updated(ij) = 1 - nstcupd = nstcupd + 1 endif + stc(ij,k) = stc(ij,k) + stc_inc_flat(ij,k)*delt endif - if ( (stc(ij,k) < con_t0c) .and. (.not. soil_freeze) .and. (k==1) ) nfrozen_upd = nfrozen_upd + 1 - - ! do not do updates if this layer or any above is frozen + ! do not do SLC updates if this layer or any above is frozen if ( (.not. soil_freeze ) .and. (.not. soil_ice ) ) then if (Land_IAU_Control%upd_slc) then - if (k==1) then - nslcupd = nslcupd + 1 - slc_updated(ij) = 1 - endif - ! apply zero limit here (higher, model-specific limits are later) - slc(ij,k) = max(slc(ij,k) + slc_inc_flat(ij,k)*delt, 0.0) - smc(ij,k) = max(smc(ij,k) + slc_inc_flat(ij,k)*delt, 0.0) + ! if soil moisture is <0.1 mm in layer, prevent DA from further reducing it + slc(ij,k) = max(slc(ij,k) + slc_inc_flat(ij,k)*delt, min(0.0001/dz(k), slc(ij,k))) + smc(ij,k) = max(smc(ij,k) + slc_inc_flat(ij,k)*delt, min(0.0001/dz(k), smc(ij,k))) endif - else - if (k==1) nfrozen = nfrozen+1 endif enddo endif ! if soil/snow point @@ -318,13 +298,10 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & errmsg = 'FATAL ERROR in noahmpdrv_timestep_init: problem in set_soilveg_noahmp' return endif - n_stc = 0 - n_slc = 0 if (Land_IAU_Control%do_stcsmc_adjustment) then if (Land_IAU_Control%upd_stc) then do i=1,lensfc if (stc_updated(i) == 1 ) then ! soil-only location - n_stc = n_stc+1 soiltype = soiltyp(i) do l = 1, lsoil_incr if (abs(stc_inc_flat(i,l)) > Land_IAU_Control%min_T_increment) then @@ -347,34 +324,10 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & enddo endif - if (Land_IAU_Control%upd_slc) then - dz(1) = -zsoil(1) - do l = 2, km - dz(l) = -zsoil(l) + zsoil(l-1) - enddo - do i=1,lensfc - if (slc_updated(i) == 1 ) then - n_slc = n_slc+1 - ! apply SM bounds (later: add upper SMC limit) - do l = 1, lsoil_incr - if (abs(slc_inc_flat(i, l)) > Land_IAU_Control%min_SLC_increment) then - ! noah-mp minimum is 1 mm per layer (in SMC) - ! no need to maintain frozen amount, would be v. small. - slc(i,l) = max( 0.001/dz(l), slc(i,l) ) - smc(i,l) = max( 0.001/dz(l), smc(i,l) ) - endif - enddo - endif - enddo - endif endif - - deallocate(stc_inc_flat, slc_inc_flat, stc_updated, slc_updated) - deallocate(mask_tile) + deallocate(stc_inc_flat, slc_inc_flat, stc_updated) + deallocate(mask_tile) - !Remove non-warning, non-error log write - !write(*,'(a,i4,a,i8)') 'noahmpdrv_timestep_init rank ', Land_IAU_Control%me, ' # of cells with stc update ', nstcupd - end subroutine noahmpdrv_timestep_init From a0f2f28669ae06137fce59dbcee5fcdf31da7da1 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Wed, 7 Jan 2026 18:10:21 +0000 Subject: [PATCH 04/19] use new cld_frac_XuRandall function from module_radiation_clouds in RRTMGP scheme --- .../UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 | 72 +++++-------------- 1 file changed, 16 insertions(+), 56 deletions(-) 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..730e1238e 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 @@ -346,7 +346,8 @@ subroutine cloud_mp_GF(nCol, nLev, lsmask, t_lay, p_lev, p_lay, qs_lay, relhum, cld_cnv_frac !< Convective cloud-fraction (1) ! Local integer :: iCol, iLay - real(kind_phys) :: tem1, deltaP, clwc, qc, qi + real(kind_phys) :: tem1, deltaP, clwc, qc, qi, play_pa + real(kind_phys) :: lambda = 0.50 tem1 = 1.0e5/con_g do iLay = 1, nLev @@ -372,8 +373,9 @@ subroutine cloud_mp_GF(nCol, nLev, lsmask, t_lay, p_lev, p_lay, qs_lay, relhum, 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) + play_pa = p_lay(iCol,iLay) *0.01 + cld_cnv_frac(iCol,iLay) = cld_frac_XuRandall(play_pa, qs_lay(iCol,iLay), & + relhum(iCol,iLay), qc+qi, alpha0, lambda, 1.0) endif enddo enddo @@ -489,7 +491,8 @@ subroutine cloud_mp_SAMF(nCol, nLev, t_lay, p_lev, p_lay, qs_lay, relhum, cld_cnv_frac !< Convective cloud-fraction ! Local integer :: iCol, iLay - real(kind_phys) :: tem0, tem1, deltaP, clwc + real(kind_phys) :: tem0, tem1, deltaP, clwc, play_pa + real(kind_phys) :: lambda = 0.50 tem0 = 1.0e5/con_g do iLay = 1, nLev @@ -504,8 +507,9 @@ subroutine cloud_mp_SAMF(nCol, nLev, t_lay, p_lev, p_lay, qs_lay, relhum, 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) + play_pa = p_lay(iCol,iLay) * 0.01 + cld_cnv_frac(iCol,iLay) = cld_frac_XuRandall(play_pa, qs_lay(iCol,iLay), & + relhum(iCol,iLay), cnv_mixratio(iCol,iLay), alpha0, lambda, 1.0) endif enddo enddo @@ -706,7 +710,8 @@ subroutine cloud_mp_thompson(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldice, i_c cld_rwp !< Cloud rain water path ! Local variables - real(kind_phys) :: tem1, pfac, cld_mr, deltaP, tem2 + real(kind_phys) :: tem1, pfac, cld_mr, deltaP, tem2, play_pa + real(kind_phys) :: lambda = 0.50 real(kind_phys), dimension(nCol, nLev, min(4,ncnd)) :: cld_condensate integer :: iCol,iLay,l @@ -740,8 +745,9 @@ subroutine cloud_mp_thompson(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldice, i_c ! Xu-Randall (1996) cloud-fraction. **Additionally, Conditioned on relative-humidity** 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) + play_pa = p_lay(iCol,iLay) * 0.01 + cld_frac(iCol,iLay) = cld_frac_XuRandall(play_pa, qs_lay(iCol,iLay), & + relhum(iCol,iLay), cld_mr, alpha0, lambda, 1.0, cond_cfrac_onRH) enddo enddo @@ -769,52 +775,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, & From bbf9cff2a7618b8b4aec51f27a5140cb5da3cd74 Mon Sep 17 00:00:00 2001 From: "Clara.Draoper-NOAA" Date: Wed, 7 Jan 2026 18:25:22 +0000 Subject: [PATCH 05/19] Tidy-up for PR. --- physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index 7070efb87..6b96074ce 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -219,7 +219,7 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & endif if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) then - print*, "adding land iau increments CSD-edited" + print*, "adding land iau increments" endif if (Land_IAU_Control%lsoil .ne. km) then @@ -259,9 +259,9 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & call calculate_landinc_mask(weasd, vegtype, soiltyp, lensfc, isice_table, mask_tile) dz(1) = -zsoil(1) - do k = 2, km - dz(k) = -zsoil(k) + zsoil(k-1) - enddo + do k = 2, km + dz(k) = -zsoil(k) + zsoil(k-1) + enddo !IAU increments are in units of 1/sec ij_loop : do ij = 1, lensfc ! mask: 1 - soil, 2 - snow, 0 - land-ice, -1 - not land @@ -269,14 +269,12 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & soil_freeze=.false. soil_ice=.false. - do k = 1, lsoil_incr + do k = 1, lsoil_incr if ( stc(ij,k) < con_t0c) soil_freeze=.true. if ( smc(ij,k) - slc(ij,k) > 0.001 ) soil_ice=.true. if (Land_IAU_Control%upd_stc) then - if (k==1) then - stc_updated(ij) = 1 - endif + if (k==1) stc_updated(ij) = 1 stc(ij,k) = stc(ij,k) + stc_inc_flat(ij,k)*delt endif From ca7e392e4dbf5ec77e98435804361d84901f2bd7 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Thu, 8 Jan 2026 14:38:39 +0000 Subject: [PATCH 06/19] modify the cld_frac_XuRandall function call in RRTMGP scheme --- .../UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 | 31 +++++++++---------- 1 file changed, 14 insertions(+), 17 deletions(-) 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 730e1238e..ae035b5b9 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_rrtmgp_cloud_mp.F90 @@ -346,8 +346,7 @@ subroutine cloud_mp_GF(nCol, nLev, lsmask, t_lay, p_lev, p_lay, qs_lay, relhum, cld_cnv_frac !< Convective cloud-fraction (1) ! Local integer :: iCol, iLay - real(kind_phys) :: tem1, deltaP, clwc, qc, qi, play_pa - real(kind_phys) :: lambda = 0.50 + real(kind_phys) :: tem1, deltaP, clwc, qc, qi tem1 = 1.0e5/con_g do iLay = 1, nLev @@ -372,10 +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. - play_pa = p_lay(iCol,iLay) *0.01 - cld_cnv_frac(iCol,iLay) = cld_frac_XuRandall(play_pa, qs_lay(iCol,iLay), & - relhum(iCol,iLay), qc+qi, alpha0, lambda, 1.0) + ! 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 @@ -491,8 +489,7 @@ subroutine cloud_mp_SAMF(nCol, nLev, t_lay, p_lev, p_lay, qs_lay, relhum, cld_cnv_frac !< Convective cloud-fraction ! Local integer :: iCol, iLay - real(kind_phys) :: tem0, tem1, deltaP, clwc, play_pa - real(kind_phys) :: lambda = 0.50 + real(kind_phys) :: tem0, tem1, deltaP, clwc tem0 = 1.0e5/con_g do iLay = 1, nLev @@ -506,10 +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. - play_pa = p_lay(iCol,iLay) * 0.01 - cld_cnv_frac(iCol,iLay) = cld_frac_XuRandall(play_pa, qs_lay(iCol,iLay), & - relhum(iCol,iLay), cnv_mixratio(iCol,iLay), alpha0, lambda, 1.0) + ! 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 @@ -710,8 +707,7 @@ subroutine cloud_mp_thompson(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldice, i_c cld_rwp !< Cloud rain water path ! Local variables - real(kind_phys) :: tem1, pfac, cld_mr, deltaP, tem2, play_pa - real(kind_phys) :: lambda = 0.50 + real(kind_phys) :: tem1, pfac, cld_mr, deltaP, tem2 real(kind_phys), dimension(nCol, nLev, min(4,ncnd)) :: cld_condensate integer :: iCol,iLay,l @@ -743,11 +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) - play_pa = p_lay(iCol,iLay) * 0.01 - cld_frac(iCol,iLay) = cld_frac_XuRandall(play_pa, qs_lay(iCol,iLay), & - relhum(iCol,iLay), cld_mr, alpha0, lambda, 1.0, 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 From 8415ad51f56668c40b06887805c49b2b0ce289ea Mon Sep 17 00:00:00 2001 From: "Clara.Draoper-NOAA" Date: Mon, 12 Jan 2026 20:04:56 +0000 Subject: [PATCH 07/19] Changes to the post L-IAU checks: * Change tfreeze to match that in NoahMP * Only set SLC to SMC if temp newly > tfreez --- physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index 6b96074ce..b0c539f60 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -197,6 +197,7 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & real(kind=kind_phys) :: smp !< for computing supercooled water real(kind=kind_phys) :: hc_incr + real(kind=kind_phys), parameter :: tfreez_noahmp=273.16 ! tfreez used in NoahMP to determine frozen ground ! --- Initialize CCPP error handling variables errmsg = '' @@ -270,7 +271,7 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & soil_freeze=.false. soil_ice=.false. do k = 1, lsoil_incr - if ( stc(ij,k) < con_t0c) soil_freeze=.true. + if ( stc(ij,k) < tfreez_noahmp ) soil_freeze=.true. if ( smc(ij,k) - slc(ij,k) > 0.001 ) soil_ice=.true. if (Land_IAU_Control%upd_stc) then @@ -306,14 +307,15 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & !the following if case applies when updated stc > melting point, it handles both !case 1: frz ==> frz, recalculate slc, smc remains !case 2: unfrz ==> frz, recalculate slc, smc remains - if (stc(i,l) .LT. con_t0c )then + if (stc(i,l) .LE. tfreez_noahmp )then !recompute supercool liquid water,smc_anl remain unchanged - smp = con_hfus*(con_t0c-stc(i,l))/(con_g*stc(i,l)) !(m) + smp = con_hfus*(tfreez_noahmp-stc(i,l))/(con_g*stc(i,l)) !(m) slc_new=maxsmc(soiltype)*(smp/satpsi(soiltype))**(-1./bb(soiltype)) slc(i,l) = max( min( slc_new, smc(i,l)), 0.0 ) endif - !case 3: frz ==> unfrz (or unfrz ==> unfrz), melt all soil ice (if any) - if (stc(i,l) .GT. con_t0c )then !do not rely on stc_bck + !case 3: frz ==> unfrz melt all soil ice (if any) + if ( stc(i,l) .GT. tfreez_noahmp .and. & + stc(i,l) - stc_inc_flat(i,l)*delt .LE. tfreez_noahmp ) then slc(i,l)=smc(i,l) endif endif From 565d62cd9c7959e0981aabc2674686595fc73359 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Tue, 27 Jan 2026 15:42:25 +0000 Subject: [PATCH 08/19] add the repository update to this PR --- .../Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f | 28 ++++++++++++++----- physics/Radiation/radiation_astronomy.f | 21 +++++++++++++- 2 files changed, 41 insertions(+), 8 deletions(-) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f index 749f778c1..9ec152b8d 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f @@ -265,6 +265,7 @@ subroutine dcyc2t3_run & & flxlwdn_adj real(kind=kind_phys) :: fluxlwnet_adj,fluxlwnet,dT_sfc, & &fluxlwDOWN_jac,lfnc,c1 + real(kind=kind_phys) :: pi2, xlon1, xlon2, sindh ! Length scale for flux-adjustment scaling real(kind=kind_phys), parameter :: & & L = 1. @@ -291,15 +292,17 @@ subroutine dcyc2t3_run & nstp = max(6, nint(tem1)) nstl = max(1, nint(nstp/tem1)) pid12 = con_pi / hour12 + pi2 = con_pi ! ! --- ... sw time-step adjustment for current cosine of zenith angle ! ---------------------------------------------------------- - if (nstl == 1) then - cns = pid12 * (solhr + deltim*f7200 - hour12) + slag - do i = 1, IM - xcosz(i) = sdec*sinlat(i) + cdec*coslat(i)*cos(cns+xlon(i)) - enddo - elseif (nstl == nstp) then +! if (nstl == 1) then +! cns = pid12 * (solhr + deltim*f7200 - hour12) + slag +! do i = 1, IM +! xcosz(i) = sdec*sinlat(i) + cdec*coslat(i)*cos(cns+xlon(i)) +! enddo +! elseif (nstl == nstp) then + if (nstl == nstp) then do i = 1, IM xcosz(i) = coszen(i) enddo @@ -314,7 +317,17 @@ subroutine dcyc2t3_run & do it=1,nstl cns = solang + (float(it)-0.5_kind_phys)*anginc + slag do i = 1, IM - coszn = sdec*sinlat(i) + cdec*coslat(i)*cos(cns+xlon(i)) + xlon1 = cns+xlon(i) + if(xlon1.gt.pi2)then + xlon1=xlon1-2.*pi2 + else if(xlon1.lt.(0.-pi2))then + xlon1=xlon1+2.*pi2 + end if + xlon2 = xlon1 + anginc + if(xlon2.gt.pi2)xlon2 = xlon1 + sindh = (sin(xlon2)-sin(xlon1))/anginc +! coszn = sdec*sinlat(i) + cdec*coslat(i)*cos(cns+xlon(i)) + coszn = sdec*sinlat(i) + cdec*coslat(i)*sindh xcosz(i) = xcosz(i) + max(zero, coszn) if (coszn > czlimt) istsun(i) = istsun(i) + 1 enddo @@ -367,6 +380,7 @@ subroutine dcyc2t3_run & xmu(i) = zero endif + print*,'qingfu, xmu=', i, xmu(i) !> - adjust \a sfc net and downward SW fluxes for zenith angle changes. ! note: sfc emiss effect will not be appied here diff --git a/physics/Radiation/radiation_astronomy.f b/physics/Radiation/radiation_astronomy.f index 90ed7cd45..0f1f41212 100644 --- a/physics/Radiation/radiation_astronomy.f +++ b/physics/Radiation/radiation_astronomy.f @@ -868,6 +868,8 @@ subroutine coszmn & ! --- locals: real (kind=kind_phys) :: coszn, cns, solang, rstp + real (kind=kind_phys) :: alat1, alon1, xlon1, xlon2 + real (kind=kind_phys) :: pi2, rad2deg, sindh integer :: istsun(IM), i, it, j, lat @@ -875,7 +877,10 @@ subroutine coszmn & solang = pid12 * (solhr - f12) ! solar angle at present time rstp = 1.0 / float(nstp) + rad2deg = 15. / pid12 + pi2 = pid12 * 12. +! print*,'qliu test=',nstp,IM,anginc*rad2deg do i = 1, IM coszen(i) = 0.0 istsun(i) = 0 @@ -884,8 +889,22 @@ subroutine coszmn & do it = 1, nstp cns = solang + (float(it)-0.5)*anginc + sollag do i = 1, IM + xlon1 = cns+xlon(i) + if(xlon1.gt.pi2)then + xlon1=xlon1-2.*pi2 + else if(xlon1.lt.(0.-pi2))then + xlon1=xlon1+2.*pi2 + end if + xlon2 = xlon1 + anginc + if(xlon2.gt.pi2)xlon2 = xlon1 + sindh = (sin(xlon2)-sin(xlon1))/anginc +! if(abs(alon1).gt.175.)then +! print*,'Qingfu,lat,lon=',alat1,alon1,xlon(i)*rad2deg,i,it +! endif + coszn = sindec * sinlat(i) + cosdec * coslat(i) & - & * cos(cns+xlon(i)) + & * sindh +! & * cos(cns+xlon(i)) coszen(i) = coszen(i) + max(0.0, coszn) if (coszn > czlimt) istsun(i) = istsun(i) + 1 enddo From c580309a8065ff44df855146142f2f10c7001c4b Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Tue, 27 Jan 2026 16:39:58 +0000 Subject: [PATCH 09/19] Revert "add the repository update to this PR" This reverts commit 565d62cd9c7959e0981aabc2674686595fc73359. --- .../Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f | 28 +++++-------------- physics/Radiation/radiation_astronomy.f | 21 +------------- 2 files changed, 8 insertions(+), 41 deletions(-) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f index 9ec152b8d..749f778c1 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f @@ -265,7 +265,6 @@ subroutine dcyc2t3_run & & flxlwdn_adj real(kind=kind_phys) :: fluxlwnet_adj,fluxlwnet,dT_sfc, & &fluxlwDOWN_jac,lfnc,c1 - real(kind=kind_phys) :: pi2, xlon1, xlon2, sindh ! Length scale for flux-adjustment scaling real(kind=kind_phys), parameter :: & & L = 1. @@ -292,17 +291,15 @@ subroutine dcyc2t3_run & nstp = max(6, nint(tem1)) nstl = max(1, nint(nstp/tem1)) pid12 = con_pi / hour12 - pi2 = con_pi ! ! --- ... sw time-step adjustment for current cosine of zenith angle ! ---------------------------------------------------------- -! if (nstl == 1) then -! cns = pid12 * (solhr + deltim*f7200 - hour12) + slag -! do i = 1, IM -! xcosz(i) = sdec*sinlat(i) + cdec*coslat(i)*cos(cns+xlon(i)) -! enddo -! elseif (nstl == nstp) then - if (nstl == nstp) then + if (nstl == 1) then + cns = pid12 * (solhr + deltim*f7200 - hour12) + slag + do i = 1, IM + xcosz(i) = sdec*sinlat(i) + cdec*coslat(i)*cos(cns+xlon(i)) + enddo + elseif (nstl == nstp) then do i = 1, IM xcosz(i) = coszen(i) enddo @@ -317,17 +314,7 @@ subroutine dcyc2t3_run & do it=1,nstl cns = solang + (float(it)-0.5_kind_phys)*anginc + slag do i = 1, IM - xlon1 = cns+xlon(i) - if(xlon1.gt.pi2)then - xlon1=xlon1-2.*pi2 - else if(xlon1.lt.(0.-pi2))then - xlon1=xlon1+2.*pi2 - end if - xlon2 = xlon1 + anginc - if(xlon2.gt.pi2)xlon2 = xlon1 - sindh = (sin(xlon2)-sin(xlon1))/anginc -! coszn = sdec*sinlat(i) + cdec*coslat(i)*cos(cns+xlon(i)) - coszn = sdec*sinlat(i) + cdec*coslat(i)*sindh + coszn = sdec*sinlat(i) + cdec*coslat(i)*cos(cns+xlon(i)) xcosz(i) = xcosz(i) + max(zero, coszn) if (coszn > czlimt) istsun(i) = istsun(i) + 1 enddo @@ -380,7 +367,6 @@ subroutine dcyc2t3_run & xmu(i) = zero endif - print*,'qingfu, xmu=', i, xmu(i) !> - adjust \a sfc net and downward SW fluxes for zenith angle changes. ! note: sfc emiss effect will not be appied here diff --git a/physics/Radiation/radiation_astronomy.f b/physics/Radiation/radiation_astronomy.f index 0f1f41212..90ed7cd45 100644 --- a/physics/Radiation/radiation_astronomy.f +++ b/physics/Radiation/radiation_astronomy.f @@ -868,8 +868,6 @@ subroutine coszmn & ! --- locals: real (kind=kind_phys) :: coszn, cns, solang, rstp - real (kind=kind_phys) :: alat1, alon1, xlon1, xlon2 - real (kind=kind_phys) :: pi2, rad2deg, sindh integer :: istsun(IM), i, it, j, lat @@ -877,10 +875,7 @@ subroutine coszmn & solang = pid12 * (solhr - f12) ! solar angle at present time rstp = 1.0 / float(nstp) - rad2deg = 15. / pid12 - pi2 = pid12 * 12. -! print*,'qliu test=',nstp,IM,anginc*rad2deg do i = 1, IM coszen(i) = 0.0 istsun(i) = 0 @@ -889,22 +884,8 @@ subroutine coszmn & do it = 1, nstp cns = solang + (float(it)-0.5)*anginc + sollag do i = 1, IM - xlon1 = cns+xlon(i) - if(xlon1.gt.pi2)then - xlon1=xlon1-2.*pi2 - else if(xlon1.lt.(0.-pi2))then - xlon1=xlon1+2.*pi2 - end if - xlon2 = xlon1 + anginc - if(xlon2.gt.pi2)xlon2 = xlon1 - sindh = (sin(xlon2)-sin(xlon1))/anginc -! if(abs(alon1).gt.175.)then -! print*,'Qingfu,lat,lon=',alat1,alon1,xlon(i)*rad2deg,i,it -! endif - coszn = sindec * sinlat(i) + cosdec * coslat(i) & - & * sindh -! & * cos(cns+xlon(i)) + & * cos(cns+xlon(i)) coszen(i) = coszen(i) + max(0.0, coszn) if (coszn > czlimt) istsun(i) = istsun(i) + 1 enddo From 954b749773efa7bbce3c6776e810b8689ecd0c47 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Wed, 28 Jan 2026 14:37:46 +0000 Subject: [PATCH 10/19] remove the printing tests --- physics/Radiation/radiation_clouds.f | 1 + 1 file changed, 1 insertion(+) diff --git a/physics/Radiation/radiation_clouds.f b/physics/Radiation/radiation_clouds.f index 54c673e91..31a033084 100644 --- a/physics/Radiation/radiation_clouds.f +++ b/physics/Radiation/radiation_clouds.f @@ -230,6 +230,7 @@ module module_radiation_clouds & find_cloudLayers, adjust_cloudIce, adjust_cloudH2O, & & adjust_cloudFinal, gethml + public cld_frac_XuRandall ! ================= contains ! ================= From 49276b5bb0a29d558dde714640d8bca23c32e887 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Mon, 2 Feb 2026 14:48:36 +0000 Subject: [PATCH 11/19] remove Zhao-Carr schemes from modules in the UFS_SCM_NEPTUNE directory --- .../GFS_PBL_generic_common.F90 | 7 ++---- .../UFS_SCM_NEPTUNE/GFS_PBL_generic_post.F90 | 14 +++-------- .../UFS_SCM_NEPTUNE/GFS_PBL_generic_post.meta | 7 ------ .../UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.F90 | 16 +++---------- .../UFS_SCM_NEPTUNE/GFS_PBL_generic_pre.meta | 7 ------ .../UFS_SCM_NEPTUNE/GFS_rad_time_vary.fv3.F90 | 15 ++---------- .../GFS_rad_time_vary.fv3.meta | 7 ------ .../UFS_SCM_NEPTUNE/GFS_rad_time_vary.scm.F90 | 15 ++---------- .../GFS_rad_time_vary.scm.meta | 7 ------ .../UFS_SCM_NEPTUNE/GFS_rrtmg_pre.F90 | 18 ++++----------- .../UFS_SCM_NEPTUNE/GFS_rrtmg_pre.meta | 14 ----------- .../UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.F90 | 8 +++---- .../UFS_SCM_NEPTUNE/GFS_rrtmgp_setup.meta | 14 ----------- .../GFS_suite_interstitial_3.F90 | 18 ++------------- .../GFS_suite_interstitial_3.meta | 14 ----------- .../GFS_suite_interstitial_4.F90 | 23 +++++++------------ .../GFS_suite_interstitial_4.meta | 16 +------------ physics/Radiation/radiation_clouds.f | 12 ++-------- 18 files changed, 32 insertions(+), 200 deletions(-) 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_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/Radiation/radiation_clouds.f b/physics/Radiation/radiation_clouds.f index 31a033084..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, ! @@ -103,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 ! @@ -338,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, & @@ -357,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, ! @@ -427,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) ! @@ -513,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) :: & From 5ce727fd82a08c98480c85d195cfec57a060afc6 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Mon, 2 Feb 2026 14:52:14 +0000 Subject: [PATCH 12/19] remove Zhao-Carr microphysics scheme from physics/MP directory --- physics/MP/Zhao_Carr/zhaocarr_gscond.f | 530 ---------------- physics/MP/Zhao_Carr/zhaocarr_gscond.meta | 300 --------- physics/MP/Zhao_Carr/zhaocarr_precpd.f | 741 ---------------------- physics/MP/Zhao_Carr/zhaocarr_precpd.meta | 270 -------- 4 files changed, 1841 deletions(-) delete mode 100644 physics/MP/Zhao_Carr/zhaocarr_gscond.f delete mode 100644 physics/MP/Zhao_Carr/zhaocarr_gscond.meta delete mode 100644 physics/MP/Zhao_Carr/zhaocarr_precpd.f delete mode 100644 physics/MP/Zhao_Carr/zhaocarr_precpd.meta 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 - From 579eb9725a26b55e2a9afffc7ab507c5eabd249b Mon Sep 17 00:00:00 2001 From: "Michael Kavulich, Jr" Date: Tue, 3 Mar 2026 14:58:49 -0600 Subject: [PATCH 13/19] Update for 10m wind option in fire_behavior module --- .../Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 index 648f6bf81..24bc6f86c 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 @@ -117,7 +117,7 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplaqm, cplchm, cplwav, cpl v1(i) = vgrs_1(i) enddo - if (cplflx .or. cplchm .or. cplwav) then + if (cplflx .or. cplchm .or. cplwav .or. cpl_fire) then do i=1,im u10mi_cpl(i) = u10m(i) v10mi_cpl(i) = v10m(i) From 90f266d4926a85f8ff201addee16a8e63b9a5069 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Wed, 4 Mar 2026 14:33:24 -0500 Subject: [PATCH 14/19] remove non-ASCII character in canopy_levs.F90 (from previous PR) --- physics/PBL/SATMEDMF/canopy_levs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/PBL/SATMEDMF/canopy_levs.F90 b/physics/PBL/SATMEDMF/canopy_levs.F90 index 87a139d48..18cfa27d8 100644 --- a/physics/PBL/SATMEDMF/canopy_levs.F90 +++ b/physics/PBL/SATMEDMF/canopy_levs.F90 @@ -780,7 +780,7 @@ subroutine canopy_levs_run(im, km, nkc, nkt, & ! Note that these changes only exist inside the chemistry part of GEM-MACH and do not affect the model physics !!! !!! Create the momentum height (layer interface) array. The original momentum layers are used above the canopy height. -!!! Below the canopy height, the "momentum" layers are assumed to be ½ way between the thermodynamic layers. +!!! Below the canopy height, the "momentum" layers are assumed to be 1/2 way between the thermodynamic layers. ! Default case: all added canopy thermodynamic layers are below the lowest resolved model thermodynamic layer ! kcan_top is either 2nd or 3rd (63 or 62) resolved model layer From e673808397a1317cec46dfde6141b1349697e8cf Mon Sep 17 00:00:00 2001 From: Lisa Date: Fri, 6 Mar 2026 16:58:20 +0000 Subject: [PATCH 15/19] Updates to progomega option of saSAS convection scheme --- physics/CONV/SAMF/samfdeepcnv.f | 72 +++++++++++++++++++-------------- physics/CONV/SAMF/samfshalcnv.f | 69 +++++++++++++++++++------------ physics/CONV/progomega_calc.f90 | 34 +++++++++------- 3 files changed, 105 insertions(+), 70 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 18be662f0..a0d5a3a14 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -221,9 +221,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km), - & tentr(im,km) - real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins - parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01) + & wc_ref(im) + real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins, + & wc_min logical flag_shallow, flag_mid c physical parameters ! parameter(grav=grav,asolfac=0.958) @@ -347,6 +347,21 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !************************************************************************ ! ! +! - Initialize parameters related to prognostic closure + if (progsigma) then + if (progomega) then + sigmind = 0.03 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + else + sigmind = 0.01 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + endif + endif + km1 = km - 1 !> - Initialize column-integrated and other single-value-per-column variable arrays. c @@ -1137,7 +1152,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do i=1,im if(cnvflg(i) .and. & (k > kbcon(i) .and. k < kmax(i))) then - tentr(i,k)=xlamue(i,k)*fent1(i,k) tem = cxlamet(i) * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem tem1 = cxlamdt(i) * frh(i,k) @@ -1788,11 +1802,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if (progomega) then call progomega_calc(first_time_step,restart,im,km, & kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout, - & grav,buo,drag,wush,tentr,bb1,bb2) + & grav,buo,drag,wush,bb1,bb2) do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then omega_u(i,k)=omegaout(i,k) omega_u(i,k)=MAX(omega_u(i,k),-80.) ! Convert to m/s for use in convective time-scale: @@ -1809,7 +1823,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) @@ -1827,7 +1841,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then rho = po(i,k)*100. / (rd * to(i,k)) omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav omega_u(i,k)=MAX(omega_u(i,k),-80.) @@ -1878,7 +1892,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then dp = 1000. * del(i,k) tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1)) omegac(i) = omegac(i) + tem * dp @@ -2907,27 +2921,25 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! compute convective turn-over time ! !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, the convective adjustment time (dtconv) is set to be proportional to the convective turnover time, which is computed using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx). - if(hwrf_samfdeep) then - do i= 1, im - if(cnvflg(i)) then - tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) - dtconv(i) = tem / wc(i) - dtconv(i) = max(dtconv(i),dtmin) - dtconv(i) = min(dtconv(i),dtmax) - endif - enddo - else - do i= 1, im - if(cnvflg(i)) then - tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) - dtconv(i) = tem / wc(i) - tfac = 1. + gdx(i) / 75000. - dtconv(i) = tfac * dtconv(i) - dtconv(i) = max(dtconv(i),dtmin) - dtconv(i) = min(dtconv(i),dtmax) - endif - enddo - endif + do i = 1, im + if (cnvflg(i)) then + tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) + if (progomega) then + wc_eff = max(wc(i), wc_min) + dtconv(i) = tem / wc_eff + else + dtconv(i) = tem / wc(i) + endif + !grid spacing scaling (disabled for HWRF SAMF deep) + if (.not. hwrf_samfdeep) then + tfac = 1. + gdx(i) / 75000. + dtconv(i) = tfac * dtconv(i) + endif + !bounds + dtconv(i) = max(dtconv(i), dtmin) + dtconv(i) = min(dtconv(i), dtmax) + endif + enddo ! !> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. do i= 1, im diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index 1ec8747f0..f8cedc769 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -169,9 +169,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), - & sigmab(im),qadv(im,km) + & sigmab(im),qadv(im,km),wc_ref(im) real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins, - & sigminm + & sigminm,wc_min logical flag_shallow,flag_mid c physical parameters ! parameter(g=grav,asolfac=0.89) @@ -205,8 +205,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) parameter(betaw=.03,dxcrtc0=9.e3) parameter(h1=0.33333333) -! progsigma - parameter(dxcrtas=500.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01) c local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), & uo(im,km), vo(im,km), qeso(im,km), @@ -295,7 +293,21 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & prsl = prslp * 0.001 del = delp * 0.001 !************************************************************************ -! +! - Initialize parameters related to prognostic closure + if (progsigma) then + if (progomega) then + sigmind = 0.03 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + else + sigmind = 0.01 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + endif + endif +! km1 = km - 1 c c initialize arrays @@ -1520,11 +1532,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (progomega) then call progomega_calc(first_time_step,restart,im,km, & kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout, - & grav,buo,drag,wush,xlamue,bb1,bb2) + & grav,buo,drag,wush,bb1,bb2) do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then omega_u(i,k)=omegaout(i,k) omega_u(i,k)=MAX(omega_u(i,k),-80.) ! Convert to m/s for use in convective time-scale: @@ -1542,7 +1554,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) @@ -1560,7 +1572,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then rho = po(i,k)*100. / (rd * to(i,k)) omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav omega_u(i,k)=MAX(omega_u(i,k),-80.) @@ -1611,7 +1623,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then dp = 1000. * del(i,k) tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1)) omegac(i) = omegac(i) + tem * dp @@ -1636,7 +1648,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then if(omega_u(i,k) .ne. 0.)then zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) else @@ -1955,21 +1967,28 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! compute convective turn-over time ! !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, calculate the convective turnover time using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx). - do i= 1, im - if(cnvflg(i)) then - tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) - dtconv(i) = tem / wc(i) - if (.not.hwrf_samfshal) then - tfac = 1. + gdx(i) / 75000. - dtconv(i) = tfac * dtconv(i) - endif - dtconv(i) = max(dtconv(i),dtmin) - dtconv(i) = max(dtconv(i),dt2) - dtconv(i) = min(dtconv(i),dtmax) - endif + do i = 1, im + if (cnvflg(i)) then + tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) + if (progomega) then + wc_eff = max(wc(i), wc_min) + dtconv(i) = tem / wc_eff + else + dtconv(i) = tem / wc(i) + endif +! - grid spacing scaling (disabled for HWRF shallow option) + if (.not. hwrf_samfshal) then + tfac = 1. + gdx(i) / 75000. + dtconv(i) = tfac * dtconv(i) + endif +! - limits + dtconv(i) = max(dtconv(i), dtmin) + dtconv(i) = max(dtconv(i),dt2) + dtconv(i) = min(dtconv(i), dtmax) + endif enddo -! -!> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. +! +! > - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. do i= 1, im if(cnvflg(i)) then sumx(i) = 0. diff --git a/physics/CONV/progomega_calc.f90 b/physics/CONV/progomega_calc.f90 index 52d99d801..d22659156 100644 --- a/physics/CONV/progomega_calc.f90 +++ b/physics/CONV/progomega_calc.f90 @@ -20,7 +20,7 @@ module progomega !!\section gen_progomega progomega_calc General Algorithm subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegain,delt,del, & - zi,cnvflg,omegaout,grav,buo,drag,wush,tentr,bb1,bb2) + zi,cnvflg,omegaout,grav,buo,drag,wush,bb1,bb2) use machine, only : kind_phys use funcphys, only : fpvs @@ -30,20 +30,17 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai integer, intent(in) :: kbcon1(im),ktcon(im) real(kind=kind_phys), intent(in) :: delt,grav,bb1,bb2 real(kind=kind_phys), intent(in) :: omegain(im,km), del(im,km),zi(im,km) - real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km),tentr(im,km) + real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km) real(kind=kind_phys), intent(inout) :: omegaout(im,km) logical, intent(in) :: cnvflg(im),first_time_step,flag_restart real(kind=kind_phys) :: termA(im,km),termB(im,km),termC(im,km),omega(im,km) real(kind=kind_phys) :: RHS(im,km),Kd(im,km) - real(kind=kind_phys) :: dp,dz,entrn,Kdn,discr,wush_pa,lbb1,lbb2,lbb3 + real(kind=kind_phys) :: dp,dz,discr,wush_pa,lbb1,lbb2,lbb3 integer :: i,k - entrn = 0.8E-4 !0.5E-4 !m^-1 - Kdn = 0.5E-4 !2.9E-4 !m^-1 - lbb1 = 0.5 !1.0 - lbb2 = 3.2 !3.0 - lbb3 = 0.5 !0.5 - + lbb1 = 1.5 + lbb2 = 0.6 + lbb3 = 1.2 !Initialization 2D do k = 1,km @@ -56,6 +53,16 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai enddo enddo + do k = 1,km + do i = 1,im + if(cnvflg(i))then + if(omega(i,k) < 1.0E-5) then + omega(i,k) = 0. + endif + endif + enddo + enddo + if(first_time_step .and. .not. flag_restart)then do k = 1,km do i = 1,im @@ -76,10 +83,7 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai do k = 2, km do i = 1, im if (cnvflg(i)) then - if (k > kbcon1(i) .and. k < ktcon(i)) then - - ! Aerodynamic drag parameter - Kd(i,k) = (tentr(i,k)/entrn)*Kdn + if (k >= kbcon1(i) .and. k < ktcon(i)) then ! Scale by dp/dz to have equation in Pa/s !(dp/dz > 0) @@ -91,8 +95,8 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai !termC - Adds buoyancy forcing !Coefficients for the quadratic equation - termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz)) + (Kd(i,k) * (dp/dz))) - termB(i,k) = -1.0 - delt * lbb3 * wush(i,k) * dp/dz + termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz))) + termB(i,k) = 1.0 - delt * lbb3 * wush(i,k) * dp/dz termC(i,k) = omega(i,k) - delt * lbb2 * buo(i,k) * (dp/dz) & - delt * omega(i,k) * (omega(i,k-1) - omega(i,k)) / dp !Compute the discriminant From cd6a29de4e992d1f21e14c8773f25c9c5c2eeadc Mon Sep 17 00:00:00 2001 From: Lisa Date: Fri, 6 Mar 2026 17:54:42 +0000 Subject: [PATCH 16/19] Correction --- physics/CONV/SAMF/samfdeepcnv.f | 2 +- physics/CONV/SAMF/samfshalcnv.f | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index a0d5a3a14..4f6176a6d 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -221,7 +221,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km), - & wc_ref(im) + & wc_eff(im) real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins, & wc_min logical flag_shallow, flag_mid diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index f8cedc769..ba2924297 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -169,7 +169,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), - & sigmab(im),qadv(im,km),wc_ref(im) + & sigmab(im),qadv(im,km),wc_eff(im) real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins, & sigminm,wc_min logical flag_shallow,flag_mid From c2991ab0daa1919b71942a2ae4230816289e0537 Mon Sep 17 00:00:00 2001 From: Lisa Date: Fri, 6 Mar 2026 18:02:53 +0000 Subject: [PATCH 17/19] correction 2 --- physics/CONV/SAMF/samfdeepcnv.f | 5 ++--- physics/CONV/SAMF/samfshalcnv.f | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 4f6176a6d..651829379 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -220,10 +220,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), - & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km), - & wc_eff(im) + & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km) real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins, - & wc_min + & wc_min, wc_eff logical flag_shallow, flag_mid c physical parameters ! parameter(grav=grav,asolfac=0.958) diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index ba2924297..e26ccacb1 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -169,9 +169,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), - & sigmab(im),qadv(im,km),wc_eff(im) + & sigmab(im),qadv(im,km) real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins, - & sigminm,wc_min + & sigminm,wc_min,wc_eff logical flag_shallow,flag_mid c physical parameters ! parameter(g=grav,asolfac=0.89) From f2fbb17a330938f51a3d7cc1a54a184661dc5570 Mon Sep 17 00:00:00 2001 From: Lisa Date: Fri, 20 Mar 2026 18:45:44 +0000 Subject: [PATCH 18/19] Correct parameter setting in samfshalcnv to pass debug runs --- physics/CONV/SAMF/samfshalcnv.f | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index e26ccacb1..55cf6f7f4 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -252,7 +252,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) tf, tcr, tcrf parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) - c----------------------------------------------------------------------- ! ! Initialize CCPP error handling variables @@ -274,8 +273,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (progsigma) then dxcrt=10.e3 + dxcrtas=500.e3 else dxcrt=15.e3 + dxcrtas=500.e3 endif c----------------------------------------------------------------------- From e5ebf607cd342b1125941ea4e843a064503008f6 Mon Sep 17 00:00:00 2001 From: Lisa Date: Mon, 23 Mar 2026 16:10:45 +0000 Subject: [PATCH 19/19] address saturation vapor pressure out of bounds error in hafs_4_nest RT --- physics/CONV/SAMF/samfdeepcnv.f | 4 ++-- physics/CONV/SAMF/samfshalcnv.f | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 651829379..2acf3b55b 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -1822,7 +1822,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) @@ -1840,7 +1840,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then rho = po(i,k)*100. / (rd * to(i,k)) omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav omega_u(i,k)=MAX(omega_u(i,k),-80.) diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index 55cf6f7f4..43d7eeec2 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -1555,7 +1555,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) @@ -1573,7 +1573,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then rho = po(i,k)*100. / (rd * to(i,k)) omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav omega_u(i,k)=MAX(omega_u(i,k),-80.)