From 54a57baa8d8874a09081cd3b51fc9dc4d53ad4e8 Mon Sep 17 00:00:00 2001 From: "Chunxi.Zhang-NOAA" Date: Tue, 15 Mar 2022 14:28:24 +0000 Subject: [PATCH] P8C updates: the TKE-EDMF PBL scheme and the saSAS cumulus scheme --- physics/mfpbltq.f | 28 ++++++++----- physics/mfscuq.f | 28 ++++++++----- physics/samfdeepcnv.f | 88 +++++++++++++++++++++++++++++++++++------ physics/samfshalcnv.f | 83 ++++++++++++++++++++++++++++++-------- physics/satmedmfvdifq.F | 78 ++++++++++++++++++++++++++++++++++-- 5 files changed, 253 insertions(+), 52 deletions(-) diff --git a/physics/mfpbltq.f b/physics/mfpbltq.f index b906052cd..a0788d5b7 100644 --- a/physics/mfpbltq.f +++ b/physics/mfpbltq.f @@ -11,7 +11,7 @@ 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,xmf, - & tcko,qcko,ucko,vcko,xlamue,a1) + & tcko,qcko,ucko,vcko,xlamueq,a1) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -35,14 +35,15 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, & buo(im,km), xmf(im,km), & tcko(im,km),qcko(im,km,ntrac1), & ucko(im,km),vcko(im,km), - & xlamue(im,km-1) + & xlamueq(im,km-1) ! c local variables and arrays ! integer i, j, k, n, ndc integer kpblx(im), kpbly(im) ! - real(kind=kind_phys) dt2, dz, ce0, cm, + real(kind=kind_phys) dt2, dz, ce0, + & cm, cq, & factor, gocp, & g, b1, f1, & bb1, bb2, @@ -56,7 +57,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, & thup, thvu, dq ! real(kind=kind_phys) rbdn(im), rbup(im), hpblx(im), - & xlamuem(im,km-1) + & xlamue(im,km-1), xlamuem(im,km-1) real(kind=kind_phys) delz(im), xlamax(im) ! real(kind=kind_phys) wu2(im,km), thlu(im,km), @@ -71,7 +72,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, parameter(g=grav) parameter(gocp=g/cp) parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) - parameter(ce0=0.4,cm=1.0) + parameter(ce0=0.4,cm=1.0,cq=1.3) parameter(qmin=1.e-8,qlmin=1.e-12) parameter(alp=1.5,vpertmax=3.0,pgcon=0.55) parameter(b1=0.5,f1=0.15) @@ -132,6 +133,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, xlamue(i,k) = xlamax(i) endif ! + xlamueq(i,k) = cq * xlamue(i,k) xlamuem(i,k) = cm * xlamue(i,k) endif enddo @@ -148,6 +150,9 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, ! thlu(i,k) = ((1.-tem)*thlu(i,k-1)+tem* & (thlx(i,k-1)+thlx(i,k)))/factor +! + tem = 0.5 * xlamueq(i,k-1) * dz + factor = 1. + tem qtu(i,k) = ((1.-tem)*qtu(i,k-1)+tem* & (qtx(i,k-1)+qtx(i,k)))/factor ! @@ -282,6 +287,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, xlamue(i,k) = xlamax(i) endif ! + xlamueq(i,k) = cq * xlamue(i,k) xlamuem(i,k) = cm * xlamue(i,k) endif enddo @@ -313,7 +319,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, do k = 1, kmpbl do i = 1, im if (cnvflg(i) .and. k < kpbl(i)) then - xmf(i,k) = a1 * sqrt(wu2(i,k)) + xmf(i,k) = sqrt(wu2(i,k)) endif enddo enddo @@ -350,7 +356,8 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, do k = 1, kmpbl do i = 1, im if (cnvflg(i) .and. k < kpbl(i)) then - xmf(i,k) = scaldfunc(i) * xmf(i,k) + tem = max(a1, sigma(i)) + xmf(i,k) = scaldfunc(i) * tem * xmf(i,k) dz = zl(i,k+1) - zl(i,k) xmmx = dz / dt2 xmf(i,k) = min(xmf(i,k),xmmx) @@ -384,6 +391,9 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, ! thlu(i,k) = ((1.-tem)*thlu(i,k-1)+tem* & (thlx(i,k-1)+thlx(i,k)))/factor +! + tem = 0.5 * xlamueq(i,k-1) * dz + factor = 1. + tem qtu(i,k) = ((1.-tem)*qtu(i,k-1)+tem* & (qtx(i,k-1)+qtx(i,k)))/factor ! @@ -432,7 +442,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, do i = 1, im if (cnvflg(i) .and. k <= kpbl(i)) then dz = zl(i,k) - zl(i,k-1) - tem = 0.5 * xlamue(i,k-1) * dz + tem = 0.5 * xlamueq(i,k-1) * dz factor = 1. + tem ! qcko(i,k,n) = ((1.-tem)*qcko(i,k-1,n)+tem* @@ -453,7 +463,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, do i = 1, im if (cnvflg(i) .and. k <= kpbl(i)) then dz = zl(i,k) - zl(i,k-1) - tem = 0.5 * xlamue(i,k-1) * dz + tem = 0.5 * xlamueq(i,k-1) * dz factor = 1. + tem ! qcko(i,k,n) = ((1.-tem)*qcko(i,k-1,n)+tem* diff --git a/physics/mfscuq.f b/physics/mfscuq.f index 3390c3e58..b41ffd13e 100644 --- a/physics/mfscuq.f +++ b/physics/mfscuq.f @@ -11,7 +11,7 @@ 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,xmfd, - & tcdo,qcdo,ucdo,vcdo,xlamde,a1) + & tcdo,qcdo,ucdo,vcdo,xlamdeq,a1) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -39,7 +39,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, & buo(im,km), xmfd(im,km), & tcdo(im,km), qcdo(im,km,ntrac1), & ucdo(im,km), vcdo(im,km), - & xlamde(im,km-1) + & xlamdeq(im,km-1) ! ! local variables and arrays ! @@ -47,7 +47,8 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, integer i,j,indx, k, n, kk, ndc integer krad1(im) ! - real(kind=kind_phys) dt2, dz, ce0, cm, + real(kind=kind_phys) dt2, dz, ce0, + & cm, cq, & gocp, factor, g, tau, & b1, f1, bb1, bb2, & a1, a2, @@ -62,7 +63,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, ! real(kind=kind_phys) wd2(im,km), thld(im,km), & qtx(im,km), qtd(im,km), - & thlvd(im), hrad(im), + & thlvd(im), hrad(im), xlamde(im,km-1), & xlamdem(im,km-1), ra1(im) real(kind=kind_phys) delz(im), xlamax(im) ! @@ -77,7 +78,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, parameter(g=grav) parameter(gocp=g/cp) parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) - parameter(ce0=0.4,cm=1.0,pgcon=0.55) + parameter(ce0=0.4,cm=1.0,cq=1.3,pgcon=0.55) parameter(qmin=1.e-8,qlmin=1.e-12) parameter(b1=0.45,f1=0.15) parameter(a2=0.5) @@ -208,6 +209,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, xlamde(i,k) = xlamax(i) endif ! + xlamdeq(i,k) = cq * xlamde(i,k) xlamdem(i,k) = cm * xlamde(i,k) endif enddo @@ -224,6 +226,9 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, ! thld(i,k) = ((1.-tem)*thld(i,k+1)+tem* & (thlx(i,k)+thlx(i,k+1)))/factor +! + tem = 0.5 * xlamdeq(i,k) * dz + factor = 1. + tem qtd(i,k) = ((1.-tem)*qtd(i,k+1)+tem* & (qtx(i,k)+qtx(i,k+1)))/factor ! @@ -347,6 +352,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, xlamde(i,k) = xlamax(i) endif ! + xlamdeq(i,k) = cq * xlamde(i,k) xlamdem(i,k) = cm * xlamde(i,k) endif enddo @@ -380,7 +386,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, do i = 1, im if(cnvflg(i) .and. & (k >= mrad(i) .and. k < krad(i))) then - xmfd(i,k) = ra1(i) * sqrt(wd2(i,k)) + xmfd(i,k) = sqrt(wd2(i,k)) endif enddo enddo @@ -418,7 +424,8 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, do i = 1, im if(cnvflg(i) .and. & (k >= mrad(i) .and. k < krad(i))) then - xmfd(i,k) = scaldfunc(i) * xmfd(i,k) + tem = max(ra1(i), sigma(i)) + xmfd(i,k) = scaldfunc(i) * tem * xmfd(i,k) dz = zl(i,k+1) - zl(i,k) xmmx = dz / dt2 xmfd(i,k) = min(xmfd(i,k),xmmx) @@ -457,6 +464,9 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, ! thld(i,k) = ((1.-tem)*thld(i,k+1)+tem* & (thlx(i,k)+thlx(i,k+1)))/factor +! + tem = 0.5 * xlamdeq(i,k) * dz + factor = 1. + tem qtd(i,k) = ((1.-tem)*qtd(i,k+1)+tem* & (qtx(i,k)+qtx(i,k+1)))/factor ! @@ -509,7 +519,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, if (cnvflg(i) .and. k < krad(i)) then if(k >= mrad(i)) then dz = zl(i,k+1) - zl(i,k) - tem = 0.5 * xlamde(i,k) * dz + tem = 0.5 * xlamdeq(i,k) * dz factor = 1. + tem ! qcdo(i,k,n) = ((1.-tem)*qcdo(i,k+1,n)+tem* @@ -532,7 +542,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, if (cnvflg(i) .and. k < krad(i)) then if(k >= mrad(i)) then dz = zl(i,k+1) - zl(i,k) - tem = 0.5 * xlamde(i,k) * dz + tem = 0.5 * xlamdeq(i,k) * dz factor = 1. + tem ! qcdo(i,k,n) = ((1.-tem)*qcdo(i,k+1,n)+tem* diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 3801e684f..0420fa1d2 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -149,16 +149,16 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & dh, dhh, dp, & dq, dqsdp, dqsdt, dt, & dt2, dtmax, dtmin, - & dxcrtas, dxcrtuf, +! & dxcrtas, dxcrtuf, dxcrtc0, + & dxcrtas, dxcrtuf, & dv1h, dv2h, dv3h, - & dv2q, & dz, dz1, e1, edtmax, & edtmaxl, edtmaxs, el2orc, elocp, & es, etah, & cthk, dthk, ! & evfact, evfactl, & fact1, fact2, factor, - & gamma, pprime, cm, + & gamma, pprime, cm, cq, & qlk, qrch, qs, & rain, rfact, shear, tfac, & val, val1, val2, @@ -225,7 +225,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! Until a realistic Nccn is provided, Nccns are assumed ! as Nccn=100 for sea and Nccn=1000 for land ! - parameter(cm=1.0) + parameter(cm=1.0,cq=1.3) ! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c) parameter(clamd=0.03,tkemx=0.65,tkemn=0.05) parameter(clamca=0.03) @@ -236,6 +236,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & parameter(cinacrmx=-120.,cinacrmn=-80.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) parameter(betaw=.03,dxcrtas=8.e3,dxcrtuf=15.e3) +! parameter(dxcrtc0=9.e3) ! ! local variables and arrays @@ -249,8 +250,12 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & wet_dep ! ! for updraft velocity calculation - real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km) - real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im) + real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km), + & wc(im) +! +! for updraft fraction & scale-aware function +! real(kind=kind_phys) scaldfunc(im), sigmagfm(im), xlamumean(im) + real(kind=kind_phys) scaldfunc(im), sigmagfm(im) ! c cloud water ! real(kind=kind_phys) tvo(im,km) @@ -392,6 +397,16 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & c0(i) = c0s endif enddo +! +!> - determine scale-aware rain conversion parameter decreasing with decreasing grid size +! do i=1,im +! if(gdx(i) < dxcrtc0) then +! tem = gdx(i) / dxcrtc0 +! tem1 = tem**2 +! c0(i) = c0(i) * tem1 +! endif +! enddo +! !> - determine rain conversion parameter above the freezing level which exponentially decreases with decreasing temperature from Han et al.'s (2017) \cite han_et_al_2017 equation 8. do k = 1, km do i = 1, im @@ -1013,6 +1028,33 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & enddo enddo endif +! +! compute mean entrainment rate in subcloud layers below cloud base +! +! do i= 1, im +! if(cnvflg(i)) then +! sumx(i) = 0. +! xlamumean(i) = 0. +! endif +! enddo +! do k = 1, km1 +! do i = 1, im +! if(cnvflg(i)) then +! if(k >= kb(i) .and. k < kbcon(i)) then +! dz = zi(i,k+1) - zi(i,k) +! tem = 0.5 * (xlamue(i,k)+xlamue(i,k+1)) +! xlamumean(i) = xlamumean(i) + tem * dz +! sumx(i) = sumx(i) + dz +! endif +! endif +! enddo +! enddo +! +! do i= 1, im +! if(cnvflg(i)) then +! xlamumean(i) = xlamumean(i) / sumx(i) +! endif +! enddo c c specify detrainment rate for the updrafts c @@ -1192,6 +1234,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if(k > kb(i) .and. k < kmax(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz + tem = cq * tem factor = 1. + tem ecko(i,k,n) = ((1.-tem)*ecko(i,k-1,n)+tem* & (ctro(i,k,n)+ctro(i,k-1,n)))/factor @@ -1209,6 +1252,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if(k > kb(i) .and. k < kmax(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz + tem = cq * tem factor = 1. + tem ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem * & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor @@ -1461,6 +1505,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & cj tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz + tem = cq * tem + tem1 = cq * tem1 factor = 1. + tem - tem1 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* & (qo(i,k)+qo(i,k-1)))/factor @@ -1636,6 +1682,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & cj tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz + tem = cq * tem + tem1 = cq * tem1 factor = 1. + tem - tem1 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* & (qo(i,k)+qo(i,k-1)))/factor @@ -1926,6 +1974,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if (cnvflg(i) .and. k < jmin(i)) then dz = zi(i,k+1) - zi(i,k) tem = 0.5 * xlamde * dz + tem = cq * tem factor = 1. + tem ecdo(i,k,n) = ((1.-tem)*ecdo(i,k+1,n)+tem* & (ctro(i,k,n)+ctro(i,k+1,n)))/factor @@ -1952,6 +2001,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & tem = xlamde * dz tem1 = 0.5 * (xlamd(i)+xlamdd) * dz endif + tem = cq * tem + tem1 = cq * tem1 factor = 1. + tem - tem1 qcdo(i,k) = ((1.-tem1)*qrcdo(i,k+1)+tem*0.5* & (qo(i,k)+qo(i,k+1)))/factor @@ -2084,7 +2135,6 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & dv1h = heo(i,k) dv2h = .5 * (heo(i,k) + heo(i,k-1)) dv3h = heo(i,k-1) - dv2q = .5 * (qo(i,k) + qo(i,k-1)) c tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) tem1 = 0.5 * (xlamud(i,k)+xlamud(i,k-1)) @@ -2107,11 +2157,12 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz & ) * factor cj + tem1 = -eta(i,k) * qrcko(i,k) + tem2 = -eta(i,k-1) * qcko(i,k-1) + ptem1 = -etad(i,k) * qrcdo(i,k) + ptem2 = -etad(i,k-1) * qcdo(i,k-1) dellaq(i,k) = dellaq(i,k) + - & (- (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz - & + aup*tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz - & + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qcdo(i,k-1))*dz - & ) * factor + & (aup*(tem1-tem2)-adw*edto(i)*(ptem1-ptem2))*factor cj tem1=eta(i,k)*(uo(i,k)-ucko(i,k)) tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1)) @@ -2502,6 +2553,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & cj tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz tem1 = 0.25 * (xlamud(i,k)+xlamud(i,k-1)) * dz + tem = cq * tem + tem1 = cq * tem1 factor = 1. + tem - tem1 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* & (qo(i,k)+qo(i,k-1)))/factor @@ -2596,6 +2649,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & tem = xlamde * dz tem1 = 0.5 * (xlamd(i)+xlamdd) * dz endif + tem = cq * tem + tem1 = cq * tem1 factor = 1. + tem - tem1 qcdo(i,k) = ((1.-tem1)*qrcd(i,k+1)+tem*0.5* & (qo(i,k)+qo(i,k+1)))/factor @@ -2775,7 +2830,8 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! tfac = tauadv(i) / dtconv(i) ! tfac = min(tfac, 1.) ! xmb(i) = tfac*betaw*rho*wc(i) - xmb(i) = betaw*rho*wc(i) +! xmb(i) = betaw*rho*wc(i) + xmb(i) = rho*wc(i) endif enddo !> - For the cases where the quasi-equilibrium assumption of Arakawa-Schubert is valid, first calculate the large scale destabilization as in equation 5 of Pan and Wu (1995) \cite pan_and_wu_1995 : @@ -2836,6 +2892,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then tem = min(max(xlamx(i), 7.e-5), 3.e-4) +! tem = min(max(xlamumean(i), 1.e-4), 1.e-3) tem = 0.2 / tem tem1 = 3.14 * tem * tem sigmagfm(i) = tem1 / garea(i) @@ -2865,7 +2922,12 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & else scaldfunc(i) = 1.0 endif - xmb(i) = xmb(i) * scaldfunc(i) + if(asqecflg(i)) then + xmb(i) = xmb(i) * scaldfunc(i) + else + tem = max(betaw, sigmagfm(i)) + xmb(i) = tem * xmb(i) * scaldfunc(i) + endif xmb(i) = min(xmb(i),xmbmax(i)) endif enddo diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 0e11ed49c..68b12d169 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -102,12 +102,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & c0l, d0, & desdt, dp, & dq, dqsdp, dqsdt, dt, - & dt2, dtmax, dtmin, dxcrt, + & dt2, dtmax, dtmin, + & dxcrt, dxcrtc0, & dv1h, dv2h, dv3h, - & dv2q, & dz, dz1, e1, - & el2orc, elocp, aafac, cm, - & es, etah, h1, + & el2orc, elocp, aafac, + & cm, cq, + & es, etah, h1, shevf, ! & evfact, evfactl, & fact1, fact2, factor, dthk, & gamma, pprime, betaw, @@ -172,16 +173,17 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! Until a realistic Nccn is provided, Nccns are assumed ! as Nccn=100 for sea and Nccn=1000 for land ! - parameter(cm=1.0) + parameter(cm=1.0,cq=1.3) ! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c) parameter(clamd=0.1,tkemx=0.65,tkemn=0.05) parameter(dtke=tkemx-tkemn) parameter(dthk=25.,sfclfac=0.2,rhcrt=0.75) parameter(cinpcrmx=180.,cinpcrmn=120.) - parameter(cinacrmx=-120.) +! shevf is an enhancing evaporation factor for shallow convection + parameter(cinacrmx=-120.,shevf=1.0) parameter(dtmax=10800.,dtmin=600.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) - parameter(betaw=.03,dxcrt=15.e3) + parameter(betaw=.03,dxcrt=15.e3,dxcrtc0=9.e3) parameter(h1=0.33333333) c local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), @@ -195,8 +197,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys), parameter :: escav = 0.8 ! wet scavenging efficiency ! ! for updraft velocity calculation - real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km) - real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im) + real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km), + & wc(im) +! +! for updraft fraction & scale-aware function +! real(kind=kind_phys) scaldfunc(im), sigmagfm(im), xlamumean(im) + real(kind=kind_phys) scaldfunc(im), sigmagfm(im) ! c cloud water ! real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km), @@ -337,6 +343,15 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif enddo ! +!> - determine scale-aware rain conversion parameter decreasing with decreasing grid size + do i=1,im + if(gdx(i) < dxcrtc0) then + tem = gdx(i) / dxcrtc0 + tem1 = tem**3 + c0(i) = c0(i) * tem1 + endif + enddo +! !> - determine rain conversion parameter above the freezing level which exponentially decreases with decreasing temperature from Han et al.'s (2017) \cite han_et_al_2017 equation 8. do k = 1, km do i = 1, im @@ -889,6 +904,33 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif enddo endif ! hwrf_samfshal +! +! compute mean entrainment rate in subcloud layers below cloud base +! +! do i= 1, im +! if(cnvflg(i)) then +! sumx(i) = 0. +! xlamumean(i) = 0. +! endif +! enddo +! do k = 1, km1 +! do i = 1, im +! if(cnvflg(i)) then +! if(k >= kb(i) .and. k < kbcon(i)) then +! dz = zi(i,k+1) - zi(i,k) +! tem = 0.5 * (xlamue(i,k)+xlamue(i,k+1)) +! xlamumean(i) = xlamumean(i) + tem * dz +! sumx(i) = sumx(i) + dz +! endif +! endif +! enddo +! enddo +! +! do i= 1, im +! if(cnvflg(i)) then +! xlamumean(i) = xlamumean(i) / sumx(i) +! endif +! enddo c c determine updraft mass flux for the subcloud layers c @@ -996,6 +1038,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(k > kb(i) .and. k < kmax(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz + tem = cq * tem factor = 1. + tem ecko(i,k,n) = ((1.-tem)*ecko(i,k-1,n)+tem* & (ctro(i,k,n)+ctro(i,k-1,n)))/factor @@ -1013,6 +1056,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(k > kb(i) .and. k < kmax(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz + tem = cq * tem factor = 1. + tem ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem * & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor @@ -1194,6 +1238,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & cj tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz tem1 = 0.5 * xlamud(i) * dz + tem = cq * tem + tem1 = cq * tem1 factor = 1. + tem - tem1 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* & (qo(i,k)+qo(i,k-1)))/factor @@ -1360,6 +1406,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & cj tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz tem1 = 0.5 * xlamud(i) * dz + tem = cq * tem + tem1 = cq * tem1 factor = 1. + tem - tem1 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* & (qo(i,k)+qo(i,k-1)))/factor @@ -1565,7 +1613,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & dv1h = heo(i,k) dv2h = .5 * (heo(i,k) + heo(i,k-1)) dv3h = heo(i,k-1) - dv2q = .5 * (qo(i,k) + qo(i,k-1)) c tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) tem1 = xlamud(i) @@ -1578,10 +1625,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz & ) * factor cj - dellaq(i,k) = dellaq(i,k) + - & ( - tem*eta(i,k-1)*dv2q*dz - & + tem1*eta(i,k-1)*.5*(qrcko(i,k)+qcko(i,k-1))*dz - & ) * factor + tem1 = -eta(i,k) * qrcko(i,k) + tem2 = -eta(i,k-1) * qcko(i,k-1) + dellaq(i,k) = dellaq(i,k) + (tem1-tem2) * factor cj tem1=eta(i,k)*(uo(i,k)-ucko(i,k)) tem2=eta(i,k-1)*(uo(i,k-1)-ucko(i,k-1)) @@ -1813,7 +1859,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! tfac = tauadv(i) / dtconv(i) ! tfac = min(tfac, 1.) ! xmb(i) = tfac*betaw*rho*wc(i) - xmb(i) = betaw*rho*wc(i) +! xmb(i) = betaw*rho*wc(i) + xmb(i) = rho*wc(i) endif enddo ! @@ -1821,6 +1868,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then tem = min(max(xlamue(i,kbcon(i)), 2.e-4), 6.e-4) +! tem = min(max(xlamumean(i), 2.e-4), 2.e-3) tem = 0.2 / tem tem1 = 3.14 * tem * tem sigmagfm(i) = tem1 / garea(i) @@ -1838,7 +1886,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & else scaldfunc(i) = 1.0 endif - xmb(i) = xmb(i) * scaldfunc(i) + tem = max(betaw, sigmagfm(i)) + xmb(i) = tem * xmb(i) * scaldfunc(i) xmb(i) = min(xmb(i),xmbmax(i)) endif enddo @@ -2145,7 +2194,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! evef = edt(i) * evfact ! if(islimsk(i) == 1) evef=edt(i) * evfactl ! if(islimsk(i) == 1) evef=.07 - qcond(i) = evef * (q1(i,k) - qeso(i,k)) + qcond(i) = shevf * evef * (q1(i,k) - qeso(i,k)) & / (1. + el2orc * qeso(i,k) / t1(i,k)**2) dp = 1000. * del(i,k) factor = dp / grav diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index be54675b0..eb2b7ad1c 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -138,7 +138,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & integer i,is,k,n,ndt,km1,kmpbl,kmscu,ntrac1,idtend integer kps,kbx,kmx integer lcld(im),kcld(im),krad(im),mrad(im) - integer kx1(im), kpblx(im) + integer kx1(im), kb1(im), kpblx(im) ! real(kind=kind_phys) tke(im,km), tkeh(im,km-1) ! @@ -198,6 +198,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & & q_diff(im,0:km-1,ntrac-1) real(kind=kind_phys) rrkp, phkp real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im) + real(kind=kind_phys) sfcpbl(im) ! logical pblflg(im), sfcflg(im), flg(im) logical scuflg(im), pcnvflg(im) @@ -233,6 +234,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & & zlup, zldn, bsum, cs0, & tem, tem1, tem2, tem3, & ptem, ptem0, ptem1, ptem2 +! + real(kind=kind_phys) slfac ! real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck ! @@ -242,7 +245,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & !! parameter(wfac=7.0,cfac=4.5) parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1) - parameter(vk=0.4,rimin=-100.) + parameter(vk=0.4,rimin=-100.,slfac=0.1) parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3) parameter(rlmn=30.,rlmn0=5.,rlmn1=5.,rlmn2=10.) parameter(prmin=0.25,prmax=4.0) @@ -573,7 +576,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & enddo enddo ! -! Find pbl height based on bulk richardson number (mrf pbl scheme) +! Find first quess pbl height based on bulk richardson number (mrf pbl scheme) ! and also for diagnostic purpose ! do i=1,im @@ -623,6 +626,73 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & if(kpbl(i) <= 1) pblflg(i)=.false. enddo ! +! update thermal at a level of slfac*hpbl for unstable pbl +! + do i=1,im + sfcpbl(i) = slfac * hpbl(i) + kb1(i) = 1 + flg(i) = .false. + if(pblflg(i)) then + flg(i) = .true. + endif + enddo + do k = 2, kmpbl + do i=1,im + if (flg(i) .and. zl(i,k) <= sfcpbl(i)) then + kb1(i) = k + else + flg(i) = .false. + endif + enddo + enddo + do i=1,im + if(pblflg(i)) kb1(i)=min(kb1(i),kpbl(i)) + enddo +! +! re-compute pbl height with the updated thermal +! + do i=1,im + flg(i) = .true. + if(pblflg(i) .and. kb1(i) > 1) then + flg(i) = .false. + rbup(i) = rbsoil(i) +! thermal(i) = thvx(i,kb1(i)) + thermal(i) = thlvx(i,kb1(i)) + kpblx(i) = kb1(i) + hpblx(i) = zl(i,kb1(i)) + endif + enddo + do k = 2, kmpbl + 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 + kpblx(i) = k + flg(i) = rbup(i) > crb(i) + endif + enddo + enddo + do i = 1,im + if(pblflg(i) .and. kb1(i) > 1) then + k = kpblx(i) + if(rbdn(i) >= crb(i)) then + rbint = 0. + elseif(rbup(i) <= crb(i)) then + rbint = 1. + else + rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i)) + endif + hpblx(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1)) + if(hpblx(i) < zi(i,kpblx(i))) kpblx(i)=kpblx(i)-1 + hpbl(i) = hpblx(i) + kpbl(i) = kpblx(i) + endif + enddo +! !> ## Compute Monin-Obukhov similarity parameters !! - Calculate the Monin-Obukhov nondimensional stability paramter, commonly !! referred to as \f$\zeta\f$ using the following equation from Businger et al.(1971) \cite businger_et_al_1971 @@ -716,7 +786,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & enddo do k = 2, kmpbl do i = 1, im - if(.not.flg(i)) then + 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))*