Skip to content

Commit

Permalink
3D-TKE EMDF GFS PBL scheme related changes from FIU (Ping Zhu, Ping.Z…
Browse files Browse the repository at this point in the history 
…[email protected]).
  • Loading branch information
@BinLiu-NOAA
BinLiu-NOAA committed Jul 18, 2024
1 parent 07b0f1b commit 9a943e2
Show file tree
Hide file tree
Showing 2 changed files with 219 additions and 3 deletions.
199 changes: 196 additions & 3 deletions physics/PBL/SATMEDMF/satmedmfvdifq.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ end subroutine satmedmfvdifq_init
!! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm
subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
& ntiw,ntke,grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, &
& def_1,def_2,sa3dtke, &
& dv,du,tdt,rtg,u1,v1,t1,q1,usfco,vsfco,icplocn2atm, &
& swh,hlw,xmu,garea,zvfun,sigmaf, &
& psk,rbsoil,zorl,u10m,v10m,fm,fh, &
Expand Down Expand Up @@ -135,6 +136,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
real(kind=kind_phys), intent(out) :: &
& dkt(:,:), dku(:,:)
!
real(kind=kind_phys), intent(in) :: &
& def_1(:,:), def_2(:,:)
logical, intent(in) :: sa3dtke !flag for scale aware 3d TKE scheme

logical, intent(in) :: dspheat
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
Expand Down Expand Up @@ -247,6 +252,19 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
& tem, tem1, tem2, tem3,
& ptem, ptem0, ptem1, ptem2
!
integer kk
real(kind=kind_phys) thetal(im,km),dku1(im,km),dkt1(im,km),
& elm_les(im,km),ele_les(im,km),pftke(im),
& dkq1(im,km),pfl(im),cpl1,cpl2,cpl3,cpl4,
& cpl5,cpl6,pfnl(im),cpnl1,cpnl2,cpnl3,cpnl4,
& cpnl5,cpnl6,cptke1,cptke2,cptke3,shrp3d
real(kind=kind_phys) ak(im,km),bk(im,km),ck(im,km),pk(im,km),
& f3(im,km),rizt(im,km)
real(kind=kind_phys) aa1,aa2,bb1,bb2,cc1,alpha1, alpha2, alpha3,
& alpha4,alpha5,ssh,ssm,gh,aa,bb,cc,dd
integer l_tkemax,kscl
real(kind=kind_phys) tkemax, scl, kmaxles, esmax

real(kind=kind_phys) slfac
!
real(kind=kind_phys) vegflo, vegfup, z0lo, z0up, vc0, zc0
Expand Down Expand Up @@ -281,6 +299,13 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
parameter(ck1=0.15,ch1=0.15)
parameter(cs0=0.4,csmf=0.5)
parameter(rchck=1.5,ndt=20)
parameter(cpl1=0.280,cpl2=0.870,cpl3=0.913)
parameter(cpl4=0.153,cpl5=0.278,cpl6=0.720)
parameter(cpnl1=0.243,cpnl2=0.936,cpnl3=1.11)
parameter(cpnl4=0.312,cpnl5=0.329,cpnl6=0.757)
parameter(cptke1=0.07,cptke2=0.142,cptke3=0.071)
parameter(aa1=0.92,aa2=0.74,bb1=16.6,bb2=10.1,cc1=0.08)
parameter(kmaxles=300.0,esmax=500.0)

if (tc_pbl == 0) then
ck0 = 0.4
Expand Down Expand Up @@ -467,7 +492,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
theta(i,k) = t1(i,k) * pix(i,k)
if(ntiw > 0) then
tem = max(q1(i,k,ntcw),qlmin)
if(sa3dtke) then
! include snow in the buoyancy calculation
tem1 = max(q1(i,k,ntiw)+q1(i,k,5),qlmin)
else
tem1 = max(q1(i,k,ntiw),qlmin)
endif
qlx(i,k) = tem + tem1
ptem = hvap*tem + (hvap+hfus)*tem1
slx(i,k) = cp * t1(i,k) + phil(i,k) - ptem
Expand Down Expand Up @@ -1292,6 +1322,110 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
enddo
enddo
! compute km, kh, kq for LES 3D TKE scheme (Deardorff 1980)
if(sa3dtke) then
! calculate thetal
do k=1,km
do i=1,im
pix(i,k) = psk(i) / prslk(i,k)
theta(i,k) = t1(i,k) * pix(i,k)
tem=theta(i,k)/t1(i,k)
if(ntiw > 0) then
tem1=max(q1(i,k,ntcw),qlmin)+
& max(q1(i,k,ntiw)+q1(i,k,5),qlmin)
thetal(i,k)=theta(i,k)-(hvap+hfus)/cp*tem*tem1
else
tem1=max(q1(i,k,ntcw),qlmin)
thetal(i,k)=theta(i,k)-hvap/cp*tem*tem1
endif
enddo
enddo
do k=1,km
do i=1,im
dku1(i,k) = 0.
dkt1(i,k) = 0.
dkq1(i,k) = 0.
enddo
enddo
do k=1,km
do i=1,im
dku1(i,k) = 0.
dkt1(i,k) = 0.
dkq1(i,k) = 0.
enddo
enddo
do k = 1, km1
do i = 1, im
dz = zi(i,k+1) - zi(i,k)
tem=grav/theta(i,1)*(thetal(i,k+1)-thetal(i,k))/dz
tem1=(garea(i)*dz)**(1.0/3.0)
! calculate LES mixing length
if(tem > 0.0) then
elm_les(i,k)=0.76*sqrt(tkeh(i,k))*tem**(-0.5)
elm_les(i,k)=min(elm_les(i,k),tem1)
else
elm_les(i,k)=tem1
endif
ele_les(i,k)=elm_les(i,k)
! calculate km, kh, and kq for LES
dku1(i,k)=0.1*elm_les(i,k)*sqrt(tkeh(i,k))
dkt1(i,k)=(1.0+2.0*elm_les(i,k)/tem1)*dku1(i,k)
dkq1(i,k)=dkt1(i,k)
dku1(i,k) = min(dku1(i,k),kmaxles)
dkt1(i,k) = min(dkt1(i,k),kmaxles)
dkq1(i,k) = min(dkq1(i,k),kmaxles)
enddo
enddo
do i = 1, im
! d/zi
scl=1000.
l_tkemax=10
kscl=10
tkemax=0.0
do k=1,km
tkemax=max(tkemax,tke(i,k))
enddo
do k=1,km
if (abs(tke(i,k)-tkemax)/tkemax .lt. 1.0e-9) then
l_tkemax=k
endif
enddo
do k=l_tkemax,km
if (tke(i,k)-0.5*tkemax .gt. 0.0) then
kscl=k
endif
enddo
kscl=min(kscl,km-10)
scl=zi(i,kscl+1)
! tem=gdx(i)/max(hpbl(i),scl)
tem=gdx(i)/max(hpbl(i),esmax)
! partition function for local fluxes
pfl(i)=cpl1*(tem**2+cpl2*tem**0.5-cpl3)/
& (tem**2+cpl4*tem**0.5+cpl5)+cpl6
pfl(i)=min(max(pfl(i),0.0),1.0)
! partition function for nonlocal fluxes
pfnl(i)=cpnl1*(tem**2+cpnl2*tem**(7./8.)-cpnl3)/
& (tem**2+cpnl4*tem**(7./8.)+cpnl5)+cpnl6
pfnl(i)=min(max(pfnl(i),0.0),1.0)
pftke(i)=(tem**2+cptke1*tem**(2./3.))/
& (tem**2+cptke2*tem**(2./3.)+cptke3)
pftke(i)=min(max(pftke(i),0.0),1.0)
enddo
! blending kq,kt, and km from LES and MS
do k = 1,km
do i=1,im
dkq(i,k)=(1.0-pfl(i))*dkq1(i,k)+pfl(i)*dkq(i,k)
dkt(i,k)=(1.0-pfl(i))*dkt1(i,k)+pfl(i)*dkt(i,k)
dku(i,k)=(1.0-pfl(i))*dku1(i,k)+pfl(i)*dku(i,k)
enddo
enddo
endif
!> ## Compute TKE.
!! - Compute a minimum TKE deduced from background diffusivity for momentum.
!
Expand Down Expand Up @@ -1430,6 +1564,17 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
endif
shrp = shrp + ptem1 + ptem2
endif
if(sa3dtke) then
! obtaining 3d shear production from dycore
if (k ==1) then
shrp3d=dku(i,k)*def_1(i,k)
else
tem1 = dku(i,k-1) * def_1(i,k-1)
tem2 = dku(i,k) * def_1(i,k)
shrp3d=0.5*(tem1+tem2)
endif
shrp=shrp3d
endif
prod(i,k) = buop + shrp
enddo
enddo
Expand All @@ -1441,14 +1586,29 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
do n = 1, ndt
do k = 1,km1
do i=1,im
tem = sqrt(tke(i,k))
tem = sqrt(tke(i,k))
if(sa3dtke) then
! calculating 3D TKE transport and pressure correlation
ptem1 = ce0 / ele(i,k)
tem1=(garea(i)*dz)**(1.0/3.0)
tem2=0.19+0.51*ele_les(i,k)/tem1
ptem2= tem2 / ele_les(i,k)
ptem=(1.0-pftke(i))*ptem2+pftke(i)*ptem1
diss(i,k) = ptem * tke(i,k) * tem
tem1 = prod(i,k) + tke(i,k) / dtn
diss(i,k)=max(min(diss(i,k), tem1), 0.)
tem=2.0*dkq(i,k)*def_2(i,k)
tem=min(tem,1.0)
tke(i,k) = tke(i,k) + dtn * (prod(i,k)-diss(i,k)+tem)
else
ptem = ce0 / ele(i,k)
diss(i,k) = ptem * tke(i,k) * tem
tem1 = prod(i,k) + tke(i,k) / dtn
diss(i,k)=max(min(diss(i,k), tem1), 0.)
tke(i,k) = tke(i,k) + dtn * (prod(i,k)-diss(i,k))
! tke(i,k) = max(tke(i,k), tkmin)
tke(i,k) = max(tke(i,k), tkmnz(i,k))
endif
! tke(i,k) = max(tke(i,k), tkmin)
tke(i,k) = max(tke(i,k), tkmnz(i,k))
enddo
enddo
enddo
Expand Down Expand Up @@ -1634,6 +1794,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
if(pcnvflg(i) .and. k < kpbl(i)) then
ptem = 0.5 * tem2 * xmf(i,k)
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke)
Expand All @@ -1646,6 +1809,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
if(scuflg(i)) then
if(k >= mrad(i) .and. k < krad(i)) then
ptem = 0.5 * tem2 * xmfd(i,k)
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke)
Expand All @@ -1657,6 +1823,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
kmx = max(kpbl(i), krad(i))
if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
ptem = tem2 * (xmf(i,k) - xmfd(i,k))
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
f1(i,k) = f1(i,k) + e_half(i,k) * ptem1
Expand Down Expand Up @@ -1825,6 +1994,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
if(pcnvflg(i) .and. k < kpbl(i)) then
ptem = 0.5 * tem2 * xmf(i,k)
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
tem = t1(i,k) + t1(i,k+1)
Expand All @@ -1843,6 +2015,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
if(scuflg(i)) then
if(k >= mrad(i) .and. k < krad(i)) then
ptem = 0.5 * tem2 * xmfd(i,k)
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
ptem = tcdo(i,k) + tcdo(i,k+1)
Expand All @@ -1858,6 +2033,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
kmx = max(kpbl(i), krad(i))
if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
ptem = tem2 * (xmf(i,k) - xmfd(i,k))
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
f2(i,k) = f2(i,k) + q_half(i,k,1) * ptem1
Expand All @@ -1879,6 +2057,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
if(pcnvflg(i) .and. k < kpbl(i)) then
ptem = 0.5 * tem2 * xmf(i,k)
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
ptem = qcko(i,k,n) + qcko(i,k+1,n)
Expand All @@ -1891,6 +2072,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
if(scuflg(i)) then
if(k >= mrad(i) .and. k < krad(i)) then
ptem = 0.5 * tem2 * xmfd(i,k)
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
ptem = qcdo(i,k,n) + qcdo(i,k+1,n)
Expand All @@ -1902,6 +2086,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
kmx = max(kpbl(i), krad(i))
if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then
ptem = tem2 * (xmf(i,k) - xmfd(i,k))
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
f2(i,k+is) = f2(i,k+is) + q_half(i,k,n) * ptem1
Expand Down Expand Up @@ -2328,6 +2515,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
!
if(pcnvflg(i) .and. k < kpbl(i)) then
ptem = 0.5 * tem2 * xmf(i,k)
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
tem = u1(i,k) + u1(i,k+1)
Expand All @@ -2346,6 +2536,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
if(scuflg(i)) then
if(k >= mrad(i) .and. k < krad(i)) then
ptem = 0.5 * tem2 * xmfd(i,k)
if(sa3dtke) then
ptem=ptem*pfnl(i)
endif
ptem1 = dtodsd * ptem
ptem2 = dtodsu * ptem
tem = u1(i,k) + u1(i,k+1)
Expand Down
23 changes: 23 additions & 0 deletions physics/PBL/SATMEDMF/satmedmfvdifq.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -256,6 +256,22 @@
type = real
kind = kind_phys
intent = in
[def_1]
standard_name = shear_prod
long_name = 3D shear production
units = m2 s-2
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
type = real
kind = kind_phys
intent = in
[def_2]
standard_name = TKE_transfer
long_name = TKE transfer and pressure correlation
units = m2 s-2
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
type = real
kind = kind_phys
intent = in
[swh]
standard_name = tendency_of_air_temperature_due_to_shortwave_heating_on_radiation_timestep
long_name = total sky shortwave heating rate
Expand Down Expand Up @@ -470,6 +486,13 @@
dimensions = ()
type = logical
intent = in
[sa3dtke]
standard_name = flag_sa_3d_tke_scheme
long_name = flag for scale-aware 3d tke scheme
units = flag
dimensions = ()
type = logical
intent = in
[dusfc]
standard_name = instantaneous_surface_x_momentum_flux
long_name = x momentum flux
Expand Down

0 comments on commit 9a943e2

Please sign in to comment.