Skip to content

Add namelist options for Tiedtke CP and introduce wind and PBLH dependent mixing-length limits in GFS PBL #359

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

Draft
BinLiu-NOAA wants to merge 6 commits into
base: ufs/dev
Choose a base branch
from
Draft

Add namelist options for Tiedtke CP and introduce wind and PBLH dependent mixing-length limits in GFS PBL #359

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
cdabdb6
From @AndrewHazelton: Add a namelist option to control scale-awarenes…
BinLiu-NOAA Jan 16, 2026
9f13565
Modifies GFS_typedefs.F90, GFS_typedefs.meta, cu_ntiedtke.F90, and cu…
JunghoonShin-NOAA Jan 22, 2026
cb19d18
Introduced a wind-shear and PBL-height dependent mixing length limit
Jan 27, 2026
8f8c7a2
Modified the Tiedtke CP code to use a different the entrainment equation
JunghoonShin-NOAA Jul 29, 2025
e2ff6f5
Sync the format cleanups from feature/hafsv2p2_final branch for
BinLiu-NOAA Feb 23, 2026
41cbca7
Improved cu_ntiedtke.F90 further based on the reviewer’s comments and…
JunghoonShin-NOAA Feb 26, 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
69 changes: 47 additions & 22 deletions physics/CONV/nTiedtke/cu_ntiedtke.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,8 +169,8 @@ end subroutine cu_ntiedtke_init
!=================================================================================================================
! level 1 subroutine 'cu_ntiedkte_run'
subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, &
evap,hfx,zprecc,lmask,lq,km,dt,dx,kbot,ktop,kcnv, &
ktrac,ud_mf,dd_mf,dt_mf,cnvw,cnvc,errmsg,errflg)
evap,hfx,zprecc,lmask,scale_fac_opt,lq,km,dt,dx,kbot,ktop,kcnv, &
ktrac,ud_mf,dd_mf,dt_mf,cnvw,cnvc,icu_zoentr,errmsg,errflg)
!=================================================================================================================
! this is the interface between the model and the mass flux convection module
! m.tiedtke e.c.m.w.f. 1989
Expand Down Expand Up @@ -200,6 +200,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,
!
implicit none
!--- input arguments:
integer, intent(in) :: scale_fac_opt,icu_zoentr

@climbfuji climbfuji Feb 23, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I think it might be useful to add links to the two issues that describe the need for these parameters and what they do as comments above line 203.

@JunghoonShin-NOAA JunghoonShin-NOAA Feb 23, 2026 •
edited
Loading

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

The option scale_fac_opt (alternative scale-awareness method) is related to the issue #357 from AndrewHazelton while the option icu_zoentr is related to the issue #358. I can add these two "links" as comments in the code near line 203.
#357 (scale_fac_opt)
#358 (icu_zoentr)

integer, intent(in) :: lq, km, ktrac
integer, intent(in), dimension(:) :: lmask

Expand Down Expand Up @@ -256,13 +257,23 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,
!
dxref = 15000.
do j=1,lq
if (dx(j).lt.dxref) then
scale_fac(j) = (1.06133+log(dxref/dx(j)))**3
scale_fac2(j) = scale_fac(j)**0.5
else
scale_fac(j) = 1.+1.33e-5*dx(j)
scale_fac2(j) = 1.
end if
if (scale_fac_opt == 1) then
if (dx(j).lt.dxref) then
scale_fac(j) = (1.06133+log(dxref/dx(j)))**2
scale_fac2(j) = scale_fac(j)
else
scale_fac(j) = 1.+1.33e-5*dx(j)
scale_fac2(j) = 1.
end if
else
if (dx(j).lt.dxref) then
scale_fac(j) = (1.06133+log(dxref/dx(j)))**3
scale_fac2(j) = scale_fac(j)**0.5
else
scale_fac(j) = 1.+1.33e-5*dx(j)
scale_fac2(j) = 1.
end if
end if
end do
!
! masv flux diagnostics.
Expand Down Expand Up @@ -370,7 +381,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,
& ktype, icbot, ictop, ztu, zqu, &
& zlu, zlude, zmfu, zmfd, zrain, &
& pcte, phhfl, lndj, pgeoh, zmfude_rate, dx, &
& scale_fac, scale_fac2)
& scale_fac, scale_fac2, icu_zoentr)
!
! to include the cloud water and cloud ice detrained from convection
!
Expand Down Expand Up @@ -446,7 +457,7 @@ subroutine cumastrn &
& ktype, kcbot, kctop, ptu, pqu, &
& plu, plude, pmfu, pmfd, prain, &
& pcte, phhfl, lndj, zgeoh, pmfude_rate, dx, &
& scale_fac, scale_fac2)
& scale_fac, scale_fac2, icu_zoentr)
implicit none
!
!***cumastrn* master routine for cumulus massflux-scheme
Expand Down Expand Up @@ -509,7 +520,7 @@ subroutine cumastrn &

!--- input arguments:
integer,intent(in):: klev,klon,klevp1,klevm1
integer,intent(in):: ktrac
integer,intent(in):: ktrac,icu_zoentr
integer,intent(in),dimension(klon):: lndj

real(kind=kind_phys),intent(in):: ztmst
Expand Down Expand Up @@ -689,7 +700,7 @@ subroutine cumastrn &
& zmfus, zmfuq, zmful, plude, zdmfup, &
& kcbot, kctop, ictop0, icum, ztmst, &
& zqsenh, zlglac, lndj, wup, wbase, &
& kdpl, pmfude_rate)
& kdpl, pmfude_rate, icu_zoentr)

!* (b) check cloud depth and change entrainment rate accordingly
! calculate precipitation rate (for downdraft calculation)
Expand Down Expand Up @@ -2010,7 +2021,7 @@ subroutine cuascn &
& pmfus, pmfuq, pmful, plude, pdmfup, &
& kcbot, kctop, kctop0, kcum, ztmst, &
& pqsenh, plglac, lndj, wup, wbase, &
& kdpl, pmfude_rate)
& kdpl, pmfude_rate, icu_zoentr)

implicit none
! this routine does the calculations for cloud ascents
Expand Down Expand Up @@ -2082,9 +2093,10 @@ subroutine cuascn &
! kctop - cloud top level
! kctop0 [ictop0] - estimate of cloud top. (cumastr)
! kcum [icum] - flag to control the call
! icu_zoentr 1: original entrainment equation, 2: Bechtold 2008 equation

!--- input arguments:
integer,intent(in):: klev,klon,klevp1,klevm1
integer,intent(in):: klev,klon,klevp1,klevm1,icu_zoentr
integer,intent(in),dimension(klon):: lndj
integer,intent(in),dimension(klon):: klwmin
integer,intent(in),dimension(klon):: kdpl
Expand Down Expand Up @@ -2135,6 +2147,8 @@ subroutine cuascn &
real(kind=kind_phys),dimension(klon):: zph,zdmfen,zdmfde,zmfuu,zmfuv,zpbase,zqold,zluold,zprecip
real(kind=kind_phys),dimension(klon,klev):: zlrain,zbuo,kup,zodetr,pdmfen

real(kind=kind_phys),parameter:: c1 = 5.0e-4 !shin
real(kind=kind_phys),parameter:: d1 = 1.0e-3 !shin
!--------------------------------
!* 1. specify parameters
!--------------------------------
Expand Down Expand Up @@ -2271,9 +2285,14 @@ subroutine cuascn &
! Why is it negative?
!---------------------------------------
if ( jk == kcbot(jl) ) then

zoentr(jl) = -entorg*(min(1.,pqen(jl,jk)/pqsen(jl,jk)) - &
1.)*(pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
if ( icu_zoentr .eq. 1 ) then
zoentr(jl) = -entorg*(min(1.,pqen(jl,jk)/pqsen(jl,jk)) - &
1.)*(pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
end if
if ( icu_zoentr .eq. 2 ) then
zoentr(jl) = (c1+d1*(1.0-min(1.,pqen(jl,jk)/pqsen(jl,jk))))* &
(pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
end if

@climbfuji climbfuji Feb 23, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

There is no catchall guard here for icu_zoentry not in [1,2].

Same for line 2498 below.

@JunghoonShin-NOAA JunghoonShin-NOAA Feb 23, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Thank you for your comments. Regarding the catchall guard for the option icu_zoentr, I could modify the code so that the code can stop with warning message if the icu_zoentr is neither 1 nor 2 as shown in the below example. If you think it would be a good idea to modify the code as shown below, please let me know.

If( icu_zoentr .eq. 1) then
original entrainment method
elseif ( icu_zoentr .eq.2 ) then
new entrainment method
else
print('Error: unsupported icu_zoentry of', icu_zeontry)
stop
endif

@climbfuji climbfuji Feb 23, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

With CCPP, the correct way is to write the error message to the errmsg string, set errflg to a value other than zero (1 is fine), and return immediately. We shouldn't stop the model with stop.

@JunghoonShin-NOAA JunghoonShin-NOAA Feb 23, 2026 •
edited
Loading

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Ok.. I checked the cu_ntiedtke.F90 code to understand your instruction. Then, do you think the below example satisfies the CCPP standard? If you think this modification looks fine, shall I change to the code as shown below?

If( icu_zoentr .eq. 1) then
original entrainment method
elseif ( icu_zoentr .eq.2 ) then
new entrainment method
else
write(errmsg,'(*(a))') 'Error: unsupported icu_zoentr’
errflg = 1
return
endif

@climbfuji climbfuji Feb 23, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Yes, that looks correct to me. Thanks!

@JunghoonShin-NOAA JunghoonShin-NOAA Feb 24, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Thank you. Then, I will modify and will test the codes. I will push the changes and will come back to you after I complete the tests.

@JunghoonShin-NOAA JunghoonShin-NOAA Feb 25, 2026 •
edited
Loading

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I tested cu_ntiedtke.F90 after I changed the code as the example shown above to add the catchall guard for the option icu_zoentr. If I used icu_zoentr=2, which is a supported option, the code worked properly. I also tested the code with the unsupported option icu_zoentr=3, and the forecast job crashed (as expected) with the error message shown below. This error message will tell users that the code doesn't work except for options 1 and 2.

An error occurred in ccpp_physics_run for group phys_ts, block/chunk 1 and thread 1 (nt= 1):
An error occured in cu_ntiedtke_run: Error: unsupported icu_zoentr

Log file location: /scratch3/HFIP/hwrfv3/scrub/Junghoon.Shin/hafsv2p2_final_change_test/2024100512/14L/hafs_atm_init.log

If you think this looks fine, I will push this change to the feature/hafs_tiedtke_varml branch for the code review.
Thank you.

@JunghoonShin-NOAA JunghoonShin-NOAA Feb 26, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

I just pushed the changed cu_ntiedtke.F90 code to the feature/hafs_tiedtke_varml branch.
Please review the code and let me know if you have any further comments. Thank you.

zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk+1)
end if
!---------------------------------------
Expand Down Expand Up @@ -2467,10 +2486,16 @@ subroutine cuascn &
ikb = kcbot(jl)
! zoentr is overwritten, but not used until
! the next jk level in the loop (ICON comment)
zoentr(jl) = entorg*(0.3-(min(1.,pqen(jl,jk-1) / &
pqsen(jl,jk-1))-1.))*(pgeoh(jl,jk-1)-pgeoh(jl,jk)) * &
zrg*min(1.,pqsen(jl,jk)/pqsen(jl,ikb))**3

if ( icu_zoentr .eq. 1 ) then
zoentr(jl) = entorg*(0.3-(min(1.,pqen(jl,jk-1) / &
pqsen(jl,jk-1))-1.))*(pgeoh(jl,jk-1)-pgeoh(jl,jk)) * &
zrg*min(1.,pqsen(jl,jk)/pqsen(jl,ikb))**3
end if
if ( icu_zoentr .eq. 2 ) then
zoentr(jl) = ( c1*(min(1.,pqsen(jl,jk)/pqsen(jl,ikb))**2) + &
d1*(1.0-min(1.,pqen(jl,jk-1)/pqsen(jl,jk-1)))*(min(1.,pqsen(jl,jk)/pqsen(jl,ikb))**3))* &
(pgeoh(jl,jk-1)-pgeoh(jl,jk))*zrg
end if
zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk)
else
zoentr(jl) = 0.
Expand Down
14 changes: 14 additions & 0 deletions physics/CONV/nTiedtke/cu_ntiedtke.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -305,6 +305,20 @@
type = real
kind = kind_phys
intent = out
[scale_fac_opt]
standard_name = control_for_scale_awareness_options_in_Tiedtke
long_name = control for scale awareness options in the Tiedtke scheme
units = none
dimensions = ()
type = integer
intent = in
[icu_zoentr]
standard_name = option_for_different_entrainment_in_Tiedtke
long_name = option for different entrainment in Tiedtke
units = none
dimensions = ()
type = integer
intent = in
[cnvc]
standard_name = convective_cloud_cover
long_name = convective cloud cover
Expand Down
33 changes: 31 additions & 2 deletions physics/PBL/SATMEDMF/satmedmfvdifq.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -320,6 +320,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &

real(kind=kind_phys) qice(im,km),qliq(im,km)

real(kind=kind_phys) vlmax(im,km),wind10(im),vlm

!PCC_CANOPY------------------------------------
integer COUNTCAN,KCAN
integer kount !IVAI
Expand Down Expand Up @@ -751,6 +753,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
do i=1,im
ustar(i) = sqrt(stress(i))
wind10(i) = sqrt(u10m(i)**2+v10m(i)**2)
enddo
!
!> - Compute buoyancy \f$\frac{\partial \theta_v}{\partial z}\f$ (bf)
Expand All @@ -763,6 +766,16 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
dw2 = (u1(i,k)-u1(i,k+1))**2
& + (v1(i,k)-v1(i,k+1))**2
shr2(i,k) = max(dw2,dw2min)*rdz*rdz
vlmax(i,k)=rlmx
if(rlmx < 0) then
tem=sqrt(u1(i,k)**2+v1(i,k)**2)
tem1=sqrt(u1(i,k+1)**2+v1(i,k+1)**2)
tem=0.5*(tem+tem1) ! mean wind speed
vlmax(i,k)=abs(rlmx)*tem/sqrt(shr2(i,k))
tem = (wind10(i)/17.5)**0.1
tem = max(tem,0.5)
vlmax(i,k)=vlmax(i,k)*tem
endif
enddo
enddo
!
Expand Down Expand Up @@ -1398,17 +1411,33 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
! Following Rodier et. al (2017), environmental wind shear effect on
! mixing length was included.
!
vlm=rlmx
if (elmx < 0) then
vlm=abs(elmx)*min(hpbl(i),2000.0) ! below PBLH
vlm=max(vLm,10.0)
if (k > kpbl(i)) vlm=min(vlmax(i,k),1.0/rdzt(i,k)) ! above PBLH
endif
!
ptem2 = min(zlup,zldn)
rlam(i,k) = elmfac * ptem2
rlam(i,k) = max(rlam(i,k), tem1)
rlam(i,k) = min(rlam(i,k), rlmx)
if (rlmx > 0) then
rlam(i,k) = min(rlam(i,k), rlmx)
else
rlam(i,k) = min(rlam(i,k), vlm)
endif
!
ptem2 = sqrt(zlup*zldn)
ele(i,k) = elefac * ptem2
ele(i,k) = max(ele(i,k), tem1)
elmh(i,k)= elmhfac * ele(i,k)
ele(i,k) = min(ele(i,k), elmx)
if (elmx > 0) then
ele(i,k) = min(ele(i,k), elmx)
else
ele(i,k) = min(ele(i,k), vlm)
endif

elmh(i,k)= min(elmh(i,k), elmhmx)
!
enddo
Expand Down