diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 18be662f0..2acf3b55b 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), - & tentr(im,km) - real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins - parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01) + & 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_eff logical flag_shallow, flag_mid c physical parameters ! parameter(grav=grav,asolfac=0.958) @@ -347,6 +346,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 +1151,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 +1801,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: @@ -1878,7 +1891,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 +2920,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..43d7eeec2 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -171,7 +171,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & omegac(im),zeta(im,km),dbyo1(im,km), & sigmab(im),qadv(im,km) real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins, - & sigminm + & sigminm,wc_min,wc_eff 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), @@ -254,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 @@ -276,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----------------------------------------------------------------------- @@ -295,7 +294,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 +1533,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: @@ -1611,7 +1624,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 +1649,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 +1968,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