Skip to content

Commit

Permalink
Update progsigma physics (see PR 18 ufs-community#18) for ccpp-physics
Browse files Browse the repository at this point in the history 
	modified:   progsigma_calc.f90
	modified:   samfdeepcnv.f
	modified:   samfshalcnv.f
  • Loading branch information
@bluefinweiwei
bluefinweiwei committed Feb 22, 2023
1 parent 8bcf91d commit 3fe7b64
Show file tree
Hide file tree
Showing 3 changed files with 71 additions and 55 deletions.
25 changes: 8 additions & 17 deletions physics/progsigma_calc.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,17 +2,15 @@
!! This file contains the subroutine that calculates the prognostic
!! updraft area fraction that is used for closure computations in
!! saSAS deep and shallow convection, based on a moisture budget
!! as described in Bengtsson et al. 2022.
!! as described in Bengtsson et al. 2022 \cite Bengtsson_2022.

!>\ingroup samfdeepcnv
!! This subroutine computes a prognostic updraft area fraction
!>\ingroup SAMFdeep
!>\ingroup SAMF_shal
!> This subroutine computes a prognostic updraft area fraction
!! used in the closure computations in the samfdeepcnv.f scheme
!>\ingroup samfshalcnv
!! This subroutine computes a prognostic updraft area fracftion
!! used in the closure computations in the samfshalcnv. scheme
!!\section progsigma General Algorithm
!> @{

!!\section gen_progsigma progsigma_calc General Algorithm
subroutine progsigma_calc (im,km,flag_init,flag_restart, &
flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, &
delt,prevsq,q,kbcon1,ktcon,cnvflg,sigmain,sigmaout, &
Expand Down Expand Up @@ -55,8 +53,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, &
rmulacvg=10.
epsilon=1.E-11
km1=km-1
betadcu = 1.5
betascu = 5.2
betadcu = 2.0
betascu = 8.0
invdelt = 1./delt

!Initialization 2D
Expand Down Expand Up @@ -212,12 +210,11 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, &
enddo

!Reduce area fraction before coupling back to mass-flux computation.
!This tuning could be addressed in updraft velocity equation instead.
if(flag_shallow)then
do i= 1, im
if(cnvflg(i)) then
sigmab(i)=sigmab(i)/betascu
sigmab(i)=MAX(0.01,sigmab(i))
sigmab(i)=MAX(0.03,sigmab(i))
endif
enddo
else
Expand All @@ -229,10 +226,4 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, &
enddo
endif


end subroutine progsigma_calc
!> @}
!! @}



50 changes: 32 additions & 18 deletions physics/samfdeepcnv.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,7 +238,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
! parameter(cinacrmx=-120.,cinacrmn=-120.)
parameter(cinacrmx=-120.,cinacrmn=-80.)
parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
parameter(betaw=.03,dxcrtuf=15.e3)
parameter(betaw=.03)

!
! local variables and arrays
Expand Down Expand Up @@ -307,7 +307,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &

gravinv = 1./grav


elocp = hvap/cp
el2orc = hvap*hvap/(rv*cp)

Expand Down Expand Up @@ -1620,7 +1619,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
c
do i = 1, im
flg(i) = cnvflg(i)
ktcon1(i) = kmax(i)
ktcon1(i) = ktcon(i)
enddo
do k = 2, km1
do i = 1, im
Expand All @@ -1639,8 +1638,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
! aa2(i) = aa2(i) +
!! & dz1 * eta(i,k) * grav * fv *
! & dz1 * grav * fv *
! & max(val,(qeso(i,k) - qo(i,k)))
if(aa2(i) < 0.) then
! & max(val,(qeso(i,k) - qo(i,k)))
!NRL MNM: Limit overshooting not to be deeper than half the actual cloud
tem = 0.5 * (zi(i,ktcon(i))-zi(i,kbcon(i)))
tem1 = zi(i,k)-zi(i,ktcon(i))
if(aa2(i) < 0. .or. tem1 >= tem) then
ktcon1(i) = k
flg(i) = .false.
endif
Expand Down Expand Up @@ -1731,8 +1733,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
enddo
enddo


if(progsigma)then
if(progsigma)then
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
Expand Down Expand Up @@ -1777,8 +1778,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
endif
enddo
c
!> - Calculate the mean updraft velocity within the cloud (wc),cast in pressure coordinates.
if(progsigma)then

!> - For progsigma = T, calculate the mean updraft velocity within the cloud (omegac),cast in pressure coordinates.
if(progsigma)then
do i = 1, im
omegac(i) = 0.
sumx(i) = 0.
Expand Down Expand Up @@ -1807,7 +1809,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
endif
enddo

!> - Calculate the xi term in Bengtsson et al. 2022 (equation 8)
!> - For progsigma = T, calculate the xi term in Bengtsson et al. 2022 \cite Bengtsson_2022 (equation 8)
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
Expand All @@ -1823,7 +1825,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
endif
enddo
enddo


endif !if progsigma

Expand Down Expand Up @@ -2466,8 +2468,10 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
!
if(progsigma)then
dxcrtas=30.e3
dxcrtuf=10.e3
else
dxcrtas=8.e3
dxcrtuf=15.e3
endif


Expand All @@ -2477,6 +2481,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
asqecflg(i) = .false.
endif
enddo

!
!> - If grid size is larger than the threshold value (i.e., asqecflg=.true.), the quasi-equilibrium assumption is used to obtain the cloud base mass flux. To begin with, calculate the change in the temperature and moisture profiles per unit cloud base mass flux.
do k = 1, km
Expand Down Expand Up @@ -2878,7 +2883,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
endif
enddo

!> - From Bengtsson et al. (2022) Prognostic closure scheme, equation 8, compute updraft area fraction based on a moisture budget
!> - 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
flag_shallow = .false.
call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
Expand All @@ -2889,6 +2894,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &

!> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer.
!! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv.
do i= 1, im
if(cnvflg(i) .and. .not.asqecflg(i)) then
k = kbcon(i)
Expand Down Expand Up @@ -3515,9 +3521,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
if(k > kb(i) .and. k < ktop(i)) then
tem = 0.5 * (eta(i,k-1) + eta(i,k)) * xmb(i)
tem1 = pfld(i,k) * 100. / (rd * t1(i,k))
sigmagfm(i) = max(sigmagfm(i), betaw)
ptem = tem / (sigmagfm(i) * tem1)
qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*sigmagfm(i)*ptem*ptem
if(progsigma)then
tem2 = sigmab(i)
else
tem2 = max(sigmagfm(i), betaw)
endif
ptem = tem / (tem2 * tem1)
qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*tem2*ptem*ptem
endif
endif
enddo
Expand All @@ -3529,9 +3539,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
if(k > 1 .and. k <= jmin(i)) then
tem = 0.5*edto(i)*(etad(i,k-1)+etad(i,k))*xmb(i)
tem1 = pfld(i,k) * 100. / (rd * t1(i,k))
sigmagfm(i) = max(sigmagfm(i), betaw)
ptem = tem / (sigmagfm(i) * tem1)
qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*sigmagfm(i)*ptem*ptem
if(progsigma)then
tem2 = sigmab(i)
else
tem2 = max(sigmagfm(i), betaw)
endif
ptem = tem / (tem2 * tem1)
qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*tem2*ptem*ptem
endif
endif
enddo
Expand Down
51 changes: 31 additions & 20 deletions physics/samfshalcnv.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,13 +50,14 @@ end subroutine samfshalcnv_init
!! \section det_samfshalcnv GFS samfshalcnv Detailed Algorithm
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, &
& 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,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,errmsg,errflg)!
& clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, &
& sigmain,sigmaout,errmsg,errflg)
!
use machine , only : kind_phys
use funcphys , only : fpvs

Expand All @@ -68,11 +69,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
& eps, epsm1, 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(:,:), &
& qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:)
& prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), &
& qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:)

real(kind=kind_phys), intent(in) :: sigmain(:,:)

!
real(kind=kind_phys), dimension(:), intent(in) :: fscav
integer, intent(inout) :: kcnv(:)
Expand Down Expand Up @@ -160,9 +160,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
& omegac(im),zeta(im,km),dbyo1(im,km),
& sigmab(im)
real(kind=kind_phys) gravinv,dxcrtas

logical flag_shallow

c physical parameters
! parameter(g=grav,asolfac=0.89)
! parameter(g=grav)
Expand Down Expand Up @@ -190,7 +188,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
parameter(cinacrmx=-120.,shevf=2.0)
parameter(dtmax=10800.,dtmin=600.)
parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5)
parameter(betaw=.03,dxcrt=15.e3,dxcrtc0=9.e3)
parameter(betaw=.03,dxcrtc0=9.e3)
parameter(h1=0.33333333)
! progsigma
parameter(dxcrtas=30.e3)
Expand Down Expand Up @@ -260,6 +258,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
cinacrmn=-80.
endif

if (progsigma) then
dxcrt=10.e3
else
dxcrt=15.e3
endif

c-----------------------------------------------------------------------
if (.not.hwrf_samfshal) then
!> ## Determine whether to perform aerosol transport
Expand Down Expand Up @@ -338,6 +342,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
enddo
endif
!!

!> - Return to the calling routine if deep convection is present or the surface buoyancy flux is negative.
totflg = .true.
do i=1,im
Expand Down Expand Up @@ -514,6 +519,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
do i = 1,im
omegac(i)=0.
enddo
do k = 1, km
do i = 1, im
dbyo1(i,k)=0.
Expand Down Expand Up @@ -1472,7 +1478,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
enddo
enddo
!
if(progsigma)then
if(progsigma)then
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
Expand Down Expand Up @@ -1517,8 +1523,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
endif
enddo
c
!> - Calculate the mean updraft velocity in pressure coordinates within the cloud (wc).
if(progsigma)then
!> - For progsigma =T, calculate the mean updraft velocity in pressure coordinates within the cloud (wc).
if(progsigma)then
do i = 1, im
omegac(i) = 0.
sumx(i) = 0.
Expand Down Expand Up @@ -1546,8 +1552,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
if (omegac(i) > val) cnvflg(i)=.false.
endif
enddo
c
c Compute zeta for prog closure

!> - For progsigma = T, calculate the xi term in Bengtsson et al. 2022 \cite Bengtsson_2022 (equation 8)
do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
Expand Down Expand Up @@ -1601,6 +1607,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
enddo
endif
c

do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
Expand Down Expand Up @@ -1926,7 +1933,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
c compute cloud base mass flux as a function of the mean
c updraft velcoity
c
c Prognostic closure
!> - 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
flag_shallow = .true.
call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
Expand All @@ -1937,6 +1944,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &

!> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity.
!! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv.
do i= 1, im
if(cnvflg(i)) then
k = kbcon(i)
Expand Down Expand Up @@ -1965,7 +1973,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
do i = 1, im
if(cnvflg(i)) then
if (gdx(i) < dxcrt) then
scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i))
if(progsigma)then
scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i))
else
Expand Down Expand Up @@ -2448,9 +2455,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
if(k > kb(i) .and. k < ktop(i)) then
tem = 0.5 * (eta(i,k-1) + eta(i,k)) * xmb(i)
tem1 = pfld(i,k) * 100. / (rd * t1(i,k))
sigmagfm(i) = max(sigmagfm(i), betaw)
ptem = tem / (sigmagfm(i) * tem1)
qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*sigmagfm(i)*ptem*ptem
if(progsigma)then
tem2 = sigmab(i)
else
tem2 = max(sigmagfm(i), betaw)
endif
ptem = tem / (tem2 * tem1)
qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*tem2*ptem*ptem
endif
endif
enddo
Expand Down

0 comments on commit 3fe7b64

Please sign in to comment.