From be146a8ebc073412b1e0364ad6225dbdef2f26fe Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Fri, 24 Mar 2023 01:36:23 +0000 Subject: [PATCH 1/3] Correction to convective cloud condensate term in prog closure for HR2 --- physics/progsigma_calc.f90 | 28 ++++-------------- physics/samfdeepcnv.f | 59 +++++++++++++++++++++----------------- physics/samfdeepcnv.meta | 6 ++-- physics/samfshalcnv.f | 54 ++++++++++++++++++++-------------- physics/samfshalcnv.meta | 6 ++-- physics/satmedmfvdifq.F | 24 +++------------- physics/satmedmfvdifq.meta | 15 ---------- 7 files changed, 81 insertions(+), 111 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index eaa1d3fda..4bbd305ae 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -13,8 +13,8 @@ !!\section gen_progsigma progsigma_calc General Algorithm 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,sigmain,sigmaout, & - sigmab,errmsg,errflg) + delt,qadv,kbcon1,ktcon,cnvflg,sigmain,sigmaout, & + sigmab) ! ! use machine, only : kind_phys @@ -25,7 +25,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! intent in integer, intent(in) :: im,km,kbcon1(im),ktcon(im) real(kind=kind_phys), intent(in) :: hvap,delt - real(kind=kind_phys), intent(in) :: prevsq(im,km), q(im,km),del(im,km), & + real(kind=kind_phys), intent(in) :: qadv(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km) logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow @@ -34,14 +34,13 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! 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) real(kind=kind_phys) :: mcons(im),fdqa(im),form(im,km), & - qadv(im,km),dp(im,km),inbu(im,km) + dp(im,km),inbu(im,km) real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & @@ -77,21 +76,6 @@ 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 @@ -133,7 +117,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & 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 + if(dbyo1(i,k)>0 .and. buy2 > 0.)then inbu(i,k)=1. endif inbu(i,k-1)=MAX(inbu(i,k-1),inbu(i,k)) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 2a3c256a9..cd130dfd0 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -102,7 +102,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(:,:), prevsq(:,:) + & 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 @@ -209,9 +209,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & 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) - real(kind=kind_phys) gravinv + real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), + & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km) + real(kind=kind_phys) gravinv,invdelt logical flag_shallow c physical parameters ! parameter(grav=grav,asolfac=0.958) @@ -306,6 +306,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & errflg = 0 gravinv = 1./grav + invdelt = 1./delt elocp = hvap/cp el2orc = hvap*hvap/(rv*cp) @@ -585,7 +586,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & 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 @@ -1515,7 +1515,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & 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 + zdqca(i,k)=dq/eta(i,k) endif ! ! compute buoyancy and drag for updraft velocity @@ -1690,7 +1690,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & 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 + zdqca(i,k)=dq/eta(i,k) endif endif endif @@ -1860,28 +1860,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if(dq > 0.) then qlko_ktcon(i) = dq qcko(i,k) = qrch - qlks(i,k) = qlko_ktcon(i) + zdqca(i,k) = dq endif endif enddo endif c -c store term needed for "termC" in prognostic area fraction closure - 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 - enddo - enddo - endif - ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then ccccc print *, ' aa1(i) before dwndrft =', aa1(i) ccccc endif @@ -2885,11 +2870,33 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - From Bengtsson et al. (2022) \cite Bengtsson_2022 prognostic closure scheme, equation 8, call progsigma_calc() to compute updraft area fraction based on a moisture budget if(progsigma)then + +!Initial computations, dynamic q-tendency + if(first_time_step .and. .not.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 = 1,km + do i = 1,im + tmfq(i,k)=tmf(i,k,1) + enddo + enddo + flag_shallow = .false. call progsigma_calc(im,km,first_time_step,restart,flag_shallow, - & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & prevsq,q,kbcon1,ktcon,cnvflg, - & sigmain,sigmaout,sigmab,errmsg,errflg) + & del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, + & qadv,kbcon1,ktcon,cnvflg, + & sigmain,sigmaout,sigmab) 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 3f28035b6..bed4d655d 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -70,10 +70,10 @@ type = logical intent = in [tmf] - standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL - long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL + standard_name = tendency_of_vertically_diffused_tracer_concentration + long_name = updated tendency of the tracers due to vertical diffusion in PBL scheme units = kg kg-1 s-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_vertical_diffusion_tracers) type = real kind = kind_phys intent = in diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 645024536..ab25e9922 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -70,7 +70,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(:,:),prevsq(:,:),q(:,:) + & qmicro(:,:),tmf(:,:,:),prevsq(:,:),q(:,:) real(kind=kind_phys), intent(in) :: sigmain(:,:) ! @@ -156,10 +156,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & cc ! parameters for prognostic sigma closure - real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), + real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), - & sigmab(im) - real(kind=kind_phys) gravinv,dxcrtas + & sigmab(im),qadv(im,km) + real(kind=kind_phys) gravinv,dxcrtas,invdelt logical flag_shallow c physical parameters ! parameter(g=grav,asolfac=0.89) @@ -247,6 +247,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & errflg = 0 gravinv = 1./grav + invdelt = 1./delt elocp = hvap/cp el2orc = hvap*hvap/(rv*cp) @@ -524,7 +525,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & 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 @@ -1270,7 +1270,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 + zdqca(i,k)=dq/eta(i,k) endif ! ! compute buoyancy and drag for updraft velocity @@ -1435,7 +1435,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 + zdqca(i,k)=dq/eta(i,k) endif endif endif @@ -1601,24 +1601,13 @@ 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) + zdqca(i,k) = dq endif endif enddo endif c - do k = 2, km1 - do i = 1, im - 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 : @@ -1935,11 +1924,32 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c !> - From Bengtsson et al. (2022) \cite Bengtsson_2022 prognostic closure scheme, equation 8, call progsigma_calc() to compute updraft area fraction based on a moisture budget if(progsigma)then +! Initial computations, dynamic q-tendency + if(first_time_step .and. .not.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 = 1,km + do i = 1,im + tmfq(i,k)=tmf(i,k,1) + enddo + enddo + flag_shallow = .true. call progsigma_calc(im,km,first_time_step,restart,flag_shallow, - & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & prevsq,q,kbcon1,ktcon,cnvflg, - & sigmain,sigmaout,sigmab,errmsg,errflg) + & del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, + & qadv,kbcon1,ktcon,cnvflg, + & sigmain,sigmaout,sigmab) 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 8c9735c32..c1fffef58 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -70,10 +70,10 @@ type = logical intent = in [tmf] - standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL - long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL + standard_name = tendency_of_vertically_diffused_tracer_concentration + long_name = updated tendency of the tracers due to vertical diffusion in PBL scheme units = kg kg-1 s-1 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_vertical_diffusion_tracers) type = real kind = kind_phys intent = in diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 08876f8f0..0387185e4 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -73,9 +73,9 @@ end subroutine satmedmfvdifq_init !! -# A mass-flux approach is also used to represent the stratocumulus-top-induced turbulence !! (mfscuq.f). !! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm - subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & + subroutine satmedmfvdifq_run(im,km,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, & + & dv,du,tdt,rtg,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, & @@ -98,7 +98,7 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & integer, intent(in) :: tc_pbl integer, intent(in) :: kinver(:) integer, intent(out) :: kpbl(:) - logical, intent(in) :: gen_tend,ldiag3d,progsigma + logical, intent(in) :: gen_tend,ldiag3d ! real(kind=kind_phys), intent(in) :: grav,rd,cp,rv,hvap,hfus,fv, & & eps,epsm1 @@ -106,7 +106,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(:,:,:), tmf(:,:) + & tdt(:,:), rtg(:,:,:) real(kind=kind_phys), intent(in) :: & & u1(:,:), v1(:,:), & & t1(:,:), q1(:,:,:), & @@ -331,14 +331,6 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & 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)) @@ -2206,14 +2198,6 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & 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 d9ab8c859..d0b11656a 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -62,13 +62,6 @@ 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 @@ -208,14 +201,6 @@ type = real kind = kind_phys intent = inout -[tmf] - 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 - kind = kind_phys - intent = inout [u1] standard_name = x_wind long_name = x component of layer wind From 6424df169766cc8c7916a6023149f50f7478db8f Mon Sep 17 00:00:00 2001 From: Michael Toy Date: Tue, 4 Apr 2023 01:42:04 +0000 Subject: [PATCH 2/3] Fixed issue with level of dividing streamline (rdxzb) being overwritten and affecting stochastic physics in ugwp --- physics/drag_suite.F90 | 11 +++-------- physics/drag_suite.meta | 2 +- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 5cb49acff..eb51a30c5 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -9,12 +9,6 @@ module drag_suite !> \defgroup gfs_drag_suite_mod GSL drag_suite Module !> This module contains the CCPP-compliant GSL orographic gravity wave drag scheme. !> @{ -!! -!> \brief This subroutine initializes the orographic gravity wave drag scheme. -!! -!> \section arg_table_drag_suite_init Argument Table -!! \htmlinclude drag_suite_init.html -!! subroutine drag_suite_init(gwd_opt, errmsg, errflg) integer, intent(in) :: gwd_opt @@ -342,7 +336,7 @@ subroutine drag_suite_run( & real(kind=kind_phys), intent(inout) :: & & dudt(:,:),dvdt(:,:), & & dtdt(:,:) - real(kind=kind_phys), intent(out) :: rdxzb(:) + real(kind=kind_phys), intent(inout) :: rdxzb(:) real(kind=kind_phys), intent(in) :: & & u1(:,:),v1(:,:), & & t1(:,:),q1(:,:), & @@ -605,7 +599,6 @@ subroutine drag_suite_run( & else xland(i)=2.0 endif - RDXZB(i) = 0.0 enddo !--- calculate scale-aware tapering factors @@ -818,6 +811,8 @@ subroutine drag_suite_run( & do i=its,im + RDXZB(i) = 0.0 + if ( ls_taper(i).GT.1.E-02 ) then ! diff --git a/physics/drag_suite.meta b/physics/drag_suite.meta index ff60290ae..66f320b98 100644 --- a/physics/drag_suite.meta +++ b/physics/drag_suite.meta @@ -536,7 +536,7 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys - intent = out + intent = inout [dx] standard_name = characteristic_grid_lengthscale long_name = size of the grid cell From eb816e3c49ea94cb45dd0da00f099864bec8c66e Mon Sep 17 00:00:00 2001 From: Michael Toy Date: Wed, 5 Apr 2023 03:31:09 +0000 Subject: [PATCH 3/3] Added back CCPP header --- physics/drag_suite.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index eb51a30c5..22f122e71 100644 --- a/physics/drag_suite.F90 +++ b/physics/drag_suite.F90 @@ -9,6 +9,12 @@ module drag_suite !> \defgroup gfs_drag_suite_mod GSL drag_suite Module !> This module contains the CCPP-compliant GSL orographic gravity wave drag scheme. !> @{ +!! +!> \brief This subroutine initializes the orographic gravity wave drag scheme. +!! +!> \section arg_table_drag_suite_init Argument Table +!! \htmlinclude drag_suite_init.html +!! subroutine drag_suite_init(gwd_opt, errmsg, errflg) integer, intent(in) :: gwd_opt