diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 1c524b800..08876f8f0 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -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 !! @@ -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) @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -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 @@ -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)) @@ -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 diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 8538c6aa7..d9ab8c859 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -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)