diff --git a/physics/mfpbltq.f b/physics/mfpbltq.f index b906052cd..c4333290b 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 @@ -384,6 +390,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 +441,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 +462,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..3c54b0bda 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 @@ -457,6 +463,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 +518,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 +541,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/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 3801e684f..ea92fda7f 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -149,23 +149,22 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & & dh, dhh, dp, & dq, dqsdp, dqsdt, dt, & dt2, dtmax, dtmin, - & dxcrtas, dxcrtuf, + & 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, & 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), @@ -225,7 +223,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) @@ -249,8 +247,11 @@ 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) ! c cloud water ! real(kind=kind_phys) tvo(im,km) @@ -365,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 @@ -392,6 +394,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & c0(i) = c0s 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 @@ -1192,6 +1195,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 +1213,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 +1466,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 +1643,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 +1935,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 +1962,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 +2096,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 +2118,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 +2514,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 +2610,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 @@ -2741,41 +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) = 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 : @@ -2815,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 !! @@ -2832,18 +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 = 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) @@ -2854,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 diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 0e11ed49c..24e01b040 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -102,15 +102,16 @@ 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, + & gamma, pprime, betaw, tauadv, & qlk, qrch, qs, & rfact, shear, tfac, & val, val1, val2, @@ -127,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), @@ -172,16 +172,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=2.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 +196,11 @@ 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) ! c cloud water ! real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), tvo(im,km), @@ -290,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 @@ -337,6 +342,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 @@ -996,6 +1010,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 +1028,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 +1210,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 +1378,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 +1585,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 +1597,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)) @@ -1775,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 @@ -1810,10 +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) = advfac(i)*betaw*rho*wc(i) endif enddo ! @@ -2145,7 +2162,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))*