From e673808397a1317cec46dfde6141b1349697e8cf Mon Sep 17 00:00:00 2001 From: Lisa Date: Fri, 6 Mar 2026 16:58:20 +0000 Subject: [PATCH 1/5] Updates to progomega option of saSAS convection scheme --- physics/CONV/SAMF/samfdeepcnv.f | 72 +++++++++++++++++++-------------- physics/CONV/SAMF/samfshalcnv.f | 69 +++++++++++++++++++------------ physics/CONV/progomega_calc.f90 | 34 +++++++++------- 3 files changed, 105 insertions(+), 70 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 18be662f0..a0d5a3a14 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -221,9 +221,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km), - & tentr(im,km) - real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins - parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01) + & wc_ref(im) + real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins, + & wc_min logical flag_shallow, flag_mid c physical parameters ! parameter(grav=grav,asolfac=0.958) @@ -347,6 +347,21 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !************************************************************************ ! ! +! - Initialize parameters related to prognostic closure + if (progsigma) then + if (progomega) then + sigmind = 0.03 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + else + sigmind = 0.01 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + endif + endif + km1 = km - 1 !> - Initialize column-integrated and other single-value-per-column variable arrays. c @@ -1137,7 +1152,6 @@ 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 - tentr(i,k)=xlamue(i,k)*fent1(i,k) tem = cxlamet(i) * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem tem1 = cxlamdt(i) * frh(i,k) @@ -1788,11 +1802,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if (progomega) then call progomega_calc(first_time_step,restart,im,km, & kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout, - & grav,buo,drag,wush,tentr,bb1,bb2) + & grav,buo,drag,wush,bb1,bb2) do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then omega_u(i,k)=omegaout(i,k) omega_u(i,k)=MAX(omega_u(i,k),-80.) ! Convert to m/s for use in convective time-scale: @@ -1809,7 +1823,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(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-1)+drag(i,k)) * dz tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) @@ -1827,7 +1841,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then rho = po(i,k)*100. / (rd * to(i,k)) omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav omega_u(i,k)=MAX(omega_u(i,k),-80.) @@ -1878,7 +1892,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then dp = 1000. * del(i,k) tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1)) omegac(i) = omegac(i) + tem * dp @@ -2907,27 +2921,25 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! compute convective turn-over time ! !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, the convective adjustment time (dtconv) is set to be proportional to the convective turnover time, which is computed using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx). - if(hwrf_samfdeep) then - do i= 1, im - if(cnvflg(i)) then - tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) - dtconv(i) = tem / wc(i) - dtconv(i) = max(dtconv(i),dtmin) - dtconv(i) = min(dtconv(i),dtmax) - endif - enddo - else - do i= 1, im - if(cnvflg(i)) then - tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) - dtconv(i) = tem / wc(i) - tfac = 1. + gdx(i) / 75000. - dtconv(i) = tfac * dtconv(i) - dtconv(i) = max(dtconv(i),dtmin) - dtconv(i) = min(dtconv(i),dtmax) - endif - enddo - endif + do i = 1, im + if (cnvflg(i)) then + tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) + if (progomega) then + wc_eff = max(wc(i), wc_min) + dtconv(i) = tem / wc_eff + else + dtconv(i) = tem / wc(i) + endif + !grid spacing scaling (disabled for HWRF SAMF deep) + if (.not. hwrf_samfdeep) then + tfac = 1. + gdx(i) / 75000. + dtconv(i) = tfac * dtconv(i) + endif + !bounds + dtconv(i) = max(dtconv(i), dtmin) + dtconv(i) = min(dtconv(i), dtmax) + endif + enddo ! !> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. do i= 1, im diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index 1ec8747f0..f8cedc769 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -169,9 +169,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), - & sigmab(im),qadv(im,km) + & sigmab(im),qadv(im,km),wc_ref(im) real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins, - & sigminm + & sigminm,wc_min logical flag_shallow,flag_mid c physical parameters ! parameter(g=grav,asolfac=0.89) @@ -205,8 +205,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) parameter(betaw=.03,dxcrtc0=9.e3) parameter(h1=0.33333333) -! progsigma - parameter(dxcrtas=500.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01) c local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), & uo(im,km), vo(im,km), qeso(im,km), @@ -295,7 +293,21 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & prsl = prslp * 0.001 del = delp * 0.001 !************************************************************************ -! +! - Initialize parameters related to prognostic closure + if (progsigma) then + if (progomega) then + sigmind = 0.03 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + else + sigmind = 0.01 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + endif + endif +! km1 = km - 1 c c initialize arrays @@ -1520,11 +1532,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (progomega) then call progomega_calc(first_time_step,restart,im,km, & kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout, - & grav,buo,drag,wush,xlamue,bb1,bb2) + & grav,buo,drag,wush,bb1,bb2) do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then omega_u(i,k)=omegaout(i,k) omega_u(i,k)=MAX(omega_u(i,k),-80.) ! Convert to m/s for use in convective time-scale: @@ -1542,7 +1554,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(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-1)+drag(i,k)) * dz tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) @@ -1560,7 +1572,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then rho = po(i,k)*100. / (rd * to(i,k)) omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav omega_u(i,k)=MAX(omega_u(i,k),-80.) @@ -1611,7 +1623,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then dp = 1000. * del(i,k) tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1)) omegac(i) = omegac(i) + tem * dp @@ -1636,7 +1648,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then if(omega_u(i,k) .ne. 0.)then zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) else @@ -1955,21 +1967,28 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! compute convective turn-over time ! !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, calculate the convective turnover time using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx). - do i= 1, im - if(cnvflg(i)) then - tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) - dtconv(i) = tem / wc(i) - if (.not.hwrf_samfshal) then - tfac = 1. + gdx(i) / 75000. - dtconv(i) = tfac * dtconv(i) - endif - dtconv(i) = max(dtconv(i),dtmin) - dtconv(i) = max(dtconv(i),dt2) - dtconv(i) = min(dtconv(i),dtmax) - endif + do i = 1, im + if (cnvflg(i)) then + tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) + if (progomega) then + wc_eff = max(wc(i), wc_min) + dtconv(i) = tem / wc_eff + else + dtconv(i) = tem / wc(i) + endif +! - grid spacing scaling (disabled for HWRF shallow option) + if (.not. hwrf_samfshal) then + tfac = 1. + gdx(i) / 75000. + dtconv(i) = tfac * dtconv(i) + endif +! - limits + dtconv(i) = max(dtconv(i), dtmin) + dtconv(i) = max(dtconv(i),dt2) + dtconv(i) = min(dtconv(i), dtmax) + endif enddo -! -!> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. +! +! > - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. do i= 1, im if(cnvflg(i)) then sumx(i) = 0. diff --git a/physics/CONV/progomega_calc.f90 b/physics/CONV/progomega_calc.f90 index 52d99d801..d22659156 100644 --- a/physics/CONV/progomega_calc.f90 +++ b/physics/CONV/progomega_calc.f90 @@ -20,7 +20,7 @@ module progomega !!\section gen_progomega progomega_calc General Algorithm subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegain,delt,del, & - zi,cnvflg,omegaout,grav,buo,drag,wush,tentr,bb1,bb2) + zi,cnvflg,omegaout,grav,buo,drag,wush,bb1,bb2) use machine, only : kind_phys use funcphys, only : fpvs @@ -30,20 +30,17 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai integer, intent(in) :: kbcon1(im),ktcon(im) real(kind=kind_phys), intent(in) :: delt,grav,bb1,bb2 real(kind=kind_phys), intent(in) :: omegain(im,km), del(im,km),zi(im,km) - real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km),tentr(im,km) + real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km) real(kind=kind_phys), intent(inout) :: omegaout(im,km) logical, intent(in) :: cnvflg(im),first_time_step,flag_restart real(kind=kind_phys) :: termA(im,km),termB(im,km),termC(im,km),omega(im,km) real(kind=kind_phys) :: RHS(im,km),Kd(im,km) - real(kind=kind_phys) :: dp,dz,entrn,Kdn,discr,wush_pa,lbb1,lbb2,lbb3 + real(kind=kind_phys) :: dp,dz,discr,wush_pa,lbb1,lbb2,lbb3 integer :: i,k - entrn = 0.8E-4 !0.5E-4 !m^-1 - Kdn = 0.5E-4 !2.9E-4 !m^-1 - lbb1 = 0.5 !1.0 - lbb2 = 3.2 !3.0 - lbb3 = 0.5 !0.5 - + lbb1 = 1.5 + lbb2 = 0.6 + lbb3 = 1.2 !Initialization 2D do k = 1,km @@ -56,6 +53,16 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai enddo enddo + do k = 1,km + do i = 1,im + if(cnvflg(i))then + if(omega(i,k) < 1.0E-5) then + omega(i,k) = 0. + endif + endif + enddo + enddo + if(first_time_step .and. .not. flag_restart)then do k = 1,km do i = 1,im @@ -76,10 +83,7 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai do k = 2, km do i = 1, im if (cnvflg(i)) then - if (k > kbcon1(i) .and. k < ktcon(i)) then - - ! Aerodynamic drag parameter - Kd(i,k) = (tentr(i,k)/entrn)*Kdn + if (k >= kbcon1(i) .and. k < ktcon(i)) then ! Scale by dp/dz to have equation in Pa/s !(dp/dz > 0) @@ -91,8 +95,8 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai !termC - Adds buoyancy forcing !Coefficients for the quadratic equation - termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz)) + (Kd(i,k) * (dp/dz))) - termB(i,k) = -1.0 - delt * lbb3 * wush(i,k) * dp/dz + termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz))) + termB(i,k) = 1.0 - delt * lbb3 * wush(i,k) * dp/dz termC(i,k) = omega(i,k) - delt * lbb2 * buo(i,k) * (dp/dz) & - delt * omega(i,k) * (omega(i,k-1) - omega(i,k)) / dp !Compute the discriminant From cd6a29de4e992d1f21e14c8773f25c9c5c2eeadc Mon Sep 17 00:00:00 2001 From: Lisa Date: Fri, 6 Mar 2026 17:54:42 +0000 Subject: [PATCH 2/5] Correction --- physics/CONV/SAMF/samfdeepcnv.f | 2 +- physics/CONV/SAMF/samfshalcnv.f | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index a0d5a3a14..4f6176a6d 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -221,7 +221,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km), - & wc_ref(im) + & wc_eff(im) real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins, & wc_min logical flag_shallow, flag_mid diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index f8cedc769..ba2924297 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -169,7 +169,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), - & sigmab(im),qadv(im,km),wc_ref(im) + & sigmab(im),qadv(im,km),wc_eff(im) real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins, & sigminm,wc_min logical flag_shallow,flag_mid From c2991ab0daa1919b71942a2ae4230816289e0537 Mon Sep 17 00:00:00 2001 From: Lisa Date: Fri, 6 Mar 2026 18:02:53 +0000 Subject: [PATCH 3/5] correction 2 --- physics/CONV/SAMF/samfdeepcnv.f | 5 ++--- physics/CONV/SAMF/samfshalcnv.f | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 4f6176a6d..651829379 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -220,10 +220,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), - & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km), - & wc_eff(im) + & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km) real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins, - & wc_min + & wc_min, wc_eff logical flag_shallow, flag_mid c physical parameters ! parameter(grav=grav,asolfac=0.958) diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index ba2924297..e26ccacb1 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -169,9 +169,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), - & sigmab(im),qadv(im,km),wc_eff(im) + & sigmab(im),qadv(im,km) real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins, - & sigminm,wc_min + & sigminm,wc_min,wc_eff logical flag_shallow,flag_mid c physical parameters ! parameter(g=grav,asolfac=0.89) From f2fbb17a330938f51a3d7cc1a54a184661dc5570 Mon Sep 17 00:00:00 2001 From: Lisa Date: Fri, 20 Mar 2026 18:45:44 +0000 Subject: [PATCH 4/5] Correct parameter setting in samfshalcnv to pass debug runs --- physics/CONV/SAMF/samfshalcnv.f | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index e26ccacb1..55cf6f7f4 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -252,7 +252,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) tf, tcr, tcrf parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) - c----------------------------------------------------------------------- ! ! Initialize CCPP error handling variables @@ -274,8 +273,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (progsigma) then dxcrt=10.e3 + dxcrtas=500.e3 else dxcrt=15.e3 + dxcrtas=500.e3 endif c----------------------------------------------------------------------- From e5ebf607cd342b1125941ea4e843a064503008f6 Mon Sep 17 00:00:00 2001 From: Lisa Date: Mon, 23 Mar 2026 16:10:45 +0000 Subject: [PATCH 5/5] address saturation vapor pressure out of bounds error in hafs_4_nest RT --- physics/CONV/SAMF/samfdeepcnv.f | 4 ++-- physics/CONV/SAMF/samfshalcnv.f | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 651829379..2acf3b55b 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -1822,7 +1822,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon(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-1)+drag(i,k)) * dz tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) @@ -1840,7 +1840,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then rho = po(i,k)*100. / (rd * to(i,k)) omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav omega_u(i,k)=MAX(omega_u(i,k),-80.) diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index 55cf6f7f4..43d7eeec2 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -1555,7 +1555,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon(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-1)+drag(i,k)) * dz tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) @@ -1573,7 +1573,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k >= kbcon1(i) .and. k < ktcon(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then rho = po(i,k)*100. / (rd * to(i,k)) omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav omega_u(i,k)=MAX(omega_u(i,k),-80.)