From 9c4545fe79fff4e7a8f7814ecb23c88fa2af17ce Mon Sep 17 00:00:00 2001 From: "Bin.Liu" Date: Fri, 9 May 2025 09:37:32 -0500 Subject: [PATCH 1/6] 3D-TKE EMDF GFS PBL scheme related changes from FIU (Ping Zhu, Ping.Zhu@noaa.gov, and Kwun Y. Fung). --- physics/PBL/SATMEDMF/satmedmfvdifq.F | 726 +++++++++++++++++++++++- physics/PBL/SATMEDMF/satmedmfvdifq.meta | 47 ++ 2 files changed, 749 insertions(+), 24 deletions(-) diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.F b/physics/PBL/SATMEDMF/satmedmfvdifq.F index 64e7bc32d..ef71bae63 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.F +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.F @@ -75,6 +75,8 @@ 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, & +!The following three variables are for SA-3D-TKE + & def_1,def_2,def_3,sa3dtke,dku3d_h,dku3d_e, & & dv,du,tdt,rtg,u1,v1,t1,q1,usfco,vsfco,icplocn2atm, & & swh,hlw,xmu,garea,zvfun,sigmaf, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & @@ -112,6 +114,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & u1(:,:), v1(:,:), & & usfco(:), vsfco(:), & & t1(:,:), q1(:,:,:), & +!The following two variables are for SA-3D-TKE + & def_1(:,:), def_2(:,:), def_3(:,:), & & swh(:,:), hlw(:,:), & & xmu(:), garea(:), & & zvfun(:), sigmaf(:), & @@ -136,10 +140,17 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & real(kind=kind_phys), intent(out) :: & & dkt(:,:), dku(:,:) ! + logical, intent(in) :: sa3dtke !flag for SA-3D-TKE scheme + logical, intent(in) :: dspheat character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg +!For passing dku to the dyn_core (SA-3D-TKE scheme) + real(kind=kind_phys), intent(out) :: + & dku3d_h(:,:),dku3d_e(:,:) + + ! flag for tke dissipative heating ! !---------------------------------------------------------------------- @@ -248,6 +259,23 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & tem, tem1, tem2, tem3, & ptem, ptem0, ptem1, ptem2 ! +!The following variables are for SA-3D-TKE + integer kk + real(kind=kind_phys) thetal(im,km),dku_les(im,km),dkt_les(im,km), + & elm_les(im,km),ele_les(im,km),pftke(im), + & dkq_les(im,km),pfl(im), + & dku_h(im,km),dkq_h(im,km),dddmp, + & cpl1,cpl2,cpl3,cpl4, + & cpl5,cpl6,pfnl(im),cpnl1,cpnl2,cpnl3,cpnl4, + & cpnl5,cpnl6,cptke1,cptke2,cptke3 + + 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 @@ -282,6 +310,15 @@ 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) +!The following variables are for SA-3D-TKE + 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,dddmp=0.1) +! parameter(aa1=0.92,aa2=0.649,bb1=17.7,bb2=9.5,cc1=0.08) if (tc_pbl == 0) then ck0 = 0.4 @@ -434,6 +471,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & hpbl(i) = 0. kpblx(i) = 1 hpblx(i) = 0. + pfl(i)=1.0 + pfnl(i)=1.0 + pftke(i)=1.0 pblflg(i)= .true. sfcflg(i)= .true. if(rbsoil(i) > 0.) sfcflg(i) = .false. @@ -462,8 +502,38 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo ! !> - Compute \f$\theta\f$(theta), and \f$q_l\f$(qlx), \f$\theta_e\f$(thetae), -!! \f$\theta_v\f$(thvx),\f$\theta_{l,v}\f$ (thlvx) including ice water - do k=1,km +!! \f$\theta_v\f$(thvx),\f$\theta_{l,v}\f$ (thlvx) including ice water +!The following is for SA-3D-TKE + if(sa3dtke) then + 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) + if(ntiw > 0) then + tem = max(q1(i,k,ntcw),qlmin) + tem1 = max(q1(i,k,ntiw)+q1(i,k,5)+q1(i,k,6),qlmin) + qlx(i,k) = tem + tem1 + ptem = hvap*tem + (hvap+hfus)*tem1 + slx(i,k) = cp * t1(i,k) + phil(i,k) - ptem + else + qlx(i,k) = max(q1(i,k,ntcw),qlmin) + slx(i,k) = cp * t1(i,k) + phil(i,k) - hvap*qlx(i,k) + endif + tem2 = 1.+fv*max(q1(i,k,1),qmin)-qlx(i,k) + thvx(i,k) = theta(i,k) * tem2 + tvx(i,k) = t1(i,k) * tem2 + qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k) + thlx(i,k) = theta(i,k) - pix(i,k)*elocp*qlx(i,k) + thlvx(i,k) = thlx(i,k) * (1. + fv * qtx(i,k)) + svx(i,k) = cp * tvx(i,k) + ptem1 = elocp * pix(i,k) * max(q1(i,k,1),qmin) + thetae(i,k)= theta(i,k) + ptem1 +! gotvx(i,k) = g / tvx(i,k) + gotvx(i,k) = g / thvx(i,k) + enddo + enddo + else + 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) @@ -489,7 +559,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! gotvx(i,k) = g / tvx(i,k) gotvx(i,k) = g / thvx(i,k) enddo - enddo + enddo + endif !sa3dtke !> - Compute an empirical cloud fraction based on !! Xu and Randall (1996) \cite xu_and_randall_1996 @@ -1115,7 +1186,87 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> ## Compute an asymtotic mixing length ! - do k = 1, km1 + if(sa3dtke) then + do k = 1, km1 + do i = 1, im + zlup = 0.0 + mlenflg = .true. + e2(i,k) = max(2.*tke(i,k), 0.001) + do n = k, km1 + if(mlenflg) then + dz = zl(i,n+1) - zl(i,n) + tem1 = 2.*gotvx(i,n+1)*(thvx(i,k)-thvx(i,n+1)) + tem2 = cs0*sqrt(e2(i,n))*sqrt(def_1(i,n)+def_2(i,n)) + e2(i,n+1) = e2(i,n) + (tem1 - tem2) * dz + zlup = zlup + dz + if(e2(i,n+1) < 0.) then + ptem = e2(i,n+1) / (e2(i,n+1) - e2(i,n)) + zlup = zlup - ptem * dz + zlup = max(zlup, 0.) + mlenflg = .false. + endif + endif + enddo + zldn = 0.0 + mlenflg = .true. + do n = k, 1, -1 + if(mlenflg) then + if(n == 1) then + dz = zl(i,1) + tem = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + tem1 = 2.*gotvx(i,n)*(tem-thvx(i,k)) + tem2 = ustar(i)*phims(i)/(vk*dz) + tem2 = cs0*sqrt(e2(i,n))*tem2 + e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz + else + dz = zl(i,n) - zl(i,n-1) + tem1 = 2.*gotvx(i,n-1)*(thvx(i,n-1)-thvx(i,k)) + tem2 = cs0*sqrt(e2(i,n))* + & sqrt(def_1(i,n-1)+def_2(i,n-1)) + e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz + endif + zldn = zldn + dz + if(e2(i,n-1) < 0.) then + ptem = e2(i,n-1) / (e2(i,n-1) - e2(i,n)) + zldn = zldn - ptem * dz + zldn = max(zldn, 0.) + mlenflg = .false. + endif + endif + enddo +! + tem = 0.5 * (zi(i,k+1)-zi(i,k)) + tem1 = min(tem, rlmnz(i,k)) +!> - Following Bougeault and Lacarrere(1989), the characteristic length +!! scale (\f$l_2\f$) (eqn 10 in Han et al.(2019) \cite Han_2019) is given by: +!!\f[ +!! l_2=min(l_{up},l_{down}) +!!\f] +!! and dissipation length scale \f$l_d\f$ is given by: +!!\f[ +!! l_d=(l_{up}l_{down})^{1/2} +!!\f] +!! where \f$l_{up}\f$ and \f$l_{down}\f$ are the distances that a parcel +!! having an initial TKE can travel upward and downward before being stopped +!! by buoyancy effects. +! +! Following Rodier et. al (2017), environmental wind shear effect on +! mixing length was included. +! + ptem2 = min(zlup,zldn) + rlam(i,k) = elmfac * ptem2 + rlam(i,k) = max(rlam(i,k), tem1) + rlam(i,k) = min(rlam(i,k), rlmx) +! + ptem2 = sqrt(zlup*zldn) + ele(i,k) = elefac * ptem2 + ele(i,k) = max(ele(i,k), tem1) + ele(i,k) = min(ele(i,k), elmx) +! + enddo + enddo + else + do k = 1, km1 do i = 1, im zlup = 0.0 mlenflg = .true. @@ -1170,7 +1321,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !! l_2=min(l_{up},l_{down}) !!\f] !! and dissipation length scale \f$l_d\f$ is given by: -!!\f[ +!!\f[ !! l_d=(l_{up}l_{down})^{1/2} !!\f] !! where \f$l_{up}\f$ and \f$l_{down}\f$ are the distances that a parcel @@ -1191,7 +1342,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ele(i,k) = min(ele(i,k), elmx) ! enddo - enddo + enddo + endif !sa3dtke + !> - Compute the surface layer length scale (\f$l_1\f$) following !! Nakanishi (2001) \cite Nakanish_2001 (eqn 9 of Han et al.(2019) \cite Han_2019) do k = 1, km1 @@ -1244,7 +1397,59 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !> ## Compute eddy diffusivities !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - do k = 1, km1 + if(sa3dtke) then + do k = 1, km1 + do i = 1, im + tem = 0.5 * (elm(i,k) + elm(i,k+1)) + tem = tem * sqrt(tkeh(i,k)) + ri = max(bf(i,k)/(def_1(i,k)+def_2(i,k)),rimin) + if(k < kpbl(i)) then + if(pcnvflg(i)) then + dku(i,k) = ckz(i,k) * tem + dkt(i,k) = dku(i,k) / prn(i,k) + else + if(ri < 0.) then ! unstable regime + dku(i,k) = ckz(i,k) * tem + dkt(i,k) = dku(i,k) / prn(i,k) + else ! stable regime + dkt(i,k) = chz(i,k) * tem + dku(i,k) = dkt(i,k) * prn(i,k) + endif + endif + else + if(ri < 0.) then ! unstable regime + dku(i,k) = ck1 * tem + dkt(i,k) = rchck * dku(i,k) + else ! stable regime + dkt(i,k) = ch1 * tem + prnum = 1.0 + 2.1 * ri + prnum = min(prnum,prmax) + dku(i,k) = dkt(i,k) * prnum + endif + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + tem1 = ckz(i,k) * tem + ptem1 = tem1 / prscu + dku(i,k) = max(dku(i,k), tem1) + dkt(i,k) = max(dkt(i,k), ptem1) + endif + endif +! + dkq(i,k) = prtke * dkt(i,k) +! + dkt(i,k) = min(dkt(i,k),dkmax) + dkt(i,k) = max(dkt(i,k),xkzo(i,k)) + dkq(i,k) = min(dkq(i,k),dkmax) + dkq(i,k) = max(dkq(i,k),xkzo(i,k)) + dku(i,k) = min(dku(i,k),dkmax) + dku(i,k) = max(dku(i,k),xkzmo(i,k)) +! + enddo + enddo + else + do k = 1, km1 do i = 1, im tem = 0.5 * (elm(i,k) + elm(i,k+1)) tem = tem * sqrt(tkeh(i,k)) @@ -1293,7 +1498,131 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & dku(i,k) = max(dku(i,k),xkzmo(i,k)) ! enddo + enddo + endif !sa3dtke + +!The following is for SA-3D-TKE + if(sa3dtke) then +! 1. compute LES component of km, kh, and kq (Deardorff 1980) +! 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)+q1(i,k,6),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 + dku_les(i,k) = 0. + dkt_les(i,k) = 0. + dkq_les(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 + dku_les(i,k)=0.1*elm_les(i,k)*sqrt(tkeh(i,k)) + dkt_les(i,k)=(1.0+2.0*elm_les(i,k)/tem1)*dku_les(i,k) + dkq_les(i,k)=dkt_les(i,k) + dku_les(i,k) = min(dku_les(i,k),kmaxles) + dkt_les(i,k) = min(dkt_les(i,k),kmaxles) + dkq_les(i,k) = min(dkq_les(i,k),kmaxles) + enddo + enddo + +! 2. compute MS horizontal km + do i = 1, im + do k = 1, km1 + tem1=sqrt(garea(i)) + dku_h(i,k)=dddmp*tem1*sqrt(tkeh(i,k)) + enddo + dku_h(i,km)=dku_h(i,km1) + enddo + +! calculate blending coefficients for km, kt, kq, and nonlocal mixing + do i = 1, im +! finding scale of large eddies from TKE +! d/zi + scl=1000. + l_tkemax=10 + kscl=10 + tkemax=0.0 + do k=1,km1 + tkemax=max(tkemax,tkeh(i,k)) + enddo + do k=1,km1 + if (abs(tkeh(i,k)-tkemax)/tkemax .lt. 1.0e-9) then + l_tkemax=k + endif + enddo + do k=l_tkemax,km1 + if (tkeh(i,k)-0.5*tkemax .gt. 0.0) then + kscl=k + endif + enddo + kscl=min(kscl,km1-10) + scl=zi(i,kscl+1) + scl=max(scl,esmax) + 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) +! partition function for TKE + 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 LES and MS components of km,kt, and kq from in the +! vertical and horizontal direction + do k = 1,km + do i=1,im + dkq(i,k)=(1.0-pfl(i))*dkq_les(i,k)+pfl(i)*dkq(i,k) + dkt(i,k)=(1.0-pfl(i))*dkt_les(i,k)+pfl(i)*dkt(i,k) + dku(i,k)=(1.0-pfl(i))*dku_les(i,k)+pfl(i)*dku(i,k) + ri = max(bf(i,k)/(def_1(i,k)+def_2(i,k)),rimin) + if(ri < 0.) then ! unstable regime + prnum=1./rchck + else + prnum=1.0 + 2.1 * ri + endif + dkq_h(i,k)=(1.0-pfl(i))*dkq_les(i,k)+ + & pfl(i)*dku_h(i,k)/prnum + dku_h(i,k)=(1.0-pfl(i))*dku_les(i,k)+pfl(i)*dku_h(i,k) + enddo enddo + + endif !sa3dtke + !> ## Compute TKE. !! - Compute a minimum TKE deduced from background diffusivity for momentum. ! @@ -1317,7 +1646,91 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !> - Compute buoyancy and shear productions of TKE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! - do k = 1, km1 + if(sa3dtke) then + do k = 1, km1 + do i = 1, im + if (k == 1) then + tem = -dkt(i,1) * bf(i,1) +! if(pcnvflg(i)) then +! ptem1 = xmf(i,1) * buou(i,1) +! else + ptem1 = 0. +! endif + if(scuflg(i) .and. mrad(i) == 1) then + ptem2 = xmfd(i,1) * buod(i,1) + else + ptem2 = 0. + endif + tem = tem + ptem1 + ptem2 + buop = 0.5 * (gotvx(i,1) * sflux(i) + tem) +! + tem = (u1(i,2)-u1(i,1))*rdzt(i,1) +! if(pcnvflg(i)) then +! ptem = xmf(i,1) * tem +! ptem1 = 0.5 * ptem * (u1(i,2)-ucko(i,2)) +! else + ptem1 = 0. +! endif + if(scuflg(i) .and. mrad(i) == 1) then + ptem = ucdo(i,1)+ucdo(i,2)-u1(i,1)-u1(i,2) + ptem = 0.5 * tem * xmfd(i,1) * ptem + else + ptem = 0. + endif + ptem1 = ptem1 + ptem +! + tem = (v1(i,2)-v1(i,1))*rdzt(i,1) +! if(pcnvflg(i)) then +! ptem = xmf(i,1) * tem +! ptem2 = 0.5 * ptem * (v1(i,2)-vcko(i,2)) +! else + ptem2 = 0. +! endif + if(scuflg(i) .and. mrad(i) == 1) then + ptem = vcdo(i,1)+vcdo(i,2)-v1(i,1)-v1(i,2) + ptem = 0.5 * tem * xmfd(i,1) * ptem + else + ptem = 0. + endif + ptem2 = ptem2 + ptem +! + tem2 = stress(i)*ustar(i)*phims(i)/(vk*zl(i,1)) + shrp = 0.5 * (ptem1 + ptem2 + tem2) + tem2 = dku(i,k)*def_1(i,k)+dku_h(i,k)*def_2(i,k) + shrp =0.5*(shrp+tem2) + else + tem1 = -dkt(i,k-1) * bf(i,k-1) + tem2 = -dkt(i,k) * bf(i,k) + tem = 0.5 * (tem1 + tem2) + if(pcnvflg(i) .and. k <= kpbl(i)) then + ptem = 0.5 * (xmf(i,k-1) + xmf(i,k)) + ptem1 = ptem * buou(i,k) + else + ptem1 = 0. + endif + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem0 = 0.5 * (xmfd(i,k-1) + xmfd(i,k)) + ptem2 = ptem0 * buod(i,k) + else + ptem2 = 0. + endif + else + ptem2 = 0. + endif + buop = tem + ptem1 + ptem2 +! +! obtaining 3d shear production from dycore + tem1 = dku(i,k-1)*def_1(i,k-1)+dku_h(i,k-1)*def_2(i,k-1) + tem2 = dku(i,k)*def_1(i,k)+dku_h(i,k)*def_2(i,k) + shrp=0.5*(tem1+tem2) + endif +! + prod(i,k) = buop + shrp + enddo + enddo + else + do k = 1, km1 do i = 1, im if (k == 1) then tem = -dkt(i,1) * bf(i,1) @@ -1434,14 +1847,40 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif prod(i,k) = buop + shrp enddo - enddo + enddo + endif !sa3dtke ! !---------------------------------------------------------------------- !> - First predict tke due to tke production & dissipation(diss) ! - dtn = dt2 / float(ndt) - do n = 1, ndt - do k = 1,km1 + if(sa3dtke) then +!The following is for SA-3D-TKE + dtn = dt2 / float(ndt) + do n = 1, ndt + do k = 1,km1 + do i=1,im + tem = sqrt(tke(i,k)) +! 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*def_3(i,k) + tem=min(tem,1.0) + tke(i,k) = tke(i,k) + dtn * (prod(i,k)-diss(i,k)+tem) +! tke(i,k) = max(tke(i,k), tkmin) + tke(i,k) = max(tke(i,k), tkmnz(i,k)) + enddo + enddo + enddo + else + dtn = dt2 / float(ndt) + do n = 1, ndt + do k = 1,km1 do i=1,im tem = sqrt(tke(i,k)) ptem = ce0 / ele(i,k) @@ -1452,8 +1891,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! tke(i,k) = max(tke(i,k), tkmin) tke(i,k) = max(tke(i,k), tkmnz(i,k)) enddo - enddo - enddo + enddo + enddo + endif !sa3dtke ! !> - Compute updraft & downdraft properties for TKE ! @@ -1620,7 +2060,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & f1(i,1) = tke(i,1) enddo ! - do k = 1,km1 +!The following is for SA-3D-TKE + if(sa3dtke) then + do k = 1,km1 do i=1,im dtodsd = dt2/del(i,k) dtodsu = dt2/del(i,k+1) @@ -1636,6 +2078,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! if(pcnvflg(i) .and. k < kpbl(i)) then ptem = 0.5 * tem2 * xmf(i,k) + ptem=ptem*pfnl(i) ! blending nolocal mixing ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke) @@ -1648,6 +2091,7 @@ 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) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke) @@ -1659,6 +2103,7 @@ 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)) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem f1(i,k) = f1(i,k) + e_half(i,k) * ptem1 @@ -1666,7 +2111,57 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif ! enddo - enddo + enddo + else + do k = 1,km1 + do i=1,im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + rdz = rdzt(i,k) + tem1 = dsig * dkq(i,k) * rdz + dsdz2 = tem1 * rdz + au(i,k) = -dtodsd*dsdz2 + al(i,k) = -dtodsu*dsdz2 + ad(i,k) = ad(i,k)-au(i,k) + ad(i,k+1)= 1.-al(i,k) + tem2 = dsig * rdz +! + if(pcnvflg(i) .and. k < kpbl(i)) then + ptem = 0.5 * tem2 * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke) + f1(i,k) = f1(i,k) - ptem * ptem1 + f1(i,k+1) = tke(i,k+1) + ptem * ptem2 + else + f1(i,k+1) = tke(i,k+1) + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem = 0.5 * tem2 * xmfd(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke) + f1(i,k) = f1(i,k) + ptem * ptem1 + f1(i,k+1) = f1(i,k+1) - ptem * ptem2 + endif + endif +! + kmx = max(kpbl(i), krad(i)) + if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then + ptem = tem2 * (xmf(i,k) - xmfd(i,k)) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + f1(i,k) = f1(i,k) + e_half(i,k) * ptem1 + f1(i,k+1) = f1(i,k+1) - e_half(i,k) * ptem2 + endif +! + enddo + enddo + endif !sa3dtke + c !> - Call tridit() to solve tridiagonal problem for TKE c @@ -1810,7 +2305,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo endif c - do k = 1,km1 +!The following is for SA-3D-TKE + if(sa3dtke) then + do k = 1,km1 do i = 1,im dtodsd = dt2/del(i,k) dtodsu = dt2/del(i,k+1) @@ -1827,6 +2324,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! if(pcnvflg(i) .and. k < kpbl(i)) then ptem = 0.5 * tem2 * xmf(i,k) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem tem = t1(i,k) + t1(i,k+1) @@ -1845,6 +2343,7 @@ 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) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem ptem = tcdo(i,k) + tcdo(i,k+1) @@ -1860,6 +2359,7 @@ 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)) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem f2(i,k) = f2(i,k) + q_half(i,k,1) * ptem1 @@ -1867,9 +2367,71 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif ! enddo - enddo + enddo + else + do k = 1,km1 + do i = 1,im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + rdz = rdzt(i,k) + tem1 = dsig * dkt(i,k) * rdz + dsdzt = tem1 * gocp + dsdz2 = tem1 * rdz + au(i,k) = -dtodsd*dsdz2 + al(i,k) = -dtodsu*dsdz2 + ad(i,k) = ad(i,k)-au(i,k) + ad(i,k+1)= 1.-al(i,k) + tem2 = dsig * rdz ! - if(ntrac1 >= 2) then + if(pcnvflg(i) .and. k < kpbl(i)) then + ptem = 0.5 * tem2 * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem = t1(i,k) + t1(i,k+1) + ptem = tcko(i,k) + tcko(i,k+1) + f1(i,k) = f1(i,k)+dtodsd*dsdzt-(ptem-tem)*ptem1 + f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+(ptem-tem)*ptem2 + ptem = qcko(i,k,1) + qcko(i,k+1,1) + f2(i,k) = f2(i,k) - ptem * ptem1 + f2(i,k+1) = q1(i,k+1,1) + ptem * ptem2 + else + f1(i,k) = f1(i,k)+dtodsd*dsdzt + f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt + f2(i,k+1) = q1(i,k+1,1) + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem = 0.5 * tem2 * xmfd(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + ptem = tcdo(i,k) + tcdo(i,k+1) + tem = t1(i,k) + t1(i,k+1) + f1(i,k) = f1(i,k) + (ptem - tem) * ptem1 + f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2 + ptem = qcdo(i,k,1) + qcdo(i,k+1,1) + f2(i,k) = f2(i,k) + ptem * ptem1 + f2(i,k+1) = f2(i,k+1) - ptem * ptem2 + endif + endif +! + kmx = max(kpbl(i), krad(i)) + if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then + ptem = tem2 * (xmf(i,k) - xmfd(i,k)) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + f2(i,k) = f2(i,k) + q_half(i,k,1) * ptem1 + f2(i,k+1) = f2(i,k+1) - q_half(i,k,1) * ptem2 + endif +! + enddo + enddo + endif !sa3dtke +! +!The following is for SA-3D-TKE + if(sa3dtke) then + if(ntrac1 >= 2) then do n = 2, ntrac1 is = (n-1) * km do k = 1, km1 @@ -1881,6 +2443,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! if(pcnvflg(i) .and. k < kpbl(i)) then ptem = 0.5 * tem2 * xmf(i,k) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem ptem = qcko(i,k,n) + qcko(i,k+1,n) @@ -1893,6 +2456,7 @@ 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) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem ptem = qcdo(i,k,n) + qcdo(i,k+1,n) @@ -1904,6 +2468,7 @@ 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)) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem f2(i,k+is) = f2(i,k+is) + q_half(i,k,n) * ptem1 @@ -1913,7 +2478,54 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo enddo enddo - endif + endif + else + if(ntrac1 >= 2) then + do n = 2, ntrac1 + is = (n-1) * km + do k = 1, km1 + do i = 1, im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + tem2 = dsig * rdzt(i,k) +! + if(pcnvflg(i) .and. k < kpbl(i)) then + ptem = 0.5 * tem2 * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + ptem = qcko(i,k,n) + qcko(i,k+1,n) + f2(i,k+is) = f2(i,k+is) - ptem * ptem1 + f2(i,k+1+is)= q1(i,k+1,n) + ptem * ptem2 + else + f2(i,k+1+is) = q1(i,k+1,n) + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem = 0.5 * tem2 * xmfd(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + ptem = qcdo(i,k,n) + qcdo(i,k+1,n) + f2(i,k+is) = f2(i,k+is) + ptem * ptem1 + f2(i,k+1+is)= f2(i,k+1+is) - ptem * ptem2 + endif + endif +! + kmx = max(kpbl(i), krad(i)) + if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then + ptem = tem2 * (xmf(i,k) - xmfd(i,k)) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + f2(i,k+is) = f2(i,k+is) + q_half(i,k,n) * ptem1 + f2(i,k+1+is) = f2(i,k+1+is) - q_half(i,k,n) * ptem2 + endif +! + enddo + enddo + enddo + endif + endif !sa3dtke c !> - Call tridin() to solve tridiagonal problem for heat and moisture c @@ -2314,7 +2926,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & f2(i,1) = v1(i,1) enddo c - do k = 1,km1 +!The following is for SA-3D-TKE + if(sa3dtke) then + do k = 1,km1 do i=1,im dtodsd = dt2/del(i,k) dtodsu = dt2/del(i,k+1) @@ -2330,6 +2944,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! if(pcnvflg(i) .and. k < kpbl(i)) then ptem = 0.5 * tem2 * xmf(i,k) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem tem = u1(i,k) + u1(i,k+1) @@ -2348,6 +2963,7 @@ 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) + ptem=ptem*pfnl(i) ptem1 = dtodsd * ptem ptem2 = dtodsu * ptem tem = u1(i,k) + u1(i,k+1) @@ -2362,7 +2978,58 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif ! enddo - enddo + enddo + else + do k = 1,km1 + do i=1,im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + rdz = rdzt(i,k) + tem1 = dsig * dku(i,k) * rdz + dsdz2 = tem1*rdz + au(i,k) = -dtodsd*dsdz2 + al(i,k) = -dtodsu*dsdz2 + ad(i,k) = ad(i,k)-au(i,k) + ad(i,k+1)= 1.-al(i,k) + tem2 = dsig * rdz +! + if(pcnvflg(i) .and. k < kpbl(i)) then + ptem = 0.5 * tem2 * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem = u1(i,k) + u1(i,k+1) + ptem = ucko(i,k) + ucko(i,k+1) + f1(i,k) = f1(i,k) - (ptem - tem) * ptem1 + f1(i,k+1) = u1(i,k+1) + (ptem - tem) * ptem2 + tem = v1(i,k) + v1(i,k+1) + ptem = vcko(i,k) + vcko(i,k+1) + f2(i,k) = f2(i,k) - (ptem - tem) * ptem1 + f2(i,k+1) = v1(i,k+1) + (ptem - tem) * ptem2 + else + f1(i,k+1) = u1(i,k+1) + f2(i,k+1) = v1(i,k+1) + endif +! + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem = 0.5 * tem2 * xmfd(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem = u1(i,k) + u1(i,k+1) + ptem = ucdo(i,k) + ucdo(i,k+1) + f1(i,k) = f1(i,k) + (ptem - tem) *ptem1 + f1(i,k+1) = f1(i,k+1) - (ptem - tem) *ptem2 + tem = v1(i,k) + v1(i,k+1) + ptem = vcdo(i,k) + vcdo(i,k+1) + f2(i,k) = f2(i,k) + (ptem - tem) * ptem1 + f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2 + endif + endif +! + enddo + enddo + endif !sa3dtke c !> - Call tridi2() to solve tridiagonal problem for momentum c @@ -2416,10 +3083,21 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> ## Save PBL height for diagnostic purpose ! - do i = 1, im + if(sa3dtke) then + do i = 1, im hpbl(i) = hpblx(i) kpbl(i) = kpblx(i) - enddo + do k=1, km + dku3d_h(i,k) = dku_h(i,k) ! pass dku3d_h to dyn_core + dku3d_e(i,k) = dkq_h(i,k) ! pass dku3d_e to dyn_core + enddo + enddo + else + do i = 1, im + hpbl(i) = hpblx(i) + kpbl(i) = kpblx(i) + enddo + endif !sa3dtke ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! return diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.meta b/physics/PBL/SATMEDMF/satmedmfvdifq.meta index 13c66399e..830452093 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.meta +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.meta @@ -256,6 +256,30 @@ type = real kind = kind_phys intent = in +[def_1] + standard_name = vert_shear_square + long_name = square of vertical shear + units = m2 s-2 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[def_2] + standard_name = hori_shear_square + long_name = square of horizontal shear + units = m2 s-2 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[def_3] + standard_name = TKE_transfer_rate + long_name = rate of TKE transfer and pressure correlation + units = m2 s-3 + 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 @@ -470,6 +494,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 @@ -534,6 +565,22 @@ type = real kind = kind_phys intent = out +[dku3d_h] + standard_name = horizontal_atmosphere_momentum_diffusivity_dyncore + long_name = horizontal atmospheric momentum diffusivity + units = m2 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out +[dku3d_e] + standard_name = atmosphere_momentum_diffusivity_tke_dyncore + long_name = atmospheric momentum diffusivity for tke + units = m2 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [kinver] standard_name = index_of_highest_temperature_inversion long_name = index of highest temperature inversion From dfa39f6baca358b6348e247a63455c64beb543a3 Mon Sep 17 00:00:00 2001 From: "Bin.Liu" Date: Tue, 4 Feb 2025 14:40:55 -0600 Subject: [PATCH 2/6] Add an option to use liquid potential temperature in temperature equation of GFS PBL scheme. (from @WeiguoWang-NOAA) --- physics/PBL/SATMEDMF/satmedmfvdifq.F | 14 +++++++++++++- physics/PBL/SATMEDMF/satmedmfvdifq.meta | 7 +++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.F b/physics/PBL/SATMEDMF/satmedmfvdifq.F index 64e7bc32d..4de3bfc73 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.F +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.F @@ -82,7 +82,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & prsi,del,prsl,prslk,phii,phil,delt, & & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku,tkeh, & & kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr, & - & rlmx,elmx,sfc_rlm,tc_pbl, & + & rlmx,elmx,sfc_rlm,tc_pbl,use_lpt, & & ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind, & & index_of_y_wind,index_of_process_pbl,gen_tend,ldiag3d, & & errmsg,errflg) @@ -97,6 +97,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & ntke, ntqv integer, intent(in) :: sfc_rlm integer, intent(in) :: tc_pbl + integer, intent(in) :: use_lpt integer, intent(in) :: kinver(:) integer, intent(out) :: kpbl(:) logical, intent(in) :: gen_tend,ldiag3d @@ -259,6 +260,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & real(kind=kind_phys) h1 real(kind=kind_phys) bfac, mffac + + real(kind=kind_phys) qice(im,km),qliq(im,km) !! parameter(bfac=100.) parameter(wfac=7.0,cfac=4.5) @@ -467,15 +470,20 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do i=1,im pix(i,k) = psk(i) / prslk(i,k) theta(i,k) = t1(i,k) * pix(i,k) + qice(i,k) = 0.0 + qliq(i,k) = 0.0 if(ntiw > 0) then tem = max(q1(i,k,ntcw),qlmin) tem1 = max(q1(i,k,ntiw),qlmin) + qice(i,k) = tem1 + qliq(i,k) = tem qlx(i,k) = tem + tem1 ptem = hvap*tem + (hvap+hfus)*tem1 slx(i,k) = cp * t1(i,k) + phil(i,k) - ptem else qlx(i,k) = max(q1(i,k,ntcw),qlmin) slx(i,k) = cp * t1(i,k) + phil(i,k) - hvap*qlx(i,k) + qliq(i,k) = qlx(i,k) endif tem2 = 1.+fv*max(q1(i,k,1),qmin)-qlx(i,k) thvx(i,k) = theta(i,k) * tem2 @@ -1818,6 +1826,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & rdz = rdzt(i,k) tem1 = dsig * dkt(i,k) * rdz dsdzt = tem1 * gocp + if (use_lpt > 0) then + dsdzt = dsdzt-tem1*elocp*(qliq(i,k+1)-qliq(i,k))*rdz + & -(1+0.33/2.5)*tem1*elocp*(qice(i,k+1)-qice(i,k))*rdz + endif dsdz2 = tem1 * rdz au(i,k) = -dtodsd*dsdz2 al(i,k) = -dtodsu*dsdz2 diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.meta b/physics/PBL/SATMEDMF/satmedmfvdifq.meta index 13c66399e..2c857be07 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.meta +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.meta @@ -619,6 +619,13 @@ dimensions = () type = integer intent = in +[use_lpt] + standard_name = control_for_using_LPT_for_TC_applications_in_the_PBL_scheme + long_name = control for using LPT in 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) From 9f259053ac220e7ebce9d305ca832c06bbc6ba4d Mon Sep 17 00:00:00 2001 From: "Kwun Y. Fung" Date: Mon, 2 Jun 2025 23:55:02 -0500 Subject: [PATCH 3/6] Update sa3dtke variables standard names in SATMEDMF --- physics/PBL/SATMEDMF/satmedmfvdifq.F | 8 +++++--- physics/PBL/SATMEDMF/satmedmfvdifq.meta | 22 +++++++++++----------- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.F b/physics/PBL/SATMEDMF/satmedmfvdifq.F index ef71bae63..523bdf424 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.F +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.F @@ -267,7 +267,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & dku_h(im,km),dkq_h(im,km),dddmp, & cpl1,cpl2,cpl3,cpl4, & cpl5,cpl6,pfnl(im),cpnl1,cpnl2,cpnl3,cpnl4, - & cpnl5,cpnl6,cptke1,cptke2,cptke3 + & cpnl5,cpnl6,cptke1,cptke2,cptke3, + & inv3 real(kind=kind_phys) ak(im,km),bk(im,km),ck(im,km),pk(im,km), & f3(im,km),rizt(im,km) @@ -318,6 +319,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & 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,dddmp=0.1) + parameter(inv3=0.33333333) ! parameter(aa1=0.92,aa2=0.649,bb1=17.7,bb2=9.5,cc1=0.08) if (tc_pbl == 0) then @@ -1533,7 +1535,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & 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) + tem1=(garea(i)*dz)**(inv3) ! calculate LES mixing length if(tem > 0.0) then elm_les(i,k)=0.76*sqrt(tkeh(i,k))*tem**(-0.5) @@ -1862,7 +1864,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & tem = sqrt(tke(i,k)) ! calculating 3D TKE transport and pressure correlation ptem1 = ce0 / ele(i,k) - tem1=(garea(i)*dz)**(1.0/3.0) + tem1=(garea(i)*dz)**(inv3) 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 diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.meta b/physics/PBL/SATMEDMF/satmedmfvdifq.meta index 830452093..25b7817b8 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.meta +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.meta @@ -257,24 +257,24 @@ kind = kind_phys intent = in [def_1] - standard_name = vert_shear_square - long_name = square of vertical shear + standard_name = square_of_vertical_shear_due_to_dynamics + long_name = square of vertical shear calculated from dynamics units = m2 s-2 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in [def_2] - standard_name = hori_shear_square - long_name = square of horizontal shear + standard_name = square_of_horizontal_shear_due_to_dynamics + long_name = square of horizontal shear calculated from dynamics units = m2 s-2 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in [def_3] - standard_name = TKE_transfer_rate - long_name = rate of TKE transfer and pressure correlation + standard_name = transfer_rate_of_tke_due_to_dynamics + long_name = rate of TKE transfer and pressure correlation calculated from dynamics units = m2 s-3 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -495,7 +495,7 @@ type = logical intent = in [sa3dtke] - standard_name = flag_sa_3d_tke_scheme + standard_name = do_scale_aware_3d_tke long_name = flag for scale-aware 3d tke scheme units = flag dimensions = () @@ -566,16 +566,16 @@ kind = kind_phys intent = out [dku3d_h] - standard_name = horizontal_atmosphere_momentum_diffusivity_dyncore - long_name = horizontal atmospheric momentum diffusivity + standard_name = horizontal_atmosphere_momentum_diffusivity_for_dyanmics + long_name = horizontal atmospheric momentum diffusivity for dynamics units = m2 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = out [dku3d_e] - standard_name = atmosphere_momentum_diffusivity_tke_dyncore - long_name = atmospheric momentum diffusivity for tke + standard_name = horizontal_atmosphere_tke_diffusivity_for_dyanmics + long_name = horizontal atmospheric tke diffusivity for dynamics units = m2 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real From a6c6bbdc829937721da9bef82f84b2162bdfb581 Mon Sep 17 00:00:00 2001 From: "Kwun Y. Fung" Date: Wed, 4 Jun 2025 16:58:37 -0500 Subject: [PATCH 4/6] Update nonlocal mixing and coding style for sa3dkte in ccpp physics --- physics/PBL/SATMEDMF/mfscuq.f | 15 +++- physics/PBL/SATMEDMF/satmedmfvdifq.F | 128 ++++++--------------------- physics/PBL/mfpbltq.f | 15 +++- 3 files changed, 54 insertions(+), 104 deletions(-) diff --git a/physics/PBL/SATMEDMF/mfscuq.f b/physics/PBL/SATMEDMF/mfscuq.f index b2a48d1b6..e9be00026 100644 --- a/physics/PBL/SATMEDMF/mfscuq.f +++ b/physics/PBL/SATMEDMF/mfscuq.f @@ -15,7 +15,9 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, & cnvflg,zl,zm,q1,t1,u1,v1,plyr,pix, & thlx,thvx,thlvx,gdx,thetae, & krad,mrad,radmin,buo,wush,tkemean,vez0fun,xmfd, - & tcdo,qcdo,ucdo,vcdo,xlamdeq,a1) + & tcdo,qcdo,ucdo,vcdo,xlamdeq,a1, +!The following flag is for SA-3D-TKE (kyf) + & sa3dtke) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -31,6 +33,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, integer krad(im), mrad(im) ! logical cnvflg(im) + logical sa3dtke !flag for SA-3D-TKE scheme (kyf) real(kind=kind_phys) delt real(kind=kind_phys) q1(ix,km,ntrac1),t1(ix,km), & u1(ix,km), v1(ix,km), @@ -427,6 +430,16 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, endif enddo ! +!> - Set updraft fraction to 0 when using SA-3D-TKE scheme (kyf) +!! Scale-aware capability is done with pfnl in satmedmfvdifq.F +!! Zhu et al. (2025) +! + if (sa3dtke) then + do i = 1, im + sigma(i) = 0. + enddo + endif +! !> - Compute scale-aware function based on !! Arakawa and Wu (2013) \cite arakawa_and_wu_2013 ! diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.F b/physics/PBL/SATMEDMF/satmedmfvdifq.F index 523bdf424..311913bab 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.F +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.F @@ -1089,14 +1089,18 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & call mfpbltq(im,im,km,kmpbl,ntcw,ntrac1,dt2, & pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx, & gdx,hpbl,kpbl,vpert,buou,wush,tkemean,vez0fun,xmf, - & tcko,qcko,ucko,vcko,xlamue,bl_upfr) + & tcko,qcko,ucko,vcko,xlamue,bl_upfr, +!The following flag is for SA-3D-TKE (kyf) + & sa3dtke) !> - Call mfscuq(), which is a new mass-flux parameterization for !! stratocumulus-top-induced turbulence mixing. For details of the mfscuq subroutine, step into its documentation ::mfscuq call mfscuq(im,im,km,kmscu,ntcw,ntrac1,dt2, & scuflg,zl,zm,q1,t1,u1,v1,plyr,pix, & thlx,thvx,thlvx,gdx,thetae, & krad,mrad,radmin,buod,wush,tkemean,vez0fun,xmfd, - & tcdo,qcdo,ucdo,vcdo,xlamde,bl_dnfr) + & tcdo,qcdo,ucdo,vcdo,xlamde,bl_dnfr, +!The following flag is for SA-3D-TKE (kyf) + & sa3dtke) if (tc_pbl == 1) then !> - unify mass fluxes with Cu @@ -1188,87 +1192,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> ## Compute an asymtotic mixing length ! - if(sa3dtke) then - do k = 1, km1 - do i = 1, im - zlup = 0.0 - mlenflg = .true. - e2(i,k) = max(2.*tke(i,k), 0.001) - do n = k, km1 - if(mlenflg) then - dz = zl(i,n+1) - zl(i,n) - tem1 = 2.*gotvx(i,n+1)*(thvx(i,k)-thvx(i,n+1)) - tem2 = cs0*sqrt(e2(i,n))*sqrt(def_1(i,n)+def_2(i,n)) - e2(i,n+1) = e2(i,n) + (tem1 - tem2) * dz - zlup = zlup + dz - if(e2(i,n+1) < 0.) then - ptem = e2(i,n+1) / (e2(i,n+1) - e2(i,n)) - zlup = zlup - ptem * dz - zlup = max(zlup, 0.) - mlenflg = .false. - endif - endif - enddo - zldn = 0.0 - mlenflg = .true. - do n = k, 1, -1 - if(mlenflg) then - if(n == 1) then - dz = zl(i,1) - tem = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) - tem1 = 2.*gotvx(i,n)*(tem-thvx(i,k)) - tem2 = ustar(i)*phims(i)/(vk*dz) - tem2 = cs0*sqrt(e2(i,n))*tem2 - e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz - else - dz = zl(i,n) - zl(i,n-1) - tem1 = 2.*gotvx(i,n-1)*(thvx(i,n-1)-thvx(i,k)) - tem2 = cs0*sqrt(e2(i,n))* - & sqrt(def_1(i,n-1)+def_2(i,n-1)) - e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz - endif - zldn = zldn + dz - if(e2(i,n-1) < 0.) then - ptem = e2(i,n-1) / (e2(i,n-1) - e2(i,n)) - zldn = zldn - ptem * dz - zldn = max(zldn, 0.) - mlenflg = .false. - endif - endif - enddo -! - tem = 0.5 * (zi(i,k+1)-zi(i,k)) - tem1 = min(tem, rlmnz(i,k)) -!> - Following Bougeault and Lacarrere(1989), the characteristic length -!! scale (\f$l_2\f$) (eqn 10 in Han et al.(2019) \cite Han_2019) is given by: -!!\f[ -!! l_2=min(l_{up},l_{down}) -!!\f] -!! and dissipation length scale \f$l_d\f$ is given by: -!!\f[ -!! l_d=(l_{up}l_{down})^{1/2} -!!\f] -!! where \f$l_{up}\f$ and \f$l_{down}\f$ are the distances that a parcel -!! having an initial TKE can travel upward and downward before being stopped -!! by buoyancy effects. -! -! Following Rodier et. al (2017), environmental wind shear effect on -! mixing length was included. -! - ptem2 = min(zlup,zldn) - rlam(i,k) = elmfac * ptem2 - rlam(i,k) = max(rlam(i,k), tem1) - rlam(i,k) = min(rlam(i,k), rlmx) -! - ptem2 = sqrt(zlup*zldn) - ele(i,k) = elefac * ptem2 - ele(i,k) = max(ele(i,k), tem1) - ele(i,k) = min(ele(i,k), elmx) -! - enddo - enddo - else - do k = 1, km1 + do k = 1, km1 do i = 1, im zlup = 0.0 mlenflg = .true. @@ -1344,8 +1268,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ele(i,k) = min(ele(i,k), elmx) ! enddo - enddo - endif !sa3dtke + enddo !> - Compute the surface layer length scale (\f$l_1\f$) following !! Nakanishi (2001) \cite Nakanish_2001 (eqn 9 of Han et al.(2019) \cite Han_2019) @@ -1533,8 +1456,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & 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 + tem=gotvx(i,1)*(thetal(i,k+1)-thetal(i,k))*rdzt(i,k) tem1=(garea(i)*dz)**(inv3) ! calculate LES mixing length if(tem > 0.0) then @@ -1555,8 +1477,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo ! 2. compute MS horizontal km - do i = 1, im - do k = 1, km1 + do k = 1, km1 + do i = 1, im tem1=sqrt(garea(i)) dku_h(i,k)=dddmp*tem1*sqrt(tkeh(i,k)) enddo @@ -1575,12 +1497,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & tkemax=max(tkemax,tkeh(i,k)) enddo do k=1,km1 - if (abs(tkeh(i,k)-tkemax)/tkemax .lt. 1.0e-9) then + if (abs(tkeh(i,k)-tkemax)/tkemax < 1.0e-9) then l_tkemax=k endif enddo do k=l_tkemax,km1 - if (tkeh(i,k)-0.5*tkemax .gt. 0.0) then + if (tkeh(i,k)-0.5*tkemax > 0.0) then kscl=k endif enddo @@ -3086,19 +3008,21 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !> ## Save PBL height for diagnostic purpose ! if(sa3dtke) then - do i = 1, im - hpbl(i) = hpblx(i) - kpbl(i) = kpblx(i) - do k=1, km - dku3d_h(i,k) = dku_h(i,k) ! pass dku3d_h to dyn_core - dku3d_e(i,k) = dkq_h(i,k) ! pass dku3d_e to dyn_core + do i = 1, im + hpbl(i) = hpblx(i) + kpbl(i) = kpblx(i) + enddo + do k = 1, km + do i = 1, im + dku3d_h(i,k) = dku_h(i,k) ! pass dku3d_h to dyn_core + dku3d_e(i,k) = dkq_h(i,k) ! pass dku3d_e to dyn_core + enddo enddo - enddo else - do i = 1, im - hpbl(i) = hpblx(i) - kpbl(i) = kpblx(i) - enddo + do i = 1, im + hpbl(i) = hpblx(i) + kpbl(i) = kpblx(i) + enddo endif !sa3dtke ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/physics/PBL/mfpbltq.f b/physics/PBL/mfpbltq.f index 2812a7eab..52cca19ef 100644 --- a/physics/PBL/mfpbltq.f +++ b/physics/PBL/mfpbltq.f @@ -16,7 +16,9 @@ module mfpbltq_mod subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, & cnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx, & gdx,hpbl,kpbl,vpert,buo,wush,tkemean,vez0fun,xmf, - & tcko,qcko,ucko,vcko,xlamueq,a1) + & tcko,qcko,ucko,vcko,xlamueq,a1, +!The following flag is for SA-3D-TKE (kyf) + & sa3dtke) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -31,6 +33,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, ! &, me integer kpbl(im) logical cnvflg(im) + logical sa3dtke !flag for SA-3D-TKE scheme (kyf) real(kind=kind_phys) delt real(kind=kind_phys) q1(ix,km,ntrac1), & t1(ix,km), u1(ix,km), v1(ix,km), @@ -359,6 +362,16 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, endif enddo ! +!> - Set updraft fraction to 0 when using SA-3D-TKE scheme (kyf) +!! Scale-aware capability is done with pfnl in satmedmfvdifq.F +!! Zhu et al. (2025) +! + if (sa3dtke) then + do i = 1, im + sigma(i) = 0. + enddo + endif +! !> - Compute scale-aware function based on !! Arakawa and Wu (2013) \cite arakawa_and_wu_2013 ! From 0eb443ba0d48a273265f8b95dfe8f17d748ee01a Mon Sep 17 00:00:00 2001 From: "Kwun Y. Fung" Date: Mon, 16 Jun 2025 09:37:39 -0500 Subject: [PATCH 5/6] Updates on SA3DTKE non-local transport --- physics/PBL/SATMEDMF/satmedmfvdifq.F | 264 +++--------------------- physics/PBL/SATMEDMF/satmedmfvdifq.meta | 8 +- 2 files changed, 33 insertions(+), 239 deletions(-) diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.F b/physics/PBL/SATMEDMF/satmedmfvdifq.F index 311913bab..ece0256e2 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.F +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.F @@ -1456,6 +1456,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do k = 1, km1 do i = 1, im + dz = zl(i,k+1) - zl(i,k) tem=gotvx(i,1)*(thetal(i,k+1)-thetal(i,k))*rdzt(i,k) tem1=(garea(i)*dz)**(inv3) ! calculate LES mixing length @@ -1545,6 +1546,24 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo enddo +! perform scale-aware for non-local transport + do k = 1, kmpbl + do i = 1, im + if (pcnvflg(i) .and. k < kpbl(i)) then + xmf(i,k) = pfnl(i) * xmf(i,k) + endif + enddo + enddo + + do k = kmscu, 1, -1 + do i = 1, im + if(scuflg(i) .and. + & (k >= mrad(i) .and. k < krad(i))) then + xmfd(i,k) = pfnl(i) * xmfd(i,k) + endif + enddo + enddo + endif !sa3dtke !> ## Compute TKE. @@ -1783,8 +1802,9 @@ 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)) ! calculating 3D TKE transport and pressure correlation + dz = zl(i,k+1) - zl(i,k) ptem1 = ce0 / ele(i,k) tem1=(garea(i)*dz)**(inv3) tem2=0.19+0.51*ele_les(i,k)/tem1 @@ -1984,60 +2004,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & f1(i,1) = tke(i,1) enddo ! -!The following is for SA-3D-TKE - if(sa3dtke) then - do k = 1,km1 - do i=1,im - dtodsd = dt2/del(i,k) - dtodsu = dt2/del(i,k+1) - dsig = prsl(i,k)-prsl(i,k+1) - rdz = rdzt(i,k) - tem1 = dsig * dkq(i,k) * rdz - dsdz2 = tem1 * rdz - au(i,k) = -dtodsd*dsdz2 - al(i,k) = -dtodsu*dsdz2 - ad(i,k) = ad(i,k)-au(i,k) - ad(i,k+1)= 1.-al(i,k) - tem2 = dsig * rdz -! - if(pcnvflg(i) .and. k < kpbl(i)) then - ptem = 0.5 * tem2 * xmf(i,k) - ptem=ptem*pfnl(i) ! blending nolocal mixing - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - ptem = qcko(i,k,ntke) + qcko(i,k+1,ntke) - f1(i,k) = f1(i,k) - ptem * ptem1 - f1(i,k+1) = tke(i,k+1) + ptem * ptem2 - else - f1(i,k+1) = tke(i,k+1) - endif -! - if(scuflg(i)) then - if(k >= mrad(i) .and. k < krad(i)) then - ptem = 0.5 * tem2 * xmfd(i,k) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - ptem = qcdo(i,k,ntke) + qcdo(i,k+1,ntke) - f1(i,k) = f1(i,k) + ptem * ptem1 - f1(i,k+1) = f1(i,k+1) - ptem * ptem2 - endif - endif -! - kmx = max(kpbl(i), krad(i)) - if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then - ptem = tem2 * (xmf(i,k) - xmfd(i,k)) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - f1(i,k) = f1(i,k) + e_half(i,k) * ptem1 - f1(i,k+1) = f1(i,k+1) - e_half(i,k) * ptem2 - endif -! - enddo - enddo - else - do k = 1,km1 + do k = 1,km1 do i=1,im dtodsd = dt2/del(i,k) dtodsu = dt2/del(i,k+1) @@ -2083,8 +2050,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif ! enddo - enddo - endif !sa3dtke + enddo c !> - Call tridit() to solve tridiagonal problem for TKE @@ -2229,71 +2195,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo endif c -!The following is for SA-3D-TKE - if(sa3dtke) then - do k = 1,km1 - do i = 1,im - dtodsd = dt2/del(i,k) - dtodsu = dt2/del(i,k+1) - dsig = prsl(i,k)-prsl(i,k+1) - rdz = rdzt(i,k) - tem1 = dsig * dkt(i,k) * rdz - dsdzt = tem1 * gocp - dsdz2 = tem1 * rdz - au(i,k) = -dtodsd*dsdz2 - al(i,k) = -dtodsu*dsdz2 - ad(i,k) = ad(i,k)-au(i,k) - ad(i,k+1)= 1.-al(i,k) - tem2 = dsig * rdz -! - if(pcnvflg(i) .and. k < kpbl(i)) then - ptem = 0.5 * tem2 * xmf(i,k) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - tem = t1(i,k) + t1(i,k+1) - ptem = tcko(i,k) + tcko(i,k+1) - f1(i,k) = f1(i,k)+dtodsd*dsdzt-(ptem-tem)*ptem1 - f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+(ptem-tem)*ptem2 - ptem = qcko(i,k,1) + qcko(i,k+1,1) - f2(i,k) = f2(i,k) - ptem * ptem1 - f2(i,k+1) = q1(i,k+1,1) + ptem * ptem2 - else - f1(i,k) = f1(i,k)+dtodsd*dsdzt - f1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt - f2(i,k+1) = q1(i,k+1,1) - endif -! - if(scuflg(i)) then - if(k >= mrad(i) .and. k < krad(i)) then - ptem = 0.5 * tem2 * xmfd(i,k) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - ptem = tcdo(i,k) + tcdo(i,k+1) - tem = t1(i,k) + t1(i,k+1) - f1(i,k) = f1(i,k) + (ptem - tem) * ptem1 - f1(i,k+1) = f1(i,k+1) - (ptem - tem) * ptem2 - ptem = qcdo(i,k,1) + qcdo(i,k+1,1) - f2(i,k) = f2(i,k) + ptem * ptem1 - f2(i,k+1) = f2(i,k+1) - ptem * ptem2 - endif - endif -! - kmx = max(kpbl(i), krad(i)) - if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then - ptem = tem2 * (xmf(i,k) - xmfd(i,k)) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - f2(i,k) = f2(i,k) + q_half(i,k,1) * ptem1 - f2(i,k+1) = f2(i,k+1) - q_half(i,k,1) * ptem2 - endif -! - enddo - enddo - else - do k = 1,km1 + do k = 1,km1 do i = 1,im dtodsd = dt2/del(i,k) dtodsu = dt2/del(i,k+1) @@ -2350,61 +2252,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif ! enddo - enddo - endif !sa3dtke -! -!The following is for SA-3D-TKE - if(sa3dtke) then - if(ntrac1 >= 2) then - do n = 2, ntrac1 - is = (n-1) * km - do k = 1, km1 - do i = 1, im - dtodsd = dt2/del(i,k) - dtodsu = dt2/del(i,k+1) - dsig = prsl(i,k)-prsl(i,k+1) - tem2 = dsig * rdzt(i,k) -! - if(pcnvflg(i) .and. k < kpbl(i)) then - ptem = 0.5 * tem2 * xmf(i,k) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - ptem = qcko(i,k,n) + qcko(i,k+1,n) - f2(i,k+is) = f2(i,k+is) - ptem * ptem1 - f2(i,k+1+is)= q1(i,k+1,n) + ptem * ptem2 - else - f2(i,k+1+is) = q1(i,k+1,n) - endif -! - if(scuflg(i)) then - if(k >= mrad(i) .and. k < krad(i)) then - ptem = 0.5 * tem2 * xmfd(i,k) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - ptem = qcdo(i,k,n) + qcdo(i,k+1,n) - f2(i,k+is) = f2(i,k+is) + ptem * ptem1 - f2(i,k+1+is)= f2(i,k+1+is) - ptem * ptem2 - endif - endif -! - kmx = max(kpbl(i), krad(i)) - if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then - ptem = tem2 * (xmf(i,k) - xmfd(i,k)) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - f2(i,k+is) = f2(i,k+is) + q_half(i,k,n) * ptem1 - f2(i,k+1+is) = f2(i,k+1+is) - q_half(i,k,n) * ptem2 - endif + enddo ! - enddo - enddo - enddo - endif - else - if(ntrac1 >= 2) then + if(ntrac1 >= 2) then do n = 2, ntrac1 is = (n-1) * km do k = 1, km1 @@ -2448,8 +2298,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo enddo enddo - endif - endif !sa3dtke + endif c !> - Call tridin() to solve tridiagonal problem for heat and moisture c @@ -2850,61 +2699,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & f2(i,1) = v1(i,1) enddo c -!The following is for SA-3D-TKE - if(sa3dtke) then - do k = 1,km1 - do i=1,im - dtodsd = dt2/del(i,k) - dtodsu = dt2/del(i,k+1) - dsig = prsl(i,k)-prsl(i,k+1) - rdz = rdzt(i,k) - tem1 = dsig * dku(i,k) * rdz - dsdz2 = tem1*rdz - au(i,k) = -dtodsd*dsdz2 - al(i,k) = -dtodsu*dsdz2 - ad(i,k) = ad(i,k)-au(i,k) - ad(i,k+1)= 1.-al(i,k) - tem2 = dsig * rdz -! - if(pcnvflg(i) .and. k < kpbl(i)) then - ptem = 0.5 * tem2 * xmf(i,k) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - tem = u1(i,k) + u1(i,k+1) - ptem = ucko(i,k) + ucko(i,k+1) - f1(i,k) = f1(i,k) - (ptem - tem) * ptem1 - f1(i,k+1) = u1(i,k+1) + (ptem - tem) * ptem2 - tem = v1(i,k) + v1(i,k+1) - ptem = vcko(i,k) + vcko(i,k+1) - f2(i,k) = f2(i,k) - (ptem - tem) * ptem1 - f2(i,k+1) = v1(i,k+1) + (ptem - tem) * ptem2 - else - f1(i,k+1) = u1(i,k+1) - f2(i,k+1) = v1(i,k+1) - endif -! - if(scuflg(i)) then - if(k >= mrad(i) .and. k < krad(i)) then - ptem = 0.5 * tem2 * xmfd(i,k) - ptem=ptem*pfnl(i) - ptem1 = dtodsd * ptem - ptem2 = dtodsu * ptem - tem = u1(i,k) + u1(i,k+1) - ptem = ucdo(i,k) + ucdo(i,k+1) - f1(i,k) = f1(i,k) + (ptem - tem) *ptem1 - f1(i,k+1) = f1(i,k+1) - (ptem - tem) *ptem2 - tem = v1(i,k) + v1(i,k+1) - ptem = vcdo(i,k) + vcdo(i,k+1) - f2(i,k) = f2(i,k) + (ptem - tem) * ptem1 - f2(i,k+1) = f2(i,k+1) - (ptem - tem) * ptem2 - endif - endif -! - enddo - enddo - else - do k = 1,km1 + do k = 1,km1 do i=1,im dtodsd = dt2/del(i,k) dtodsu = dt2/del(i,k+1) @@ -2952,8 +2747,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif ! enddo - enddo - endif !sa3dtke + enddo c !> - Call tridi2() to solve tridiagonal problem for momentum c diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.meta b/physics/PBL/SATMEDMF/satmedmfvdifq.meta index 25b7817b8..ef412ec47 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.meta +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.meta @@ -273,8 +273,8 @@ kind = kind_phys intent = in [def_3] - standard_name = transfer_rate_of_tke_due_to_dynamics - long_name = rate of TKE transfer and pressure correlation calculated from dynamics + standard_name = horizontal_transfer_rate_of_tke_due_to_dynamics + long_name = rate of horizontal TKE transfer and pressure correlation calculated from dynamics units = m2 s-3 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -566,7 +566,7 @@ kind = kind_phys intent = out [dku3d_h] - standard_name = horizontal_atmosphere_momentum_diffusivity_for_dyanmics + standard_name = horizontal_atmosphere_momentum_diffusivity_for_dynamics long_name = horizontal atmospheric momentum diffusivity for dynamics units = m2 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) @@ -574,7 +574,7 @@ kind = kind_phys intent = out [dku3d_e] - standard_name = horizontal_atmosphere_tke_diffusivity_for_dyanmics + standard_name = horizontal_atmosphere_tke_diffusivity_for_dynamics long_name = horizontal atmospheric tke diffusivity for dynamics units = m2 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) From 11a126b1568d3d527c96d8cf888016abf89778a4 Mon Sep 17 00:00:00 2001 From: "Kwun Y. Fung" Date: Tue, 17 Jun 2025 12:17:16 -0500 Subject: [PATCH 6/6] Update nonlocal mass-flux transport for Shear production of TKE --- physics/PBL/SATMEDMF/satmedmfvdifq.F | 46 +++++++++++++++++++++++++--- 1 file changed, 42 insertions(+), 4 deletions(-) diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.F b/physics/PBL/SATMEDMF/satmedmfvdifq.F index ece0256e2..bb00d9c61 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.F +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.F @@ -1606,6 +1606,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif tem = tem + ptem1 + ptem2 buop = 0.5 * (gotvx(i,1) * sflux(i) + tem) +! + tem1 = dku(i,1)*def_1(i,1)+dku_h(i,1)*def_2(i,1) ! tem = (u1(i,2)-u1(i,1))*rdzt(i,1) ! if(pcnvflg(i)) then @@ -1638,9 +1640,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ptem2 = ptem2 + ptem ! tem2 = stress(i)*ustar(i)*phims(i)/(vk*zl(i,1)) - shrp = 0.5 * (ptem1 + ptem2 + tem2) - tem2 = dku(i,k)*def_1(i,k)+dku_h(i,k)*def_2(i,k) - shrp =0.5*(shrp+tem2) + shrp = 0.5 * (tem1 + ptem1 + ptem2 + tem2) else tem1 = -dkt(i,k-1) * bf(i,k-1) tem2 = -dkt(i,k) * bf(i,k) @@ -1666,7 +1666,45 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! obtaining 3d shear production from dycore tem1 = dku(i,k-1)*def_1(i,k-1)+dku_h(i,k-1)*def_2(i,k-1) tem2 = dku(i,k)*def_1(i,k)+dku_h(i,k)*def_2(i,k) - shrp=0.5*(tem1+tem2) + tem = 0.5*(tem1+tem2) + tem1 = (u1(i,k+1)-u1(i,k))*rdzt(i,k) + tem2 = (u1(i,k)-u1(i,k-1))*rdzt(i,k-1) + if(pcnvflg(i) .and. k <= kpbl(i)) then + ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2 + ptem1 = 0.5 * ptem * (u1(i,k)-ucko(i,k)) + else + ptem1 = 0. + endif + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2 + ptem2 = 0.5 * ptem0 * (ucdo(i,k)-u1(i,k)) + else + ptem2 = 0. + endif + else + ptem2 = 0. + endif + shrp = tem + ptem1 + ptem2 + tem1 = (v1(i,k+1)-v1(i,k))*rdzt(i,k) + tem2 = (v1(i,k)-v1(i,k-1))*rdzt(i,k-1) + if(pcnvflg(i) .and. k <= kpbl(i)) then + ptem = xmf(i,k) * tem1 + xmf(i,k-1) * tem2 + ptem1 = 0.5 * ptem * (v1(i,k)-vcko(i,k)) + else + ptem1 = 0. + endif + if(scuflg(i)) then + if(k >= mrad(i) .and. k < krad(i)) then + ptem0 = xmfd(i,k) * tem1 + xmfd(i,k-1) * tem2 + ptem2 = 0.5 * ptem0 * (vcdo(i,k)-v1(i,k)) + else + ptem2 = 0. + endif + else + ptem2 = 0. + endif + shrp = shrp + ptem1 + ptem2 endif ! prod(i,k) = buop + shrp