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

Correction to prog closure (convection) and ugwp stoch phys fix #57

Merged
merged 6 commits into from
Apr 6, 2023
Merged

Correction to prog closure (convection) and ugwp stoch phys fix #57

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
be146a8
Correction to convective cloud condensate term in prog closure for HR2
lisa-bengtsson Mar 24, 2023
1553636
Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physi…
lisa-bengtsson Mar 24, 2023
aa7088e
Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physi…
lisa-bengtsson Apr 3, 2023
6424df1
Fixed issue with level of dividing streamline (rdxzb) being overwritt…
mdtoyNOAA Apr 4, 2023
eb816e3
Added back CCPP header
mdtoyNOAA Apr 5, 2023
94eaf51
Merge remote-tracking branch 'mike/ufs/dev_ugwp_stoch_phys_fix' into …
lisa-bengtsson Apr 5, 2023
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
28 changes: 6 additions & 22 deletions physics/progsigma_calc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@
!!\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, &
sigmab,errmsg,errflg)
delt,qadv,kbcon1,ktcon,cnvflg,sigmain,sigmaout, &
sigmab)
!
!
use machine, only : kind_phys
Expand All @@ -25,7 +25,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, &
! intent in
integer, intent(in) :: im,km,kbcon1(im),ktcon(im)
real(kind=kind_phys), intent(in) :: hvap,delt
real(kind=kind_phys), intent(in) :: prevsq(im,km), q(im,km),del(im,km), &
real(kind=kind_phys), intent(in) :: qadv(im,km),del(im,km), &
qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), &
omega_u(im,km),zeta(im,km)
logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow
Expand All @@ -34,14 +34,13 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, &
! intent out
real(kind=kind_phys), intent(out) :: sigmaout(im,km)
real(kind=kind_phys), intent(out) :: sigmab(im)
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg


! 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), &
qadv(im,km),dp(im,km),inbu(im,km)
dp(im,km),inbu(im,km)


real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, &
Expand Down Expand Up @@ -77,21 +76,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, &
mcons(i)=0.
enddo

!Initial computations, dynamic q-tendency
if(flag_init .and. .not.flag_restart)then
do k = 1,km
do i = 1,im
qadv(i,k)=0.
enddo
enddo
else
do k = 1,km
do i = 1,im
qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt
enddo
enddo
endif

do k = 2,km1
do i = 1,im
if(cnvflg(i))then
Expand Down Expand Up @@ -133,7 +117,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, &
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(k > kbcon1(i) .and. k < ktcon(i) .and. buy2 > 0.)then
if(dbyo1(i,k)>0 .and. buy2 > 0.)then
inbu(i,k)=1.
endif
inbu(i,k-1)=MAX(inbu(i,k-1),inbu(i,k))
Expand Down
59 changes: 33 additions & 26 deletions physics/samfdeepcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
real(kind=kind_phys), intent(in) :: nthresh
real(kind=kind_phys), intent(in) :: ca_deep(:)
real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), &
& tmf(:,:),q(:,:), prevsq(:,:)
& tmf(:,:,:),q(:,:), prevsq(:,:)
real(kind=kind_phys), intent(out) :: rainevap(:)
real(kind=kind_phys), intent(out) :: sigmaout(:,:)
logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger
Expand Down Expand Up @@ -209,9 +209,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
& bb1, bb2, wucb
!
! parameters for prognostic sigma closure
real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km),
& omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im)
real(kind=kind_phys) gravinv
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)
real(kind=kind_phys) gravinv,invdelt
logical flag_shallow
c physical parameters
! parameter(grav=grav,asolfac=0.958)
Expand Down Expand Up @@ -306,6 +306,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
errflg = 0

gravinv = 1./grav
invdelt = 1./delt

elocp = hvap/cp
el2orc = hvap*hvap/(rv*cp)
Expand Down Expand Up @@ -585,7 +586,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
do i = 1, im
dbyo1(i,k)=0.
zdqca(i,k)=0.
qlks(i,k)=0.
omega_u(i,k)=0.
zeta(i,k)=1.0
enddo
Expand Down Expand Up @@ -1515,7 +1515,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
pwavo(i) = pwavo(i) + pwo(i,k)
! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp
cnvwt(i,k) = etah * qlk * grav / dp
qlks(i,k)=qlk
zdqca(i,k)=dq/eta(i,k)
endif
!
! compute buoyancy and drag for updraft velocity
Expand Down Expand Up @@ -1690,7 +1690,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
pwavo(i) = pwavo(i) + pwo(i,k)
! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp
cnvwt(i,k) = etah * qlk * grav / dp
qlks(i,k)=qlk
zdqca(i,k)=dq/eta(i,k)
endif
endif
endif
Expand Down Expand Up @@ -1860,28 +1860,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &
if(dq > 0.) then
qlko_ktcon(i) = dq
qcko(i,k) = qrch
qlks(i,k) = qlko_ktcon(i)
zdqca(i,k) = dq
endif
endif
enddo
endif
c

c store term needed for "termC" in prognostic area fraction closure
if(progsigma)then
do k = 2, km1
do i = 1, im
dp = 1000. * del(i,k)
if (cnvflg(i)) then
if(k > kbcon(i) .and. k < ktcon(i)) then
zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) +
& pwo(i,k)+dellal(i,k))
endif
endif
enddo
enddo
endif

ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then
ccccc print *, ' aa1(i) before dwndrft =', aa1(i)
ccccc endif
Expand Down Expand Up @@ -2885,11 +2870,33 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, &

!> - 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

!Initial computations, dynamic q-tendency
if(first_time_step .and. .not.restart)then
do k = 1,km
do i = 1,im
qadv(i,k)=0.
enddo
enddo
else
do k = 1,km
do i = 1,im
qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt
enddo
enddo
endif

do k = 1,km
do i = 1,im
tmfq(i,k)=tmf(i,k,1)
grantfirl marked this conversation as resolved.
Show resolved Hide resolved
enddo
enddo

flag_shallow = .false.
call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
& del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
& prevsq,q,kbcon1,ktcon,cnvflg,
& sigmain,sigmaout,sigmab,errmsg,errflg)
& del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
& qadv,kbcon1,ktcon,cnvflg,
& sigmain,sigmaout,sigmab)
endif

!> - 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.
Expand Down
6 changes: 3 additions & 3 deletions physics/samfdeepcnv.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,10 +70,10 @@
type = logical
intent = in
[tmf]
standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
standard_name = tendency_of_vertically_diffused_tracer_concentration
long_name = updated tendency of the tracers due to vertical diffusion in PBL scheme
units = kg kg-1 s-1
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_vertical_diffusion_tracers)
type = real
kind = kind_phys
intent = in
Expand Down
54 changes: 32 additions & 22 deletions physics/samfshalcnv.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
real(kind=kind_phys), intent(in) :: delt
real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), &
& prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), &
& qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:)
& qmicro(:,:),tmf(:,:,:),prevsq(:,:),q(:,:)

real(kind=kind_phys), intent(in) :: sigmain(:,:)
!
Expand Down Expand Up @@ -156,10 +156,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
cc

! parameters for prognostic sigma closure
real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km),
real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km),
& omegac(im),zeta(im,km),dbyo1(im,km),
& sigmab(im)
real(kind=kind_phys) gravinv,dxcrtas
& sigmab(im),qadv(im,km)
real(kind=kind_phys) gravinv,dxcrtas,invdelt
logical flag_shallow
c physical parameters
! parameter(g=grav,asolfac=0.89)
Expand Down Expand Up @@ -247,6 +247,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
errflg = 0

gravinv = 1./grav
invdelt = 1./delt

elocp = hvap/cp
el2orc = hvap*hvap/(rv*cp)
Expand Down Expand Up @@ -524,7 +525,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
do i = 1, im
dbyo1(i,k)=0.
zdqca(i,k)=0.
qlks(i,k)=0.
omega_u(i,k)=0.
zeta(i,k)=1.0
enddo
Expand Down Expand Up @@ -1270,7 +1270,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
qcko(i,k)= qlk + qrch
pwo(i,k) = etah * c0t(i,k) * dz * qlk
cnvwt(i,k) = etah * qlk * grav / dp
qlks(i,k)=qlk
zdqca(i,k)=dq/eta(i,k)
endif
!
! compute buoyancy and drag for updraft velocity
Expand Down Expand Up @@ -1435,7 +1435,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
qcko(i,k) = qlk + qrch
pwo(i,k) = etah * c0t(i,k) * dz * qlk
cnvwt(i,k) = etah * qlk * grav / dp
qlks(i,k)=qlk
zdqca(i,k)=dq/eta(i,k)
endif
endif
endif
Expand Down Expand Up @@ -1601,24 +1601,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
if(dq > 0.) then
qlko_ktcon(i) = dq
qcko(i,k) = qrch
qlks(i,k) = qlko_ktcon(i)
zdqca(i,k) = dq
endif
endif
enddo
endif
c

do k = 2, km1
do i = 1, im
if (cnvflg(i)) then
if(k > kbcon(i) .and. k < ktcon(i)) then
zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) +
& pwo(i,k)+dellal(i,k))
endif
endif
enddo
enddo

c--- compute precipitation efficiency in terms of windshear
c
!! - Calculate the wind shear and precipitation efficiency according to equation 58 in Fritsch and Chappell (1980) \cite fritsch_and_chappell_1980 :
Expand Down Expand Up @@ -1935,11 +1924,32 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, &
c
!> - 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
! Initial computations, dynamic q-tendency
if(first_time_step .and. .not.restart)then
do k = 1,km
do i = 1,im
qadv(i,k)=0.
enddo
enddo
else
do k = 1,km
do i = 1,im
qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt
enddo
enddo
endif

do k = 1,km
do i = 1,im
tmfq(i,k)=tmf(i,k,1)
enddo
enddo

flag_shallow = .true.
call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
& del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
& prevsq,q,kbcon1,ktcon,cnvflg,
& sigmain,sigmaout,sigmab,errmsg,errflg)
& del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
& qadv,kbcon1,ktcon,cnvflg,
& sigmain,sigmaout,sigmab)
endif

!> - 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.
Expand Down
6 changes: 3 additions & 3 deletions physics/samfshalcnv.meta
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,10 +70,10 @@
type = logical
intent = in
[tmf]
standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
standard_name = tendency_of_vertically_diffused_tracer_concentration
long_name = updated tendency of the tracers due to vertical diffusion in PBL scheme
units = kg kg-1 s-1
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_vertical_diffusion_tracers)
type = real
kind = kind_phys
intent = in
Expand Down
24 changes: 4 additions & 20 deletions physics/satmedmfvdifq.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,9 +73,9 @@ end subroutine satmedmfvdifq_init
!! -# A mass-flux approach is also used to represent the stratocumulus-top-induced turbulence
!! (mfscuq.f).
!! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm
subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
& ntiw,ntke,grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, &
& dv,du,tdt,rtg,tmf,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, &
& dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, &
& psk,rbsoil,zorl,u10m,v10m,fm,fh, &
& tsea,heat,evap,stress,spd1,kpbl, &
& prsi,del,prsl,prslk,phii,phil,delt, &
Expand All @@ -98,15 +98,15 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
integer, intent(in) :: tc_pbl
integer, intent(in) :: kinver(:)
integer, intent(out) :: kpbl(:)
logical, intent(in) :: gen_tend,ldiag3d,progsigma
logical, intent(in) :: gen_tend,ldiag3d
!
real(kind=kind_phys), intent(in) :: grav,rd,cp,rv,hvap,hfus,fv, &
& eps,epsm1
real(kind=kind_phys), intent(in) :: delt, xkzm_m, xkzm_h, xkzm_s
real(kind=kind_phys), intent(in) :: dspfac, bl_upfr, bl_dnfr
real(kind=kind_phys), intent(in) :: rlmx, elmx
real(kind=kind_phys), intent(inout) :: dv(:,:), du(:,:), &
& tdt(:,:), rtg(:,:,:), tmf(:,:)
& tdt(:,:), rtg(:,:,:)
real(kind=kind_phys), intent(in) :: &
& u1(:,:), v1(:,:), &
& t1(:,:), q1(:,:,:), &
Expand Down Expand Up @@ -331,14 +331,6 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
zm(i,k) = zi(i,k+1)
enddo
enddo
!> - Initialize variables needed for prognostic cumulus closure
if(progsigma)then
do k=1,km
do i=1,im
tmf(i,k) = 0.
enddo
enddo
endif
!> - Compute horizontal grid size (\p gdx)
do i=1,im
gdx(i) = sqrt(garea(i))
Expand Down Expand Up @@ -2206,14 +2198,6 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
enddo
enddo

if(progsigma)then
do k = 1,km
do i = 1,im
tmf(i,k)=(f2(i,k)-q1(i,k,1))*rdt
enddo
enddo
endif

do i = 1,im
dtsfc(i) = rho_a(i) * cp * heat(i)
dqsfc(i) = rho_a(i) * hvap * evap(i)
Expand Down
Loading