From 54a57baa8d8874a09081cd3b51fc9dc4d53ad4e8 Mon Sep 17 00:00:00 2001 From: "Chunxi.Zhang-NOAA" Date: Tue, 15 Mar 2022 14:28:24 +0000 Subject: [PATCH 1/2] 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))* From e8dc7233804baa920a3c044a39d29a56d9b18930 Mon Sep 17 00:00:00 2001 From: "Chunxi.Zhang-NOAA" Date: Tue, 29 Mar 2022 17:35:45 +0000 Subject: [PATCH 2/2] The second updates for the saSAS cumulus scheme to improve the TC intensity forecast and a bug fix related to SL sedimentation of graupel in the Thompson scheme --- physics/mfpbltq.f | 5 +- physics/mfscuq.f | 5 +- physics/module_mp_thompson.F90 | 9 ++- physics/samfdeepcnv.f | 129 +++++++++------------------------ physics/samfshalcnv.f | 98 +++++++++---------------- 5 files changed, 78 insertions(+), 168 deletions(-) diff --git a/physics/mfpbltq.f b/physics/mfpbltq.f index a0788d5b7..c4333290b 100644 --- a/physics/mfpbltq.f +++ b/physics/mfpbltq.f @@ -319,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) = sqrt(wu2(i,k)) + xmf(i,k) = a1 * sqrt(wu2(i,k)) endif enddo enddo @@ -356,8 +356,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 - tem = max(a1, sigma(i)) - xmf(i,k) = scaldfunc(i) * tem * xmf(i,k) + xmf(i,k) = scaldfunc(i) * xmf(i,k) dz = zl(i,k+1) - zl(i,k) xmmx = dz / dt2 xmf(i,k) = min(xmf(i,k),xmmx) diff --git a/physics/mfscuq.f b/physics/mfscuq.f index b41ffd13e..3c54b0bda 100644 --- a/physics/mfscuq.f +++ b/physics/mfscuq.f @@ -386,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) = sqrt(wd2(i,k)) + xmfd(i,k) = ra1(i) * sqrt(wd2(i,k)) endif enddo enddo @@ -424,8 +424,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 - tem = max(ra1(i), sigma(i)) - xmfd(i,k) = scaldfunc(i) * tem * xmfd(i,k) + xmfd(i,k) = scaldfunc(i) * xmfd(i,k) dz = zl(i,k+1) - zl(i,k) xmmx = dz / dt2 xmfd(i,k) = min(xmfd(i,k),xmmx) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index c23b6d1d8..6d7327e8c 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -4067,7 +4067,14 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & do k = kte, kts, -1 vtg = 0. if (rg(k).gt. R1) then - vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g + ygra1 = alog10(max(1.E-9, rg(k))) + zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 + N0_exp = 10.**(zans1) + N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) + lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 + lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + + vtg = rhof(k)*av_g*cgg(6)*ogg3 * (1./lamg)**bv_g if (temp(k).gt. T_0) then vtgk(k) = MAX(vtg, vtrk(k)) else diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 0420fa1d2..ea92fda7f 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -149,7 +149,6 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & dh, dhh, dp, & dq, dqsdp, dqsdt, dt, & dt2, dtmax, dtmin, -! & dxcrtas, dxcrtuf, dxcrtc0, & dxcrtas, dxcrtuf, & dv1h, dv2h, dv3h, & dz, dz1, e1, edtmax, @@ -165,7 +164,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & w1, w1l, w1s, w2, & w2l, w2s, w3, w3l, & w3s, w4, w4l, w4s, - & rho, betaw, + & rho, betaw, tauadv, & xdby, xpw, xpwd, ! & xqrch, mbdt, tem, & xqrch, tem, tem1, tem2, @@ -179,8 +178,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im), real(kind=kind_phys) aa1(im), tkemean(im),clamt(im), & ps(im), del(im,km), prsl(im,km), -! & umean(im), tauadv(im), gdx(im), - & gdx(im), + & umean(im), advfac(im), gdx(im), & delhbar(im), delq(im), delq2(im), & delqbar(im), delqev(im), deltbar(im), & deltv(im), dtconv(im), edt(im), @@ -236,7 +234,6 @@ 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 @@ -254,7 +251,6 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & 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 @@ -370,6 +366,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & xpwav(i)= 0. xpwev(i)= 0. vshear(i) = 0. + advfac(i) = 0. rainevap(i) = 0. gdx(i) = sqrt(garea(i)) enddo @@ -398,15 +395,6 @@ subroutine samfdeepcnv_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**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 @@ -1028,33 +1016,6 @@ 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 @@ -2796,42 +2757,40 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif ! !> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. -! do i= 1, im -! if(cnvflg(i)) then -! sumx(i) = 0. -! umean(i) = 0. -! endif -! enddo -! do k = 2, km1 -! do i = 1, im -! if(cnvflg(i)) then -! if(k >= kbcon1(i) .and. k < ktcon1(i)) then -! dz = zi(i,k) - zi(i,k-1) -! tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)) -! umean(i) = umean(i) + tem * dz -! sumx(i) = sumx(i) + dz -! endif -! endif -! enddo -! enddo -! do i= 1, im -! if(cnvflg(i)) then -! umean(i) = umean(i) / sumx(i) -! umean(i) = max(umean(i), 1.) -! tauadv(i) = gdx(i) / umean(i) -! endif -! enddo + do i= 1, im + if(cnvflg(i)) then + sumx(i) = 0. + umean(i) = 0. + endif + enddo + do k = 2, km1 + do i = 1, im + if(cnvflg(i)) then + if(k >= kbcon1(i) .and. k < ktcon1(i)) then + dz = zi(i,k) - zi(i,k-1) + tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)) + umean(i) = umean(i) + tem * dz + sumx(i) = sumx(i) + dz + endif + endif + enddo + enddo + do i= 1, im + if(cnvflg(i)) then + umean(i) = umean(i) / sumx(i) + umean(i) = max(umean(i), 1.) + tauadv = gdx(i) / umean(i) + advfac(i) = tauadv / dtconv(i) + advfac(i) = min(advfac(i), 1.) + endif + enddo !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) -! tfac = tauadv(i) / dtconv(i) -! tfac = min(tfac, 1.) -! xmb(i) = tfac*betaw*rho*wc(i) -! xmb(i) = betaw*rho*wc(i) - xmb(i) = rho*wc(i) + xmb(i) = advfac(i)*betaw*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 : @@ -2871,10 +2830,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & !! !! Again when dtconv is larger than tauadv, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. if(asqecflg(i)) then -! tfac = tauadv(i) / dtconv(i) -! tfac = min(tfac, 1.) -! xmb(i) = -tfac * fld(i) / xk(i) - xmb(i) = -fld(i) / xk(i) + xmb(i) = -advfac(i) * fld(i) / xk(i) endif enddo !! @@ -2888,19 +2844,6 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & !! ! !> - For scale-aware parameterization, the updraft fraction (sigmagfm) is first computed as a function of the lateral entrainment rate at cloud base (see Han et al.'s (2017) \cite han_et_al_2017 equation 4 and 5), following the study by Grell and Freitas (2014) \cite grell_and_freitas_2014. - if(hwrf_samfdeep) then - 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) - sigmagfm(i) = max(sigmagfm(i), 0.001) - sigmagfm(i) = min(sigmagfm(i), 0.999) - endif - enddo - else do i = 1, im if(cnvflg(i)) then tem = min(max(xlamue(i,kbcon(i)), 7.e-5), 3.e-4) @@ -2911,7 +2854,6 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & sigmagfm(i) = min(sigmagfm(i), 0.999) endif enddo - endif ! !> - Then, calculate the reduction factor (scaldfunc) of the vertical convective eddy transport of mass flux as a function of updraft fraction from the studies by Arakawa and Wu (2013) \cite arakawa_and_wu_2013 (also see Han et al.'s (2017) \cite han_et_al_2017 equation 1 and 2). The final cloud base mass flux with scale-aware parameterization is obtained from the mass flux when sigmagfm << 1, multiplied by the reduction factor (Han et al.'s (2017) \cite han_et_al_2017 equation 2). do i = 1, im @@ -2922,12 +2864,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & else scaldfunc(i) = 1.0 endif - 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) = xmb(i) * scaldfunc(i) xmb(i) = min(xmb(i),xmbmax(i)) endif enddo diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 68b12d169..24e01b040 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -111,7 +111,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & es, etah, h1, shevf, ! & evfact, evfactl, & fact1, fact2, factor, dthk, - & gamma, pprime, betaw, + & gamma, pprime, betaw, tauadv, & qlk, qrch, qs, & rfact, shear, tfac, & val, val1, val2, @@ -128,8 +128,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) aa1(im), cina(im), & tkemean(im), clamt(im), & ps(im), del(im,km), prsl(im,km), -! & umean(im), tauadv(im), gdx(im), - & gdx(im), + & umean(im), advfac(im), gdx(im), & delhbar(im), delq(im), delq2(im), & delqbar(im), delqev(im), deltbar(im), ! & deltv(im), dtconv(im), edt(im), @@ -180,7 +179,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & parameter(dthk=25.,sfclfac=0.2,rhcrt=0.75) parameter(cinpcrmx=180.,cinpcrmn=120.) ! shevf is an enhancing evaporation factor for shallow convection - parameter(cinacrmx=-120.,shevf=1.0) + parameter(cinacrmx=-120.,shevf=2.0) parameter(dtmax=10800.,dtmin=600.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) parameter(betaw=.03,dxcrt=15.e3,dxcrtc0=9.e3) @@ -201,7 +200,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & 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 @@ -296,6 +294,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & aa1(i) = 0. cina(i) = 0. ! vshear(i) = 0. + advfac(i) = 0. gdx(i) = sqrt(garea(i)) xmb(i) = 0. scaldfunc(i)=-1.0 ! wang initialized @@ -904,33 +903,6 @@ 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 @@ -1821,31 +1793,33 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo ! !> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. -! do i= 1, im -! if(cnvflg(i)) then -! sumx(i) = 0. -! umean(i) = 0. -! endif -! enddo -! do k = 2, km1 -! do i = 1, im -! if(cnvflg(i)) then -! if(k >= kbcon1(i) .and. k < ktcon1(i)) then -! dz = zi(i,k) - zi(i,k-1) -! tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)) -! umean(i) = umean(i) + tem * dz -! sumx(i) = sumx(i) + dz -! endif -! endif -! enddo -! enddo -! do i= 1, im -! if(cnvflg(i)) then -! umean(i) = umean(i) / sumx(i) -! umean(i) = max(umean(i), 1.) -! tauadv(i) = gdx(i) / umean(i) -! endif -! enddo + do i= 1, im + if(cnvflg(i)) then + sumx(i) = 0. + umean(i) = 0. + endif + enddo + do k = 2, km1 + do i = 1, im + if(cnvflg(i)) then + if(k >= kbcon1(i) .and. k < ktcon1(i)) then + dz = zi(i,k) - zi(i,k-1) + tem = sqrt(u1(i,k)*u1(i,k)+v1(i,k)*v1(i,k)) + umean(i) = umean(i) + tem * dz + sumx(i) = sumx(i) + dz + endif + endif + enddo + enddo + do i= 1, im + if(cnvflg(i)) then + umean(i) = umean(i) / sumx(i) + umean(i) = max(umean(i), 1.) + tauadv = gdx(i) / umean(i) + advfac(i) = tauadv / dtconv(i) + advfac(i) = min(advfac(i), 1.) + endif + enddo c c compute cloud base mass flux as a function of the mean c updraft velcoity @@ -1856,11 +1830,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(cnvflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) -! tfac = tauadv(i) / dtconv(i) -! tfac = min(tfac, 1.) -! xmb(i) = tfac*betaw*rho*wc(i) -! xmb(i) = betaw*rho*wc(i) - xmb(i) = rho*wc(i) + xmb(i) = advfac(i)*betaw*rho*wc(i) endif enddo ! @@ -1868,7 +1838,6 @@ 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) @@ -1886,8 +1855,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & else scaldfunc(i) = 1.0 endif - tem = max(betaw, sigmagfm(i)) - xmb(i) = tem * xmb(i) * scaldfunc(i) + xmb(i) = xmb(i) * scaldfunc(i) xmb(i) = min(xmb(i),xmbmax(i)) endif enddo