From 14c2f52273cfcea106ff7d07001dff3c04d2fe1b Mon Sep 17 00:00:00 2001 From: rhaesung Date: Fri, 25 Oct 2024 18:25:03 +0000 Subject: [PATCH 01/13] modified prognostic updraft fraction (sigmab) calcultation --- physics/CONV/progsigma_calc.f90 | 126 +++++++++++++++++--------------- 1 file changed, 66 insertions(+), 60 deletions(-) diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index 469df49f6..b6f6428b1 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.0) + 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) + 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 From 36e7b63271dffae8f3a41957eea9d4f8452a417f Mon Sep 17 00:00:00 2001 From: rhaesung Date: Fri, 25 Oct 2024 18:27:08 +0000 Subject: [PATCH 02/13] fix in missing vertical transport of turbulent kinetic energy (TKE) when aerosol transport is turned on --- physics/CONV/SAMF/samfdeepcnv.f | 25 ++++++++++++++++++++----- physics/CONV/SAMF/samfdeepcnv.meta | 8 ++++++++ physics/CONV/SAMF/samfshalcnv.f | 27 +++++++++++++++++++++------ physics/CONV/SAMF/samfshalcnv.meta | 8 ++++++++ 4 files changed, 57 insertions(+), 11 deletions(-) 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 From 1ace72d8cf168909ef9fdc7012ceb4c33a52713c Mon Sep 17 00:00:00 2001 From: rhaesung Date: Fri, 25 Oct 2024 18:31:34 +0000 Subject: [PATCH 03/13] introduce TKE at model layer interfaces (tkeh) --- physics/PBL/SATMEDMF/satmedmfvdifq.F | 14 +++++++------- physics/PBL/SATMEDMF/satmedmfvdifq.meta | 8 ++++++++ 2 files changed, 15 insertions(+), 7 deletions(-) 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 From b325da7a1937a70b40ee1fae3ca6b805e3c9d9f5 Mon Sep 17 00:00:00 2001 From: rhaesung Date: Thu, 31 Oct 2024 15:11:03 +0000 Subject: [PATCH 04/13] fix compilation error --- physics/CONV/C3/cu_c3_deep.F90 | 6 +++--- physics/CONV/C3/cu_c3_sh.F90 | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) 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 From f6f18927fed34af525753c9635f95895f4577c49 Mon Sep 17 00:00:00 2001 From: rhaesung Date: Tue, 12 Nov 2024 15:00:28 +0000 Subject: [PATCH 05/13] replace the hardcoded precision constants with the appropriate precision --- physics/CONV/progsigma_calc.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index b6f6428b1..49a925867 100644 --- a/physics/CONV/progsigma_calc.f90 +++ b/physics/CONV/progsigma_calc.f90 @@ -115,7 +115,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& sigmab(i) = sigmain(i,k) else sigmab(i) = sigmab(i) / sumx(i) - sigmab(i) = min(sigmab(i), 1.0) + sigmab(i) = min(sigmab(i), 1._kind_phys) if(sigmab(i) < 1.E-5) sigmab(i)=0. endif endif @@ -184,7 +184,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& !sigmab do i = 1,im if(cnvflg(i))then - DEN=MIN(termC(i)+termB(i),1.E8) + DEN=MIN(termC(i)+termB(i),1._kind_phys) cvg=termD(i)*delt ZZ=MAX(0.0,SIGN(1.0,termA(i))) & *MAX(0.0,SIGN(1.0,termB(i))) & From 80609a95dc465374c87ae0bb984f624bdf37a89d Mon Sep 17 00:00:00 2001 From: rhaesung Date: Mon, 9 Dec 2024 17:58:29 +0000 Subject: [PATCH 06/13] Fixing the incorrectly modified precision constant --- physics/CONV/progsigma_calc.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index 49a925867..a7bd4b257 100644 --- a/physics/CONV/progsigma_calc.f90 +++ b/physics/CONV/progsigma_calc.f90 @@ -184,7 +184,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& !sigmab do i = 1,im if(cnvflg(i))then - DEN=MIN(termC(i)+termB(i),1._kind_phys) + 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))) & From a7f6a11a8a7a916aa92b1848e76ff5fc927347aa Mon Sep 17 00:00:00 2001 From: rhaesung Date: Mon, 10 Feb 2025 19:14:11 +0000 Subject: [PATCH 07/13] Fix bug --- physics/CONV/progsigma_calc.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index a7bd4b257..fc050679b 100644 --- a/physics/CONV/progsigma_calc.f90 +++ b/physics/CONV/progsigma_calc.f90 @@ -225,7 +225,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& do i= 1, im if(cnvflg(i)) then sigmab(i)=sigmab(i)/betadcu - sigmab(i)=MAX(sigmind,sigmab(i)) + sigmab(i)=MAX(sigmind_new,sigmab(i)) endif enddo endif From a2a4f60202be1e3d0bcb065fdd6c7d8b68fff4c4 Mon Sep 17 00:00:00 2001 From: rhaesung Date: Thu, 13 Feb 2025 01:46:25 +0000 Subject: [PATCH 08/13] fix bug --- physics/CONV/progsigma_calc.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index fc050679b..8926b1084 100644 --- a/physics/CONV/progsigma_calc.f90 +++ b/physics/CONV/progsigma_calc.f90 @@ -63,7 +63,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& km1=km-1 invdelt = 1./delt - if (flag_init) then + if(flag_init .and. .not. flag_restart) then sigmind_new=0.0 else sigmind_new=sigmind From 1bf0b07a2e42eb09862be5e387cf8a2974e06c60 Mon Sep 17 00:00:00 2001 From: rhaesung Date: Fri, 14 Feb 2025 13:02:59 +0000 Subject: [PATCH 09/13] make tkeh intent(inout) --- physics/CONV/SAMF/samfdeepcnv.f | 4 ++-- physics/CONV/SAMF/samfdeepcnv.meta | 2 +- physics/CONV/SAMF/samfshalcnv.f | 4 ++-- physics/CONV/SAMF/samfshalcnv.meta | 2 +- physics/PBL/SATMEDMF/satmedmfvdifq.F | 4 ++-- physics/PBL/SATMEDMF/satmedmfvdifq.meta | 2 +- 6 files changed, 9 insertions(+), 9 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 05b801c7e..c58d95886 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -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(:,:), tkeh(:,:) + & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), real(kind=kind_phys), dimension(:), intent(in) :: fscav logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, & & progsigma,do_mynnedmf @@ -115,7 +115,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH real(kind=kind_phys), intent(inout) :: qtr(:,:,:), & & q1(:,:), t1(:,:), u1(:,:), v1(:,:), & - & cnvw(:,:), cnvc(:,:) + & cnvw(:,:), cnvc(:,:), tkeh(:,:) integer, intent(out) :: kbot(:), ktop(:) real(kind=kind_phys), intent(out) :: cldwrk(:), & diff --git a/physics/CONV/SAMF/samfdeepcnv.meta b/physics/CONV/SAMF/samfdeepcnv.meta index 756c95ba5..f44c197be 100644 --- a/physics/CONV/SAMF/samfdeepcnv.meta +++ b/physics/CONV/SAMF/samfdeepcnv.meta @@ -249,7 +249,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = in + intent = inout [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 1d3b3dbd2..eb7a6c0f2 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -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(:,:), tkeh(:,:), & + & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), & & tmf(:,:,:), q(:,:) real(kind=kind_phys), intent(in), optional :: qmicro(:,:), & & prevsq(:,:) @@ -81,7 +81,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & integer, intent(inout) :: kcnv(:) ! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH real(kind=kind_phys), intent(inout) :: qtr(:,:,:), & - & q1(:,:), t1(:,:), u1(:,:), v1(:,:) + & q1(:,:), t1(:,:), u1(:,:), v1(:,:), tkeh(:,:) ! integer, intent(out) :: kbot(:), ktop(:) real(kind=kind_phys), intent(out) :: rn(:), & diff --git a/physics/CONV/SAMF/samfshalcnv.meta b/physics/CONV/SAMF/samfshalcnv.meta index fcf964e2b..fd0f8cb10 100644 --- a/physics/CONV/SAMF/samfshalcnv.meta +++ b/physics/CONV/SAMF/samfshalcnv.meta @@ -249,7 +249,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = in + intent = inout [qtr] standard_name = convective_transportable_tracers long_name = array to contain cloud water and other convective trans. tracers diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.F b/physics/PBL/SATMEDMF/satmedmfvdifq.F index 49e9f32f7..4f772c23b 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.F +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.F @@ -107,7 +107,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & real(kind=kind_phys), intent(in) :: dspfac, bl_upfr, bl_dnfr real(kind=kind_phys), intent(in) :: rlmx, elmx real(kind=kind_phys), intent(inout) :: dv(:,:), du(:,:), & - & tdt(:,:), rtg(:,:,:) + & tdt(:,:), rtg(:,:,:), tkeh(:,:) real(kind=kind_phys), intent(in) :: & & u1(:,:), v1(:,:), & & usfco(:), vsfco(:), & @@ -134,7 +134,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & dtsfc(:), dqsfc(:), & & hpbl(:) real(kind=kind_phys), intent(out) :: & - & dkt(:,:), dku(:,:), tkeh(:,:) + & dkt(:,:), dku(:,:) ! logical, intent(in) :: dspheat character(len=*), intent(out) :: errmsg diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.meta b/physics/PBL/SATMEDMF/satmedmfvdifq.meta index c89062889..4216b2102 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.meta +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.meta @@ -517,7 +517,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [dkt] standard_name = atmosphere_heat_diffusivity long_name = atmospheric heat diffusivity From 621a04df4a0930c5b1f3d4bb85ea52b08c908e12 Mon Sep 17 00:00:00 2001 From: rhaesung Date: Fri, 14 Feb 2025 13:34:28 +0000 Subject: [PATCH 10/13] fix error --- physics/CONV/SAMF/samfdeepcnv.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index c58d95886..9da4dfaff 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -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(:,:) real(kind=kind_phys), dimension(:), intent(in) :: fscav logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, & & progsigma,do_mynnedmf From 408d8e2d96a3c6866d99ac30646cafa8c1e461f4 Mon Sep 17 00:00:00 2001 From: rhaesung Date: Fri, 14 Feb 2025 14:37:02 +0000 Subject: [PATCH 11/13] set tkeh=0 before doing a loop --- physics/PBL/SATMEDMF/satmedmfvdifq.F | 1 + 1 file changed, 1 insertion(+) diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.F b/physics/PBL/SATMEDMF/satmedmfvdifq.F index 4f772c23b..41bdcc0d8 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.F +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.F @@ -351,6 +351,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do k=1,km do i=1,im tke(i,k) = max(q1(i,k,ntke), tkmin) + tkeh(i,k) = 0 enddo enddo do k=1,km1 From 467451fe88074ad0884b713b6e337f4ac918cf0d Mon Sep 17 00:00:00 2001 From: rhaesung Date: Fri, 14 Feb 2025 17:03:22 +0000 Subject: [PATCH 12/13] address reviewer comment --- physics/CONV/SAMF/samfdeepcnv.meta | 2 +- physics/CONV/SAMF/samfshalcnv.meta | 2 +- physics/PBL/SATMEDMF/satmedmfvdifq.meta | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.meta b/physics/CONV/SAMF/samfdeepcnv.meta index f44c197be..32e734318 100644 --- a/physics/CONV/SAMF/samfdeepcnv.meta +++ b/physics/CONV/SAMF/samfdeepcnv.meta @@ -246,7 +246,7 @@ 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) + dimensions = (horizontal_loop_extent,vertical_interface_dimension) type = real kind = kind_phys intent = inout diff --git a/physics/CONV/SAMF/samfshalcnv.meta b/physics/CONV/SAMF/samfshalcnv.meta index fd0f8cb10..df1fa3320 100644 --- a/physics/CONV/SAMF/samfshalcnv.meta +++ b/physics/CONV/SAMF/samfshalcnv.meta @@ -246,7 +246,7 @@ 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) + dimensions = (horizontal_loop_extent,vertical_interface_dimension) type = real kind = kind_phys intent = inout diff --git a/physics/PBL/SATMEDMF/satmedmfvdifq.meta b/physics/PBL/SATMEDMF/satmedmfvdifq.meta index 4216b2102..13c66399e 100644 --- a/physics/PBL/SATMEDMF/satmedmfvdifq.meta +++ b/physics/PBL/SATMEDMF/satmedmfvdifq.meta @@ -514,7 +514,7 @@ 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) + dimensions = (horizontal_loop_extent,vertical_interface_dimension) type = real kind = kind_phys intent = inout From b3b8cc675b0f78b3aa053d185b2f73961f83c205 Mon Sep 17 00:00:00 2001 From: rhaesung Date: Fri, 14 Feb 2025 18:21:34 +0000 Subject: [PATCH 13/13] address reviewer comment --- physics/CONV/SAMF/samfdeepcnv.f | 28 +++++++++++++++------------- physics/CONV/SAMF/samfshalcnv.f | 28 +++++++++++++++------------- 2 files changed, 30 insertions(+), 26 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 9da4dfaff..4b6921658 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -1285,22 +1285,24 @@ 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 * + if(ntk > 2) then + 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) + ercko(i,k,kk) = ecko(i,k,kk) + endif endif - endif + enddo enddo - enddo + endif endif endif c diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index eb7a6c0f2..c5c5ddfe8 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -1092,22 +1092,24 @@ 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 * + if(ntk > 2) then + 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) + ercko(i,k,kk) = ecko(i,k,kk) + endif endif - endif + enddo enddo - enddo + endif endif endif c