Skip to content

Commit

Permalink
Merge pull request NCAR#27 from hafs-community/feature/hafs_sync_202210
Browse files Browse the repository at this point in the history 
Add the tc_pbl option used by HAFSv1 for GFS TKE EDMF PBL scheme
  • Loading branch information
@grantfirl
grantfirl authored Dec 14, 2022
2 parents befc336 + bf4ac7e commit 688c771
Show file tree
Hide file tree
Showing 2 changed files with 117 additions and 31 deletions.
141 changes: 110 additions & 31 deletions physics/satmedmfvdifq.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,12 @@ module satmedmfvdifq
!! Updated version of satmedmfvdif.f (May 2019) to have better low level
!! inversion, to reduce the cold bias in lower troposphere,
!! and to reduce the negative wind speed bias in upper troposphere
!!
!! Incorporate the LES-based changes for TC simulation
!! (Chen et al.,2022, https://doi.org/10.1175/WAF-D-21-0168.1)
!! with additional improvements on MF working with Cu schemes
!! Xiaomin Chen, 5/2/2022
!!
!> \section arg_table_satmedmfvdifq_init Argument Table
!! \htmlinclude satmedmfvdifq_init.html
!!
Expand Down Expand Up @@ -75,7 +81,7 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
& prsi,del,prsl,prslk,phii,phil,delt, &
& dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku, &
& kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr, &
& rlmx,elmx,sfc_rlm, &
& rlmx,elmx,sfc_rlm,tc_pbl, &
& ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind, &
& index_of_y_wind,index_of_process_pbl,gen_tend,ldiag3d, &
& errmsg,errflg)
Expand All @@ -89,6 +95,7 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
integer, intent(in) :: im, km, ntrac, ntcw, ntrw, ntiw, &
& ntke, ntqv
integer, intent(in) :: sfc_rlm
integer, intent(in) :: tc_pbl
integer, intent(in) :: kinver(:)
integer, intent(out) :: kpbl(:)
logical, intent(in) :: gen_tend,ldiag3d,progsigma
Expand Down Expand Up @@ -242,7 +249,10 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
real(kind=kind_phys) qlcr, zstblmax, hcrinv
!
real(kind=kind_phys) h1

real(kind=kind_phys) bfac, mffac
!!
parameter(bfac=100.)
parameter(wfac=7.0,cfac=4.5)
parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1)
parameter(vk=0.4,rimin=-100.,slfac=0.1)
Expand All @@ -259,10 +269,19 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
parameter(qlcr=3.5e-5,zstblmax=2500.)
parameter(xkinv1=0.15,xkinv2=0.3)
parameter(h1=0.33333333,hcrinv=250.)
parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15)
parameter(ce0=0.4,cs0=0.2)
parameter(ck1=0.15,ch1=0.15)
parameter(cs0=0.2)
parameter(rchck=1.5,ndt=20)

if (tc_pbl == 0) then
ck0 = 0.4
ch0 = 0.4
ce0 = 0.4
else if (tc_pbl == 1) then
ck0 = 0.55
ch0 = 0.55
ce0 = 0.12
endif
gravi = 1.0 / grav
g = grav
gocp = g / cp
Expand Down Expand Up @@ -600,11 +619,18 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
do i = 1, im
if(.not.flg(i)) then
rbdn(i) = rbup(i)
spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
! rbup(i) = (thvx(i,k)-thermal(i))*
! & (g*zl(i,k)/thvx(i,1))/spdk2
rbup(i) = (thlvx(i,k)-thermal(i))*
& (g*zl(i,k)/thlvx(i,1))/spdk2
if (tc_pbl == 0) then
spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
! rbup(i) = (thvx(i,k)-thermal(i))*
! & (g*zl(i,k)/thvx(i,1))/spdk2
rbup(i) = (thlvx(i,k)-thermal(i))*
& (g*zl(i,k)/thlvx(i,1))/spdk2
else if (tc_pbl == 1) then
spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.)
& + bfac*ustar(i)**2
rbup(i) = (thlvx(i,k)-thermal(i))*
& (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2
endif
kpblx(i) = k
flg(i) = rbup(i) > crb(i)
endif
Expand Down Expand Up @@ -674,11 +700,18 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
do i = 1, im
if(.not.flg(i) .and. k > kb1(i)) then
rbdn(i) = rbup(i)
spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
! rbup(i) = (thvx(i,k)-thermal(i))*
! & (g*zl(i,k)/thvx(i,1))/spdk2
rbup(i) = (thlvx(i,k)-thermal(i))*
& (g*zl(i,k)/thlvx(i,1))/spdk2
if (tc_pbl == 0) then
spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
! rbup(i) = (thvx(i,k)-thermal(i))*
! & (g*zl(i,k)/thvx(i,1))/spdk2
rbup(i) = (thlvx(i,k)-thermal(i))*
& (g*zl(i,k)/thlvx(i,1))/spdk2
else if (tc_pbl == 1) then
spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.)
& + bfac*ustar(i)**2
rbup(i) = (thlvx(i,k)-thermal(i))*
& (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2
endif
kpblx(i) = k
flg(i) = rbup(i) > crb(i)
endif
Expand Down Expand Up @@ -796,9 +829,16 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
do i = 1, im
if(.not.flg(i) .and. k > kb1(i)) then
rbdn(i) = rbup(i)
spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
rbup(i) = (thlvx(i,k)-thermal(i))*
& (g*zl(i,k)/thlvx(i,1))/spdk2
if (tc_pbl == 0) then
spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.)
rbup(i) = (thlvx(i,k)-thermal(i))*
& (g*zl(i,k)/thlvx(i,1))/spdk2
else if (tc_pbl == 1) then
spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.)
& + bfac*ustar(i)**2
rbup(i) = (thlvx(i,k)-thermal(i))*
& (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2
endif
kpbl(i) = k
flg(i) = rbup(i) > crb(i)
endif
Expand Down Expand Up @@ -929,6 +969,26 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
& thlx,thvx,thlvx,gdx,thetae,
& krad,mrad,radmin,buod,xmfd,
& tcdo,qcdo,ucdo,vcdo,xlamde,bl_dnfr)
if (tc_pbl == 1) then
!> - unify mass fluxes with Cu
do i=1,im
if(zol(i) > -0.5) then
do k = 1, km
xmf(i,k) = 0.0
end do
end if
end do
!> - taper off MF in high-wind conditions
do i = 1,im
tem = sqrt(u10m(i)**2+v10m(i)**2)
mffac = (1. - MIN(MAX(tem - 20.0, 0.0), 10.0)/10.)
do k = 1, km
xmf(i,k) = xmf(i,k)*mffac
xmfd(i,k) = xmfd(i,k)*mffac
enddo
enddo
endif
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Expand Down Expand Up @@ -1013,10 +1073,15 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
tem3=((u1(i,n+1)-u1(i,n))/dz)**2
tem3=tem3+((v1(i,n+1)-v1(i,n))/dz)**2
tem3=cs0*sqrt(tem3)*sqrt(tke(i,k))
ptem = (gotvx(i,n)*(thvx(i,n+1)-thvx(i,k))+tem3)*dz
! ptem = (gotvx(i,n)*(thlvx(i,n+1)-thlvx(i,k)+tem3)*dz
! ptem = (gotvx(i,n)*(tem1-thvx(i,k))+tem3)*dz
!! ptem = (gotvx(i,n)*(tem1-thlvx(i,k)+tem3)*dz
if (tc_pbl == 0) then
ptem = (gotvx(i,n)*(thvx(i,n+1)-thvx(i,k))+tem3)*dz
! ptem = (gotvx(i,n)*(thlvx(i,n+1)-thlvx(i,k)+tem3)*dz
! ptem = (gotvx(i,n)*(tem1-thvx(i,k))+tem3)*dz
!! ptem = (gotvx(i,n)*(tem1-thlvx(i,k)+tem3)*dz
else if (tc_pbl == 1) then
tem1 = 0.5*(thvx(i,n+1)+thvx(i,n))
ptem = (gotvx(i,n)*(tem1-thvx(i,k))+tem3)*dz
endif
bsum = bsum + ptem
zlup = zlup + dz
if(bsum >= tke(i,k)) then
Expand Down Expand Up @@ -1047,10 +1112,14 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
tem3 = cs0*sqrt(tem3)*sqrt(tke(i,1))
else
dz = zl(i,n) - zl(i,n-1)
tem1 = thvx(i,n-1)
! tem1 = thlvx(i,n-1)
! tem1 = 0.5 * (thvx(i,n-1) + thvx(i,n))
!! tem1 = 0.5 * (thlvx(i,n-1) + thlvx(i,n))
if (tc_pbl == 0) then
tem1 = thvx(i,n-1)
! tem1 = thlvx(i,n-1)
! tem1 = 0.5 * (thvx(i,n-1) + thvx(i,n))
!! tem1 = 0.5 * (thlvx(i,n-1) + thlvx(i,n))
else if (tc_pbl == 1) then
tem1 = 0.5*(thvx(i,n)+thvx(i,n-1))
endif
tem3 = ((u1(i,n)-u1(i,n-1))/dz)**2
tem3 = tem3+((v1(i,n)-v1(i,n-1))/dz)**2
tem3 = cs0*sqrt(tem3)*sqrt(tke(i,k))
Expand Down Expand Up @@ -1117,19 +1186,29 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, &
else
ptem = 1. + 2.7 * zol(i)
zk = tem / ptem
endif
elm(i,k) = zk*rlam(i,k)/(rlam(i,k)+zk)
endif
if (tc_pbl == 0) then
elm(i,k) = zk*rlam(i,k)/(rlam(i,k)+zk)
!> - If sfc_rlm=1, use zk for elm within surface layer
if ( sfc_rlm == 1 ) then
if ( sfcflg(i) .and.
& zl(i,k) < min(100.0,hpbl(i)*0.05) ) elm(i,k)=zk
if ( sfc_rlm == 1 ) then
if ( sfcflg(i) .and.
& zl(i,k) < min(100.0,hpbl(i)*0.05) ) elm(i,k)=zk
endif
else if (tc_pbl == 1) then
! new blending method for mixing length
elm(i,k) = sqrt( 1.0/( 1.0/(zk**2)+1.0/(rlam(i,k)**2) ) )
endif
!
dz = zi(i,k+1) - zi(i,k)
tem = max(gdx(i),dz)
elm(i,k) = min(elm(i,k), tem)
ele(i,k) = min(ele(i,k), tem)
if (tc_pbl == 0) then
ele(i,k) = min(ele(i,k), tem)
else if (tc_pbl == 1) then
ele(i,k) = elm(i,k)
endif
!
enddo
enddo
Expand Down
7 changes: 7 additions & 0 deletions physics/satmedmfvdifq.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -588,6 +588,13 @@
dimensions = ()
type = integer
intent = in
[tc_pbl]
standard_name = control_for_TC_applications_in_the_PBL_scheme
long_name = control for TC applications in the PBL scheme
units = none
dimensions = ()
type = integer
intent = in
[ntqv]
standard_name = index_of_specific_humidity_in_tracer_concentration_array
long_name = tracer index for water vapor (specific humidity)
Expand Down

0 comments on commit 688c771

Please sign in to comment.