Skip to content
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

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

Open
wants to merge 7 commits into
base: ufs/dev
Choose a base branch
from
Open

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 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
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
124 changes: 73 additions & 51 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(:,:), &
Copy link
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 @@ -1821,11 +1845,10 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
val = 1.e-4
if (wc(i) < val) cnvflg(i)=.false.
endif
enddo
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 +2935,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(:,:)
Copy link
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