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

Jump to bottom

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

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

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

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
7 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
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 @@ -2016,9 +2016,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
25 changes: 20 additions & 5 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(:,:), tkeh(:,:)
real(kind=kind_phys), dimension(:), intent(in) :: fscav
logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, &
& progsigma,do_mynnedmf
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,22 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
enddo
enddo
enddo
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
c
Expand Down Expand Up @@ -2941,7 +2956,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_layer_dimension)
type = real
kind = kind_phys
intent = in
[qtr]
standard_name = convective_transportable_tracers
long_name = array to contain cloud water and other convective trans. tracers
Expand Down
27 changes: 21 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 @@ -53,7 +53,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 All @@ -71,7 +71,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
& betamcu
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(:,:), tkeh(:,:), &
& tmf(:,:,:), q(:,:)
real(kind=kind_phys), intent(in), optional :: qmicro(:,:), &
& prevsq(:,:)
Expand Down Expand Up @@ -872,14 +872,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 @@ -1093,6 +1092,22 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
enddo
enddo
enddo
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
c
Expand Down Expand Up @@ -1982,7 +1997,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_layer_dimension)
type = real
kind = kind_phys
intent = in
[qtr]
standard_name = convective_transportable_tracers
long_name = array to contain cloud water and other convective trans. tracers
Expand Down
126 changes: 66 additions & 60 deletions physics/CONV/progsigma_calc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,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 @@ -31,7 +31,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 @@ -48,13 +48,13 @@ 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, &
DEN,dp1,invdelt
DEN,dp1,invdelt,sigmind_new

!Parameters
gcvalmx = 0.1
Expand All @@ -63,6 +63,12 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,&
km1=km-1
invdelt = 1./delt

if (flag_init) then
sigmind_new=0.0
else
sigmind_new=sigmind
end if

!Initialization 2D
do k = 1,km
do i = 1,im
Expand All @@ -80,7 +86,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 @@ -89,47 +95,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.0)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Change 1.0 to 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 @@ -138,8 +146,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 @@ -148,8 +158,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 @@ -158,39 +170,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
if(flag_init .and. .not. flag_restart)then
do i = 1,im
if(cnvflg(i))then
sigmab(i)=0.03
do i = 1,im
if(cnvflg(i))then
DEN=MIN(termC(i)+termB(i),1.E8)
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 as above for all the fortran intrinsics in this block. (Explicitly setting the precision for these calls removes any ambiguity and makes to code more portable)

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@dustinswales Thanks! I've updated the code to replace the hardcoded precision constants with the appropriate precision.

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
enddo
else
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),0.01)
endif
endif!cnvflg
enddo
endif
endif!cnvflg
enddo

do k=1,km
do i=1,im
Expand Down
Loading
Loading