Skip to content

Updates to progomega option of saSAS convection scheme #364

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 8 commits into from
Mar 25, 2026
Merged

Updates to progomega option of saSAS convection scheme #364

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
e673808
Updates to progomega option of saSAS convection scheme
lisa-bengtsson Mar 6, 2026
cd6a29d
Correction
lisa-bengtsson Mar 6, 2026
c2991ab
correction 2
lisa-bengtsson Mar 6, 2026
68e8b6d
Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physi…
lisa-bengtsson Mar 16, 2026
223d093
Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physi…
lisa-bengtsson Mar 18, 2026
e841285
Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physi…
lisa-bengtsson Mar 20, 2026
f2fbb17
Correct parameter setting in samfshalcnv to pass debug runs
lisa-bengtsson Mar 20, 2026
e5ebf60
address saturation vapor pressure out of bounds error in hafs_4_nest RT
lisa-bengtsson Mar 23, 2026
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
73 changes: 42 additions & 31 deletions physics/CONV/SAMF/samfdeepcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -220,10 +220,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
!
! parameters for prognostic sigma closure
real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
& omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km),
& tentr(im,km)
real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins
parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01)
& omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km)
real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins,
& wc_min, wc_eff
logical flag_shallow, flag_mid
c physical parameters
! parameter(grav=grav,asolfac=0.958)
Expand Down Expand Up @@ -347,6 +346,21 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
!************************************************************************
!
!
! - Initialize parameters related to prognostic closure
if (progsigma) then
if (progomega) then
sigmind = 0.03
sigmins = 0.03
sigminm = 0.03
wc_min = 0.2
else
sigmind = 0.01
sigmins = 0.03
sigminm = 0.03
wc_min = 0.2
endif
endif

km1 = km - 1
!> - Initialize column-integrated and other single-value-per-column variable arrays.
c
Expand Down Expand Up @@ -1137,7 +1151,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
do i=1,im
if(cnvflg(i) .and.
& (k > kbcon(i) .and. k < kmax(i))) then
tentr(i,k)=xlamue(i,k)*fent1(i,k)
tem = cxlamet(i) * frh(i,k) * fent2(i,k)
xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
tem1 = cxlamdt(i) * frh(i,k)
Expand Down Expand Up @@ -1788,11 +1801,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
if (progomega) then
call progomega_calc(first_time_step,restart,im,km,
& kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout,
& grav,buo,drag,wush,tentr,bb1,bb2)
& grav,buo,drag,wush,bb1,bb2)
do k = 1, km
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
omega_u(i,k)=omegaout(i,k)
omega_u(i,k)=MAX(omega_u(i,k),-80.)
! Convert to m/s for use in convective time-scale:
Expand All @@ -1809,7 +1822,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz
tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k))
Expand All @@ -1827,7 +1840,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
do k = 1, km
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
rho = po(i,k)*100. / (rd * to(i,k))
omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav
omega_u(i,k)=MAX(omega_u(i,k),-80.)
Expand Down Expand Up @@ -1878,7 +1891,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
dp = 1000. * del(i,k)
tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1))
omegac(i) = omegac(i) + tem * dp
Expand Down Expand Up @@ -2907,27 +2920,25 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
! compute convective turn-over time
!
!> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, the convective adjustment time (dtconv) is set to be proportional to the convective turnover time, which is computed using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx).
if(hwrf_samfdeep) then
do i= 1, im
if(cnvflg(i)) then
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
dtconv(i) = tem / wc(i)
dtconv(i) = max(dtconv(i),dtmin)
dtconv(i) = min(dtconv(i),dtmax)
endif
enddo
else
do i= 1, im
if(cnvflg(i)) then
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
dtconv(i) = tem / wc(i)
tfac = 1. + gdx(i) / 75000.
dtconv(i) = tfac * dtconv(i)
dtconv(i) = max(dtconv(i),dtmin)
dtconv(i) = min(dtconv(i),dtmax)
endif
enddo
endif
do i = 1, im
if (cnvflg(i)) then
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
if (progomega) then
wc_eff = max(wc(i), wc_min)
dtconv(i) = tem / wc_eff
else
dtconv(i) = tem / wc(i)
endif
!grid spacing scaling (disabled for HWRF SAMF deep)
if (.not. hwrf_samfdeep) then
tfac = 1. + gdx(i) / 75000.
dtconv(i) = tfac * dtconv(i)
endif
!bounds
dtconv(i) = max(dtconv(i), dtmin)
dtconv(i) = min(dtconv(i), dtmax)
endif
enddo
!
!> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed.
do i= 1, im
Expand Down
67 changes: 43 additions & 24 deletions physics/CONV/SAMF/samfshalcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
& omegac(im),zeta(im,km),dbyo1(im,km),
& sigmab(im),qadv(im,km)
real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins,
& sigminm
& sigminm,wc_min,wc_eff
logical flag_shallow,flag_mid
c physical parameters
! parameter(g=grav,asolfac=0.89)
Expand Down Expand Up @@ -205,8 +205,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
! parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
parameter(betaw=.03,dxcrtc0=9.e3)
parameter(h1=0.33333333)
! progsigma
parameter(dxcrtas=500.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01)
c local variables and arrays
real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km),
& uo(im,km), vo(im,km), qeso(im,km),
Expand Down Expand Up @@ -295,7 +293,21 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
prsl = prslp * 0.001
del = delp * 0.001
!************************************************************************
!
! - Initialize parameters related to prognostic closure
if (progsigma) then
if (progomega) then
sigmind = 0.03
sigmins = 0.03
sigminm = 0.03
wc_min = 0.2
else
sigmind = 0.01
sigmins = 0.03
sigminm = 0.03
wc_min = 0.2
endif
endif
!
km1 = km - 1
c
c initialize arrays
Expand Down Expand Up @@ -1520,11 +1532,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
if (progomega) then
call progomega_calc(first_time_step,restart,im,km,
& kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout,
& grav,buo,drag,wush,xlamue,bb1,bb2)
& grav,buo,drag,wush,bb1,bb2)
do k = 1, km
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
omega_u(i,k)=omegaout(i,k)
omega_u(i,k)=MAX(omega_u(i,k),-80.)
! Convert to m/s for use in convective time-scale:
Expand All @@ -1542,7 +1554,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz
tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k))
Expand All @@ -1560,7 +1572,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
rho = po(i,k)*100. / (rd * to(i,k))
omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav
omega_u(i,k)=MAX(omega_u(i,k),-80.)
Expand Down Expand Up @@ -1611,7 +1623,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
dp = 1000. * del(i,k)
tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1))
omegac(i) = omegac(i) + tem * dp
Expand All @@ -1636,7 +1648,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon1(i) .and. k < ktcon(i)) then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
if(omega_u(i,k) .ne. 0.)then
zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k))
else
Expand Down Expand Up @@ -1955,21 +1967,28 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
! compute convective turn-over time
!
!> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, calculate the convective turnover time using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx).
do i= 1, im
if(cnvflg(i)) then
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
dtconv(i) = tem / wc(i)
if (.not.hwrf_samfshal) then
tfac = 1. + gdx(i) / 75000.
dtconv(i) = tfac * dtconv(i)
endif
dtconv(i) = max(dtconv(i),dtmin)
dtconv(i) = max(dtconv(i),dt2)
dtconv(i) = min(dtconv(i),dtmax)
endif
do i = 1, im
if (cnvflg(i)) then
tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i))
if (progomega) then
wc_eff = max(wc(i), wc_min)
dtconv(i) = tem / wc_eff
else
dtconv(i) = tem / wc(i)
endif
! - grid spacing scaling (disabled for HWRF shallow option)
if (.not. hwrf_samfshal) then
tfac = 1. + gdx(i) / 75000.
dtconv(i) = tfac * dtconv(i)
endif
! - limits
dtconv(i) = max(dtconv(i), dtmin)
dtconv(i) = max(dtconv(i),dt2)
dtconv(i) = min(dtconv(i), dtmax)
endif
enddo
!
!> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed.
!
! > - Calculate advective time scale (tauadv) using a mean cloud layer wind speed.
do i= 1, im
if(cnvflg(i)) then
sumx(i) = 0.
Expand Down
34 changes: 19 additions & 15 deletions physics/CONV/progomega_calc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ module progomega
!!\section gen_progomega progomega_calc General Algorithm

subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegain,delt,del, &
zi,cnvflg,omegaout,grav,buo,drag,wush,tentr,bb1,bb2)
zi,cnvflg,omegaout,grav,buo,drag,wush,bb1,bb2)

use machine, only : kind_phys
use funcphys, only : fpvs
Expand All @@ -30,20 +30,17 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
integer, intent(in) :: kbcon1(im),ktcon(im)
real(kind=kind_phys), intent(in) :: delt,grav,bb1,bb2
real(kind=kind_phys), intent(in) :: omegain(im,km), del(im,km),zi(im,km)
real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km),tentr(im,km)
real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km)
real(kind=kind_phys), intent(inout) :: omegaout(im,km)
logical, intent(in) :: cnvflg(im),first_time_step,flag_restart
real(kind=kind_phys) :: termA(im,km),termB(im,km),termC(im,km),omega(im,km)
real(kind=kind_phys) :: RHS(im,km),Kd(im,km)
real(kind=kind_phys) :: dp,dz,entrn,Kdn,discr,wush_pa,lbb1,lbb2,lbb3
real(kind=kind_phys) :: dp,dz,discr,wush_pa,lbb1,lbb2,lbb3
integer :: i,k

entrn = 0.8E-4 !0.5E-4 !m^-1
Kdn = 0.5E-4 !2.9E-4 !m^-1
lbb1 = 0.5 !1.0
lbb2 = 3.2 !3.0
lbb3 = 0.5 !0.5

lbb1 = 1.5
lbb2 = 0.6
lbb3 = 1.2

!Initialization 2D
do k = 1,km
Expand All @@ -56,6 +53,16 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
enddo
enddo

do k = 1,km
do i = 1,im
if(cnvflg(i))then
if(omega(i,k) < 1.0E-5) then
omega(i,k) = 0.
endif
endif
enddo
enddo

if(first_time_step .and. .not. flag_restart)then
do k = 1,km
do i = 1,im
Expand All @@ -76,10 +83,7 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
do k = 2, km
do i = 1, im
if (cnvflg(i)) then
if (k > kbcon1(i) .and. k < ktcon(i)) then

! Aerodynamic drag parameter
Kd(i,k) = (tentr(i,k)/entrn)*Kdn
if (k >= kbcon1(i) .and. k < ktcon(i)) then

! Scale by dp/dz to have equation in Pa/s
!(dp/dz > 0)
Expand All @@ -91,8 +95,8 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai
!termC - Adds buoyancy forcing

!Coefficients for the quadratic equation
termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz)) + (Kd(i,k) * (dp/dz)))
termB(i,k) = -1.0 - delt * lbb3 * wush(i,k) * dp/dz
termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz)))
termB(i,k) = 1.0 - delt * lbb3 * wush(i,k) * dp/dz
termC(i,k) = omega(i,k) - delt * lbb2 * buo(i,k) * (dp/dz) &
- delt * omega(i,k) * (omega(i,k-1) - omega(i,k)) / dp
!Compute the discriminant
Expand Down
Loading