Skip to content

Introduce prognostic updraft velocity to saSAS and C3 cumulus convection schemes #246

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rhaesung merged 13 commits into from
Mar 24, 2025
Merged

Introduce prognostic updraft velocity to saSAS and C3 cumulus convection schemes #246

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
a69aa8b
Introduce prognostic updraft velocity
lisa-bengtsson Jan 23, 2025
d333da7
Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physi…
lisa-bengtsson Jan 28, 2025
3629c45
Address review comment from Jongil Han for a more stable solution of …
lisa-bengtsson Jan 29, 2025
177371a
Fix ascii charachter in French name not liked by CCPP
lisa-bengtsson Jan 29, 2025
0f51c08
After review from Jongil Han, update initialization value in case of …
lisa-bengtsson Jan 29, 2025
3b2bedb
Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physi…
lisa-bengtsson Feb 4, 2025
02d73b6
Minor cleanup of white spaces
lisa-bengtsson Feb 6, 2025
e920f4a
sync branch upstream with ccpp/physics ufs/dev
lisa-bengtsson Feb 19, 2025
375aa85
update code to address debug mode error
lisa-bengtsson Feb 25, 2025
43807e6
Address review comments from Grant Firl on inten inout
lisa-bengtsson Feb 26, 2025
2c1887d
minor clean up in convection routines
lisa-bengtsson Feb 27, 2025
bdd8baf
pull in the latest ufs/dev branch
lisa-bengtsson Mar 7, 2025
b977174
Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physi…
lisa-bengtsson Mar 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
123 changes: 74 additions & 49 deletions physics/CONV/SAMF/samfdeepcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,8 @@ module samfdeepcnv

use samfcnv_aerosols, only : samfdeepcnv_aerosols
use progsigma, only : progsigma_calc

use progomega, only : progomega_calc

contains

subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, &
Expand Down Expand Up @@ -77,14 +78,14 @@ 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,cldwrk,rn,kbot,ktop,kcnv, &
& hwrf_samfdeep,progsigma,progomega,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,betadcu,betamcu,betascu, &
& maxMF, do_mynnedmf,errmsg,errflg)
& rainevap,sigmain,sigmaout,omegain,omegaout, &
& betadcu,betamcu,betascu,maxMF,do_mynnedmf,errmsg,errflg)
!
use machine , only : kind_phys
use funcphys , only : fpvs
Expand All @@ -100,16 +101,17 @@ 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,do_mynnedmf
& progsigma,progomega,do_mynnedmf
real(kind=kind_phys), intent(in) :: nthresh,betadcu,betamcu, &
& betascu
real(kind=kind_phys), intent(in), optional :: ca_deep(:)
real(kind=kind_phys), intent(in), optional :: sigmain(:,:), &
& qmicro(:,:), prevsq(:,:)
& qmicro(:,:), prevsq(:,:), omegain(:,:)
real(kind=kind_phys), intent(in) :: tmf(:,:,:),q(:,:)
real(kind=kind_phys), dimension (:), intent(in), optional :: maxMF
real(kind=kind_phys), intent(out) :: rainevap(:)
real(kind=kind_phys), intent(out), optional :: sigmaout(:,:)
real(kind=kind_phys), intent(out), optional :: sigmaout(:,:), &

@grantfirl grantfirl Feb 12, 2025

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sigmaout and omegaout are intent(out) but are not always give a value if present, which is a CCPP rule. The fact that they're optional maybe makes that rule cloudy, though. I'm thinking that it is probably best to add if present(sigmaout) sigmaout = 0.0_kind_phys (same for omegaout). @dustinswales You're more involved with optional arguments, do you think this is necessary?

& omegaout(:,:)
logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger
integer, intent(inout) :: kcnv(:)
! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH
Expand Down Expand Up @@ -216,7 +218,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
! parameters for prognostic sigma closure
real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
& omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km),
& sigmaoutx(im)
& sigmaoutx(im),tentr(im,km)
real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins
parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01)
logical flag_shallow, flag_mid
Expand Down Expand Up @@ -333,6 +335,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
c-----------------------------------------------------------------------
!> ## Compute preliminary quantities needed for static, dynamic, and feedback control portions of the algorithm.
!> - Convert input pressure terms to centibar units.

!************************************************************************
! convert input Pa terms to Cb terms -- Moorthi
ps = psp * 0.001
Expand Down Expand Up @@ -1131,7 +1134,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
do k = 2, km1
do i=1,im
if(cnvflg(i) .and.
& (k > kbcon(i) .and. k < kmax(i))) then
& (k > kbcon(i) .and. k < kmax(i))) then
tentr(i,k)=xlamue(i,k)*fent1(i,k)
tem = cxlamet(i) * frh(i,k) * fent2(i,k)
xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
tem1 = cxlamdt(i) * frh(i,k)
Expand Down Expand Up @@ -1743,43 +1747,62 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
enddo
!
! compute updraft velocity square(wu2)
!> - Calculate updraft velocity square(wu2) according to Han et al.'s (2017) \cite han_et_al_2017 equation 7.
!
!> - Calculate diagnostic updraft velocity square(wu2) according to Han et al.'s (2017) \cite han_et_al_2017 equation 7.
!> - if progomega = true, calculate prognostic updraft velocity (Pa/s) according to progomega routine.

if (hwrf_samfdeep) then
do i = 1, im
if (cnvflg(i)) then
k = kbcon1(i)
tem = po(i,k) / (rd * to(i,k))
wucb = -0.01 * dot(i,k) / (tem * grav)
if(wucb.gt.0.) then
wu2(i,k) = wucb * wucb
else
wu2(i,k) = 0.
endif
endif
enddo
endif
!
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz
tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k))
tem2 = wush(i,k) * sqrt(wu2(i,k-1))
tem2 = (tem1 - tem2) * dz
ptem = (1. - tem) * wu2(i,k-1)
ptem1 = 1. + tem
wu2(i,k) = (ptem + tem2) / ptem1
wu2(i,k) = max(wu2(i,k), 0.)
do i = 1, im
if (cnvflg(i)) then
k = kbcon1(i)
tem = po(i,k) / (rd * to(i,k))
wucb = -0.01 * dot(i,k) / (tem * grav)
if(wucb.gt.0.) then
wu2(i,k) = wucb * wucb
else
wu2(i,k) = 0.
endif
endif
endif
enddo
enddo

if(progsigma)then
do k = 2, km1
enddo
endif
!
if (progomega) then
call progomega_calc(first_time_step,restart,im,km,
& kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout,
& grav,buo,drag,wush,tentr,bb1,bb2)
do k = 1, km
do i = 1, im
if (cnvflg(i)) then
omega_u(i,k)=omegaout(i,k)
omega_u(i,k)=MAX(omega_u(i,k),-80.)
! Convert to m/s for use in convective time-scale:
rho = po(i,k)*100. / (rd * to(i,k))
tem = (-omega_u(i,k)) / ((rho * grav))
wu2(i,k) = tem**2
wu2(i,k) = max(wu2(i,k), 0.)
endif
enddo
enddo
else
! diagnostic method:
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz
tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k))
tem2 = wush(i,k) * sqrt(wu2(i,k-1))
tem2 = (tem1 - tem2) * dz
ptem = (1. - tem) * wu2(i,k-1)
ptem1 = 1. + tem
wu2(i,k) = (ptem + tem2) / ptem1
wu2(i,k) = max(wu2(i,k), 0.)
endif
endif
enddo
enddo
! convert to Pa/s for use in closure
do k = 1, km
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
Expand All @@ -1790,10 +1813,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
endif
enddo
enddo
endif

endif !progomega

!
! compute updraft velocity average over the whole cumulus
!
!> - Calculate the mean updraft velocity within the cloud (wc).
do i = 1, im
wc(i) = 0.
Expand Down Expand Up @@ -1822,10 +1846,12 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
if (wc(i) < val) cnvflg(i)=.false.
endif
enddo


c

!> - For progsigma = T, calculate the mean updraft velocity within the cloud (omegac),cast in pressure coordinates.
if(progsigma)then
if(progsigma)then
do i = 1, im
omegac(i) = 0.
sumx(i) = 0.
Expand Down Expand Up @@ -2912,10 +2938,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
advfac(i) = min(advfac(i), 1.)
endif
enddo

!> - 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
Expand Down
27 changes: 26 additions & 1 deletion physics/CONV/SAMF/samfdeepcnv.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
[ccpp-table-properties]
name = samfdeepcnv
type = scheme
dependencies = ../../tools/funcphys.f90,../../hooks/machine.F,samfaerosols.F,../progsigma_calc.f90
dependencies = ../../tools/funcphys.f90,../../hooks/machine.F,samfaerosols.F,../progsigma_calc.f90,../progomega_calc.f90

########################################################################
[ccpp-arg-table]
Expand Down Expand Up @@ -321,6 +321,13 @@
dimensions = ()
type = logical
intent = in
[progomega]
standard_name = do_prognostic_updraft_velocity
long_name = flag for prognostic omega in cumuls scheme
units = flag
dimensions = ()
type = logical
intent = in
[cldwrk]
standard_name = cloud_work_function
long_name = cloud work function
Expand Down Expand Up @@ -455,6 +462,24 @@
kind = kind_phys
intent = out
optional = True
[omegain]
standard_name = prognostic_updraft_velocity_in_convection
long_name = convective updraft velocity
units = frac
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
type = real
kind = kind_phys
intent = in
optional = True
[omegaout]
standard_name = updraft_velocity_updated_by_physics
long_name = convective updraft velocity updated by physics
units = frac
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
type = real
kind = kind_phys
intent = out
optional = True
[betascu]
standard_name = tuning_param_for_shallow_cu
long_name = tuning param for shallow cu in case prognostic closure is used
Expand Down
77 changes: 50 additions & 27 deletions physics/CONV/SAMF/samfshalcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ module samfshalcnv

use samfcnv_aerosols, only : samfshalcnv_aerosols
use progsigma, only : progsigma_calc

use progomega, only : progomega_calc
contains

subroutine samfshalcnv_init(imfshalcnv, imfshalcnv_samf, &
Expand Down Expand Up @@ -52,12 +52,13 @@ end subroutine samfshalcnv_init
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, &
& tmf,qmicro,progsigma,progomega, &
& 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, &
& sigmain,sigmaout,betadcu,betamcu,betascu,errmsg,errflg)
& sigmain,sigmaout,omegain,omegaout,betadcu,betamcu,betascu, &
& errmsg,errflg)
!
use machine , only : kind_phys
use funcphys , only : fpvs
Expand All @@ -75,7 +76,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
& tmf(:,:,:), q(:,:)
real(kind=kind_phys), intent(in), optional :: qmicro(:,:), &
& prevsq(:,:)
real(kind=kind_phys), intent(in), optional :: sigmain(:,:)
real(kind=kind_phys), intent(in), optional :: sigmain(:,:), &
& omegain(:,:)
!
real(kind=kind_phys), dimension(:), intent(in) :: fscav
integer, intent(inout) :: kcnv(:)
Expand All @@ -88,11 +90,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
& cnvw(:,:), cnvc(:,:), dt_mf(:,:)
!
real(kind=kind_phys), intent(out), optional :: ud_mf(:,:), &
& sigmaout(:,:)
& sigmaout(:,:), omegaout(:,:)

@grantfirl grantfirl Feb 12, 2025

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment about initializing intent(out) variables.

real(kind=kind_phys), intent(in) :: clam, c0s, c1, &
& asolfac, evef, pgcon
logical, intent(in) :: hwrf_samfshal,first_time_step, &
& restart,progsigma
& restart,progsigma,progomega
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
!
Expand Down Expand Up @@ -1478,7 +1480,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
!
! compute updraft velocity square(wu2)
!> - Calculate updraft velocity square(wu2) according to Han et al.'s (2017) \cite han_et_al_2017 equation 7.
!
!!> - if progomega = true, calculate prognostic updraft velocity (Pa/s) according to progomega routine.
if (hwrf_samfshal) then
do i = 1, im
if (cnvflg(i)) then
Expand All @@ -1493,26 +1495,46 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
endif
enddo
endif
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz
tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k))
tem2 = wush(i,k) * sqrt(wu2(i,k-1))
tem2 = (tem1 - tem2) * dz
ptem = (1. - tem) * wu2(i,k-1)
ptem1 = 1. + tem
wu2(i,k) = (ptem + tem2) / ptem1
wu2(i,k) = max(wu2(i,k), 0.)
endif
endif
enddo
enddo
!
if(progsigma)then
do k = 2, km1
if (progomega) then
call progomega_calc(first_time_step,restart,im,km,
& kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout,
& grav,buo,drag,wush,xlamue,bb1,bb2)
do k = 1, km
do i = 1, im
if (cnvflg(i)) then
omega_u(i,k)=omegaout(i,k)
omega_u(i,k)=MAX(omega_u(i,k),-80.)
! Convert to m/s for use in convective time-scale:
rho = po(i,k)*100. / (rd * to(i,k))
tem = (-omega_u(i,k)) / ((rho * grav))
wu2(i,k) = tem**2
wu2(i,k) = max(wu2(i,k), 0.)
endif
enddo
enddo

else
! diagnostic updraft velocity
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz
tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k))
tem2 = wush(i,k) * sqrt(wu2(i,k-1))
tem2 = (tem1 - tem2) * dz
ptem = (1. - tem) * wu2(i,k-1)
ptem1 = 1. + tem
wu2(i,k) = (ptem + tem2) / ptem1
wu2(i,k) = max(wu2(i,k), 0.)
endif
endif
enddo
enddo
!convert to Pa/s for use in closure
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
Expand All @@ -1523,8 +1545,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
endif
enddo
enddo
endif

endif !progomega

! compute updraft velocity averaged over the whole cumulus
!
!> - Calculate the mean updraft velocity within the cloud (wc).
Expand Down
Loading