diff --git a/physics/CONV/C3/cu_c3_deep.F90 b/physics/CONV/C3/cu_c3_deep.F90 index 9bcbe5910..f3ec0a5e4 100644 --- a/physics/CONV/C3/cu_c3_deep.F90 +++ b/physics/CONV/C3/cu_c3_deep.F90 @@ -2016,9 +2016,9 @@ subroutine cu_c3_deep_run( & endif enddo call progsigma_calc(itf,ktf,flag_init,flag_restart,flag_shallow, & - flag_mid,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,xlv,dtime, & - forceqv_spechum,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, & - sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab) + flag_mid,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,xlv,dtime, & + forceqv_spechum,k22,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, & + sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab) endif !$acc end kernels diff --git a/physics/CONV/C3/cu_c3_sh.F90 b/physics/CONV/C3/cu_c3_sh.F90 index 0625b2949..f927d1fbd 100644 --- a/physics/CONV/C3/cu_c3_sh.F90 +++ b/physics/CONV/C3/cu_c3_sh.F90 @@ -983,9 +983,9 @@ subroutine cu_c3_sh_run ( & endif enddo call progsigma_calc(itf,ktf,flag_init,flag_restart,flag_shallow, & - flag_mid,del,tmf,qmicro,dbyo,zdqca,omega_u,zeta,xlv,dtime, & - forceqv_spechum,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, & - sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab) + flag_mid,del,tmf,qmicro,dbyo,zdqca,omega_u,zeta,xlv,dtime, & + forceqv_spechum,k22,kbcon,ktop,cnvflg,betascu,betamcu,betadcu, & + sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab) endif diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 3ad067657..05b801c7e 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -76,7 +76,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & tmf,qmicro,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,tkeh,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & hwrf_samfdeep,progsigma,cldwrk,rn,kbot,ktop,kcnv, & & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & @@ -97,7 +97,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & fv, grav, hvap, rd, rv, t0c real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & - & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:) + & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), tkeh(:,:) real(kind=kind_phys), dimension(:), intent(in) :: fscav logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, & & progsigma,do_mynnedmf @@ -953,8 +953,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if(cnvflg(i)) then if(k >= kb(i) .and. k < kbcon(i)) then dz = zo(i,k+1) - zo(i,k) - tem = 0.5 * (qtr(i,k,ntk)+qtr(i,k+1,ntk)) - tkemean(i) = tkemean(i) + tem * dz + tkemean(i) = tkemean(i) + tkeh(i,k) * dz sumx(i) = sumx(i) + dz endif endif @@ -1286,6 +1285,22 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & enddo enddo enddo + kk = ntk -2 + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + 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 + ercko(i,k,kk) = ecko(i,k,kk) + endif + endif + enddo + enddo endif endif c @@ -2941,7 +2956,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & flag_mid = .false. call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, - & delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu, + & delt,qadv,kb,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu, & sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab) endif diff --git a/physics/CONV/SAMF/samfdeepcnv.meta b/physics/CONV/SAMF/samfdeepcnv.meta index 3652a0d27..756c95ba5 100644 --- a/physics/CONV/SAMF/samfdeepcnv.meta +++ b/physics/CONV/SAMF/samfdeepcnv.meta @@ -242,6 +242,14 @@ type = real kind = kind_phys intent = in +[tkeh] + standard_name = vertical_turbulent_kinetic_energy_at_interface + long_name = vertical turbulent kinetic energy at model layer interfaces + units = m2 s-2 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [qtr] standard_name = convective_transportable_tracers long_name = array to contain cloud water and other convective trans. tracers diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index ce783ea15..1d3b3dbd2 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -53,7 +53,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp,first_time_step,restart, & & tmf,qmicro,progsigma, & - & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,tkeh,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, & @@ -71,7 +71,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & betamcu real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & - & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), & + & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), tkeh(:,:), & & tmf(:,:,:), q(:,:) real(kind=kind_phys), intent(in), optional :: qmicro(:,:), & & prevsq(:,:) @@ -872,14 +872,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & tkemean(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 = zo(i,k+1) - zo(i,k) - tem = 0.5 * (qtr(i,k,ntk)+qtr(i,k+1,ntk)) - tkemean(i) = tkemean(i) + tem * dz + tkemean(i) = tkemean(i) + tkeh(i,k) * dz sumx(i) = sumx(i) + dz endif endif @@ -1093,6 +1092,22 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo + kk = ntk -2 + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + 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 + ercko(i,k,kk) = ecko(i,k,kk) + endif + endif + enddo + enddo endif endif c @@ -1982,7 +1997,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & flag_mid = .false. call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, - & delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu, + & delt,qadv,kb,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu, & sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab) endif diff --git a/physics/CONV/SAMF/samfshalcnv.meta b/physics/CONV/SAMF/samfshalcnv.meta index 4dfa8ac20..fcf964e2b 100644 --- a/physics/CONV/SAMF/samfshalcnv.meta +++ b/physics/CONV/SAMF/samfshalcnv.meta @@ -242,6 +242,14 @@ type = real kind = kind_phys intent = in +[tkeh] + standard_name = vertical_turbulent_kinetic_energy_at_interface + long_name = vertical turbulent kinetic energy at model layer interfaces + units = m2 s-2 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [qtr] standard_name = convective_transportable_tracers long_name = array to contain cloud water and other convective trans. tracers diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index 469df49f6..a7bd4b257 100644 --- a/physics/CONV/progsigma_calc.f90 +++ b/physics/CONV/progsigma_calc.f90 @@ -21,7 +21,7 @@ module progsigma !!\section gen_progsigma progsigma_calc General Algorithm subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& flag_mid,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & - delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu, & + delt,qadv,kb,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu, & sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab) ! ! @@ -31,7 +31,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& implicit none ! intent in - integer, intent(in) :: im,km,kbcon1(im),ktcon(im) + integer, intent(in) :: im,km,kb(im),kbcon1(im),ktcon(im) real(kind=kind_phys), intent(in) :: hvap,delt,betascu,betamcu,betadcu, & sigmind,sigminm,sigmins real(kind=kind_phys), intent(in) :: qadv(im,km),del(im,km), & @@ -48,13 +48,13 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& ! Local variables integer :: i,k,km1 real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im) - real(kind=kind_phys) :: mcons(im),fdqa(im),form(im,km), & + real(kind=kind_phys) :: fdqa(im),form(im,km), & dp(im,km),inbu(im,km) - + real(kind=kind_phys) :: sumx(im) real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & fdqb,dtdyn,dxlim,rmulacvg,tem, & - DEN,dp1,invdelt + DEN,dp1,invdelt,sigmind_new !Parameters gcvalmx = 0.1 @@ -63,6 +63,12 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& km1=km-1 invdelt = 1./delt + if (flag_init) then + sigmind_new=0.0 + else + sigmind_new=sigmind + end if + !Initialization 2D do k = 1,km do i = 1,im @@ -80,7 +86,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& termC(i)=0. termD(i)=0. fdqa(i)=0. - mcons(i)=0. + sumx(i)=0. enddo do k = 2,km1 @@ -89,47 +95,49 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& dp(i,k) = 1000. * del(i,k) endif enddo - enddo + enddo - !Initial computations, place maximum sigmain in sigmab - do i=1,im +!compute sigmain averaged over cloud layers after advection and place it in sigmab + do k=2,km1 + do i=1,im if(cnvflg(i))then - do k=2,km - if(sigmain(i,k)>sigmab(i))then - sigmab(i)=sigmain(i,k) - endif - enddo - endif - enddo - - do i=1,im - if(cnvflg(i))then - if(sigmab(i) < 1.E-5)then !after advection - sigmab(i)=0. + if(k > kbcon1(i) .and. k < ktcon(i)) then + sigmab(i) = sigmab(i) + sigmain(i,k) * dp(i,k) + sumx(i) = sumx(i) + dp(i,k) endif - endif + endif + enddo + enddo + do i = 1, im + if(cnvflg(i)) then + if(sumx(i) == 0.) then + k = kbcon1(i) + sigmab(i) = sigmain(i,k) + else + sigmab(i) = sigmab(i) / sumx(i) + sigmab(i) = min(sigmab(i), 1._kind_phys) + if(sigmab(i) < 1.E-5) sigmab(i)=0. + endif + endif enddo !compute termD "The vertical integral of the latent heat convergence is limited to the - !buoyant layers with positive moisture convergence (accumulated from the surface). - !Lowest level: - do i = 1,im - dp1 = 1000. * del(i,1) - mcons(i)=(hvap*(qadv(i,1)+tmf(i,1)+qmicro(i,1))*dp1) - enddo - !Levels above: - do k = 2,km1 + ! layers with positive moisture convergence (accumulated from the updraft starting level). + do k = 1,km1 do i = 1,im if(cnvflg(i))then - mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i,k)) - buy2 = termD(i)+mcon+mcons(i) -! Do the integral over buoyant layers with positive mcon acc from surface - if(dbyo1(i,k)>0 .and. buy2 > 0.)then + if(k >= kb(i) .and. k < ktcon(i)) then + mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i,k)) + buy2 = termD(i)+mcon +! +! Do the integral over buoyant layers with positive mcon acc from +! updraft starting level +! + if(buy2 > 0.)then inbu(i,k)=1. - endif - inbu(i,k-1)=MAX(inbu(i,k-1),inbu(i,k)) - termD(i) = termD(i) + inbu(i,k-1)*mcons(i) - mcons(i)=mcon + termD(i) = termD(i) + mcon + endif + endif endif enddo enddo @@ -138,8 +146,10 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& do k = 2,km1 do i = 1,im if(cnvflg(i))then + if(k >= kbcon1(i) .and. k < ktcon(i)) then tem=sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k)*dp(i,k) termA(i)=termA(i)+tem + endif endif enddo enddo @@ -148,8 +158,10 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& do k = 2,km1 do i = 1,im if(cnvflg(i))then + if(k >= kbcon1(i) .and. k < ktcon(i)) then tem=zeta(i,k)*dbyo1(i,k)*inbu(i,k)*dp(i,k) termB(i)=termB(i)+tem + endif endif enddo enddo @@ -158,39 +170,33 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& do k = 2,km1 do i = 1,im if(cnvflg(i))then + if(k >= kbcon1(i) .and. k < ktcon(i)) then form(i,k)=-1.0*inbu(i,k)*(omega_u(i,k)*delt) fdqb=0.5*((form(i,k)*zdqca(i,k))) termC(i)=termC(i)+inbu(i,k)* & (fdqb+fdqa(i))*hvap*zeta(i,k) fdqa(i)=fdqb + endif endif enddo enddo !sigmab - if(flag_init .and. .not. flag_restart)then - do i = 1,im - if(cnvflg(i))then - sigmab(i)=0.03 + do i = 1,im + if(cnvflg(i))then + DEN=MIN(termC(i)+termB(i),1.e8_kind_phys) + cvg=termD(i)*delt + ZZ=MAX(0.0,SIGN(1.0,termA(i))) & + *MAX(0.0,SIGN(1.0,termB(i))) & + *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) + cvg=MAX(0.0,cvg) + sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) + if(sigmab(i)>0.)then + sigmab(i)=MIN(sigmab(i),0.95) + sigmab(i)=MAX(sigmab(i),sigmind_new) endif - enddo - else - do i = 1,im - if(cnvflg(i))then - DEN=MIN(termC(i)+termB(i),1.E8) - cvg=termD(i)*delt - ZZ=MAX(0.0,SIGN(1.0,termA(i))) & - *MAX(0.0,SIGN(1.0,termB(i))) & - *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) - cvg=MAX(0.0,cvg) - sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) - if(sigmab(i)>0.)then - sigmab(i)=MIN(sigmab(i),0.95) - sigmab(i)=MAX(sigmab(i),0.01) - endif - endif!cnvflg - enddo - endif + endif!cnvflg + enddo do k=1,km do i=1,im diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.F b/physics/PBL/SATMEDMF/satmedmfvdifq.F index c0e43e809..49e9f32f7 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.F +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.F @@ -80,7 +80,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & & tsea,heat,evap,stress,spd1,kpbl, & & prsi,del,prsl,prslk,phii,phil,delt, & - & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku, & + & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku,tkeh, & & kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr, & & rlmx,elmx,sfc_rlm,tc_pbl, & & ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind, & @@ -134,7 +134,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & dtsfc(:), dqsfc(:), & & hpbl(:) real(kind=kind_phys), intent(out) :: & - & dkt(:,:), dku(:,:) + & dkt(:,:), dku(:,:), tkeh(:,:) ! logical, intent(in) :: dspheat character(len=*), intent(out) :: errmsg @@ -152,7 +152,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), e2(im,0:km) + real(kind=kind_phys) tke(im,km), tkei(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), @@ -1561,7 +1561,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! do k=1,kps do i=1,im - tkeh(i,k) = 0.5 * (tke(i,k)+tke(i,k+1)) + tkei(i,k) = 0.5 * (tke(i,k)+tke(i,k+1)) enddo enddo do k=1,kps @@ -1582,7 +1582,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do k = 1, kps do i = 1, im kmx = max(kpbl(i), krad(i)) - e_half(i,k) = tkeh(i,k) + e_half(i,k) = tkei(i,k) if((pcnvflg(i) .or. scuflg(i)) .and. (k < kmx)) then tem = 0. if(pcnvflg(i) .and. k < kpbl(i)) then @@ -1598,14 +1598,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & rrkp = e_diff(i,k+1) / e_diff(i,k) phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) e_half(i,k) = tke(i,k+1) + - & phkp*(tkeh(i,k)-tke(i,k+1)) + & phkp*(tkei(i,k)-tke(i,k+1)) elseif (tem < 0.) then rrkp = 0. if(abs(e_diff(i,k)) > 1.e-22) & rrkp = e_diff(i,k-1) / e_diff(i,k) phkp = (rrkp+abs(rrkp)) / (1.+abs(rrkp)) e_half(i,k) = tke(i,k) + - & phkp*(tkeh(i,k)-tke(i,k)) + & phkp*(tkei(i,k)-tke(i,k)) endif endif enddo diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.meta b/physics/PBL/SATMEDMF/satmedmfvdifq.meta index cdbfa67b7..c89062889 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.meta +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.meta @@ -510,6 +510,14 @@ type = real kind = kind_phys intent = out +[tkeh] + standard_name = vertical_turbulent_kinetic_energy_at_interface + long_name = vertical turbulent kinetic energy at model layer interfaces + units = m2 s-2 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [dkt] standard_name = atmosphere_heat_diffusivity long_name = atmospheric heat diffusivity