Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
141 changes: 110 additions & 31 deletions physics/satmedmfvdifq.F
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
Comment thread
BinLiu-NOAA marked this conversation as resolved.
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
Comment thread
BinLiu-NOAA marked this conversation as resolved.
!!
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
Comment thread
BinLiu-NOAA marked this conversation as resolved.
!> - 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
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