diff --git a/physics/mfpbltq.f b/physics/mfpbltq.f index fb775e2e1..a93862a41 100644 --- a/physics/mfpbltq.f +++ b/physics/mfpbltq.f @@ -11,7 +11,7 @@ module mfpbltq_mod !> @{ 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, + & gdx,hpbl,kpbl,vpert,buo,wush,tkemean,vez0fun,xmf, & tcko,qcko,ucko,vcko,xlamueq,a1) ! use machine , only : kind_phys @@ -33,7 +33,8 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, & plyr(im,km),pix(im,km),thlx(im,km), & thvx(im,km),zl(im,km), zm(im,km), & gdx(im), hpbl(im), vpert(im), - & buo(im,km), xmf(im,km), + & buo(im,km), wush(im,km), + & tkemean(im),vez0fun(im),xmf(im,km), & tcko(im,km),qcko(im,km,ntrac1), & ucko(im,km),vcko(im,km), & xlamueq(im,km-1) @@ -44,8 +45,8 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, integer kpblx(im), kpbly(im) ! real(kind=kind_phys) dt2, dz, ce0, - & cm, cq, - & factor, gocp, + & cm, cq, tkcrt, + & factor, gocp, cmxfac, & g, b1, f1, & bb1, bb2, & alp, vpertmax,a1, pgcon, @@ -59,7 +60,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, ! real(kind=kind_phys) rbdn(im), rbup(im), hpblx(im), & xlamue(im,km-1), xlamuem(im,km-1) - real(kind=kind_phys) delz(im), xlamax(im) + real(kind=kind_phys) delz(im), xlamax(im), ce0t(im) ! real(kind=kind_phys) wu2(im,km), thlu(im,km), & qtx(im,km), qtu(im,km) @@ -73,7 +74,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,cq=1.3) + parameter(ce0=0.4,cm=1.0,cq=1.0,tkcrt=2.,cmxfac=5.) 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) @@ -112,13 +113,27 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, enddo ! !> - Compute entrainment rate +! +! if tkemean>tkcrt, ce0t=sqrt(tkemean/tkcrt)*ce0 +! + do i=1,im + if(cnvflg(i)) then + ce0t(i) = ce0 * vez0fun(i) + if(tkemean(i) > tkcrt) then + tem = sqrt(tkemean(i)/tkcrt) + tem1 = min(tem, cmxfac) + tem2 = tem1 * ce0 + ce0t(i) = max(ce0t(i), tem2) + endif + endif + enddo ! do i=1,im if(cnvflg(i)) then k = kpbl(i) / 2 k = max(k, 1) delz(i) = zl(i,k+1) - zl(i,k) - xlamax(i) = ce0 / delz(i) + xlamax(i) = ce0t(i) / delz(i) endif enddo ! @@ -129,7 +144,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, ptem = 1./(zm(i,k)+delz(i)) tem = max((hpbl(i)-zm(i,k)+delz(i)) ,delz(i)) ptem1 = 1./tem - xlamue(i,k) = ce0 * (ptem+ptem1) + xlamue(i,k) = ce0t(i) * (ptem+ptem1) else xlamue(i,k) = xlamax(i) endif @@ -210,11 +225,13 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, do i = 1, im if(cnvflg(i)) then dz = zm(i,k) - zm(i,k-1) - tem = 0.25*bb1*(xlamue(i,k)+xlamue(i,k-1))*dz - tem1 = bb2 * buo(i,k) * dz + tem = 0.25*bb1*(xlamue(i,k-1)+xlamue(i,k))*dz + tem1 = max(wu2(i,k-1), 0.) + tem1 = bb2 * buo(i,k) - wush(i,k) * sqrt(tem1) + tem2 = tem1 * dz ptem = (1. - tem) * wu2(i,k-1) ptem1 = 1. + tem - wu2(i,k) = (ptem + tem1) / ptem1 + wu2(i,k) = (ptem + tem2) / ptem1 endif enddo enddo @@ -271,7 +288,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, k = kpbl(i) / 2 k = max(k, 1) delz(i) = zl(i,k+1) - zl(i,k) - xlamax(i) = ce0 / delz(i) + xlamax(i) = ce0t(i) / delz(i) endif enddo ! @@ -283,7 +300,7 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, ptem = 1./(zm(i,k)+delz(i)) tem = max((hpbl(i)-zm(i,k)+delz(i)) ,delz(i)) ptem1 = 1./tem - xlamue(i,k) = ce0 * (ptem+ptem1) + xlamue(i,k) = ce0t(i) * (ptem+ptem1) else xlamue(i,k) = xlamax(i) endif diff --git a/physics/mfscuq.f b/physics/mfscuq.f index cafa61b55..d690dce05 100644 --- a/physics/mfscuq.f +++ b/physics/mfscuq.f @@ -12,7 +12,7 @@ module mfscuq_mod 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, + & krad,mrad,radmin,buo,wush,tkemean,vez0fun,xmfd, & tcdo,qcdo,ucdo,vcdo,xlamdeq,a1) ! use machine , only : kind_phys @@ -38,9 +38,10 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, & gdx(im), & zl(im,km), zm(im,km), & thetae(im,km), radmin(im), - & buo(im,km), xmfd(im,km), - & tcdo(im,km), qcdo(im,km,ntrac1), - & ucdo(im,km), vcdo(im,km), + & buo(im,km), wush(im,km), + & tkemean(im),vez0fun(im),xmfd(im,km), + & tcdo(im,km),qcdo(im,km,ntrac1), + & ucdo(im,km),vcdo(im,km), & xlamdeq(im,km-1) ! ! local variables and arrays @@ -51,6 +52,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, ! real(kind=kind_phys) dt2, dz, ce0, & cm, cq, + & tkcrt, cmxfac, & gocp, factor, g, tau, & b1, f1, bb1, bb2, & a1, a2, @@ -67,7 +69,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, & qtx(im,km), qtd(im,km), & thlvd(im), hrad(im), xlamde(im,km-1), & xlamdem(im,km-1), ra1(im) - real(kind=kind_phys) delz(im), xlamax(im) + real(kind=kind_phys) delz(im), xlamax(im), ce0t(im) ! real(kind=kind_phys) xlamavg(im), sigma(im), & scaldfunc(im), sumx(im) @@ -80,7 +82,8 @@ 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,cq=1.3,pgcon=0.55) + parameter(ce0=0.4,cm=1.0,cq=1.0,pgcon=0.55) + parameter(tkcrt=2.,cmxfac=5.) parameter(qmin=1.e-8,qlmin=1.e-12) parameter(b1=0.45,f1=0.15) parameter(a2=0.5) @@ -185,13 +188,27 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, !! ! !> - Compute entrainment rate +! +! if tkemean>tkcrt, ce0t=sqrt(tkemean/tkcrt)*ce0 +! + do i=1,im + if(cnvflg(i)) then + ce0t(i) = ce0 * vez0fun(i) + if(tkemean(i) > tkcrt) then + tem = sqrt(tkemean(i)/tkcrt) + tem1 = min(tem, cmxfac) + tem2 = tem1 * ce0 + ce0t(i) = max(ce0t(i), tem2) + endif + endif + enddo ! do i=1,im if(cnvflg(i)) then k = mrad(i) + (krad(i)-mrad(i)) / 2 k = max(k, mrad(i)) delz(i) = zl(i,k+1) - zl(i,k) - xlamax(i) = ce0 / delz(i) + xlamax(i) = ce0t(i) / delz(i) endif enddo ! @@ -206,7 +223,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, endif tem = max((hrad(i)-zm(i,k)+delz(i)) ,delz(i)) ptem1 = 1./tem - xlamde(i,k) = ce0 * (ptem+ptem1) + xlamde(i,k) = ce0t(i) * (ptem+ptem1) else xlamde(i,k) = xlamax(i) endif @@ -289,10 +306,12 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, if(cnvflg(i) .and. k < krad1(i)) then dz = zm(i,k+1) - zm(i,k) tem = 0.25*bb1*(xlamde(i,k)+xlamde(i,k+1))*dz - tem1 = bb2 * buo(i,k+1) * dz + tem1 = max(wd2(i,k+1), 0.) + tem1 = bb2*buo(i,k+1) - wush(i,k+1)*sqrt(tem1) + tem2 = tem1 * dz ptem = (1. - tem) * wd2(i,k+1) ptem1 = 1. + tem - wd2(i,k) = (ptem + tem1) / ptem1 + wd2(i,k) = (ptem + tem2) / ptem1 endif enddo enddo @@ -334,7 +353,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, k = mrad(i) + (krad(i)-mrad(i)) / 2 k = max(k, mrad(i)) delz(i) = zl(i,k+1) - zl(i,k) - xlamax(i) = ce0 / delz(i) + xlamax(i) = ce0t(i) / delz(i) endif enddo ! @@ -349,7 +368,7 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, endif tem = max((hrad(i)-zm(i,k)+delz(i)) ,delz(i)) ptem1 = 1./tem - xlamde(i,k) = ce0 * (ptem+ptem1) + xlamde(i,k) = ce0t(i) * (ptem+ptem1) else xlamde(i,k) = xlamax(i) endif diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index b828c9ab0..6a4ef5e02 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -469,11 +469,11 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, & end if if (mpirank==mpiroot) then if (is_aerosol_aware) then - write (0,'(a)') 'Using aerosol-aware version of Thompson microphysics' + write (*,'(a)') 'Using aerosol-aware version of Thompson microphysics' else if(merra2_aerosol_aware) then - write (0,'(a)') 'Using merra2 aerosol-aware version of Thompson microphysics' + write (*,'(a)') 'Using merra2 aerosol-aware version of Thompson microphysics' else - write (0,'(a)') 'Using non-aerosol-aware version of Thompson microphysics' + write (*,'(a)') 'Using non-aerosol-aware version of Thompson microphysics' end if end if @@ -896,22 +896,22 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, & !> - Call table_ccnact() to read a static file containing CCN activation of aerosols. The !! data were created from a parcel model by Feingold & Heymsfield with !! further changes by Eidhammer and Kriedenweis - if (mpirank==mpiroot) write(0,*) ' calling table_ccnAct routine' + if (mpirank==mpiroot) write(*,*) ' calling table_ccnAct routine' call table_ccnAct(errmsg,errflg) if (.not. errflg==0) return !> - Call table_efrw() and table_efsw() to creat collision efficiency table !! between rain/snow and cloud water - if (mpirank==mpiroot) write(0,*) ' creating qc collision eff tables' + if (mpirank==mpiroot) write(*,*) ' creating qc collision eff tables' call table_Efrw call table_Efsw !> - Call table_dropevap() to creat rain drop evaporation table - if (mpirank==mpiroot) write(0,*) ' creating rain evap table' + if (mpirank==mpiroot) write(*,*) ' creating rain evap table' call table_dropEvap !> - Call qi_aut_qs() to create conversion of some ice mass into snow category - if (mpirank==mpiroot) write(0,*) ' creating ice converting to snow table' + if (mpirank==mpiroot) write(*,*) ' creating ice converting to snow table' call qi_aut_qs call cpu_time(etime) @@ -942,7 +942,7 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, & call cpu_time(stime) !> - Call qr_acr_qg() to create rain collecting graupel & graupel collecting rain table - if (mpirank==mpiroot) write(0,*) ' creating rain collecting graupel table' + if (mpirank==mpiroot) write(*,*) ' creating rain collecting graupel table' call cpu_time(stime) call qr_acr_qg call cpu_time(etime) @@ -956,7 +956,7 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, & if (mpirank==mpiroot) print '("Computing rain collecting snow table took ",f10.3," seconds.")', etime-stime !> - Call freezeh2o() to create cloud water and rain freezing (Bigg, 1953) table - if (mpirank==mpiroot) write(0,*) ' creating freezing of water drops table' + if (mpirank==mpiroot) write(*,*) ' creating freezing of water drops table' call cpu_time(stime) call freezeH2O(threads) call cpu_time(etime) @@ -969,7 +969,7 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, & endif if_not_iiwarm - if (mpirank==mpiroot) write(0,*) ' ... DONE microphysical lookup tables' + if (mpirank==mpiroot) write(*,*) ' ... DONE microphysical lookup tables' endif if_micro_init diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 3029398e9..079218f5a 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -2081,7 +2081,8 @@ subroutine progcld_thompson_wsm6 & ! --- constant values real (kind=kind_phys), parameter :: xrc3 = 100. - + real (kind=kind_phys), parameter :: snow2ice = 0.25 + real (kind=kind_phys), parameter :: coef_t = 0.025 ! !===> ... begin here @@ -2097,7 +2098,7 @@ subroutine progcld_thompson_wsm6 & rei (i,k) = re_ice(i,k) rer (i,k) = rrain_def ! default rain radius to 1000 micron res (i,k) = re_snow(i,K) - tem2d (i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) ) + tem2d (i,k) = min(1.0, max( 0.0, (con_ttp-tlyr(i,k))*coef_t)) clwf(i,k) = 0.0 enddo enddo @@ -2138,12 +2139,14 @@ subroutine progcld_thompson_wsm6 & if(tem1 > 1.e-12 .and. clw(i,k,ntcw) < 1.e-12) & rew(i,k)=reliq_def tem2 = cnvw(i,k)*tem2d(i,k) - cip(i,k) = max(0.0, (clw(i,k,ntiw) + tem2 ) - & *gfac * delp(i,k)) + cip(i,k) = max(0.0, (clw(i,k,ntiw) + + & snow2ice*clw(i,k,ntsw) + tem2) * + & gfac * delp(i,k)) if(tem2 > 1.e-12 .and. clw(i,k,ntiw) < 1.e-12) & rei(i,k)=reice_def crp(i,k) = max(0.0, clw(i,k,ntrw) * gfac * delp(i,k)) - csp(i,k) = max(0.0, clw(i,k,ntsw) * gfac * delp(i,k)) + csp(i,k) = max(0.0, (1.-snow2ice)*clw(i,k,ntsw) * + & gfac * delp(i,k)) enddo enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 0d4f9fd0f..8a36fe34c 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -191,7 +191,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & pwavo(im), pwevo(im), mbdt(im), & qcdo(im,km), qcond(im), qevap(im), & rntot(im), vshear(im), xaa0(im), - & xlamd(im), xk(im), cina(im), + & xlamd(im), xlamdet(im),xlamddt(im), + & cxlamet(im), cxlamdt(im), + & xk(im), cina(im), & xmb(im), xmbmax(im), xpwav(im), & xpwev(im), xlamx(im), delebar(im,ntr), ! & xpwev(im), delebar(im,ntr), @@ -201,13 +203,12 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & cj real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn, & cinacr, cinacrmx, cinacrmn, - & sfclfac, rhcrt + & sfclfac, rhcrt, + & tkcrt, cmxfac cj ! ! parameters for updraft velocity calculation - real(kind=kind_phys) bet1, cd1, f1, gam1, -! & bb1, bb2 - & bb1, bb2, wucb + real(kind=kind_phys) bb1, bb2, csmf, wucb ! ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), @@ -229,7 +230,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! Until a realistic Nccn is provided, Nccns are assumed ! as Nccn=100 for sea and Nccn=1000 for land ! - parameter(cm=1.0,cq=1.3) + parameter(cm=1.0,cq=1.0) ! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c) parameter(clamd=0.03,tkemx=0.65,tkemn=0.05) parameter(clamca=0.03) @@ -238,7 +239,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & parameter(cinpcrmx=180.,cinpcrmn=120.) ! parameter(cinacrmx=-120.,cinacrmn=-120.) parameter(cinacrmx=-120.,cinacrmn=-80.) - parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) + parameter(bb1=4.0,bb2=0.8,csmf=0.2) + parameter(tkcrt=2.,cmxfac=15.) parameter(betaw=.03) ! @@ -254,7 +256,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! ! for updraft velocity calculation real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km), - & wc(im) + & wush(im,km), wc(im) ! ! for updraft fraction & scale-aware function real(kind=kind_phys) scaldfunc(im), sigmagfm(im) @@ -577,6 +579,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! vo(i,k) = v1(i,k) * rcs(i) wu2(i,k) = 0. buo(i,k) = 0. + wush(i,k) = 0. drag(i,k) = 0. cnvwt(i,k)= 0. endif @@ -805,7 +808,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ptem1= .5*(cinpcrmx-cinpcrmn) cinpcr = cinpcrmx - ptem * ptem1 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i)) - if(tem1 > cinpcr) then + if(tem1 > cinpcr .and. + & zi(i,kbcon(i)) > hpbl(i)) then cnvflg(i) = .false. endif endif @@ -975,7 +979,25 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & endif enddo endif - +! + do i=1,im + if(cnvflg(i)) then + xlamdet(i) = xlamde + xlamddt(i) = xlamdd + cxlamet(i) = cxlame + cxlamdt(i) = cxlamd + if(tkemean(i) > tkcrt) then + tem = 1. + tkemean(i)/tkcrt + tem1 = min(tem, cmxfac) + clamt(i) = tem1 * clam + xlamdet(i) = tem1 * xlamdet(i) + xlamddt(i) = tem1 * xlamddt(i) + cxlamet(i) = tem1 * cxlamet(i) + cxlamdt(i) = tem1 * cxlamdt(i) + endif + endif + enddo +! else ! if(do_ca .and. ca_entr)then @@ -995,6 +1017,15 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & endif enddo endif +! + do i=1,im + if(cnvflg(i)) then + xlamdet(i) = xlamde + xlamddt(i) = xlamdd + cxlamet(i) = cxlame + cxlamdt(i) = cxlamd + endif + enddo ! endif !(.not. hwrf_samfdeep .and. ntk > 0) ! @@ -1083,7 +1114,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do i=1,im if(cnvflg(i) .and. & (k > kbcon(i) .and. k < kmax(i))) then - tem = cxlame * frh(i,k) * fent2(i,k) + tem = cxlamet(i) * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem endif enddo @@ -1093,9 +1124,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do i=1,im if(cnvflg(i) .and. & (k > kbcon(i) .and. k < kmax(i))) then - tem = cxlame * frh(i,k) * fent2(i,k) + tem = cxlamet(i) * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem - tem1 = cxlamd * frh(i,k) + tem1 = cxlamdt(i) * frh(i,k) xlamud(i,k) = xlamud(i,k) + tem1 endif enddo @@ -1531,6 +1562,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & buo(i,k) = buo(i,k) + grav * fv * & max(val,(qeso(i,k) - qo(i,k))) drag(i,k) = max(xlamue(i,k),xlamud(i,k)) +! + tem = ((uo(i,k)-uo(i,k-1))/dz)**2 + tem = tem+((vo(i,k)-vo(i,k-1))/dz)**2 + wush(i,k) = csmf * sqrt(tem) +! endif ! endif @@ -1701,8 +1737,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! compute updraft velocity square(wu2) !> - Calculate updraft velocity square(wu2) according to Han et al.'s (2017) \cite han_et_al_2017 equation 7. ! - bb1 = 4.0 - bb2 = 0.8 if (hwrf_samfdeep) then do i = 1, im if (cnvflg(i)) then @@ -1723,11 +1757,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if (cnvflg(i)) then if(k > kbcon1(i) .and. k < ktcon(i)) then dz = zi(i,k) - zi(i,k-1) - tem = 0.25 * bb1 * (drag(i,k)+drag(i,k-1)) * dz - tem1 = 0.5 * bb2 * (buo(i,k)+buo(i,k-1)) * dz + tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz + tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) + tem2 = wush(i,k) * sqrt(wu2(i,k-1)) + tem2 = (tem1 - tem2) * dz ptem = (1. - tem) * wu2(i,k-1) ptem1 = 1. + tem - wu2(i,k) = (ptem + tem1) / ptem1 + wu2(i,k) = (ptem + tem2) / ptem1 wu2(i,k) = max(wu2(i,k), 0.) endif endif @@ -1953,11 +1989,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if (cnvflg(i) .and. k <= kmax(i)-1) then if(k < jmin(i) .and. k >= kd94(i)) then dz = zi(i,k+1) - zi(i,k) - ptem = xlamdd - xlamde + ptem = xlamddt(i) - xlamdet(i) etad(i,k) = etad(i,k+1) * (1. - ptem * dz) else if(k < kd94(i)) then dz = zi(i,k+1) - zi(i,k) - ptem = xlamd(i) + xlamdd - xlamde + ptem = xlamd(i) + xlamddt(i) - xlamdet(i) etad(i,k) = etad(i,k+1) * (1. - ptem * dz) endif endif @@ -1996,11 +2032,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if (cnvflg(i) .and. k < jmin(i)) then dz = zi(i,k+1) - zi(i,k) if(k >= kd94(i)) then - tem = xlamde * dz - tem1 = 0.5 * xlamdd * dz + tem = xlamdet(i) * dz + tem1 = 0.5 * xlamddt(i) * dz else - tem = xlamde * dz - tem1 = 0.5 * (xlamd(i)+xlamdd) * dz + tem = xlamdet(i) * dz + tem1 = 0.5 * (xlamd(i)+xlamddt(i)) * dz endif factor = 1. + tem - tem1 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* @@ -2024,7 +2060,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do i = 1, im if (cnvflg(i) .and. k < jmin(i)) then dz = zi(i,k+1) - zi(i,k) - tem = 0.5 * xlamde * dz + tem = 0.5 * xlamdet(i) * dz tem = cq * tem factor = 1. + tem ecdo(i,k,n) = ((1.-tem)*ecdo(i,k+1,n)+tem* @@ -2045,7 +2081,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! detad = etad(i,k+1) - etad(i,k) cj dz = zi(i,k+1) - zi(i,k) - tem = 0.5 * xlamde * dz + tem = 0.5 * xlamdet(i) * dz tem = cq * tem factor = 1. + tem qcdo(i,k) = ((1.-tem)*qrcdo(i,k+1)+tem* @@ -2184,11 +2220,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & tem1 = 0.5 * (xlamud(i,k)+xlamud(i,k-1)) c if(k <= kd94(i)) then - ptem = xlamde - ptem1 = xlamd(i)+xlamdd + ptem = xlamdet(i) + ptem1 = xlamd(i)+xlamddt(i) else - ptem = xlamde - ptem1 = xlamdd + ptem = xlamdet(i) + ptem1 = xlamddt(i) endif factor = grav / dp @@ -2670,11 +2706,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if (asqecflg(i) .and. k < jmin(i)) then dz = zi(i,k+1) - zi(i,k) if(k >= kd94(i)) then - tem = xlamde * dz - tem1 = 0.5 * xlamdd * dz + tem = xlamdet(i) * dz + tem1 = 0.5 * xlamddt(i) * dz else - tem = xlamde * dz - tem1 = 0.5 * (xlamd(i)+xlamdd) * dz + tem = xlamdet(i) * dz + tem1 = 0.5 * (xlamd(i)+xlamddt(i)) * dz endif factor = 1. + tem - tem1 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* @@ -2694,7 +2730,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! detad = etad(i,k+1) - etad(i,k) cj dz = zi(i,k+1) - zi(i,k) - tem = 0.5 * xlamde * dz + tem = 0.5 * xlamdet(i) * dz tem = cq * tem factor = 1. + tem qcdo(i,k) = ((1.-tem)*qrcd(i,k+1)+tem* diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 2f8b188a5..a7682342f 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -110,7 +110,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & cm, cq, & es, etah, h1, shevf, ! & evfact, evfactl, - & fact1, fact2, factor, dthk, + & fact1, fact2, factor, + & cthk, cthkmn, dthk, & gamma, pprime, betaw, tauadv, & qlk, qrch, qs, & rfact, shear, tfac, @@ -122,7 +123,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & ptem, ptem1 ! integer kb(im), kb1(im), kbcon(im), kbcon1(im), - & ktcon(im), ktcon1(im), ktconn(im), + & ktcon(im), ktcon1(im), & kbm(im), kmax(im) ! real(kind=kind_phys) aa1(im), cina(im), @@ -147,13 +148,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn, & cinacr, cinacrmx, cinacrmn, - & sfclfac, rhcrt + & sfclfac, rhcrt, + & tkcrt, cmxfac ! ! parameters for updraft velocity calculation - real(kind=kind_phys) bet1, cd1, f1, gam1, -! & bb1, bb2 - & bb1, bb2, wucb - + real(kind=kind_phys) bb1, bb2, csmf, wucb cc ! parameters for prognostic sigma closure @@ -179,16 +178,19 @@ 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,cq=1.3) + parameter(cm=1.0,cq=1.0) ! 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(cthk=200.,cthkmn=0.,dthk=25.) + parameter(sfclfac=0.2,rhcrt=0.75) parameter(cinpcrmx=180.,cinpcrmn=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(bb1=4.0,bb2=0.8,csmf=0.2) + parameter(tkcrt=2.,cmxfac=15.) +! parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) parameter(betaw=.03,dxcrtc0=9.e3) parameter(h1=0.33333333) ! progsigma @@ -206,7 +208,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! ! for updraft velocity calculation real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km), - & wc(im) + & wush(im,km), wc(im) ! ! for updraft fraction & scale-aware function real(kind=kind_phys) scaldfunc(im), sigmagfm(im) @@ -304,7 +306,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & rn(i)=0. kbcon(i)=km ktcon(i)=1 - ktconn(i)=1 kb(i)=km pdot(i) = 0. qlko_ktcon(i) = 0. @@ -331,7 +332,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & rn(i)=0. kbcon(i)=km ktcon(i)=1 - ktconn(i)=1 kb(i)=km pdot(i) = 0. qlko_ktcon(i) = 0. @@ -512,6 +512,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! vo(i,k) = v1(i,k) * rcs(i) wu2(i,k) = 0. buo(i,k) = 0. + wush(i,k) = 0. drag(i,k) = 0. cnvwt(i,k) = 0. endif @@ -746,7 +747,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & cinpcr = cinpcrmx - ptem * ptem1 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i)) - if(tem1 > cinpcr) then + if(tem1 > cinpcr .and. + & zi(i,kbcon(i)) > hpbl(i)) then cnvflg(i) = .false. endif endif @@ -893,6 +895,16 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif endif enddo +! + do i=1,im + if(cnvflg(i)) then + if(tkemean(i) > tkcrt) then + tem = 1. + tkemean(i)/tkcrt + tem1 = min(tem, cmxfac) + clamt(i) = tem1 * clam + endif + endif + enddo ! else ! @@ -970,7 +982,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & eta(i,k) = eta(i,k-1) * (1 + ptem * dz) if(eta(i,k) <= 0.) then kmax(i) = k - ktconn(i) = k kbm(i) = min(kbm(i),kmax(i)) flg(i) = .false. endif @@ -1200,7 +1211,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & !> - Calculate the cloud top as the first level where parcel buoyancy becomes negative; the maximum possible value is at \f$p=0.7p_{sfc}\f$. do i = 1, im flg(i) = cnvflg(i) - if(flg(i)) ktcon(i) = kbm(i) + if(flg(i)) ktcon(i) = 1 enddo do k = 2, km1 do i=1,im @@ -1213,6 +1224,18 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo enddo c +c turn off shallow convection if cloud depth is larger than cthk or less than cthkmn +c + do i = 1, im + if(cnvflg(i)) then + tem = pfld(i,kbcon(i))-pfld(i,ktcon(i)) + if(tem > cthk .or. tem < cthkmn) then + cnvflg(i) = .false. + endif + endif + enddo + +c c specify upper limit of mass flux at cloud base c !> - Calculate the maximum value of the cloud base mass flux using the CFL-criterion-based formula of Han and Pan (2011) \cite han_and_pan_2011, equation 7. @@ -1286,6 +1309,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & buo(i,k) = buo(i,k) + grav * fv * & max(val,(qeso(i,k) - qo(i,k))) drag(i,k) = max(xlamue(i,k),xlamud(i)) +! + tem = ((uo(i,k)-uo(i,k-1))/dz)**2 + tem = tem+((vo(i,k)-vo(i,k-1))/dz)**2 + wush(i,k) = csmf * sqrt(tem) +! endif ! endif @@ -1445,9 +1473,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! ! compute updraft velocity square(wu2) !> - Calculate updraft velocity square(wu2) according to Han et al.'s (2017) \cite han_et_al_2017 equation 7. -! - bb1 = 4.0 - bb2 = 0.8 ! if (hwrf_samfshal) then do i = 1, im @@ -1468,11 +1493,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (cnvflg(i)) then if(k > kbcon1(i) .and. k < ktcon(i)) then dz = zi(i,k) - zi(i,k-1) - tem = 0.25 * bb1 * (drag(i,k)+drag(i,k-1)) * dz - tem1 = 0.5 * bb2 * (buo(i,k)+buo(i,k-1)) * dz + tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz + tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) + tem2 = wush(i,k) * sqrt(wu2(i,k-1)) + tem2 = (tem1 - tem2) * dz ptem = (1. - tem) * wu2(i,k-1) ptem1 = 1. + tem - wu2(i,k) = (ptem + tem1) / ptem1 + wu2(i,k) = (ptem + tem2) / ptem1 wu2(i,k) = max(wu2(i,k), 0.) endif endif diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 0387185e4..75925a5f3 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -75,7 +75,8 @@ end subroutine satmedmfvdifq_init !! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & ntiw,ntke,grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & - & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, & + & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu, & + & garea,zvfun,sigmaf, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & & tsea,heat,evap,stress,spd1,kpbl, & & prsi,del,prsl,prslk,phii,phil,delt, & @@ -112,7 +113,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & t1(:,:), q1(:,:,:), & & swh(:,:), hlw(:,:), & & xmu(:), garea(:), & - & zvfun(:), & + & zvfun(:), sigmaf(:), & & psk(:), rbsoil(:), & & zorl(:), tsea(:), & & u10m(:), v10m(:), & @@ -147,7 +148,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & integer lcld(im),kcld(im),krad(im),mrad(im) integer kx1(im), kb1(im), kpblx(im) ! - real(kind=kind_phys) tke(im,km), tkeh(im,km-1) + real(kind=kind_phys) tke(im,km), tkeh(im,km-1), e2(im,0:km) ! real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km), & qlx(im,km), thetae(im,km),thlx(im,km), @@ -159,15 +160,16 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & qstl(im,km) ! real(kind=kind_phys) dtdz1(im), gdx(im), - & phih(im), phim(im), prn(im,km-1), + & phih(im), phim(im), + & phims(im), prn(im,km-1), & rbdn(im), rbup(im), thermal(im), & ustar(im), wstar(im), hpblx(im), & ust3(im), wst3(im), rho_a(im), - & z0(im), crb(im), + & z0(im), crb(im), tkemean(im), & hgamt(im), hgamq(im), & wscale(im),vpert(im), & zol(im), sflux(im), - & tx1(im), tx2(im) + & sumx(im), tx1(im), tx2(im) ! real(kind=kind_phys) radmin(im) ! @@ -181,7 +183,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & real(kind=kind_phys) elm(im,km), ele(im,km), & ckz(im,km), chz(im,km), & diss(im,km-1),prod(im,km-1), - & bf(im,km-1), shr2(im,km-1), + & bf(im,km-1), shr2(im,km-1), wush(im,km), & xlamue(im,km-1), xlamde(im,km-1), & gotvx(im,km), rlam(im,km-1) ! @@ -205,7 +207,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & 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) + real(kind=kind_phys) sfcpbl(im), vez0fun(im) ! logical pblflg(im), sfcflg(im), flg(im) logical scuflg(im), pcnvflg(im) @@ -238,11 +240,13 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & zfac, zfmin, vk, spdk2, & tkmin, tkbmx, xkgdx, & xkinv1, xkinv2, - & zlup, zldn, bsum, cs0, + & zlup, zldn, cs0, csmf, & tem, tem1, tem2, tem3, & ptem, ptem0, ptem1, ptem2 ! real(kind=kind_phys) slfac +! + real(kind=kind_phys) vegflo, vegfup, z0lo, z0up, vc0, zc0 ! real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck ! @@ -269,8 +273,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & parameter(qlcr=3.5e-5,zstblmax=2500.) parameter(xkinv1=0.15,xkinv2=0.3) parameter(h1=0.33333333,hcrinv=250.) + parameter(vegflo=0.1,vegfup=1.0,z0lo=0.1,z0up=1.0) + parameter(vc0=1.0,zc0=1.0) parameter(ck1=0.15,ch1=0.15) - parameter(cs0=0.2) + parameter(cs0=0.4,csmf=0.5) parameter(rchck=1.5,ndt=20) if (tc_pbl == 0) then @@ -318,6 +324,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & xmfd(i,k) = 0. buou(i,k) = 0. buod(i,k) = 0. + wush(i,k) = 0. ckz(i,k) = ck1 chz(i,k) = ch1 rlmnz(i,k) = rlmn0 @@ -435,7 +442,18 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & kcld(i) = km1 endif enddo - +! +! compute a function for green vegetation fraction and surface roughness +! + do i = 1,im + tem = (sigmaf(i) - vegflo) / (vegfup - vegflo) + tem = min(max(tem, 0.), 1.) + tem1 = sqrt(tem) + ptem = (z0(i) - z0lo) / (z0up - z0lo) + ptem = min(max(ptem, 0.), 1.) + vez0fun(i) = (1. + vc0 * tem1) * (1. + zc0 * ptem) + enddo +! !> - Compute \f$\theta\f$(theta), and \f$q_l\f$(qlx), \f$\theta_e\f$(thetae), !! \f$\theta_v\f$(thvx),\f$\theta_{l,v}\f$ (thlvx) including ice water do k=1,km @@ -461,7 +479,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & svx(i,k) = cp * tvx(i,k) ptem1 = elocp * pix(i,k) * max(q1(i,k,1),qmin) thetae(i,k)= theta(i,k) + ptem1 - gotvx(i,k) = g / tvx(i,k) +! gotvx(i,k) = g / tvx(i,k) + gotvx(i,k) = g / thvx(i,k) enddo enddo @@ -726,6 +745,40 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif enddo ! +! compute mean tke within pbl +! + do i = 1, im + sumx(i) = 0. + tkemean(i) = 0. + enddo + do k = 1, kmpbl + do i = 1, im + if(k < kpbl(i)) then + dz = zi(i,k+1) - zi(i,k) + tkemean(i) = tkemean(i) + tke(i,k) * dz + sumx(i) = sumx(i) + dz + endif + enddo + enddo + do i = 1, im + if(tkemean(i) > 0. .and. sumx(i) > 0.) then + tkemean(i) = tkemean(i) / sumx(i) + endif + enddo +! +! compute wind shear term as a sink term for updraft and downdraft +! velocity +! + kps = max(kmpbl, kmscu) + do k = 2, kps + do i = 1, im + dz = zi(i,k+1) - zi(i,k) + tem = (0.5*(u1(i,k-1)-u1(i,k+1))/dz)**2 + tem1 = tem+(0.5*(v1(i,k-1)-v1(i,k+1))/dz)**2 + wush(i,k) = csmf * sqrt(tem1) + enddo + 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 @@ -758,9 +811,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & tem = 1.0 / (1. - aphi16*zol1) phih(i) = sqrt(tem) phim(i) = sqrt(phih(i)) + tem1 = 1.0 / (1. - aphi16*zol(i)) + phims(i) = sqrt(sqrt(tem1)) else phim(i) = 1. + aphi5*zol1 phih(i) = phim(i) + phims(i) = 1. + aphi5*zol(i) endif enddo ! @@ -952,14 +1008,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !! to take into account nonlocal transport by large eddies. For details of the mfpbltq subroutine, step into its documentation ::mfpbltq call mfpbltq(im,im,km,kmpbl,ntcw,ntrac1,dt2, & pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx, - & gdx,hpbl,kpbl,vpert,buou,xmf, + & gdx,hpbl,kpbl,vpert,buou,wush,tkemean,vez0fun,xmf, & tcko,qcko,ucko,vcko,xlamue,bl_upfr) !> - Call mfscuq(), which is a new mass-flux parameterization for !! stratocumulus-top-induced turbulence mixing. For details of the mfscuq subroutine, step into its documentation ::mfscuq call mfscuq(im,im,km,kmscu,ntcw,ntrac1,dt2, & scuflg,zl,zm,q1,t1,u1,v1,plyr,pix, & thlx,thvx,thlvx,gdx,thetae, - & krad,mrad,radmin,buod,xmfd, + & krad,mrad,radmin,buod,wush,tkemean,vez0fun,xmfd, & tcdo,qcdo,ucdo,vcdo,xlamde,bl_dnfr) if (tc_pbl == 1) then @@ -1055,79 +1111,44 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do k = 1, km1 do i = 1, im zlup = 0.0 - bsum = 0.0 mlenflg = .true. + e2(i,k) = max(2.*tke(i,k), 0.001) do n = k, km1 if(mlenflg) then dz = zl(i,n+1) - zl(i,n) -! tem1 = 0.5 * (thvx(i,n) + thvx(i,n+1)) -!! tem1 = 0.5 * (thlvx(i,n) + thlvx(i,n+1)) - 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)) - 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 + tem1 = 2.*gotvx(i,n+1)*(thvx(i,k)-thvx(i,n+1)) + tem2 = cs0*sqrt(e2(i,n))*sqrt(shr2(i,n)) + e2(i,n+1) = e2(i,n) + (tem1 - tem2) * dz zlup = zlup + dz - if(bsum >= tke(i,k)) then - if(ptem >= 0.) then - tem2 = max(ptem, zfmin) - else - tem2 = min(ptem, -zfmin) - endif - ptem1 = (bsum - tke(i,k)) / tem2 - zlup = zlup - ptem1 * dz + if(e2(i,n+1) < 0.) then + ptem = e2(i,n+1) / (e2(i,n+1) - e2(i,n)) + zlup = zlup - ptem * dz zlup = max(zlup, 0.) mlenflg = .false. endif endif enddo zldn = 0.0 - bsum = 0.0 mlenflg = .true. do n = k, 1, -1 if(mlenflg) then if(n == 1) then dz = zl(i,1) - tem1 = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) -! tem1 = 0.5 * (tem1 + thvx(i,n)) -!! tem1 = 0.5 * (tem1 + thlvx(i,n)) - tem3 = (u1(i,1)/dz)**2 - tem3 = tem3+(v1(i,1)/dz)**2 - tem3 = cs0*sqrt(tem3)*sqrt(tke(i,1)) + tem = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + tem1 = 2.*gotvx(i,n)*(tem-thvx(i,k)) + tem2 = ustar(i)*phims(i)/(vk*dz) + tem2 = cs0*sqrt(e2(i,n))*tem2 + e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz else dz = zl(i,n) - zl(i,n-1) - 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)) + tem1 = 2.*gotvx(i,n-1)*(thvx(i,n-1)-thvx(i,k)) + tem2 = cs0*sqrt(e2(i,n))*sqrt(shr2(i,n-1)) + e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz endif - ptem = (gotvx(i,n)*(thvx(i,k)-tem1)+tem3)*dz -! ptem = (gotvx(i,n)*(thlvx(i,k)-tem1)+tem3)*dz - bsum = bsum + ptem zldn = zldn + dz - if(bsum >= tke(i,k)) then - if(ptem >= 0.) then - tem2 = max(ptem, zfmin) - else - tem2 = min(ptem, -zfmin) - endif - ptem1 = (bsum - tke(i,k)) / tem2 - zldn = zldn - ptem1 * dz + if(e2(i,n-1) < 0.) then + ptem = e2(i,n-1) / (e2(i,n-1) - e2(i,n)) + zldn = zldn - ptem * dz zldn = max(zldn, 0.) mlenflg = .false. endif @@ -1192,6 +1213,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & elm(i,k) = sqrt( 1.0/( 1.0/(zk**2)+1.0/(rlam(i,k)**2) ) ) endif +! + if(k == 1) elm(i,k)=zk +! dz = zi(i,k+1) - zi(i,k) tem = max(gdx(i),dz) elm(i,k) = min(elm(i,k), tem) @@ -1335,8 +1359,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif ptem2 = ptem2 + ptem ! -! tem2 = stress(i)*spd1(i)/zl(i,1) - tem2 = stress(i)*ustar(i)*phim(i)/(vk*zl(i,1)) + tem2 = stress(i)*ustar(i)*phims(i)/(vk*zl(i,1)) shrp = 0.5 * (tem1 + ptem1 + ptem2 + tem2) else tem1 = -dkt(i,k-1) * bf(i,k-1) @@ -1664,8 +1687,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif if((pcnvflg(i) .or. scuflg(i)) .and. & (k >= kbx .and. k <= kmx)) then - dz = zi(i,k+1) - zi(i,k) - tem = f1(i,k) * dz + tem = f1(i,k) * del(i,k) * gravi if(f1(i,k) < 0.) tsumn(i) = tsumn(i) + tem if(f1(i,k) > 0.) tsump(i) = tsump(i) + tem endif @@ -1720,8 +1742,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo do k = 1,km do i = 1,im - dz = zi(i,k+1) - zi(i,k) - tem = f1(i,k) * dz + tem = f1(i,k) * del(i,k) * gravi if(f1(i,k) < 0.) tsumn(i) = tsumn(i) + tem if(f1(i,k) > 0.) tsump(i) = tsump(i) + tem enddo @@ -1913,7 +1934,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif if((pcnvflg(i) .or. scuflg(i)) .and. & (k >= kbx .and. k <= kmx)) then - tem = f2(i,k) * del(i,k) / grav + tem = f2(i,k) * del(i,k) * gravi if(f2(i,k) < 0.) tsumn(i) = tsumn(i) + tem if(f2(i,k) > 0.) tsump(i) = tsump(i) + tem endif @@ -1969,7 +1990,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo do k = 1,km do i = 1,im - tem = f2(i,k) * del(i,k) / grav + tem = f2(i,k) * del(i,k) * gravi if(f2(i,k) < 0.) tsumn(i) = tsumn(i) + tem if(f2(i,k) > 0.) tsump(i) = tsump(i) + tem enddo @@ -2098,7 +2119,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif if((pcnvflg(i) .or. scuflg(i)) .and. & (k >= kbx .and. k <= kmx)) then - tem = f2(i,k+is) * del(i,k) / grav + tem = f2(i,k+is) * del(i,k) * gravi if(f2(i,k+is) < 0.) tsumn(i) = tsumn(i) + tem if(f2(i,k+is) > 0.) tsump(i) = tsump(i) + tem endif @@ -2154,7 +2175,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo do k = 1,km do i = 1,im - tem = f2(i,k+is) * del(i,k) / grav + tem = f2(i,k+is) * del(i,k) * gravi if(f2(i,k+is) < 0.) tsumn(i) = tsumn(i) + tem if(f2(i,k+is) > 0.) tsump(i) = tsump(i) + tem enddo @@ -2197,7 +2218,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend enddo enddo - +! do i = 1,im dtsfc(i) = rho_a(i) * cp * heat(i) dqsfc(i) = rho_a(i) * hvap * evap(i) diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index d0b11656a..b6680dccb 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -273,6 +273,14 @@ type = real kind = kind_phys intent = in +[sigmaf] + standard_name = bounded_vegetation_area_fraction + long_name = areal fractional cover of green vegetation bounded on the bottom + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in [psk] standard_name = surface_dimensionless_exner_function long_name = dimensionless Exner function at the surface interface