From 1590813076b29c1e76acd531cfd747bc14936cf1 Mon Sep 17 00:00:00 2001 From: AndrewHazelton <50676372+AndrewHazelton@users.noreply.github.com> Date: Fri, 19 Aug 2022 20:58:24 -0400 Subject: [PATCH 1/4] Add files via upload Updated Physics Files --- physics/satmedmfvdifq.F | 142 +++++++++++++++++++++++++++++-------- physics/satmedmfvdifq.meta | 7 ++ 2 files changed, 121 insertions(+), 28 deletions(-) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 9c5ad4029..c281b28b2 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -17,7 +17,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 !! @@ -77,7 +82,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & & 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) @@ -91,6 +96,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & 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 @@ -244,6 +250,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & real(kind=kind_phys) qlcr, zstblmax, hcrinv ! real(kind=kind_phys) h1 + + real(kind=kind_phys) bfac, mffac !! parameter(wfac=7.0,cfac=4.5) parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1) @@ -261,10 +269,19 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & 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 @@ -594,11 +611,19 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & 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 + bfac = 100.0 + 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 @@ -668,11 +693,19 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & 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 + bfac = 100.0 + 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 @@ -790,9 +823,17 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & 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 + bfac = 100. + 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 @@ -923,6 +964,31 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & & 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 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1007,10 +1073,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & 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 @@ -1041,10 +1112,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & 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)) @@ -1111,8 +1186,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & 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) + 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 !> - If sfc_rlm=1, use zk for elm within surface layer if ( sfc_rlm == 1 ) then @@ -1123,7 +1204,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & 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 db89f488d..e0ff383d5 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -573,6 +573,13 @@ dimensions = () type = integer intent = in +[tc_pbl] + standard_name = option_of_tc_application_in_boundary_layer_mass_flux_scheme + long_name = option of tc application in boundary layer mass flux 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 3dea70eeebbef01053f61d3a5506786cebd235ef Mon Sep 17 00:00:00 2001 From: Bin Liu Date: Wed, 7 Dec 2022 19:48:29 +0000 Subject: [PATCH 2/4] Update satmedmfvdifq.F to address some PR review comments related to the tc_pbl option. --- physics/satmedmfvdifq.F | 56 ++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 31 deletions(-) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 3b7a72b3c..a314b1a27 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -273,13 +273,13 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & parameter(rchck=1.5,ndt=20) if (tc_pbl == 0) then - ck0=0.4 - ch0=0.4 - ce0=0.4 + ck0 = 0.4 + ch0 = 0.4 + ce0 = 0.4 else if (tc_pbl == 1) then - ck0=0.55 - ch0=0.55 - ce0=0.12 + ck0 = 0.55 + ch0 = 0.55 + ce0 = 0.12 endif gravi = 1.0 / grav g = grav @@ -973,28 +973,23 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & & 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 -! + 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 + 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 - enddo - endif ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1197,17 +1192,16 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & 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 + 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 -!> - 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 - endif -! dz = zi(i,k+1) - zi(i,k) tem = max(gdx(i),dz) elm(i,k) = min(elm(i,k), tem) From cee2bbcc265548098053d19cbccedf2fdb943a24 Mon Sep 17 00:00:00 2001 From: Bin Liu Date: Fri, 9 Dec 2022 11:08:32 +0000 Subject: [PATCH 3/4] For the tc_pbl related changes in satmedmfvdifq.F, set bfac as a parameter (according to @grantfirl's suggestion). --- physics/satmedmfvdifq.F | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index a314b1a27..08876f8f0 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -252,6 +252,7 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & 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) @@ -625,7 +626,6 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & rbup(i) = (thlvx(i,k)-thermal(i))* & (g*zl(i,k)/thlvx(i,1))/spdk2 else if (tc_pbl == 1) then - bfac = 100.0 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))* @@ -707,7 +707,6 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & rbup(i) = (thlvx(i,k)-thermal(i))* & (g*zl(i,k)/thlvx(i,1))/spdk2 else if (tc_pbl == 1) then - bfac = 100.0 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))* @@ -835,7 +834,6 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & rbup(i) = (thlvx(i,k)-thermal(i))* & (g*zl(i,k)/thlvx(i,1))/spdk2 else if (tc_pbl == 1) then - bfac = 100. 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))* From a3c7ac03e1cbb4c2fba2972d595c36991246d2d1 Mon Sep 17 00:00:00 2001 From: Bin Liu Date: Fri, 9 Dec 2022 19:21:11 +0000 Subject: [PATCH 4/4] Update tc_pbl standard and long names in satmedmfvdifq.meta according to PR review comments/suggestions. --- physics/satmedmfvdifq.meta | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 124e52945..d9ab8c859 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -589,8 +589,8 @@ type = integer intent = in [tc_pbl] - standard_name = option_of_tc_application_in_boundary_layer_mass_flux_scheme - long_name = option of tc application in boundary layer mass flux scheme + 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