Skip to content

Minor modification and bug fix in GFS cumulus convection schemes #232

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
grantfirl merged 22 commits into from
Mar 6, 2025
Merged

Minor modification and bug fix in GFS cumulus convection schemes #232

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
14c2f52
modified prognostic updraft fraction (sigmab) calcultation
rhaesung Oct 25, 2024
36e7b63
fix in missing vertical transport of turbulent kinetic energy (TKE) w…
rhaesung Oct 25, 2024
1ace72d
introduce TKE at model layer interfaces (tkeh)
rhaesung Oct 25, 2024
b325da7
fix compilation error
rhaesung Oct 31, 2024
f6f1892
replace the hardcoded precision constants with the appropriate precision
rhaesung Nov 12, 2024
a1aaf1a
Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update
rhaesung Nov 20, 2024
39e00e2
Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update
rhaesung Nov 26, 2024
80609a9
Fixing the incorrectly modified precision constant
rhaesung Dec 9, 2024
c241b99
Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update
rhaesung Jan 13, 2025
2ba9a65
Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update
rhaesung Jan 28, 2025
daeabfa
Merge remote-tracking branch 'upstream' into conv_hr5update
rhaesung Feb 4, 2025
d63e3f1
Merge remote-tracking branch 'upstream' into conv_hr5update
rhaesung Feb 4, 2025
485a7c7
Merge branch 'conv_hr5update' of https://github.com/rhaesung/ccpp-phy…
rhaesung Feb 5, 2025
a7f6a11
Fix bug
rhaesung Feb 10, 2025
a2a4f60
fix bug
rhaesung Feb 13, 2025
1bf0b07
make tkeh intent(inout)
rhaesung Feb 14, 2025
621a04d
fix error
rhaesung Feb 14, 2025
408d8e2
set tkeh=0 before doing a loop
rhaesung Feb 14, 2025
467451f
address reviewer comment
rhaesung Feb 14, 2025
b3b8cc6
address reviewer comment
rhaesung Feb 14, 2025
84a31db
Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update
rhaesung Feb 18, 2025
52331b2
Merge remote-tracking branch 'upstream/ufs/dev' into conv_hr5update
rhaesung Feb 21, 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
6 changes: 3 additions & 3 deletions physics/CONV/C3/cu_c3_deep.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2033,9 +2033,9 @@ subroutine cu_c3_deep_run( &
endif
enddo
call progsigma_calc(itf,ktf,flag_init,flag_restart,flag_shallow, &
flag_mid,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,xlv,dtime, &
forceqv_spechum,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, &
sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
flag_mid,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,xlv,dtime, &
forceqv_spechum,k22,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, &
sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
endif

!$acc end kernels
Expand Down
6 changes: 3 additions & 3 deletions physics/CONV/C3/cu_c3_sh.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -983,9 +983,9 @@ subroutine cu_c3_sh_run ( &
endif
enddo
call progsigma_calc(itf,ktf,flag_init,flag_restart,flag_shallow, &
flag_mid,del,tmf,qmicro,dbyo,zdqca,omega_u,zeta,xlv,dtime, &
forceqv_spechum,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, &
sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
flag_mid,del,tmf,qmicro,dbyo,zdqca,omega_u,zeta,xlv,dtime, &
forceqv_spechum,k22,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, &
sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)

endif

Expand Down
29 changes: 23 additions & 6 deletions physics/CONV/SAMF/samfdeepcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
& tmf,qmicro,itc,ntc,cliq,cp,cvap, &
& eps,epsm1,fv,grav,hvap,rd,rv, &
& t0c,delt,ntk,ntr,delp, &
& prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, &
& prslp,psp,phil,tkeh,qtr,prevsq,q,q1,t1,u1,v1,fscav, &
& hwrf_samfdeep,progsigma,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, &
Expand All @@ -97,7 +97,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
& fv, grav, hvap, rd, rv, t0c
real(kind=kind_phys), intent(in) :: delt
real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), &
& prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:)
& 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,sigmab_coldstart
Expand All @@ -115,7 +115,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH
real(kind=kind_phys), intent(inout) :: qtr(:,:,:), &
& q1(:,:), t1(:,:), u1(:,:), v1(:,:), &
& cnvw(:,:), cnvc(:,:)
& cnvw(:,:), cnvc(:,:), tkeh(:,:)

integer, intent(out) :: kbot(:), ktop(:)
real(kind=kind_phys), intent(out) :: cldwrk(:), &
Expand Down Expand Up @@ -953,8 +953,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
if(cnvflg(i)) then
if(k >= kb(i) .and. k < kbcon(i)) then
dz = zo(i,k+1) - zo(i,k)
tem = 0.5 * (qtr(i,k,ntk)+qtr(i,k+1,ntk))
tkemean(i) = tkemean(i) + tem * dz
tkemean(i) = tkemean(i) + tkeh(i,k) * dz
sumx(i) = sumx(i) + dz
endif
endif
Expand Down Expand Up @@ -1286,6 +1285,24 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
enddo
enddo
enddo
if(ntk > 2) then
kk = ntk -2
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kb(i) .and. k < kmax(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
tem = cq * tem
factor = 1. + tem
ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem *
& (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor
ercko(i,k,kk) = ecko(i,k,kk)
endif
endif
enddo
enddo
endif
endif
endif
c
Expand Down Expand Up @@ -2942,7 +2959,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
flag_mid = .false.
call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
& flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,
& delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu,
& delt,qadv,kb,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu,
& sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
endif

Expand Down
8 changes: 8 additions & 0 deletions physics/CONV/SAMF/samfdeepcnv.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,6 +242,14 @@
type = real
kind = kind_phys
intent = in
[tkeh]
standard_name = vertical_turbulent_kinetic_energy_at_interface
long_name = vertical turbulent kinetic energy at model layer interfaces
units = m2 s-2
dimensions = (horizontal_loop_extent,vertical_interface_dimension)
type = real
kind = kind_phys
intent = inout
[qtr]
standard_name = convective_transportable_tracers
long_name = array to contain cloud water and other convective trans. tracers
Expand Down
29 changes: 23 additions & 6 deletions physics/CONV/SAMF/samfshalcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@ 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, &
& prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, &
& prslp,psp,phil,tkeh,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, &
Expand Down Expand Up @@ -82,7 +82,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
integer, intent(inout) :: kcnv(:)
! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH
real(kind=kind_phys), intent(inout) :: qtr(:,:,:), &
& q1(:,:), t1(:,:), u1(:,:), v1(:,:)
& q1(:,:), t1(:,:), u1(:,:), v1(:,:), tkeh(:,:)
!
integer, intent(out) :: kbot(:), ktop(:)
real(kind=kind_phys), intent(out) :: rn(:), &
Expand Down Expand Up @@ -873,14 +873,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
tkemean(i) = 0.
endif
enddo

!
do k = 1, km1
do i = 1, im
if(cnvflg(i)) then
if(k >= kb(i) .and. k < kbcon(i)) then
dz = zo(i,k+1) - zo(i,k)
tem = 0.5 * (qtr(i,k,ntk)+qtr(i,k+1,ntk))
tkemean(i) = tkemean(i) + tem * dz
tkemean(i) = tkemean(i) + tkeh(i,k) * dz
sumx(i) = sumx(i) + dz
endif
endif
Expand Down Expand Up @@ -1094,6 +1093,24 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
enddo
enddo
enddo
if(ntk > 2) then
kk = ntk -2
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kb(i) .and. k < kmax(i)) then
dz = zi(i,k) - zi(i,k-1)
tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz
tem = cq * tem
factor = 1. + tem
ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem *
& (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor
ercko(i,k,kk) = ecko(i,k,kk)
endif
endif
enddo
enddo
endif
endif
endif
c
Expand Down Expand Up @@ -1983,7 +2000,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
flag_mid = .false.
call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
& flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,
& delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu,
& delt,qadv,kb,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu,
& sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
endif

Expand Down
8 changes: 8 additions & 0 deletions physics/CONV/SAMF/samfshalcnv.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,6 +242,14 @@
type = real
kind = kind_phys
intent = in
[tkeh]
standard_name = vertical_turbulent_kinetic_energy_at_interface
long_name = vertical turbulent kinetic energy at model layer interfaces
units = m2 s-2
dimensions = (horizontal_loop_extent,vertical_interface_dimension)
type = real
kind = kind_phys
intent = inout
[qtr]
standard_name = convective_transportable_tracers
long_name = array to contain cloud water and other convective trans. tracers
Expand Down
112 changes: 60 additions & 52 deletions physics/CONV/progsigma_calc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ module progsigma
!!\section gen_progsigma progsigma_calc General Algorithm
subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,&
flag_mid,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, &
delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu, &
delt,qadv,kb,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu, &
sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
!
!
Expand All @@ -29,7 +29,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,&
implicit none

! intent in
integer, intent(in) :: im,km,kbcon1(im),ktcon(im)
integer, intent(in) :: im,km,kb(im),kbcon1(im),ktcon(im)
real(kind=kind_phys), intent(in) :: hvap,delt,betascu,betamcu,betadcu, &
sigmind,sigminm,sigmins
real(kind=kind_phys), intent(in) :: qadv(im,km),del(im,km), &
Expand All @@ -46,9 +46,9 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,&
! 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), &
real(kind=kind_phys) :: fdqa(im),form(im,km), &
dp(im,km),inbu(im,km)
real(kind=kind_phys) :: sumx(im)

real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, &
fdqb,dtdyn,dxlim,rmulacvg,tem, &
Expand Down Expand Up @@ -84,7 +84,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,&
termC(i)=0.
termD(i)=0.
fdqa(i)=0.
mcons(i)=0.
sumx(i)=0.
enddo

do k = 2,km1
Expand All @@ -93,47 +93,49 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,&
dp(i,k) = 1000. * del(i,k)
endif
enddo
enddo
enddo

!Initial computations, place maximum sigmain in sigmab
do i=1,im
!compute sigmain averaged over cloud layers after advection and place it in sigmab
do k=2,km1
do i=1,im
if(cnvflg(i))then
do k=2,km
if(sigmain(i,k)>sigmab(i))then
sigmab(i)=sigmain(i,k)
endif
enddo
endif
enddo

do i=1,im
if(cnvflg(i))then
if(sigmab(i) < 1.E-5)then !after advection
sigmab(i)=0.
if(k > kbcon1(i) .and. k < ktcon(i)) then
sigmab(i) = sigmab(i) + sigmain(i,k) * dp(i,k)
sumx(i) = sumx(i) + dp(i,k)
endif
endif
endif
enddo
enddo
do i = 1, im
if(cnvflg(i)) then
if(sumx(i) == 0.) then
k = kbcon1(i)
sigmab(i) = sigmain(i,k)
else
sigmab(i) = sigmab(i) / sumx(i)
sigmab(i) = min(sigmab(i), 1._kind_phys)
if(sigmab(i) < 1.E-5) sigmab(i)=0.
endif
endif
enddo

!compute termD "The vertical integral of the latent heat convergence is limited to the
!buoyant layers with positive moisture convergence (accumulated from the surface).
!Lowest level:
do i = 1,im
dp1 = 1000. * del(i,1)
mcons(i)=(hvap*(qadv(i,1)+tmf(i,1)+qmicro(i,1))*dp1)
enddo
!Levels above:
do k = 2,km1
! layers with positive moisture convergence (accumulated from the updraft starting level).
do k = 1,km1
do i = 1,im
if(cnvflg(i))then
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(dbyo1(i,k)>0 .and. buy2 > 0.)then
if(k >= kb(i) .and. k < ktcon(i)) then
mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i,k))
buy2 = termD(i)+mcon
!
! Do the integral over buoyant layers with positive mcon acc from
! updraft starting level
!
if(buy2 > 0.)then
inbu(i,k)=1.
endif
inbu(i,k-1)=MAX(inbu(i,k-1),inbu(i,k))
termD(i) = termD(i) + inbu(i,k-1)*mcons(i)
mcons(i)=mcon
termD(i) = termD(i) + mcon
endif
endif
endif
enddo
enddo
Expand All @@ -142,8 +144,10 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,&
do k = 2,km1
do i = 1,im
if(cnvflg(i))then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
tem=sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k)*dp(i,k)
termA(i)=termA(i)+tem
endif
endif
enddo
enddo
Expand All @@ -152,8 +156,10 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,&
do k = 2,km1
do i = 1,im
if(cnvflg(i))then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
tem=zeta(i,k)*dbyo1(i,k)*inbu(i,k)*dp(i,k)
termB(i)=termB(i)+tem
endif
endif
enddo
enddo
Expand All @@ -162,31 +168,33 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,&
do k = 2,km1
do i = 1,im
if(cnvflg(i))then
if(k >= kbcon1(i) .and. k < ktcon(i)) then
form(i,k)=-1.0*inbu(i,k)*(omega_u(i,k)*delt)
fdqb=0.5*((form(i,k)*zdqca(i,k)))
termC(i)=termC(i)+inbu(i,k)* &
(fdqb+fdqa(i))*hvap*zeta(i,k)
fdqa(i)=fdqb
endif
endif
enddo
enddo

!sigmab
do i = 1,im
if(cnvflg(i))then
DEN=MIN(termC(i)+termB(i),1.E8)
cvg=termD(i)*delt
ZZ=MAX(0.0,SIGN(1.0,termA(i))) &
*MAX(0.0,SIGN(1.0,termB(i))) &
*MAX(0.0,SIGN(1.0,termC(i)-epsilon))
cvg=MAX(0.0,cvg)
sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ))
if(sigmab(i)>0.)then
sigmab(i)=MIN(sigmab(i),0.95)
sigmab(i)=MAX(sigmab(i),sigmind_new)
endif
endif!cnvflg
enddo
do i = 1,im
if(cnvflg(i))then
DEN=MIN(termC(i)+termB(i),1.e8_kind_phys)
cvg=termD(i)*delt
ZZ=MAX(0.0,SIGN(1.0,termA(i))) &
*MAX(0.0,SIGN(1.0,termB(i))) &
*MAX(0.0,SIGN(1.0,termC(i)-epsilon))
cvg=MAX(0.0,cvg)
sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ))
if(sigmab(i)>0.)then
sigmab(i)=MIN(sigmab(i),0.95)
sigmab(i)=MAX(sigmab(i),sigmind_new)
endif
endif!cnvflg
enddo

do k=1,km
do i=1,im
Comment thread
rhaesung marked this conversation as resolved.
Expand Down
Loading