From 235ec3825715dfb646afe00dc71b0c010ff9b661 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 13 Apr 2022 17:44:27 +0000 Subject: [PATCH 01/11] add progsigma_calc --- physics/progsigma_calc.f90 | 260 +++++++++++++++++++++++++++++++++++++ 1 file changed, 260 insertions(+) create mode 100644 physics/progsigma_calc.f90 diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 new file mode 100644 index 000000000..378d43ef4 --- /dev/null +++ b/physics/progsigma_calc.f90 @@ -0,0 +1,260 @@ +!>\file progsigma +!! This file contains the subroutine that calculates the prognostic +!! updraft area fraction that is used for closure computations in +!! saSAS deep and shallow convection. + +!>\ingroup samfdeepcnv +!! This subroutine computes a prognostic updraft area fraction +!! used in the closure computations in the samfdeepcnv.f scheme +!>\ingroup samfshalcnv +!! This subroutine computes a prognostic updraft area fracftion +!! used in the closure computations in the samfshalcnv. scheme +!!\section progsigma General Algorithm +!> @{ + + subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & + del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & + qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & + do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & + ca_turb,ca_micro,ca_shal,ca_rad,convcount,ca1,ca2,ca3,ca4, & + sigmain,sigmaout,sigmab,errmsg,errflg) +! +! + use machine, only : kind_phys + use funcphys, only : fpvs + + implicit none + +! intent in + integer, intent(in) :: im,km,kbcon1(im),ktcon(im) + real, intent(in) :: hvap,delt + real, intent(in) :: qgrs_dsave(im,km), q(im,km),del(im,km), & + qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & + omega_u(im,km),zeta(im,km),gdx(im) + logical, intent(in) :: flag_init,flag_restart,flag_deep,cnvflg(im) + real(kind=kind_phys), intent(in) :: nthresh + real(kind=kind_phys), intent(in) :: ca_deep(im) + real(kind=kind_phys), intent(out):: ca_turb(im), & + ca_micro(im),ca_rad(im),ca_shal(im),convcount(im),ca1(im), & + ca2(im),ca3(im),ca4(im) + logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger + + real(kind=kind_phys), intent(in) :: sigmain(im,km) + +! intent out + real(kind=kind_phys), intent(out) :: sigmaout(im,km) + real(kind=kind_phys), intent(out) :: sigmab(im) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + +! Local variables + integer :: i,k,km1 + real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & + mcons(im),zfdqa(im),zform(im,km), & + qadv(im,km),sigmamax(im) + + + real(kind=kind_phys) :: gcvalmx,ZEPS7,ZZ,ZCVG,mcon,buy2, & + zfdqb,dtdyn,dxlim,rmulacvg,dp,tem, & + alpha,DEN + integer :: inbu(im,km) + + !Parameters + gcvalmx = 0.1 + rmulacvg=10. + ZEPS7=1.E-11 + km1=km-1 + alpha=7000. + + !Initialization 2D + do k = 1,km + do i = 1,im + sigmaout(i,k)=0. + inbu(i,k)=0 + zform(i,k)=0. + enddo + enddo + + !Initialization 1D + do i=1,im + sigmab(i)=0. + sigmamax(i)=0.95 + termA(i)=0. + termB(i)=0. + termC(i)=0. + termD(i)=0. + zfdqa(i)=0. + mcons(i)=0. + enddo + + !Temporary Initialization output: + do i = 1,im + if(flag_deep)then + !ca_turb(i)=0. + ca_shal(i)=0. + endif + if(.not. flag_deep)then + ca_rad(i)=0. + convcount(i)=0. + ca1(i)=0. + endif + enddo + + !Initial computations, place maximum sigmain in sigmab + + do k=2,km + do i=1,im + if(flag_init .and. .not. flag_restart)then + if(cnvflg(i))then + sigmab(i)=0.03 + endif + else + if(cnvflg(i))then + !if(sigmain(i,k)<1.E-5)then + ! sigmain(i,k)=0. + !endif + if(sigmain(i,k)>sigmab(i))then + sigmab(i)=sigmain(i,k) + endif + endif + endif + enddo + enddo + + do i=1,im + if(sigmab(i) < 1.E-5)then !after advection + sigmab(i)=0. + endif + enddo + + !Initial computations, sigmamax + do i=1,im + sigmamax(i)=alpha/gdx(i) + sigmamax(i)=MIN(0.95,sigmamax(i)) + enddo + + !Initial computations, dynamic q-tendency + do k = 1,km + do i = 1,im + if(flag_init .and. .not.flag_restart)then + qadv(i,k)=0. + else + qadv(i,k)=(q(i,k) - qgrs_dsave(i,k))/delt + endif + enddo + enddo + + !compute termD "The vertical integral of the latent heat convergence is limited to the + !buoyant layers with positive moisture convergence (accumulated from the surface). + !Lowest level: + do i = 1,im + dp = 1000. * del(i,1) + mcons(i)=(hvap*(qadv(i,1)+tmf(i,1)+qmicro(i,1))*dp) + enddo + !Levels above: + do k = 2,km1 + do i = 1,im + dp = 1000. * del(i,k) + if(cnvflg(i))then + mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp) + buy2 = termD(i)+mcon+mcons(i) +! Do the integral over buoyant layers with positive mcon acc from surface + if(k > kbcon1(i) .and. k < ktcon(i) .and. buy2 > 0.)then + inbu(i,k)=1 + endif + inbu(i,k-1)=MAX(inbu(i,k-1),inbu(i,k)) + termD(i) = termD(i) + float(inbu(i,k-1))*mcons(i) + mcons(i)=mcon + endif + enddo + enddo + + !termA + do k = 2,km1 + do i = 1,im + dp = 1000. * del(i,k) + if(cnvflg(i))then + tem=(sigmab(i)*zeta(i,k)*float(inbu(i,k))*dbyo1(i,k))*dp + termA(i)=termA(i)+tem + endif + enddo + enddo + + !termB + do k = 2,km1 + do i = 1,im + dp = 1000. * del(i,k) + if(cnvflg(i))then + tem=(dbyo1(i,k)*float(inbu(i,k)))*dp + termB(i)=termB(i)+tem + endif + enddo + enddo + + !termC + do k = 2,km1 + do i = 1,im + if(cnvflg(i))then + dp = 1000. * del(i,k) + zform(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt) + zfdqb=0.5*((zform(i,k)*zdqca(i,k))) + termC(i)=termC(i)+(float(inbu(i,k))* & + (zfdqb+zfdqa(i))*hvap*zeta(i,k)) + zfdqa(i)=zfdqb + endif + enddo + enddo + + !sigmab + do i = 1,im + if(cnvflg(i))then + + DEN=MIN(termC(i)+termB(i),1.E8) !1.E8 + !DEN=MAX(termC(i)+termB(i),1.E7) !1.E7 + + ZCVG=termD(i)*delt + + ZZ=MAX(0.0,SIGN(1.0,termA(i))) & + *MAX(0.0,SIGN(1.0,termB(i))) & + *MAX(0.0,SIGN(1.0,termC(i)-ZEPS7)) + + + ZCVG=MAX(0.0,ZCVG) + + if(flag_init)then + sigmab(i)=0.03 + else + sigmab(i)=(ZZ*(termA(i)+ZCVG))/(DEN+(1.0-ZZ)) + endif + + if(sigmab(i)>0.)then + sigmab(i)=MIN(sigmab(i),sigmamax(i)) + sigmab(i)=MAX(sigmab(i),0.01) + endif + + if(flag_deep)then + !ca_turb(i)=ZCVG + ca_shal(i)=termC(i) + else + ca_rad(i)=ZCVG + ca1(i)=termC(i) + endif + !ca3(i)=sigmab(i) + + endif!cnvflg + enddo + + do k=1,km + do i=1,im + if(cnvflg(i))then + sigmaout(i,k)=sigmab(i) + endif + enddo + enddo + + end subroutine progsigma_calc +!> @} +!! @} + + + From 842eae33944aafb5100ccb7660a04027a7452147 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Mon, 18 Apr 2022 15:05:40 +0000 Subject: [PATCH 02/11] Ensuring the moisture budget is correct via PBL, microphysics coupling --- physics/GFS_MP_generic.F90 | 26 +++++-- physics/GFS_MP_generic.meta | 23 ++++++ physics/progsigma_calc.f90 | 33 +------- physics/samfdeepcnv.f | 149 +++++++++++++++++++++++++++++++----- physics/samfdeepcnv.meta | 86 ++++++++++++++++++++- physics/satmedmfvdifq.F | 5 +- physics/satmedmfvdifq.meta | 8 ++ 7 files changed, 274 insertions(+), 56 deletions(-) diff --git a/physics/GFS_MP_generic.F90 b/physics/GFS_MP_generic.F90 index e106cb908..dbf2d15fa 100644 --- a/physics/GFS_MP_generic.F90 +++ b/physics/GFS_MP_generic.F90 @@ -88,14 +88,15 @@ end subroutine GFS_MP_generic_post_init !> @{ subroutine GFS_MP_generic_post_run( & im, levs, kdt, nrcm, nncl, ntcw, ntrac, imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_nssl, & - imp_physics_mg, imp_physics_fer_hires, cal_pre, cplflx, cplchm, con_g, rainmin, dtf, frain, rainc, & + imp_physics_mg, imp_physics_fer_hires, cal_pre, cplflx, cplchm, progsigma, con_g, rainmin, dtf, frain, rainc, & rain1, rann, xlat, xlon, gt0, gq0, prsl, prsi, phii, tsfc, ice, snow, graupel, save_t, save_q, rain0, ice0, snow0,& graupel0, del, rain, domr_diag, domzr_diag, domip_diag, doms_diag, tprcp, srflag, sr, cnvprcp, totprcp, totice, & totsnw, totgrp, cnvprcpb, totprcpb, toticeb, totsnwb, totgrpb, rain_cpl, rainc_cpl, snow_cpl, pwat, & drain_cpl, dsnow_cpl, lsm, lsm_ruc, lsm_noahmp, raincprv, rainncprv, iceprv, snowprv, & graupelprv, draincprv, drainncprv, diceprv, dsnowprv, dgraupelprv, dtp, dfi_radar_max_intervals, & - dtend, dtidx, index_of_temperature, index_of_process_mp,ldiag3d, qdiag3d, lssav, num_dfi_radar, fh_dfi_radar, & - index_of_process_dfi_radar, ix_dfi_radar, dfi_radar_tten, radar_tten_limits, fhour, errmsg, errflg) + dtend, dtidx, index_of_temperature, index_of_process_mp,ldiag3d, qdiag3d,dqdt_qmicro, lssav, num_dfi_radar, & + fh_dfi_radar,index_of_process_dfi_radar, ix_dfi_radar, dfi_radar_tten, radar_tten_limits, fhour, qgrs_dsave, & + errmsg, errflg) ! use machine, only: kind_phys @@ -104,7 +105,7 @@ subroutine GFS_MP_generic_post_run( integer, intent(in) :: im, levs, kdt, nrcm, nncl, ntcw, ntrac, num_dfi_radar, index_of_process_dfi_radar integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_mg, imp_physics_fer_hires integer, intent(in) :: imp_physics_nssl - logical, intent(in) :: cal_pre, lssav, ldiag3d, qdiag3d, cplflx, cplchm + logical, intent(in) :: cal_pre, lssav, ldiag3d, qdiag3d, cplflx, cplchm, progsigma integer, intent(in) :: index_of_temperature,index_of_process_mp integer :: dfi_radar_max_intervals @@ -148,7 +149,8 @@ subroutine GFS_MP_generic_post_run( real(kind=kind_phys), dimension(:), intent(inout) :: diceprv real(kind=kind_phys), dimension(:), intent(inout) :: dsnowprv real(kind=kind_phys), dimension(:), intent(inout) :: dgraupelprv - + real(kind=kind_phys), dimension(:,:), intent(out) :: dqdt_qmicro + real(kind=kind_phys), dimension(:,:), intent(out) :: qgrs_dsave real(kind=kind_phys), intent(in) :: dtp ! CCPP error handling @@ -420,6 +422,15 @@ subroutine GFS_MP_generic_post_run( endif if_tendency_diagnostics endif if_save_fields + !If prognostic updraft area fraction is used in saSAS + if(progsigma)then + do k=1,levs + do i=1,im + dqdt_qmicro(i,k)=(gq0(i,k,1)-save_q(i,k,1))/dtp + enddo + enddo + endif + if (cplflx .or. cplchm) then do i = 1, im dsnow_cpl(i)= max(zero, rain(i) * srflag(i)) @@ -455,6 +466,11 @@ subroutine GFS_MP_generic_post_run( pwat(i) = pwat(i) * onebg enddo + do k = 1, levs + do i=1, im + qgrs_dsave(i,k) = gq0(i,k,1) + enddo + enddo end subroutine GFS_MP_generic_post_run !> @} diff --git a/physics/GFS_MP_generic.meta b/physics/GFS_MP_generic.meta index 6177b1344..f8cb0acae 100644 --- a/physics/GFS_MP_generic.meta +++ b/physics/GFS_MP_generic.meta @@ -248,6 +248,13 @@ dimensions = () type = logical intent = in +[progsigma] + standard_name = flag_for_prognostic_sigma + long_name = flag for prognostic sigma + units = flag + dimensions = () + type = logical + intent = in [con_g] standard_name = gravitational_acceleration long_name = gravitational acceleration @@ -844,6 +851,22 @@ dimensions = () type = logical intent = in +[dqdt_qmicro] + standard_name = instantanious_moisture_tendency_due_to_microphysics + long_name = moisture tendency due to microphysics + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out +[qgrs_dsave] + standard_name = tracer_concentration_dsave + long_name = model layer mean tracer concentration dsave + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [lssav] standard_name = flag_for_diagnostics long_name = logical flag for storing diagnostics diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 378d43ef4..6772bc8d4 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -12,12 +12,11 @@ !!\section progsigma General Algorithm !> @{ - subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & + subroutine progsigma_calc (im,km,flag_init,flag_restart, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - ca_turb,ca_micro,ca_shal,ca_rad,convcount,ca1,ca2,ca3,ca4, & - sigmain,sigmaout,sigmab,errmsg,errflg) + ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) ! ! use machine, only : kind_phys @@ -31,12 +30,10 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & real, intent(in) :: qgrs_dsave(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) - logical, intent(in) :: flag_init,flag_restart,flag_deep,cnvflg(im) + logical, intent(in) :: flag_init,flag_restart,cnvflg(im) real(kind=kind_phys), intent(in) :: nthresh real(kind=kind_phys), intent(in) :: ca_deep(im) - real(kind=kind_phys), intent(out):: ca_turb(im), & - ca_micro(im),ca_rad(im),ca_shal(im),convcount(im),ca1(im), & - ca2(im),ca3(im),ca4(im) + real(kind=kind_phys), intent(out):: ca_micro(im) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger real(kind=kind_phys), intent(in) :: sigmain(im,km) @@ -87,19 +84,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & mcons(i)=0. enddo - !Temporary Initialization output: - do i = 1,im - if(flag_deep)then - !ca_turb(i)=0. - ca_shal(i)=0. - endif - if(.not. flag_deep)then - ca_rad(i)=0. - convcount(i)=0. - ca1(i)=0. - endif - enddo - !Initial computations, place maximum sigmain in sigmab do k=2,km @@ -232,15 +216,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & sigmab(i)=MAX(sigmab(i),0.01) endif - if(flag_deep)then - !ca_turb(i)=ZCVG - ca_shal(i)=termC(i) - else - ca_rad(i)=ZCVG - ca1(i)=termC(i) - endif - !ca3(i)=sigmab(i) - endif!cnvflg enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index ea92fda7f..7be31a04a 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -75,17 +75,18 @@ end subroutine samfdeepcnv_finalize !! !! \section samfdeep_detailed GFS samfdeepcnv Detailed Algorithm !! @{ - subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & + subroutine samfdeepcnv_run (im,km,first_time_step,restart, & + & tmf,qmicro,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,q1,t1,u1,v1,fscav,hwrf_samfdeep, & - & cldwrk,rn,kbot,ktop,kcnv,islimsk,garea, & - & dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & + & prslp,psp,phil,qtr,qgrs_dsave,q,q1,t1,u1,v1,fscav, & + & hwrf_samfdeep,progsigma,wclosureflg,cldwrk,rn,kbot,ktop,kcnv, & + & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, & + & rainevap, sigmain, sigmaout, ca_micro, & & errmsg,errflg) ! use machine , only : kind_phys @@ -101,10 +102,14 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:) real(kind=kind_phys), dimension(:), intent(in) :: fscav - logical, intent(in) :: hwrf_samfdeep + logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, & + & progsigma, wclosureflg real(kind=kind_phys), intent(in) :: nthresh real(kind=kind_phys), intent(in) :: ca_deep(:) - real(kind=kind_phys), intent(out) :: rainevap(:) + real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & + & tmf(:,:),q(:,:), qgrs_dsave(:,:) + real(kind=kind_phys), intent(out) :: rainevap(:),ca_micro(:) + real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger integer, intent(inout) :: kcnv(:) @@ -208,6 +213,10 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! & bb1, bb2 & bb1, bb2, wucb ! +! parameters for prognostic sigma closure + real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), + & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) + c physical parameters ! parameter(grav=grav,asolfac=0.958) ! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) @@ -368,6 +377,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & vshear(i) = 0. advfac(i) = 0. rainevap(i) = 0. + omegac(i)=0. gdx(i) = sqrt(garea(i)) enddo @@ -570,6 +580,11 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & buo(i,k) = 0. drag(i,k) = 0. cnvwt(i,k)= 0. + dbyo1(i,k)=0. + zdqca(i,k)=0. + qlks(i,k)=0. + omega_u(i,k)=0. + zeta(i,k)=1.0 endif enddo enddo @@ -1497,6 +1512,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & pwavo(i) = pwavo(i) + pwo(i,k) ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp cnvwt(i,k) = etah * qlk * grav / dp + qlks(i,k)=qlk endif ! ! compute buoyancy and drag for updraft velocity @@ -1569,6 +1585,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & dz1 = zo(i,k+1) - zo(i,k) ! aa1(i) = aa1(i) + buo(i,k) * dz1 * eta(i,k) aa1(i) = aa1(i) + buo(i,k) * dz1 + dbyo1(i,k) = hcko(i,k) - heso(i,k) endif endif enddo @@ -1669,6 +1686,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & pwavo(i) = pwavo(i) + pwo(i,k) ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp cnvwt(i,k) = etah * qlk * grav / dp + qlks(i,k)=qlk endif endif endif @@ -1710,6 +1728,20 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + + if(progsigma)then + do k = 2, km1 + do i = 1, im + if (cnvflg(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.) + endif + endif + enddo + enddo + endif ! ! compute updraft velocity average over the whole cumulus ! @@ -1742,6 +1774,54 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif enddo c + +!> - Calculate the mean updraft velocity within the cloud (wc),cast in pressure coordinates. + if(progsigma)then + + do i = 1, im + omegac(i) = 0. + sumx(i) = 0. + enddo + do k = 2, km1 + do i = 1, im + if (cnvflg(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 + sumx(i) = sumx(i) + dp + endif + endif + enddo + enddo + do i = 1, im + if(cnvflg(i)) then + if(sumx(i) == 0.) then + cnvflg(i)=.false. + else + omegac(i) = omegac(i) / sumx(i) + endif + val = -1.2 + if (omegac(i) > val) cnvflg(i)=.false. + endif + enddo + +!> - Calculate the xi term in Bengtsson et al. 2022 (equation 8) + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then + zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + zeta(i,k)=MAX(0.,zeta(i,k)) + zeta(i,k)=MIN(1.,zeta(i,k)) + endif + endif + enddo + enddo + + + endif !if progsigma + c exchange ktcon with ktcon1 c !> - Swap the indices of the convective cloud top (ktcon) and the overshooting convection top (ktcon1) to use the same cloud top level in the calculations of \f$A^+\f$ and \f$A^*\f$. @@ -1773,11 +1853,26 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if(dq > 0.) then qlko_ktcon(i) = dq qcko(i,k) = qrch + qlks(i,k) = qlko_ktcon(i) endif endif enddo endif c + +c store term needed for "termC" in prognostic area fraction closure + do k = 2, km1 + do i = 1, im + dp = 1000. * del(i,k) + if (cnvflg(i)) then + if(k > kbcon(i) .and. k < ktcon(i)) then + zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + + & pwo(i,k)+dellal(i,k)) + endif + endif + enddo + enddo + ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then ccccc print *, ' aa1(i) before dwndrft =', aa1(i) ccccc endif @@ -2375,6 +2470,14 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & asqecflg(i) = .false. endif enddo + +!> - If wclosureflg is true, then quasi-equilibrium closure of Arakawa-Schubert is not used any longer, regardless of resolution + if(wclosureflg)then + do i = 1, im + asqecflg(i) = .false. + enddo + endif + ! !> - If grid size is larger than the threshold value (i.e., asqecflg=.true.), the quasi-equilibrium assumption is used to obtain the cloud base mass flux. To begin with, calculate the change in the temperature and moisture profiles per unit cloud base mass flux. do k = 1, km @@ -2784,13 +2887,27 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & advfac(i) = min(advfac(i), 1.) endif enddo + +!> - From Bengtsson et al. (2022) Prognostic closure scheme, equation 8, compute updraft area fraction based on a moisture budget + if(progsigma)then + call progsigma_calc(im,km,first_time_step,restart, + & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, + & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, + & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, + & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + endif + !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) - xmb(i) = advfac(i)*betaw*rho*wc(i) + if(progsigma)then + xmb(i) = sigmab(i)*((-1.0*omegac(i))/grav) + else + xmb(i) = advfac(i)*betaw*rho*wc(i) + endif endif enddo !> - For the cases where the quasi-equilibrium assumption of Arakawa-Schubert is valid, first calculate the large scale destabilization as in equation 5 of Pan and Wu (1995) \cite pan_and_wu_1995 : @@ -2859,7 +2976,11 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then if (gdx(i) < dxcrtuf) then - scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) + if(progsigma)then + scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i)) + else + scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) + endif scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) else scaldfunc(i) = 1.0 @@ -2869,16 +2990,6 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif enddo ! - if (do_ca .and. ca_closure)then - do i = 1, im - if(cnvflg(i)) then - if (ca_deep(i) > nthresh) then - xmb(i) = xmb(i)*1.25 - endif - endif - enddo - endif - !> - Transport aerosols if present ! ! if (do_aerosols) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index baf01fb8e..3eb330551 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = samfdeepcnv type = scheme - dependencies = funcphys.f90,machine.F,samfaerosols.F + dependencies = funcphys.f90,machine.F,samfaerosols.F,progsigma_calc.f90 ######################################################################## [ccpp-arg-table] @@ -55,6 +55,36 @@ dimensions = () type = integer intent = in +[first_time_step] + standard_name = flag_for_first_timestep + long_name = flag for first time step for time integration loop (cold/warmstart) + units = flag + dimensions = () + type = logical + intent = in +[restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in +[tmf] + standard_name = turbulence_moisture_flux + long_name = turbulence_moisture_flux + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[qmicro] + standard_name = instantanious_moisture_tendency_due_to_microphysics + long_name = moisture tendency due to microphysics + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [itc] standard_name = index_of_first_chemical_tracer_for_convection long_name = index of first chemical tracer transported/scavenged by convection @@ -219,6 +249,22 @@ type = real kind = kind_phys intent = inout +[qgrs_dsave] + standard_name = tracer_concentration_dsave + long_name = model layer mean tracer concentration dsave + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[q] + standard_name = specific_humidity + long_name = water vapor specific humidity + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [q1] standard_name = specific_humidity_of_new_state long_name = updated vapor specific humidity @@ -266,6 +312,28 @@ dimensions = () type = logical intent = in +[wclosureflg] + standard_name = flag_for_wclosure + long_name = flag for vertical velocity closure + units = flag + dimensions = () + type = logical + intent = in +[progsigma] + standard_name = flag_for_prognostic_sigma + long_name = flag for prognostic sigma + units = flag + dimensions = () + type = logical + intent = in +[ca_micro] + standard_name = output_prognostic_sigma_two + long_name = output of prognostic area fraction two + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [cldwrk] standard_name = cloud_work_function long_name = cloud work function @@ -381,6 +449,22 @@ type = real kind = kind_phys intent = inout +[sigmain] + standard_name = updraft_area_fraction + long_name = convective updraft area fraction + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[sigmaout] + standard_name = updraft_area_fraction_updated_by_physics + long_name = convective updraft area fraction updated by physics + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [qlcn] standard_name = mass_fraction_of_convective_cloud_liquid_water long_name = mass fraction of convective cloud liquid water diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index eb2b7ad1c..b0367b4d4 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -69,7 +69,7 @@ end subroutine satmedmfvdifq_finalize !! @{ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & & grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & - & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, & + & dv,du,tdt,rtg,tmf,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & & tsea,heat,evap,stress,spd1,kpbl, & & prsi,del,prsl,prslk,phii,phil,delt, & @@ -121,7 +121,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & real(kind=kind_phys), intent(out) :: & & dusfc(:), dvsfc(:), & & dtsfc(:), dqsfc(:), & - & hpbl(:) + & hpbl(:), tmf(:,:) real(kind=kind_phys), intent(out) :: & & dkt(:,:), dku(:,:) ! @@ -2114,6 +2114,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & qtend = (f2(i,k)-q1(i,k,1))*rdt tdt(i,k) = tdt(i,k)+ttend rtg(i,k,1) = rtg(i,k,1)+qtend + tmf(i,k) = qtend ! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend ! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend enddo diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index db89f488d..211bbaec6 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -201,6 +201,14 @@ type = real kind = kind_phys intent = inout +[tmf] + standard_name = turbulence_moisture_flux + long_name = turbulence_moisture_flux + units = kg kg-1 s-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [u1] standard_name = x_wind long_name = x component of layer wind From 4f84ed749af88d8d833d78ec352c2a02dfb52589 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Tue, 19 Apr 2022 19:22:53 +0000 Subject: [PATCH 03/11] add shallow convection closure updates, add ntsigma in generic files --- physics/GFS_DCNV_generic.F90 | 14 ++-- physics/GFS_DCNV_generic.meta | 14 ++++ physics/GFS_MP_generic.meta | 2 +- physics/GFS_SCNV_generic.F90 | 12 +-- physics/GFS_SCNV_generic.meta | 14 ++++ physics/progsigma_calc.f90 | 73 +++++++++--------- physics/samfdeepcnv.f | 31 ++++---- physics/samfdeepcnv.meta | 10 +-- physics/samfshalcnv.f | 139 +++++++++++++++++++++++++++++++--- physics/samfshalcnv.meta | 79 ++++++++++++++++++- physics/satmedmfvdifq.F | 5 +- physics/satmedmfvdifq.meta | 6 +- 12 files changed, 315 insertions(+), 84 deletions(-) diff --git a/physics/GFS_DCNV_generic.F90 b/physics/GFS_DCNV_generic.F90 index a9e0ba7e0..07defcba1 100644 --- a/physics/GFS_DCNV_generic.F90 +++ b/physics/GFS_DCNV_generic.F90 @@ -18,8 +18,8 @@ end subroutine GFS_DCNV_generic_pre_finalize subroutine GFS_DCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, do_cnvgwd, cplchm, & gu0, gv0, gt0, gq0, nsamftrac, ntqv, & save_u, save_v, save_t, save_q, clw, & - ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl, & - ntgnc, nthl, nthnc, nthv, ntgv, & + ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,& + ntgl,ntgnc, nthl, nthnc, nthv, ntgv, & cscnv, satmedmf, trans_trac, ras, ntrac, & dtidx, index_of_process_dcnv, errmsg, errflg) @@ -28,7 +28,7 @@ subroutine GFS_DCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, do_cnvgwd, cplc implicit none integer, intent(in) :: im, levs, nsamftrac, ntqv, index_of_process_dcnv, dtidx(:,:), & - ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntrac,ntgnc,nthl,nthnc,nthv,ntgv + ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntrac,ntgnc,nthl,nthnc,nthv,ntgv,ntsigma logical, intent(in) :: ldiag3d, qdiag3d, do_cnvgwd, cplchm real(kind=kind_phys), dimension(:,:), intent(in) :: gu0 real(kind=kind_phys), dimension(:,:), intent(in) :: gv0 @@ -74,7 +74,7 @@ subroutine GFS_DCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, do_cnvgwd, cplc n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & n /= ntsnc .and. n /= ntgl .and. n /= ntgnc .and. & n /= nthl .and. n /= nthnc .and. n /= nthv .and. & - n /= ntgv ) then + n /= ntgv .and. n /= ntsigma) then tracers = tracers + 1 if(dtidx(100+n,index_of_process_dcnv)>0) then save_q(:,:,n) = clw(:,:,tracers) @@ -114,7 +114,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, & rainc, cldwrk, upd_mf, dwn_mf, det_mf, dtend, dtidx, index_of_process_dcnv, & index_of_temperature, index_of_x_wind, index_of_y_wind, ntqv, gq0, save_q, & cnvw, cnvc, cnvw_phy_f3d, cnvc_phy_f3d, flag_for_dcnv_generic_tend, & - ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl, & + ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl, & ntgnc, nthl, nthnc, nthv, ntgv, ntrac,clw, & satmedmf, trans_trac, errmsg, errflg) @@ -145,7 +145,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, & integer, intent(in) :: dtidx(:,:), index_of_process_dcnv, index_of_temperature, & index_of_x_wind, index_of_y_wind, ntqv integer, intent(in) :: ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl, & - ntgnc, nthl, nthnc, nthv, ntgv, ntrac + ntgnc, nthl, nthnc, nthv, ntgv, ntsigma, ntrac real(kind=kind_phys), dimension(:,:,:), intent(in) :: clw @@ -212,7 +212,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, & n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & n /= ntsnc .and. n /= ntgl .and. n /= ntgnc .and. & n /= nthl .and. n /= nthnc .and. n /= nthv .and. & - n /= ntgv ) then + n /= ntgv .and. n /= ntsigma) then tracers = tracers + 1 idtend = dtidx(100+n,index_of_process_dcnv) if(idtend>0) then diff --git a/physics/GFS_DCNV_generic.meta b/physics/GFS_DCNV_generic.meta index e15acaf1c..f095259d4 100644 --- a/physics/GFS_DCNV_generic.meta +++ b/physics/GFS_DCNV_generic.meta @@ -190,6 +190,13 @@ dimensions = () type = integer intent = in +[ntsigma] + standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array + long_name = tracer index of updraft_area_fraction + units = index + dimensions = () + type = integer + intent = in [ntrw] standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array long_name = tracer index for rain water @@ -670,6 +677,13 @@ dimensions = () type = integer intent = in +[ntsigma] + standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array + long_name = tracer index of updraft_area_fraction + units = index + dimensions = () + type = integer + intent = in [ntrw] standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array long_name = tracer index for rain water diff --git a/physics/GFS_MP_generic.meta b/physics/GFS_MP_generic.meta index f8cb0acae..763cad85a 100644 --- a/physics/GFS_MP_generic.meta +++ b/physics/GFS_MP_generic.meta @@ -852,7 +852,7 @@ type = logical intent = in [dqdt_qmicro] - standard_name = instantanious_moisture_tendency_due_to_microphysics + standard_name = instantaneous_moisture_tendency_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) diff --git a/physics/GFS_SCNV_generic.F90 b/physics/GFS_SCNV_generic.F90 index 58447f6bf..cbef02bb0 100644 --- a/physics/GFS_SCNV_generic.F90 +++ b/physics/GFS_SCNV_generic.F90 @@ -17,14 +17,14 @@ end subroutine GFS_SCNV_generic_pre_finalize subroutine GFS_SCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, gu0, gv0, gt0, gq0, & save_u, save_v, save_t, save_q, ntqv, nsamftrac, flag_for_scnv_generic_tend, & dtidx, index_of_process_scnv, ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc, & - cscnv, satmedmf, trans_trac, ras, ntrac, clw, errmsg, errflg) + ntsigma, cscnv, satmedmf, trans_trac, ras, ntrac, clw, errmsg, errflg) use machine, only: kind_phys implicit none integer, intent(in) :: im, levs, ntqv, nsamftrac, index_of_process_scnv, dtidx(:,:) - integer, intent(in) :: ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac + integer, intent(in) :: ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac logical, intent(in) :: ldiag3d, qdiag3d, flag_for_scnv_generic_tend real(kind=kind_phys), dimension(:,:), intent(in) :: gu0, gv0, gt0 real(kind=kind_phys), dimension(:,:,:), intent(in) :: gq0 @@ -55,7 +55,7 @@ subroutine GFS_SCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, gu0, gv0, gt0, do n=2,ntrac if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt .and. & n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & - n /= ntsnc .and. n /= ntgl .and. n /= ntgnc) then + n /= ntsnc .and. n /= ntgl .and. n /= ntgnc .and. n /= ntsigma) then tracers = tracers + 1 if(dtidx(100+n,index_of_process_scnv)>0) then save_q(:,:,n) = clw(:,:,tracers) @@ -97,7 +97,7 @@ subroutine GFS_SCNV_generic_post_run (im, levs, nn, lssav, ldiag3d, qdiag3d, & rainc, cnvprcp, cnvprcpb, cnvw_phy_f3d, cnvc_phy_f3d, & dtend, dtidx, index_of_temperature, index_of_x_wind, index_of_y_wind, & index_of_process_scnv, ntqv, flag_for_scnv_generic_tend, & - ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc, & + ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc, & imfshalcnv, imfshalcnv_sas, imfshalcnv_samf, ntrac, & cscnv, satmedmf, trans_trac, ras, errmsg, errflg) @@ -106,7 +106,7 @@ subroutine GFS_SCNV_generic_post_run (im, levs, nn, lssav, ldiag3d, qdiag3d, & implicit none integer, intent(in) :: im, levs, nn, ntqv, nsamftrac - integer, intent(in) :: ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac + integer, intent(in) :: ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac logical, intent(in) :: lssav, ldiag3d, qdiag3d, flag_for_scnv_generic_tend real(kind=kind_phys), intent(in) :: frain real(kind=kind_phys), dimension(:,:), intent(in) :: gu0, gv0, gt0 @@ -186,7 +186,7 @@ subroutine GFS_SCNV_generic_post_run (im, levs, nn, lssav, ldiag3d, qdiag3d, & do n=2,ntrac if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt .and. & n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & - n /= ntsnc .and. n /= ntgl .and. n /= ntgnc) then + n /= ntsnc .and. n /= ntgl .and. n /= ntgnc .and. n/= ntsigma) then tracers = tracers + 1 idtend = dtidx(100+n,index_of_process_scnv) if(idtend>0) then diff --git a/physics/GFS_SCNV_generic.meta b/physics/GFS_SCNV_generic.meta index 5cbda127c..4fd189948 100644 --- a/physics/GFS_SCNV_generic.meta +++ b/physics/GFS_SCNV_generic.meta @@ -183,6 +183,13 @@ dimensions = () type = integer intent = in +[ntsigma] + standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array + long_name = tracer index of updraft_area_fraction + units = index + dimensions = () + type = integer + intent = in [ntrw] standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array long_name = tracer index for rain water @@ -614,6 +621,13 @@ dimensions = () type = integer intent = in +[ntsigma] + standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array + long_name = tracer index of updraft_area_fraction + units = index + dimensions = () + type = integer + intent = in [ntrw] standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array long_name = tracer index for rain water diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 6772bc8d4..ca05f6778 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -1,7 +1,8 @@ !>\file progsigma !! This file contains the subroutine that calculates the prognostic !! updraft area fraction that is used for closure computations in -!! saSAS deep and shallow convection. +!! saSAS deep and shallow convection, based on a moisture budget +!! as described in Bengtsson et al. 2022. !>\ingroup samfdeepcnv !! This subroutine computes a prognostic updraft area fraction @@ -13,9 +14,8 @@ !> @{ subroutine progsigma_calc (im,km,flag_init,flag_restart, & - del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & - qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & - do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & + flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & + delt,qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) ! ! @@ -30,12 +30,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & real, intent(in) :: qgrs_dsave(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) - logical, intent(in) :: flag_init,flag_restart,cnvflg(im) - real(kind=kind_phys), intent(in) :: nthresh - real(kind=kind_phys), intent(in) :: ca_deep(im) + logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow real(kind=kind_phys), intent(out):: ca_micro(im) - logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger - real(kind=kind_phys), intent(in) :: sigmain(im,km) ! intent out @@ -47,28 +43,29 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! Local variables integer :: i,k,km1 real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & - mcons(im),zfdqa(im),zform(im,km), & + mcons(im),fdqa(im),form(im,km), & qadv(im,km),sigmamax(im) - real(kind=kind_phys) :: gcvalmx,ZEPS7,ZZ,ZCVG,mcon,buy2, & - zfdqb,dtdyn,dxlim,rmulacvg,dp,tem, & - alpha,DEN + real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & + fdqb,dtdyn,dxlim,rmulacvg,dp,tem, & + alpha,DEN,betascu integer :: inbu(im,km) !Parameters gcvalmx = 0.1 rmulacvg=10. - ZEPS7=1.E-11 + epsilon=1.E-11 km1=km-1 alpha=7000. + betascu = 3.0 !Initialization 2D do k = 1,km do i = 1,im sigmaout(i,k)=0. inbu(i,k)=0 - zform(i,k)=0. + form(i,k)=0. enddo enddo @@ -80,8 +77,9 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & termB(i)=0. termC(i)=0. termD(i)=0. - zfdqa(i)=0. + fdqa(i)=0. mcons(i)=0. + ca_micro(i)=0. enddo !Initial computations, place maximum sigmain in sigmab @@ -94,9 +92,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif else if(cnvflg(i))then - !if(sigmain(i,k)<1.E-5)then - ! sigmain(i,k)=0. - !endif if(sigmain(i,k)>sigmab(i))then sigmab(i)=sigmain(i,k) endif @@ -107,7 +102,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do i=1,im if(sigmab(i) < 1.E-5)then !after advection - sigmab(i)=0. + sigmab(i)=0. endif enddo @@ -180,11 +175,11 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do i = 1,im if(cnvflg(i))then dp = 1000. * del(i,k) - zform(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt) - zfdqb=0.5*((zform(i,k)*zdqca(i,k))) + form(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt) + fdqb=0.5*((form(i,k)*zdqca(i,k))) termC(i)=termC(i)+(float(inbu(i,k))* & - (zfdqb+zfdqa(i))*hvap*zeta(i,k)) - zfdqa(i)=zfdqb + (fdqb+fdqa(i))*hvap*zeta(i,k)) + fdqa(i)=fdqb endif enddo enddo @@ -193,29 +188,26 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do i = 1,im if(cnvflg(i))then - DEN=MIN(termC(i)+termB(i),1.E8) !1.E8 - !DEN=MAX(termC(i)+termB(i),1.E7) !1.E7 - - ZCVG=termD(i)*delt - + DEN=MIN(termC(i)+termB(i),1.E8) + cvg=termD(i)*delt ZZ=MAX(0.0,SIGN(1.0,termA(i))) & *MAX(0.0,SIGN(1.0,termB(i))) & - *MAX(0.0,SIGN(1.0,termC(i)-ZEPS7)) + *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) - ZCVG=MAX(0.0,ZCVG) + cvg=MAX(0.0,cvg) - if(flag_init)then + if(flag_init .and. .not. flag_restart)then sigmab(i)=0.03 else - sigmab(i)=(ZZ*(termA(i)+ZCVG))/(DEN+(1.0-ZZ)) + sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) endif if(sigmab(i)>0.)then sigmab(i)=MIN(sigmab(i),sigmamax(i)) sigmab(i)=MAX(sigmab(i),0.01) endif - + ca_micro(i)=sigmab(i) endif!cnvflg enddo @@ -226,7 +218,20 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo enddo + + !Since updraft velocity is much lower in shallow cu region, termC becomes small in shallow cu application, thus the area fraction + !in this regime becomes too large compared with the deep cu region. To address this simply apply a scaling factor for shallow cu + !before computing the massflux to reduce the total strength of the SC MF: + if(flag_shallow)then + do i= 1, im + if(cnvflg(i)) then + sigmab(i)=sigmab(i)/betascu + endif + enddo + endif + + end subroutine progsigma_calc !> @} !! @} diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 7be31a04a..35aea0eb1 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -216,7 +216,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),qlks(im,km), & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) - + logical flag_shallow c physical parameters ! parameter(grav=grav,asolfac=0.958) ! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) @@ -1729,7 +1729,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & enddo enddo - if(progsigma)then + if(progsigma)then do k = 2, km1 do i = 1, im if (cnvflg(i)) then @@ -1776,8 +1776,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & c !> - Calculate the mean updraft velocity within the cloud (wc),cast in pressure coordinates. - if(progsigma)then - + if(progsigma)then do i = 1, im omegac(i) = 0. sumx(i) = 0. @@ -1861,17 +1860,19 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & c c store term needed for "termC" in prognostic area fraction closure - do k = 2, km1 - do i = 1, im - dp = 1000. * del(i,k) - if (cnvflg(i)) then - if(k > kbcon(i) .and. k < ktcon(i)) then - zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + - & pwo(i,k)+dellal(i,k)) + if(progsigma)then + do k = 2, km1 + do i = 1, im + dp = 1000. * del(i,k) + if (cnvflg(i)) then + if(k > kbcon(i) .and. k < ktcon(i)) then + zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + + & pwo(i,k)+dellal(i,k)) + endif endif - endif + enddo enddo - enddo + endif ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then ccccc print *, ' aa1(i) before dwndrft =', aa1(i) @@ -2890,10 +2891,10 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - From Bengtsson et al. (2022) Prognostic closure scheme, equation 8, compute updraft area fraction based on a moisture budget if(progsigma)then - call progsigma_calc(im,km,first_time_step,restart, + flag_shallow = .false. + call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, - & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) endif diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 3eb330551..71f9b87a5 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -70,15 +70,15 @@ type = logical intent = in [tmf] - standard_name = turbulence_moisture_flux - long_name = turbulence_moisture_flux + standard_name = turbulence_moisture_flux_for_coupling_to_convection + long_name = turbulence_moisture_flux_for_coupling_to_convection units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in [qmicro] - standard_name = instantanious_moisture_tendency_due_to_microphysics + standard_name = instantaneous_moisture_tendency_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) @@ -450,13 +450,13 @@ kind = kind_phys intent = inout [sigmain] - standard_name = updraft_area_fraction + standard_name = prognostic_updraft_area_fraction_in_convection long_name = convective updraft area fraction units = frac dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = inout + intent = in [sigmaout] standard_name = updraft_area_fraction_updated_by_physics long_name = convective updraft area fraction updated by physics diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 24e01b040..6a682e9eb 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -57,11 +57,13 @@ end subroutine samfshalcnv_finalize !! @{ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & - & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,q1,t1,u1,v1,fscav, & + & t0c,delt,ntk,ntr,delp,first_time_step,restart, & + & tmf,qmicro,progsigma, & + & prslp,psp,phil,qtr,qgrs_dsave,q,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & - & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal,errmsg,errflg) + & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, + & ca_micro,sigmain,sigmaout,errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -74,7 +76,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & eps, epsm1, fv, grav, hvap, rd, rv, t0c real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & - & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:) + & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), & + & qmicro(:,:),tmf(:,:),qgrs_dsave(:,:),q(:,:),sigmain(:,:) ! real(kind=kind_phys), dimension(:), intent(in) :: fscav integer, intent(inout) :: kcnv(:) @@ -83,12 +86,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & q1(:,:), t1(:,:), u1(:,:), v1(:,:) ! integer, intent(out) :: kbot(:), ktop(:) - real(kind=kind_phys), intent(out) :: rn(:), & - & cnvw(:,:), cnvc(:,:), ud_mf(:,:), dt_mf(:,:) + real(kind=kind_phys), intent(out) :: rn(:), ca_micro(:), & + & cnvw(:,:), cnvc(:,:), ud_mf(:,:), dt_mf(:,:), sigmaout(:,:) ! real(kind=kind_phys), intent(in) :: clam, c0s, c1, & & asolfac, evef, pgcon - logical, intent(in) :: hwrf_samfshal + logical, intent(in) :: hwrf_samfshal,first_time_step, & + & restart,progsigma character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! @@ -155,6 +159,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & bb1, bb2, wucb cc + +! parameters for prognostic sigma closure + real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), + & omegac(im),zeta(im,km),dbyo1(im,km), + & sigmab(im) + logical flag_shallow + c physical parameters ! parameter(g=grav,asolfac=0.89) ! parameter(g=grav) @@ -323,6 +334,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! vshear(i) = 0. gdx(i) = sqrt(garea(i)) xmb(i) = 0. + ca_micro(i) = 0. enddo endif !! @@ -498,6 +510,21 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + + + do i = 1,im + omegac(i)=0. + enddo + + do k = 1, km + do i = 1, im + dbyo1(i,k)=0. + zdqca(i,k)=0. + qlks(i,k)=0. + omega_u(i,k)=0. + zeta(i,k)=1.0 + enddo + enddo ! ! initialize tracer variables ! @@ -1237,6 +1264,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & qcko(i,k)= qlk + qrch pwo(i,k) = etah * c0t(i,k) * dz * qlk cnvwt(i,k) = etah * qlk * grav / dp + qlks(i,k)=qlk endif ! ! compute buoyancy and drag for updraft velocity @@ -1304,6 +1332,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(k >= kbcon(i) .and. k < ktcon(i)) then dz1 = zo(i,k+1) - zo(i,k) aa1(i) = aa1(i) + buo(i,k) * dz1 + dbyo1(i,k) = hcko(i,k) - heso(i,k) endif endif enddo @@ -1402,6 +1431,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & qcko(i,k) = qlk + qrch pwo(i,k) = etah * c0t(i,k) * dz * qlk cnvwt(i,k) = etah * qlk * grav / dp + qlks(i,k)=qlk endif endif endif @@ -1444,6 +1474,20 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo enddo ! + if(progsigma)then + do k = 2, km1 + do i = 1, im + if (cnvflg(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.) + endif + endif + enddo + enddo + endif + ! compute updraft velocity averaged over the whole cumulus ! !> - Calculate the mean updraft velocity within the cloud (wc). @@ -1475,6 +1519,50 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif enddo c +!> - Calculate the mean updraft velocity in pressure coordinates within the cloud (wc). + if(progsigma)then + do i = 1, im + omegac(i) = 0. + sumx(i) = 0. + enddo + do k = 2, km1 + do i = 1, im + if (cnvflg(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 + sumx(i) = sumx(i) + dp + endif + endif + enddo + enddo + do i = 1, im + if(cnvflg(i)) then + if(sumx(i) == 0.) then + cnvflg(i)=.false. + else + omegac(i) = omegac(i) / sumx(i) + endif + val = -1.2 + if (omegac(i) > val) cnvflg(i)=.false. + endif + enddo +c +! > - Calculate the mean updraft velocity within the cloud (omega). + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then + zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + zeta(i,k)=MAX(0.,zeta(i,k)) + zeta(i,k)=MIN(1.,zeta(i,k)) + endif + endif + enddo + enddo + endif !if progsigma + c exchange ktcon with ktcon1 c do i = 1, im @@ -1505,11 +1593,25 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(dq > 0.) then qlko_ktcon(i) = dq qcko(i,k) = qrch + qlks(i,k) = qlko_ktcon(i) endif endif enddo endif c + + do k = 2, km1 + do i = 1, im + dp = 1000. * del(i,k) + if (cnvflg(i)) then + if(k > kbcon(i) .and. k < ktcon(i)) then + zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + + & pwo(i,k)+dellal(i,k)) + endif + endif + enddo + enddo + c--- compute precipitation efficiency in terms of windshear c !! - Calculate the wind shear and precipitation efficiency according to equation 58 in Fritsch and Chappell (1980) \cite fritsch_and_chappell_1980 : @@ -1824,13 +1926,26 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c compute cloud base mass flux as a function of the mean c updraft velcoity c +c Prognostic closure + if(progsigma)then + flag_shallow = .true. + call progsigma_calc(im,km,first_time_step,restart,flag_shallow, + & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, + & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, + & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + endif + !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. do i= 1, im if(cnvflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) - xmb(i) = advfac(i)*betaw*rho*wc(i) + if(progsigma)then + xmb(i) = sigmab(i)*((-1.0*omegac(i))/grav) + else + xmb(i) = advfac(i)*betaw*rho*wc(i) + endif endif enddo ! @@ -1850,8 +1965,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then if (gdx(i) < dxcrt) then - scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) - scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) + if(progsigma)then + scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i)) + else + scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) + endif + scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) else scaldfunc(i) = 1.0 endif diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index d768d4451..4383b8a67 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = samfshalcnv type = scheme - dependencies = funcphys.f90,machine.F,samfaerosols.F + dependencies = funcphys.f90,machine.F,samfaerosols.F,progsigma_calc.f90 ######################################################################## [ccpp-arg-table] @@ -55,6 +55,36 @@ dimensions = () type = integer intent = in +[first_time_step] + standard_name = flag_for_first_timestep + long_name = flag for first time step for time integration loop (cold/warmstart) + units = flag + dimensions = () + type = logical + intent = in +[restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in +[tmf] + standard_name = turbulence_moisture_flux_for_coupling_to_convection + long_name = turbulence_moisture_flux_for_coupling_to_convection + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[qmicro] + standard_name = instantaneous_moisture_tendency_due_to_microphysics + long_name = moisture tendency due to microphysics + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [itc] standard_name = index_of_first_chemical_tracer_for_convection long_name = index of first chemical tracer transported/scavenged by convection @@ -219,6 +249,22 @@ type = real kind = kind_phys intent = inout +[qgrs_dsave] + standard_name = tracer_concentration_dsave + long_name = model layer mean tracer concentration dsave + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[q] + standard_name = specific_humidity + long_name = water vapor specific humidity + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [q1] standard_name = specific_humidity_of_new_state long_name = updated vapor specific humidity @@ -413,6 +459,37 @@ dimensions = () type = logical intent = in +[progsigma] + standard_name = flag_for_prognostic_sigma + long_name = flag for prognostic sigma + units = flag + dimensions = () + type = logical + intent = in +[ca_micro] + standard_name = output_prognostic_sigma_two + long_name = output of prognostic area fraction two + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[sigmain] + standard_name = prognostic_updraft_area_fraction_in_convection + long_name = convective updraft area fraction + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[sigmaout] + standard_name = updraft_area_fraction_updated_by_physics + long_name = convective updraft area fraction updated by physics + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index b0367b4d4..0fce7dd9a 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -121,9 +121,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & real(kind=kind_phys), intent(out) :: & & dusfc(:), dvsfc(:), & & dtsfc(:), dqsfc(:), & - & hpbl(:), tmf(:,:) + & hpbl(:) real(kind=kind_phys), intent(out) :: & - & dkt(:,:), dku(:,:) + & dkt(:,:), dku(:,:), tmf(:,:) ! logical, intent(in) :: dspheat character(len=*), intent(out) :: errmsg @@ -299,6 +299,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & xmfd(i,k) = 0. buou(i,k) = 0. buod(i,k) = 0. + tmf(i,k) = 0. ckz(i,k) = ck1 chz(i,k) = ch1 rlmnz(i,k) = rlmn0 diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 211bbaec6..9b803e4a5 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -202,10 +202,10 @@ kind = kind_phys intent = inout [tmf] - standard_name = turbulence_moisture_flux - long_name = turbulence_moisture_flux + standard_name = turbulence_moisture_flux_for_coupling_to_convection + long_name = turbulence_moisture_flux_for_coupling_to_convection units = kg kg-1 s-1 - dimensions = (horizontal_dimension,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = out From b530db11801ca67a10e2757aef04ed9492a06e7e Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 20 Apr 2022 01:29:13 +0000 Subject: [PATCH 04/11] cleaning some diagnostics --- physics/progsigma_calc.f90 | 5 +---- physics/samfdeepcnv.f | 8 +++----- physics/samfdeepcnv.meta | 8 -------- physics/samfshalcnv.f | 7 +++---- physics/samfshalcnv.meta | 8 -------- 5 files changed, 7 insertions(+), 29 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index ca05f6778..7673602b6 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -16,7 +16,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & delt,qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & - ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + sigmain,sigmaout,sigmab,errmsg,errflg) ! ! use machine, only : kind_phys @@ -31,7 +31,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow - real(kind=kind_phys), intent(out):: ca_micro(im) real(kind=kind_phys), intent(in) :: sigmain(im,km) ! intent out @@ -79,7 +78,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & termD(i)=0. fdqa(i)=0. mcons(i)=0. - ca_micro(i)=0. enddo !Initial computations, place maximum sigmain in sigmab @@ -207,7 +205,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & sigmab(i)=MIN(sigmab(i),sigmamax(i)) sigmab(i)=MAX(sigmab(i),0.01) endif - ca_micro(i)=sigmab(i) endif!cnvflg enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 35aea0eb1..45cbf70e1 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -86,7 +86,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, sigmain, sigmaout, ca_micro, & + & rainevap, sigmain, sigmaout, & & errmsg,errflg) ! use machine , only : kind_phys @@ -108,7 +108,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), qgrs_dsave(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:),ca_micro(:) + real(kind=kind_phys), intent(out) :: rainevap(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -919,8 +919,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & enddo if(totflg) return !! -! -!Lisa: at this point only trigger criteria is set ! turbulent entrainment rate assumed to be proportional ! to subcloud mean TKE @@ -2895,7 +2893,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, - & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + & sigmain,sigmaout,sigmab,errmsg,errflg) endif !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 71f9b87a5..71c78036d 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -326,14 +326,6 @@ dimensions = () type = logical intent = in -[ca_micro] - standard_name = output_prognostic_sigma_two - long_name = output of prognostic area fraction two - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [cldwrk] standard_name = cloud_work_function long_name = cloud work function diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 6a682e9eb..343691279 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -63,7 +63,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, - & ca_micro,sigmain,sigmaout,errmsg,errflg) + & sigmain,sigmaout,errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -86,7 +86,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & q1(:,:), t1(:,:), u1(:,:), v1(:,:) ! integer, intent(out) :: kbot(:), ktop(:) - real(kind=kind_phys), intent(out) :: rn(:), ca_micro(:), & + real(kind=kind_phys), intent(out) :: rn(:), & & cnvw(:,:), cnvc(:,:), ud_mf(:,:), dt_mf(:,:), sigmaout(:,:) ! real(kind=kind_phys), intent(in) :: clam, c0s, c1, & @@ -334,7 +334,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! vshear(i) = 0. gdx(i) = sqrt(garea(i)) xmb(i) = 0. - ca_micro(i) = 0. enddo endif !! @@ -1932,7 +1931,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, - & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + & sigmain,sigmaout,sigmab,errmsg,errflg) endif !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity. diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index 4383b8a67..895460ffd 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -466,14 +466,6 @@ dimensions = () type = logical intent = in -[ca_micro] - standard_name = output_prognostic_sigma_two - long_name = output of prognostic area fraction two - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [sigmain] standard_name = prognostic_updraft_area_fraction_in_convection long_name = convective updraft area fraction From fc7e7a0e226663e0fec80729097a7dc69ed5551d Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Fri, 22 Apr 2022 04:00:31 +0000 Subject: [PATCH 05/11] addressing some review comments --- physics/GFS_MP_generic.F90 | 16 +++--- physics/GFS_MP_generic.meta | 12 ++--- physics/progsigma_calc.f90 | 103 ++++++++++++++++++------------------ physics/samfdeepcnv.f | 10 ++-- physics/samfdeepcnv.meta | 14 ++--- physics/samfshalcnv.f | 12 +++-- physics/samfshalcnv.meta | 14 ++--- physics/satmedmfvdifq.meta | 4 +- 8 files changed, 95 insertions(+), 90 deletions(-) diff --git a/physics/GFS_MP_generic.F90 b/physics/GFS_MP_generic.F90 index dbf2d15fa..f8dd8ab6c 100644 --- a/physics/GFS_MP_generic.F90 +++ b/physics/GFS_MP_generic.F90 @@ -95,7 +95,7 @@ subroutine GFS_MP_generic_post_run( drain_cpl, dsnow_cpl, lsm, lsm_ruc, lsm_noahmp, raincprv, rainncprv, iceprv, snowprv, & graupelprv, draincprv, drainncprv, diceprv, dsnowprv, dgraupelprv, dtp, dfi_radar_max_intervals, & dtend, dtidx, index_of_temperature, index_of_process_mp,ldiag3d, qdiag3d,dqdt_qmicro, lssav, num_dfi_radar, & - fh_dfi_radar,index_of_process_dfi_radar, ix_dfi_radar, dfi_radar_tten, radar_tten_limits, fhour, qgrs_dsave, & + fh_dfi_radar,index_of_process_dfi_radar, ix_dfi_radar, dfi_radar_tten, radar_tten_limits, fhour, prevsq, & errmsg, errflg) ! use machine, only: kind_phys @@ -150,7 +150,7 @@ subroutine GFS_MP_generic_post_run( real(kind=kind_phys), dimension(:), intent(inout) :: dsnowprv real(kind=kind_phys), dimension(:), intent(inout) :: dgraupelprv real(kind=kind_phys), dimension(:,:), intent(out) :: dqdt_qmicro - real(kind=kind_phys), dimension(:,:), intent(out) :: qgrs_dsave + real(kind=kind_phys), dimension(:,:), intent(out) :: prevsq real(kind=kind_phys), intent(in) :: dtp ! CCPP error handling @@ -466,11 +466,13 @@ subroutine GFS_MP_generic_post_run( pwat(i) = pwat(i) * onebg enddo - do k = 1, levs - do i=1, im - qgrs_dsave(i,k) = gq0(i,k,1) - enddo - enddo + if(progsigma)then + do k = 1, levs + do i=1, im + prevsq(i,k) = gq0(i,k,1) + enddo + enddo + endif end subroutine GFS_MP_generic_post_run !> @} diff --git a/physics/GFS_MP_generic.meta b/physics/GFS_MP_generic.meta index 763cad85a..eb0f17fa8 100644 --- a/physics/GFS_MP_generic.meta +++ b/physics/GFS_MP_generic.meta @@ -249,8 +249,8 @@ type = logical intent = in [progsigma] - standard_name = flag_for_prognostic_sigma - long_name = flag for prognostic sigma + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumulus scheme units = flag dimensions = () type = logical @@ -852,16 +852,16 @@ type = logical intent = in [dqdt_qmicro] - standard_name = instantaneous_moisture_tendency_due_to_microphysics + standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = out -[qgrs_dsave] - standard_name = tracer_concentration_dsave - long_name = model layer mean tracer concentration dsave +[prevsq] + standard_name = specific_humidity_on_previous_timestep + long_name = specific_humidity_on_previous_timestep units = kg kg-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 7673602b6..58a6fc0ef 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -1,4 +1,4 @@ -!>\file progsigma +!>\file progsigma_calc.f90 !! This file contains the subroutine that calculates the prognostic !! updraft area fraction that is used for closure computations in !! saSAS deep and shallow convection, based on a moisture budget @@ -15,7 +15,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & - delt,qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & + delt,prevsq,q,kbcon1,ktcon,cnvflg,gdx, & sigmain,sigmaout,sigmab,errmsg,errflg) ! ! @@ -27,7 +27,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! intent in integer, intent(in) :: im,km,kbcon1(im),ktcon(im) real, intent(in) :: hvap,delt - real, intent(in) :: qgrs_dsave(im,km), q(im,km),del(im,km), & + real, intent(in) :: prevsq(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow @@ -43,33 +43,32 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & integer :: i,k,km1 real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & mcons(im),fdqa(im),form(im,km), & - qadv(im,km),sigmamax(im) + qadv(im,km),sigmamax(im),dp(im),inbu(im,km) real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & - fdqb,dtdyn,dxlim,rmulacvg,dp,tem, & - alpha,DEN,betascu - integer :: inbu(im,km) + fdqb,dtdyn,dxlim,rmulacvg,tem, & + alpha,DEN,betascu,dp1 !Parameters - gcvalmx = 0.1 - rmulacvg=10. - epsilon=1.E-11 - km1=km-1 - alpha=7000. - betascu = 3.0 + gcvalmx = 0.1 + rmulacvg=10. + epsilon=1.E-11 + km1=km-1 + alpha=7000. + betascu = 3.0 !Initialization 2D - do k = 1,km - do i = 1,im - sigmaout(i,k)=0. - inbu(i,k)=0 - form(i,k)=0. - enddo - enddo + do k = 1,km + do i = 1,im + sigmaout(i,k)=0. + inbu(i,k)=0. + form(i,k)=0. + enddo + enddo !Initialization 1D - do i=1,im + do i=1,im sigmab(i)=0. sigmamax(i)=0.95 termA(i)=0. @@ -80,23 +79,32 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & mcons(i)=0. enddo - !Initial computations, place maximum sigmain in sigmab + do k = 2,km1 + do i = 1,im + if(cnvflg(i))then + dp(i) = 1000. * del(i,k) + endif + enddo + enddo - do k=2,km - do i=1,im - if(flag_init .and. .not. flag_restart)then + !Initial computations, place maximum sigmain in sigmab + if(flag_init .and. .not. flag_restart)then + do i=1,im if(cnvflg(i))then sigmab(i)=0.03 endif - else + enddo + else + do i=1,im if(cnvflg(i))then - if(sigmain(i,k)>sigmab(i))then - sigmab(i)=sigmain(i,k) - endif + do k=2,km + if(sigmain(i,k)>sigmab(i))then + sigmab(i)=sigmain(i,k) + endif + enddo endif - endif - enddo - enddo + enddo + endif do i=1,im if(sigmab(i) < 1.E-5)then !after advection @@ -116,7 +124,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & if(flag_init .and. .not.flag_restart)then qadv(i,k)=0. else - qadv(i,k)=(q(i,k) - qgrs_dsave(i,k))/delt + qadv(i,k)=(q(i,k) - prevsq(i,k))/delt endif enddo enddo @@ -125,22 +133,21 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !buoyant layers with positive moisture convergence (accumulated from the surface). !Lowest level: do i = 1,im - dp = 1000. * del(i,1) - mcons(i)=(hvap*(qadv(i,1)+tmf(i,1)+qmicro(i,1))*dp) + dp1 = 1000. * del(i,1) + mcons(i)=(hvap*(qadv(i,1)+tmf(i,1)+qmicro(i,1))*dp1) enddo !Levels above: do k = 2,km1 do i = 1,im - dp = 1000. * del(i,k) if(cnvflg(i))then - mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp) + mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i)) buy2 = termD(i)+mcon+mcons(i) ! Do the integral over buoyant layers with positive mcon acc from surface if(k > kbcon1(i) .and. k < ktcon(i) .and. buy2 > 0.)then - inbu(i,k)=1 + inbu(i,k)=1. endif inbu(i,k-1)=MAX(inbu(i,k-1),inbu(i,k)) - termD(i) = termD(i) + float(inbu(i,k-1))*mcons(i) + termD(i) = termD(i) + inbu(i,k-1)*mcons(i) mcons(i)=mcon endif enddo @@ -149,9 +156,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !termA do k = 2,km1 do i = 1,im - dp = 1000. * del(i,k) if(cnvflg(i))then - tem=(sigmab(i)*zeta(i,k)*float(inbu(i,k))*dbyo1(i,k))*dp + tem=(sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k))*dp(i) termA(i)=termA(i)+tem endif enddo @@ -160,9 +166,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !termB do k = 2,km1 do i = 1,im - dp = 1000. * del(i,k) if(cnvflg(i))then - tem=(dbyo1(i,k)*float(inbu(i,k)))*dp + tem=(dbyo1(i,k)*inbu(i,k))*dp(i) termB(i)=termB(i)+tem endif enddo @@ -172,10 +177,9 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - dp = 1000. * del(i,k) - form(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt) + form(i,k)=-1.0*inbu(i,k)*(omega_u(i,k)*delt) fdqb=0.5*((form(i,k)*zdqca(i,k))) - termC(i)=termC(i)+(float(inbu(i,k))* & + termC(i)=termC(i)+inbu(i,k)* & (fdqb+fdqa(i))*hvap*zeta(i,k)) fdqa(i)=fdqb endif @@ -185,22 +189,17 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !sigmab do i = 1,im if(cnvflg(i))then - DEN=MIN(termC(i)+termB(i),1.E8) cvg=termD(i)*delt ZZ=MAX(0.0,SIGN(1.0,termA(i))) & *MAX(0.0,SIGN(1.0,termB(i))) & - *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) - - + *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) cvg=MAX(0.0,cvg) - if(flag_init .and. .not. flag_restart)then sigmab(i)=0.03 else sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) endif - if(sigmab(i)>0.)then sigmab(i)=MIN(sigmab(i),sigmamax(i)) sigmab(i)=MAX(sigmab(i),0.01) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 45cbf70e1..02b2dcb83 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -79,7 +79,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & tmf,qmicro,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,qgrs_dsave,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & hwrf_samfdeep,progsigma,wclosureflg,cldwrk,rn,kbot,ktop,kcnv, & & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & @@ -107,7 +107,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: nthresh real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & - & tmf(:,:),q(:,:), qgrs_dsave(:,:) + & tmf(:,:),q(:,:), prevsq(:,:) real(kind=kind_phys), intent(out) :: rainevap(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -217,6 +217,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) logical flag_shallow + real(kind=kind_phys) gravinv c physical parameters ! parameter(grav=grav,asolfac=0.958) ! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) @@ -309,6 +310,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & errmsg = '' errflg = 0 + gravinv = 1./grav elocp = hvap/cp el2orc = hvap*hvap/(rv*cp) @@ -2892,7 +2894,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & flag_shallow = .false. call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, + & prevsq,q,kbcon1,ktcon,cnvflg,gdx, & sigmain,sigmaout,sigmab,errmsg,errflg) endif @@ -2903,7 +2905,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) if(progsigma)then - xmb(i) = sigmab(i)*((-1.0*omegac(i))/grav) + xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv) else xmb(i) = advfac(i)*betaw*rho*wc(i) endif diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 71c78036d..e956d24ed 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -70,8 +70,8 @@ type = logical intent = in [tmf] - standard_name = turbulence_moisture_flux_for_coupling_to_convection - long_name = turbulence_moisture_flux_for_coupling_to_convection + standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL + long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -249,9 +249,9 @@ type = real kind = kind_phys intent = inout -[qgrs_dsave] - standard_name = tracer_concentration_dsave - long_name = model layer mean tracer concentration dsave +[prevsq] + standard_name = specific_humidity_on_previous_timestep + long_name = specific_humidity_on_previous_timestep units = kg kg-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -320,8 +320,8 @@ type = logical intent = in [progsigma] - standard_name = flag_for_prognostic_sigma - long_name = flag for prognostic sigma + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumuls scheme units = flag dimensions = () type = logical diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 343691279..c3bb842b3 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -59,7 +59,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp,first_time_step,restart, & & tmf,qmicro,progsigma, & - & prslp,psp,phil,qtr,qgrs_dsave,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, @@ -77,7 +77,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), & - & qmicro(:,:),tmf(:,:),qgrs_dsave(:,:),q(:,:),sigmain(:,:) + & qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:),sigmain(:,:) ! real(kind=kind_phys), dimension(:), intent(in) :: fscav integer, intent(inout) :: kcnv(:) @@ -165,6 +165,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & omegac(im),zeta(im,km),dbyo1(im,km), & sigmab(im) logical flag_shallow + real(kind=kind_phys) gravinv c physical parameters ! parameter(g=grav,asolfac=0.89) @@ -249,6 +250,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & errmsg = '' errflg = 0 + gravinv = 1./grav + elocp = hvap/cp el2orc = hvap*hvap/(rv*cp) @@ -1601,7 +1604,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im - dp = 1000. * del(i,k) if (cnvflg(i)) then if(k > kbcon(i) .and. k < ktcon(i)) then zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + @@ -1930,7 +1932,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & flag_shallow = .true. call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, + & prevsq,q,kbcon1,ktcon,cnvflg,gdx, & sigmain,sigmaout,sigmab,errmsg,errflg) endif @@ -1941,7 +1943,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) if(progsigma)then - xmb(i) = sigmab(i)*((-1.0*omegac(i))/grav) + xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv) else xmb(i) = advfac(i)*betaw*rho*wc(i) endif diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index 895460ffd..a4cca64b8 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -70,8 +70,8 @@ type = logical intent = in [tmf] - standard_name = turbulence_moisture_flux_for_coupling_to_convection - long_name = turbulence_moisture_flux_for_coupling_to_convection + standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL + long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -249,9 +249,9 @@ type = real kind = kind_phys intent = inout -[qgrs_dsave] - standard_name = tracer_concentration_dsave - long_name = model layer mean tracer concentration dsave +[prevsq] + standard_name = specific_humidity_on_previous_timestep + long_name = specific_humidity_on_previous_timestep units = kg kg-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -460,8 +460,8 @@ type = logical intent = in [progsigma] - standard_name = flag_for_prognostic_sigma - long_name = flag for prognostic sigma + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumulus scheme units = flag dimensions = () type = logical diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 9b803e4a5..fa30cd9f7 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -202,8 +202,8 @@ kind = kind_phys intent = inout [tmf] - standard_name = turbulence_moisture_flux_for_coupling_to_convection - long_name = turbulence_moisture_flux_for_coupling_to_convection + standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL + long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real From 0200e2d05757af34a070e52fe8558ecbaa73a42a Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 27 Apr 2022 18:48:38 +0000 Subject: [PATCH 06/11] addressing some review comments --- physics/progsigma_calc.f90 | 2 +- physics/samfdeepcnv.f | 2 +- physics/samfdeepcnv.meta | 2 +- physics/samfshalcnv.f | 2 +- physics/samfshalcnv.meta | 2 +- physics/satmedmfvdifq.F | 25 ++++++++++++++++++++----- physics/satmedmfvdifq.meta | 7 +++++++ 7 files changed, 32 insertions(+), 10 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 58a6fc0ef..55f6b5e3a 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -180,7 +180,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & form(i,k)=-1.0*inbu(i,k)*(omega_u(i,k)*delt) fdqb=0.5*((form(i,k)*zdqca(i,k))) termC(i)=termC(i)+inbu(i,k)* & - (fdqb+fdqa(i))*hvap*zeta(i,k)) + (fdqb+fdqa(i))*hvap*zeta(i,k) fdqa(i)=fdqb endif enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index f8a60a5e3..5fd54a2ec 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -79,7 +79,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & tmf,qmicro,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & hwrf_samfdeep,progsigma,wclosureflg,cldwrk,rn,kbot,ktop,kcnv, & & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index e956d24ed..2b2942812 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -78,7 +78,7 @@ kind = kind_phys intent = in [qmicro] - standard_name = instantaneous_moisture_tendency_due_to_microphysics + standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 56571457a..325566877 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -59,7 +59,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp,first_time_step,restart, & & tmf,qmicro,progsigma, & - & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index a4cca64b8..8c9735c32 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -78,7 +78,7 @@ kind = kind_phys intent = in [qmicro] - standard_name = instantaneous_moisture_tendency_due_to_microphysics + standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 0fce7dd9a..c7a6fadc9 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -67,8 +67,8 @@ end subroutine satmedmfvdifq_finalize !! (mfscuq.f). !! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm !! @{ - subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & - & grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & + subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & + & ntiw,ntke,grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & & dv,du,tdt,rtg,tmf,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & & tsea,heat,evap,stress,spd1,kpbl, & @@ -91,7 +91,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & integer, intent(in) :: sfc_rlm integer, intent(in) :: kinver(:) integer, intent(out) :: kpbl(:) - logical, intent(in) :: gen_tend,ldiag3d + logical, intent(in) :: gen_tend,ldiag3d,progsigma ! real(kind=kind_phys), intent(in) :: grav,rd,cp,rv,hvap,hfus,fv, & & eps,epsm1 @@ -299,7 +299,6 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & xmfd(i,k) = 0. buou(i,k) = 0. buod(i,k) = 0. - tmf(i,k) = 0. ckz(i,k) = ck1 chz(i,k) = ch1 rlmnz(i,k) = rlmn0 @@ -313,6 +312,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & zm(i,k) = zi(i,k+1) enddo enddo +!> - Initialize variables needed for prognostic cumulus closure + if(progsigma)then + do k=1,km + do i=1,im + tmf(i,k) = 0. + enddo + enddo + endif !> - Compute horizontal grid size (\p gdx) do i=1,im gdx(i) = sqrt(garea(i)) @@ -2115,11 +2122,19 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & qtend = (f2(i,k)-q1(i,k,1))*rdt tdt(i,k) = tdt(i,k)+ttend rtg(i,k,1) = rtg(i,k,1)+qtend - tmf(i,k) = qtend ! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend ! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend enddo enddo + + if(progsigma)then + do k = 1,km + do i = 1,im + tmf(i,k)=(f2(i,k)-q1(i,k,1))*rdt + enddo + enddo + endif + do i = 1,im dtsfc(i) = rho_a(i) * cp * heat(i) dqsfc(i) = rho_a(i) * hvap * evap(i) diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index fa30cd9f7..88ab676b8 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -62,6 +62,13 @@ dimensions = () type = integer intent = in +[progsigma] + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumuls scheme + units = flag + dimensions = () + type = logical + intent = in [ntrac] standard_name = number_of_vertical_diffusion_tracers long_name = number of tracers to diffuse vertically From e2d5a2a6310a3aa8313ba9c65ff29dc16bf6b954 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 27 Apr 2022 22:40:40 +0000 Subject: [PATCH 07/11] cleaning out some print statements --- physics/progsigma_calc.f90 | 12 ++++++------ physics/samfdeepcnv.f | 36 +++++++++++++++++++++++------------- physics/samfdeepcnv.meta | 15 ++++++++------- 3 files changed, 37 insertions(+), 26 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 55f6b5e3a..21612fd6c 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -43,7 +43,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & integer :: i,k,km1 real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & mcons(im),fdqa(im),form(im,km), & - qadv(im,km),sigmamax(im),dp(im),inbu(im,km) + qadv(im,km),sigmamax(im),dp(im,km),inbu(im,km) real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & @@ -82,7 +82,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - dp(i) = 1000. * del(i,k) + dp(i,k) = 1000. * del(i,k) endif enddo enddo @@ -128,7 +128,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo enddo - + !compute termD "The vertical integral of the latent heat convergence is limited to the !buoyant layers with positive moisture convergence (accumulated from the surface). !Lowest level: @@ -140,7 +140,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i)) + mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i,k)) buy2 = termD(i)+mcon+mcons(i) ! Do the integral over buoyant layers with positive mcon acc from surface if(k > kbcon1(i) .and. k < ktcon(i) .and. buy2 > 0.)then @@ -157,7 +157,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - tem=(sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k))*dp(i) + tem=(sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k))*dp(i,k) termA(i)=termA(i)+tem endif enddo @@ -167,7 +167,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - tem=(dbyo1(i,k)*inbu(i,k))*dp(i) + tem=(dbyo1(i,k)*inbu(i,k))*dp(i,k) termB(i)=termB(i)+tem endif enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 5fd54a2ec..bf48d2035 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -80,13 +80,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & - & hwrf_samfdeep,progsigma,wclosureflg,cldwrk,rn,kbot,ktop,kcnv, & + & hwrf_samfdeep,progsigma,cldwrk,rn,kbot,ktop,kcnv, & & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, sigmain, sigmaout, & + & rainevap, sigmain, sigmaout, ca_micro, & & errmsg,errflg) ! use machine , only : kind_phys @@ -103,12 +103,12 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:) real(kind=kind_phys), dimension(:), intent(in) :: fscav logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, & - & progsigma, wclosureflg + & progsigma real(kind=kind_phys), intent(in) :: nthresh real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), prevsq(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:) + real(kind=kind_phys), intent(out) :: rainevap(:), ca_micro(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -243,7 +243,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! parameter(cinacrmx=-120.,cinacrmn=-120.) parameter(cinacrmx=-120.,cinacrmn=-80.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) - parameter(betaw=.03,dxcrtas=8.e3,dxcrtuf=15.e3) + parameter(betaw=.03,dxcrtuf=15.e3) ! ! local variables and arrays @@ -380,6 +380,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & advfac(i) = 0. rainevap(i) = 0. omegac(i)=0. + ca_micro(i)=0. gdx(i) = sqrt(garea(i)) enddo @@ -2456,8 +2457,15 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & c c------- final changed variable per unit mass flux c -!> - If grid size is less than a threshold value (dxcrtas: currently 8km), the quasi-equilibrium assumption of Arakawa-Schubert is not used any longer. +!> - If grid size is less than a threshold value (dxcrtas: currently 8km if progsigma is not used and 30km if progsigma is used), the quasi-equilibrium assumption of Arakawa-Schubert is not used any longer. ! + if(progsigma)then + dxcrtas=30.e3 + else + dxcrtas=8.e3 + endif + + do i = 1, im asqecflg(i) = cnvflg(i) if(asqecflg(i) .and. gdx(i) < dxcrtas) then @@ -2465,13 +2473,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & endif enddo -!> - If wclosureflg is true, then quasi-equilibrium closure of Arakawa-Schubert is not used any longer, regardless of resolution - if(wclosureflg)then - do i = 1, im - asqecflg(i) = .false. - enddo - endif - ! !> - If grid size is larger than the threshold value (i.e., asqecflg=.true.), the quasi-equilibrium assumption is used to obtain the cloud base mass flux. To begin with, calculate the change in the temperature and moisture profiles per unit cloud base mass flux. do k = 1, km @@ -2884,6 +2885,15 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. + + if(progsigma)then + do i= 1, im + if(cnvflg(i))then + ca_micro(i)=sigmab(i) + endif + enddo + endif + do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 2b2942812..5e589b318 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -312,13 +312,6 @@ dimensions = () type = logical intent = in -[wclosureflg] - standard_name = flag_for_wclosure - long_name = flag for vertical velocity closure - units = flag - dimensions = () - type = logical - intent = in [progsigma] standard_name = do_prognostic_updraft_area_fraction long_name = flag for prognostic sigma in cumuls scheme @@ -667,6 +660,14 @@ type = real kind = kind_phys intent = out +[ca_micro] + standard_name = output_prognostic_sigma_two + long_name = output of prognostic area fraction two + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 8b815e026ddcd8e660e437dc94bb89058e6f80de Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Fri, 29 Apr 2022 03:03:46 +0000 Subject: [PATCH 08/11] address some bugs caught by debug flag --- physics/progsigma_calc.f90 | 11 ++++++----- physics/samfdeepcnv.f | 34 +++++++++++++++++----------------- physics/samfdeepcnv.meta | 8 -------- physics/samfshalcnv.f | 9 ++++++--- 4 files changed, 29 insertions(+), 33 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 21612fd6c..c05af3003 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -26,8 +26,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! intent in integer, intent(in) :: im,km,kbcon1(im),ktcon(im) - real, intent(in) :: hvap,delt - real, intent(in) :: prevsq(im,km), q(im,km),del(im,km), & + real(kind=kind_phys), intent(in) :: hvap,delt + real(kind=kind_phys), intent(in) :: prevsq(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow @@ -63,7 +63,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do i = 1,im sigmaout(i,k)=0. inbu(i,k)=0. - form(i,k)=0. + form(i,k)=0. + dp(i,k)=0. enddo enddo @@ -157,7 +158,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - tem=(sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k))*dp(i,k) + tem=sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k)*dp(i,k) termA(i)=termA(i)+tem endif enddo @@ -167,7 +168,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - tem=(dbyo1(i,k)*inbu(i,k))*dp(i,k) + tem=zeta(i,k)*dbyo1(i,k)*inbu(i,k)*dp(i,k) termB(i)=termB(i)+tem endif enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index bf48d2035..071bf0557 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -86,8 +86,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, sigmain, sigmaout, ca_micro, & - & errmsg,errflg) + & rainevap, sigmain, sigmaout, errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -108,7 +107,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), prevsq(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:), ca_micro(:) + real(kind=kind_phys), intent(out) :: rainevap(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -380,7 +379,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & advfac(i) = 0. rainevap(i) = 0. omegac(i)=0. - ca_micro(i)=0. gdx(i) = sqrt(garea(i)) enddo @@ -583,14 +581,20 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & buo(i,k) = 0. drag(i,k) = 0. cnvwt(i,k)= 0. + endif + enddo + enddo + + do k = 1, km + do i = 1, im dbyo1(i,k)=0. zdqca(i,k)=0. qlks(i,k)=0. omega_u(i,k)=0. zeta(i,k)=1.0 - endif - enddo + enddo enddo + ! ! initialize tracer variables ! @@ -1811,9 +1815,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do i = 1, im if (cnvflg(i)) then if(k >= kbcon1(i) .and. k < ktcon(i)) then - zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) - zeta(i,k)=MAX(0.,zeta(i,k)) - zeta(i,k)=MIN(1.,zeta(i,k)) + if(omega_u(i,k) .ne. 0.)then + zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + else + zeta(i,k)=0. + endif + zeta(i,k)=MAX(0.,zeta(i,k)) + zeta(i,k)=MIN(1.,zeta(i,k)) endif endif enddo @@ -2886,14 +2894,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. - if(progsigma)then - do i= 1, im - if(cnvflg(i))then - ca_micro(i)=sigmab(i) - endif - enddo - endif - do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 5e589b318..3f28035b6 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -660,14 +660,6 @@ type = real kind = kind_phys intent = out -[ca_micro] - standard_name = output_prognostic_sigma_two - long_name = output of prognostic area fraction two - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 325566877..ef0366b84 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -513,7 +513,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo enddo - do i = 1,im omegac(i)=0. enddo @@ -1551,12 +1550,16 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif enddo c -! > - Calculate the mean updraft velocity within the cloud (omega). +c Compute zeta for prog closure do k = 2, km1 do i = 1, im if (cnvflg(i)) then if(k > kbcon1(i) .and. k < ktcon(i)) then - zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + if(omega_u(i,k) .ne. 0.)then + zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + else + zeta(i,k)=0. + endif zeta(i,k)=MAX(0.,zeta(i,k)) zeta(i,k)=MIN(1.,zeta(i,k)) endif From be534e7be9c335c9d8736fe0d5bf8fe83e14cce9 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Thu, 5 May 2022 04:08:06 +0000 Subject: [PATCH 09/11] addressing some review comments --- physics/progsigma_calc.f90 | 114 +++++++++++++++++++------------------ physics/samfdeepcnv.f | 16 +++++- physics/samfdeepcnv.meta | 8 +++ 3 files changed, 80 insertions(+), 58 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index c05af3003..fe74dc0c1 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -48,7 +48,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & fdqb,dtdyn,dxlim,rmulacvg,tem, & - alpha,DEN,betascu,dp1 + alpha,DEN,betascu,dp1,invdelt !Parameters gcvalmx = 0.1 @@ -57,11 +57,11 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & km1=km-1 alpha=7000. betascu = 3.0 + invdelt = 1./delt !Initialization 2D do k = 1,km do i = 1,im - sigmaout(i,k)=0. inbu(i,k)=0. form(i,k)=0. dp(i,k)=0. @@ -70,7 +70,9 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !Initialization 1D do i=1,im - sigmab(i)=0. + if(cnvflg(i))then + sigmab(i)=0. + endif sigmamax(i)=0.95 termA(i)=0. termB(i)=0. @@ -89,24 +91,16 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & enddo !Initial computations, place maximum sigmain in sigmab - if(flag_init .and. .not. flag_restart)then - do i=1,im - if(cnvflg(i))then - sigmab(i)=0.03 - endif - enddo - else - do i=1,im - if(cnvflg(i))then - do k=2,km - if(sigmain(i,k)>sigmab(i))then - sigmab(i)=sigmain(i,k) - endif - enddo - endif - enddo - endif - + do i=1,im + if(cnvflg(i))then + do k=2,km + if(sigmain(i,k)>sigmab(i))then + sigmab(i)=sigmain(i,k) + endif + enddo + endif + enddo + do i=1,im if(sigmab(i) < 1.E-5)then !after advection sigmab(i)=0. @@ -120,15 +114,20 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & enddo !Initial computations, dynamic q-tendency - do k = 1,km - do i = 1,im - if(flag_init .and. .not.flag_restart)then + if(flag_init .and. .not.flag_restart)then + do k = 1,km + do i = 1,im qadv(i,k)=0. - else - qadv(i,k)=(q(i,k) - prevsq(i,k))/delt - endif + enddo enddo - enddo + else + do k = 1,km + do i = 1,im + qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt + enddo + enddo + endif + !compute termD "The vertical integral of the latent heat convergence is limited to the !buoyant layers with positive moisture convergence (accumulated from the surface). @@ -173,7 +172,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo enddo - + !termC do k = 2,km1 do i = 1,im @@ -188,33 +187,38 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & enddo !sigmab - do i = 1,im - if(cnvflg(i))then - DEN=MIN(termC(i)+termB(i),1.E8) - cvg=termD(i)*delt - ZZ=MAX(0.0,SIGN(1.0,termA(i))) & - *MAX(0.0,SIGN(1.0,termB(i))) & - *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) - cvg=MAX(0.0,cvg) - if(flag_init .and. .not. flag_restart)then - sigmab(i)=0.03 - else - sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) - endif - if(sigmab(i)>0.)then - sigmab(i)=MIN(sigmab(i),sigmamax(i)) - sigmab(i)=MAX(sigmab(i),0.01) - endif - endif!cnvflg - enddo + if(flag_init .and. .not. flag_restart)then + do i = 1,im + if(cnvflg(i))then + sigmab(i)=0.03 + endif + enddo + else + do i = 1,im + if(cnvflg(i))then + DEN=MIN(termC(i)+termB(i),1.E8) + cvg=termD(i)*delt + ZZ=MAX(0.0,SIGN(1.0,termA(i))) & + *MAX(0.0,SIGN(1.0,termB(i))) & + *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) + cvg=MAX(0.0,cvg) + sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) + if(sigmab(i)>0.)then + sigmab(i)=MIN(sigmab(i),sigmamax(i)) + sigmab(i)=MAX(sigmab(i),0.01) + endif + endif!cnvflg + enddo + endif - do k=1,km - do i=1,im - if(cnvflg(i))then - sigmaout(i,k)=sigmab(i) - endif - enddo - enddo + do k=1,km + do i=1,im + if(cnvflg(i))then + sigmaout(i,k)=sigmab(i) + endif + enddo + enddo + !Since updraft velocity is much lower in shallow cu region, termC becomes small in shallow cu application, thus the area fraction !in this regime becomes too large compared with the deep cu region. To address this simply apply a scaling factor for shallow cu diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 071bf0557..8398af769 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -85,8 +85,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & - & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, sigmain, sigmaout, errmsg,errflg) + & do_ca, ca_closure, ca_entr, ca_trigger, nthresh,ca_deep, & + & ca_micro, rainevap, sigmain, sigmaout, errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -107,7 +107,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), prevsq(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:) + real(kind=kind_phys), intent(out) :: rainevap(:), ca_micro(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -2894,6 +2894,16 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. + do i=1,im + ca_micro(i)=0. + enddo + + do i=1,im + if(cnvflg(i))then + ca_micro(i)=sigmab(i) + endif + enddo + do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 3f28035b6..1764a74fd 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -652,6 +652,14 @@ type = real kind = kind_phys intent = in +[ca_micro] + standard_name = output_prognostic_sigma_two + long_name = output of prognostic area fraction two + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [rainevap] standard_name = physics_field_for_coupling long_name = physics_field_for_coupling From 6f38cc6e83526f907ec2c34712eb230d677f4a2a Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 18 May 2022 23:58:09 +0000 Subject: [PATCH 10/11] address some review comments, fix decomposition error, correct bug in initialization --- physics/GFS_suite_interstitial_3.F90 | 37 +++++++++-- physics/GFS_suite_interstitial_3.meta | 73 +++++++++++++++++++++ physics/progsigma_calc.f90 | 93 ++++++++++++--------------- physics/samfdeepcnv.f | 22 ++----- physics/samfdeepcnv.meta | 8 --- physics/samfshalcnv.f | 20 +++--- 6 files changed, 162 insertions(+), 91 deletions(-) diff --git a/physics/GFS_suite_interstitial_3.F90 b/physics/GFS_suite_interstitial_3.F90 index 79ab481ec..9fa7d69b7 100644 --- a/physics/GFS_suite_interstitial_3.F90 +++ b/physics/GFS_suite_interstitial_3.F90 @@ -9,10 +9,13 @@ module GFS_suite_interstitial_3 !! \htmlinclude GFS_suite_interstitial_3_run.html !! subroutine GFS_suite_interstitial_3_run (otsptflag, & - im, levs, nn, cscnv, & + im, levs, nn, cscnv,imfshalcnv, imfdeepcnv, & + imfshalcnv_samf, imfdeepcnv_samf,progsigma, & + first_time_step, restart, & satmedmf, trans_trac, do_shoc, ltaerosol, ntrac, ntcw, & ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, & - xlon, xlat, gt0, gq0, imp_physics, imp_physics_mg, & + xlon, xlat, gt0, gq0, sigmain,sigmaout,qmicro, & + 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, & @@ -33,8 +36,9 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6,imp_physics_fer_hires, & imp_physics_nssl, me, index_of_process_conv_trans integer, intent(in ), dimension(:) :: islmsk, kpbl, kinver - logical, intent(in ) :: cscnv, satmedmf, trans_trac, do_shoc, ltaerosol, ras - + logical, intent(in ) :: cscnv, satmedmf, trans_trac, do_shoc, ltaerosol, ras, progsigma + logical, intent(in ) :: first_time_step, restart + integer, intent(in ) :: imfshalcnv, imfdeepcnv, imfshalcnv_samf,imfdeepcnv_samf integer, intent(in) :: ntinc, ntlnc logical, intent(in) :: ldiag3d, qdiag3d integer, dimension(:,:), intent(in) :: dtidx @@ -48,6 +52,8 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & real(kind=kind_phys), intent(in ), dimension(:,:) :: gt0 real(kind=kind_phys), intent(in ), dimension(:,:,:) :: gq0 + real(kind=kind_phys), intent(out ), dimension(:,:) :: sigmain + real(kind=kind_phys), intent(out ), dimension(:,:) :: sigmaout,qmicro 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 @@ -73,6 +79,27 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & errmsg = '' errflg = 0 + ! In case of using prognostic updraf area fraction, initialize area fraction here + ! since progsigma_calc is called from both deep and shallow schemes. + if(((imfshalcnv == imfshalcnv_samf) .or. (imfdeepcnv == imfdeepcnv_samf)) & + .and. progsigma)then + if(first_time_step .and. .not. restart)then + do k=1,levs + do i=1,im + sigmain(i,k)=0.0 + sigmaout(i,k)=0.0 + qmicro(i,k)=0.0 + enddo + enddo + endif + do k=1,levs + do i=1,im + sigmaout(i,k)=0.0 + enddo + enddo + endif + + if (cscnv .or. satmedmf .or. trans_trac .or. ras) then tracers = 2 do n=2,ntrac @@ -192,4 +219,4 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & end subroutine GFS_suite_interstitial_3_run - end module GFS_suite_interstitial_3 \ No newline at end of file + end module GFS_suite_interstitial_3 diff --git a/physics/GFS_suite_interstitial_3.meta b/physics/GFS_suite_interstitial_3.meta index 22a11d0ea..fbeb9f03c 100644 --- a/physics/GFS_suite_interstitial_3.meta +++ b/physics/GFS_suite_interstitial_3.meta @@ -43,6 +43,55 @@ dimensions = () type = logical intent = in +[imfdeepcnv] + standard_name = control_for_deep_convection_scheme + long_name = flag for mass-flux deep convection scheme + units = flag + dimensions = () + type = integer + intent = in +[imfdeepcnv_samf] + standard_name = identifer_for_scale_aware_mass_flux_deep_convection + long_name = flag for SAMF deep convection scheme + units = flag + dimensions = () + type = integer + intent = in +[imfshalcnv] + standard_name = control_for_shallow_convection_scheme + long_name = flag for mass-flux shallow convection scheme + units = flag + dimensions = () + type = integer + intent = in +[imfshalcnv_samf] + standard_name = identifier_for_scale_aware_mass_flux_shallow_convection + long_name = flag for SAMF shallow convection scheme + units = flag + dimensions = () + type = integer + intent = in +[progsigma] + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumuls scheme + units = flag + dimensions = () + type = logical + intent = in +[first_time_step] + standard_name = flag_for_first_timestep + long_name = flag for first time step for time integration loop (cold/warmstart) + units = flag + dimensions = () + type = logical + intent = in +[restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in [satmedmf] standard_name = flag_for_scale_aware_TKE_moist_EDMF_PBL long_name = flag for scale-aware TKE moist EDMF PBL scheme @@ -173,6 +222,30 @@ type = real kind = kind_phys intent = in +[sigmain] + standard_name = prognostic_updraft_area_fraction_in_convection + long_name = convective updraft area fraction + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out +[sigmaout] + standard_name = updraft_area_fraction_updated_by_physics + long_name = convective updraft area fraction updated by physics + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out +[qmicro] + standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics + long_name = moisture tendency due to microphysics + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [imp_physics] standard_name = control_for_microphysics_scheme long_name = choice of microphysics scheme diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index fe74dc0c1..49a5e2a4f 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -14,9 +14,9 @@ !> @{ subroutine progsigma_calc (im,km,flag_init,flag_restart, & - flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & - delt,prevsq,q,kbcon1,ktcon,cnvflg,gdx, & - sigmain,sigmaout,sigmab,errmsg,errflg) + del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & + delt,prevsq,q,kbcon1,ktcon,cnvflg,sigmain,sigmaout, & + sigmab,errmsg,errflg) ! ! use machine, only : kind_phys @@ -29,8 +29,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & real(kind=kind_phys), intent(in) :: hvap,delt real(kind=kind_phys), intent(in) :: prevsq(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & - omega_u(im,km),zeta(im,km),gdx(im) - logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow + omega_u(im,km),zeta(im,km) + logical, intent(in) :: flag_init,flag_restart,cnvflg(im) real(kind=kind_phys), intent(in) :: sigmain(im,km) ! intent out @@ -41,21 +41,20 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! Local variables integer :: i,k,km1 - real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & - mcons(im),fdqa(im),form(im,km), & - qadv(im,km),sigmamax(im),dp(im,km),inbu(im,km) + real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im) + real(kind=kind_phys) :: mcons(im),fdqa(im),form(im,km), & + qadv(im,km),dp(im,km),inbu(im,km) real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & fdqb,dtdyn,dxlim,rmulacvg,tem, & - alpha,DEN,betascu,dp1,invdelt + DEN,betascu,dp1,invdelt !Parameters gcvalmx = 0.1 rmulacvg=10. epsilon=1.E-11 km1=km-1 - alpha=7000. betascu = 3.0 invdelt = 1./delt @@ -70,10 +69,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !Initialization 1D do i=1,im - if(cnvflg(i))then - sigmab(i)=0. - endif - sigmamax(i)=0.95 + sigmab(i)=0. termA(i)=0. termB(i)=0. termC(i)=0. @@ -82,6 +78,21 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & mcons(i)=0. enddo + !Initial computations, dynamic q-tendency + if(flag_init .and. .not.flag_restart)then + do k = 1,km + do i = 1,im + qadv(i,k)=0. + enddo + enddo + else + do k = 1,km + do i = 1,im + qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt + enddo + enddo + endif + do k = 2,km1 do i = 1,im if(cnvflg(i))then @@ -102,33 +113,13 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & enddo do i=1,im - if(sigmab(i) < 1.E-5)then !after advection - sigmab(i)=0. + if(cnvflg(i))then + if(sigmab(i) < 1.E-5)then !after advection + sigmab(i)=0. + endif endif enddo - - !Initial computations, sigmamax - do i=1,im - sigmamax(i)=alpha/gdx(i) - sigmamax(i)=MIN(0.95,sigmamax(i)) - enddo - - !Initial computations, dynamic q-tendency - if(flag_init .and. .not.flag_restart)then - do k = 1,km - do i = 1,im - qadv(i,k)=0. - enddo - enddo - else - do k = 1,km - do i = 1,im - qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt - enddo - enddo - endif - - + !compute termD "The vertical integral of the latent heat convergence is limited to the !buoyant layers with positive moisture convergence (accumulated from the surface). !Lowest level: @@ -194,7 +185,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo else - do i = 1,im + do i = 1,im if(cnvflg(i))then DEN=MIN(termC(i)+termB(i),1.E8) cvg=termD(i)*delt @@ -204,7 +195,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & cvg=MAX(0.0,cvg) sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) if(sigmab(i)>0.)then - sigmab(i)=MIN(sigmab(i),sigmamax(i)) + sigmab(i)=MIN(sigmab(i),0.95) sigmab(i)=MAX(sigmab(i),0.01) endif endif!cnvflg @@ -218,19 +209,15 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo enddo - - !Since updraft velocity is much lower in shallow cu region, termC becomes small in shallow cu application, thus the area fraction - !in this regime becomes too large compared with the deep cu region. To address this simply apply a scaling factor for shallow cu - !before computing the massflux to reduce the total strength of the SC MF: - - if(flag_shallow)then - do i= 1, im - if(cnvflg(i)) then - sigmab(i)=sigmab(i)/betascu - endif - enddo - endif + !Reduce area fraction before coupling back to mass-flux computation. + !This tuning could be addressed in updraft velocity equation instead. + do i= 1, im + if(cnvflg(i)) then + sigmab(i)=sigmab(i)/betascu + endif + enddo + end subroutine progsigma_calc diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 8398af769..3968061a2 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -86,7 +86,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh,ca_deep, & - & ca_micro, rainevap, sigmain, sigmaout, errmsg,errflg) + & rainevap,sigmain, sigmaout, errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -107,7 +107,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), prevsq(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:), ca_micro(:) + real(kind=kind_phys), intent(out) :: rainevap(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -214,8 +214,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),qlks(im,km), - & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) - logical flag_shallow + & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) real(kind=kind_phys) gravinv c physical parameters ! parameter(grav=grav,asolfac=0.958) @@ -2884,25 +2883,14 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - From Bengtsson et al. (2022) Prognostic closure scheme, equation 8, compute updraft area fraction based on a moisture budget if(progsigma)then - flag_shallow = .false. - call progsigma_calc(im,km,first_time_step,restart,flag_shallow, + call progsigma_calc(im,km,first_time_step,restart, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & prevsq,q,kbcon1,ktcon,cnvflg,gdx, + & prevsq,q,kbcon1,ktcon,cnvflg, & sigmain,sigmaout,sigmab,errmsg,errflg) endif !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. - - do i=1,im - ca_micro(i)=0. - enddo - - do i=1,im - if(cnvflg(i))then - ca_micro(i)=sigmab(i) - endif - enddo do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 1764a74fd..3f28035b6 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -652,14 +652,6 @@ type = real kind = kind_phys intent = in -[ca_micro] - standard_name = output_prognostic_sigma_two - long_name = output of prognostic area fraction two - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [rainevap] standard_name = physics_field_for_coupling long_name = physics_field_for_coupling diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index ef0366b84..b7943725b 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -62,7 +62,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & - & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, + & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, & & sigmain,sigmaout,errmsg,errflg) ! use machine , only : kind_phys @@ -77,7 +77,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), & - & qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:),sigmain(:,:) + & qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:) + + real(kind=kind_phys), intent(in) :: sigmain(:,:) ! real(kind=kind_phys), dimension(:), intent(in) :: fscav integer, intent(inout) :: kcnv(:) @@ -164,8 +166,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), & sigmab(im) - logical flag_shallow - real(kind=kind_phys) gravinv + real(kind=kind_phys) gravinv,dxcrtas c physical parameters ! parameter(g=grav,asolfac=0.89) @@ -196,6 +197,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) parameter(betaw=.03,dxcrt=15.e3,dxcrtc0=9.e3) parameter(h1=0.33333333) +! progsigma + parameter(dxcrtas=30.e3) 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), @@ -340,6 +343,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo endif !! + !> - Return to the calling routine if deep convection is present or the surface buoyancy flux is negative. totflg = .true. do i=1,im @@ -1932,20 +1936,20 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c c Prognostic closure if(progsigma)then - flag_shallow = .true. - call progsigma_calc(im,km,first_time_step,restart,flag_shallow, + call progsigma_calc(im,km,first_time_step,restart, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & prevsq,q,kbcon1,ktcon,cnvflg,gdx, + & prevsq,q,kbcon1,ktcon,cnvflg, & sigmain,sigmaout,sigmab,errmsg,errflg) endif !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. + do i= 1, im if(cnvflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) - if(progsigma)then + if(progsigma .and. gdx(i) < dxcrtas)then xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv) else xmb(i) = advfac(i)*betaw*rho*wc(i) From fc79cc39ecaef4bc567b052842dc560e3253f0d8 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Thu, 19 May 2022 19:59:43 +0000 Subject: [PATCH 11/11] Change intent to inout for conditional variables --- physics/GFS_MP_generic_post.F90 | 4 ++-- physics/GFS_MP_generic_post.meta | 4 ++-- physics/GFS_suite_interstitial_3.F90 | 4 ++-- physics/GFS_suite_interstitial_3.meta | 4 ++-- physics/satmedmfvdifq.F | 4 ++-- physics/satmedmfvdifq.meta | 2 +- 6 files changed, 11 insertions(+), 11 deletions(-) diff --git a/physics/GFS_MP_generic_post.F90 b/physics/GFS_MP_generic_post.F90 index 3cdfe91f3..992c9b969 100644 --- a/physics/GFS_MP_generic_post.F90 +++ b/physics/GFS_MP_generic_post.F90 @@ -81,8 +81,8 @@ subroutine GFS_MP_generic_post_run( real(kind=kind_phys), dimension(:), intent(inout) :: diceprv real(kind=kind_phys), dimension(:), intent(inout) :: dsnowprv real(kind=kind_phys), dimension(:), intent(inout) :: dgraupelprv - real(kind=kind_phys), dimension(:,:), intent(out) :: dqdt_qmicro - real(kind=kind_phys), dimension(:,:), intent(out) :: prevsq + real(kind=kind_phys), dimension(:,:), intent(inout) :: dqdt_qmicro + real(kind=kind_phys), dimension(:,:), intent(inout) :: prevsq real(kind=kind_phys), intent(in) :: dtp ! CCPP error handling diff --git a/physics/GFS_MP_generic_post.meta b/physics/GFS_MP_generic_post.meta index 20881fbb3..7ba09363a 100644 --- a/physics/GFS_MP_generic_post.meta +++ b/physics/GFS_MP_generic_post.meta @@ -738,7 +738,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [prevsq] standard_name = specific_humidity_on_previous_timestep long_name = specific_humidity_on_previous_timestep @@ -746,7 +746,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [lssav] standard_name = flag_for_diagnostics long_name = logical flag for storing diagnostics diff --git a/physics/GFS_suite_interstitial_3.F90 b/physics/GFS_suite_interstitial_3.F90 index 9fa7d69b7..0f666d4a5 100644 --- a/physics/GFS_suite_interstitial_3.F90 +++ b/physics/GFS_suite_interstitial_3.F90 @@ -52,8 +52,8 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & real(kind=kind_phys), intent(in ), dimension(:,:) :: gt0 real(kind=kind_phys), intent(in ), dimension(:,:,:) :: gq0 - real(kind=kind_phys), intent(out ), dimension(:,:) :: sigmain - real(kind=kind_phys), intent(out ), dimension(:,:) :: sigmaout,qmicro + real(kind=kind_phys), intent(inout ), dimension(:,:) :: sigmain + real(kind=kind_phys), intent(inout ), dimension(:,:) :: sigmaout,qmicro 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 diff --git a/physics/GFS_suite_interstitial_3.meta b/physics/GFS_suite_interstitial_3.meta index fbeb9f03c..fe52a1adc 100644 --- a/physics/GFS_suite_interstitial_3.meta +++ b/physics/GFS_suite_interstitial_3.meta @@ -229,7 +229,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [sigmaout] standard_name = updraft_area_fraction_updated_by_physics long_name = convective updraft area fraction updated by physics @@ -237,7 +237,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [qmicro] standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics long_name = moisture tendency due to microphysics diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index c7a6fadc9..865e4481c 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -99,7 +99,7 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & real(kind=kind_phys), intent(in) :: dspfac, bl_upfr, bl_dnfr real(kind=kind_phys), intent(in) :: rlmx, elmx real(kind=kind_phys), intent(inout) :: dv(:,:), du(:,:), & - & tdt(:,:), rtg(:,:,:) + & tdt(:,:), rtg(:,:,:), tmf(:,:) real(kind=kind_phys), intent(in) :: & & u1(:,:), v1(:,:), & & t1(:,:), q1(:,:,:), & @@ -123,7 +123,7 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & & dtsfc(:), dqsfc(:), & & hpbl(:) real(kind=kind_phys), intent(out) :: & - & dkt(:,:), dku(:,:), tmf(:,:) + & dkt(:,:), dku(:,:) ! logical, intent(in) :: dspheat character(len=*), intent(out) :: errmsg diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 88ab676b8..8538c6aa7 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -215,7 +215,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [u1] standard_name = x_wind long_name = x component of layer wind