Introduce prognostic updraft velocity to saSAS and C3 cumulus convection schemes

physics/CONV/SAMF/samfdeepcnv.f

@@ -8,7 +8,8 @@ module samfdeepcnv
 
       use samfcnv_aerosols, only : samfdeepcnv_aerosols
       use progsigma, only : progsigma_calc
-
+      use progomega, only : progomega_calc
+
       contains
 
       subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf,            &
@@ -77,14 +78,15 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
      &    eps,epsm1,fv,grav,hvap,rd,rv,                                 &
      &    t0c,delt,ntk,ntr,delp,                                        &
      &    prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav,                &
-     &    hwrf_samfdeep,progsigma,cldwrk,rn,kbot,ktop,kcnv,             &
+     &    hwrf_samfdeep,progsigma,progomega,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,                           &
      &    CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,&
      &    clam,c0s,c1,betal,betas,evef,pgcon,asolfac,                   &
      &    do_ca, ca_closure, ca_entr, ca_trigger, nthresh,ca_deep,      &
-     &    rainevap,sigmain,sigmaout,betadcu,betamcu,betascu,            &
-     &    maxMF, do_mynnedmf,sigmab_coldstart,errmsg,errflg)
+     &    rainevap,sigmain,sigmaout,omegain,omegaout,betadcu,betamcu,   &
+     &    betascu,maxMF,do_mynnedmf,sigmab_coldstart,errmsg,errflg)
+
 !
       use machine , only : kind_phys
       use funcphys , only : fpvs
@@ -100,16 +102,17 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
      &   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,sigmab_coldstart
+     &     progsigma,progomega,do_mynnedmf,sigmab_coldstart
       real(kind=kind_phys), intent(in) :: nthresh,betadcu,betamcu,      &
      &                                    betascu
       real(kind=kind_phys), intent(in), optional :: ca_deep(:)
       real(kind=kind_phys), intent(in), optional :: sigmain(:,:),       &
-     &     qmicro(:,:),  prevsq(:,:)
+     &     qmicro(:,:),  prevsq(:,:), omegain(:,:)
       real(kind=kind_phys), intent(in) :: tmf(:,:,:),q(:,:)
       real(kind=kind_phys), dimension (:), intent(in), optional :: maxMF
       real(kind=kind_phys), intent(out) :: rainevap(:)
-      real(kind=kind_phys), intent(out), optional :: sigmaout(:,:)
+      real(kind=kind_phys), intent(inout), optional :: sigmaout(:,:),     &
+     &     omegaout(:,:)
       logical, intent(in)  :: do_ca,ca_closure,ca_entr,ca_trigger
       integer, intent(inout)  :: kcnv(:)
       ! DH* TODO - check dimensions of qtr, ntr+2 correct?  *DH
@@ -216,7 +219,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),
-     &     sigmaoutx(im)
+     &     sigmaoutx(im),tentr(im,km)
       real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins
       parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01)
       logical flag_shallow, flag_mid
@@ -333,6 +336,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
 c-----------------------------------------------------------------------
 !>  ## Compute preliminary quantities needed for static, dynamic, and feedback control portions of the algorithm.
 !>  - Convert input pressure terms to centibar units.
+
 !************************************************************************
 !     convert input Pa terms to Cb terms  -- Moorthi
       ps   = psp   * 0.001
@@ -1131,7 +1135,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
         do k = 2, km1
         do i=1,im
           if(cnvflg(i) .and.
-     &      (k > kbcon(i) .and. k < kmax(i))) then
+     &       (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)
@@ -1743,43 +1748,64 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
       enddo
 !
 !  compute updraft velocity square(wu2)
-!> - Calculate updraft velocity square(wu2) according to Han et al.'s (2017) \cite han_et_al_2017 equation 7.
-!
+!> - Calculate diagnostic updraft velocity square(wu2) according to Han et al.'s (2017) \cite han_et_al_2017 equation 7.
+!> - if progomega = true, calculate prognostic updraft velocity (Pa/s) according to progomega routine.
+
       if (hwrf_samfdeep) then
-      do i = 1, im
-        if (cnvflg(i)) then
-          k = kbcon1(i)
-          tem = po(i,k) / (rd * to(i,k))
-          wucb = -0.01 * dot(i,k) / (tem * grav)
-          if(wucb.gt.0.) then
-            wu2(i,k) = wucb * wucb
-          else
-            wu2(i,k) = 0.
-          endif
-        endif
-      enddo
-      endif
-!
-      do k = 2, km1
-        do i = 1, im
-          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-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 + tem2) / ptem1
-              wu2(i,k) = max(wu2(i,k), 0.)
+         do i = 1, im
+            if (cnvflg(i)) then
+               k = kbcon1(i)
+               tem = po(i,k) / (rd * to(i,k))
+               wucb = -0.01 * dot(i,k) / (tem * grav)
+               if(wucb.gt.0.) then
+                  wu2(i,k) = wucb * wucb
+               else
+                  wu2(i,k) = 0.
+               endif
             endif
-          endif
-        enddo
-      enddo
-
-      if(progsigma)then                                                                                                                                                                   
-          do k = 2, km1
+         enddo
+      endif
+!                  
+      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)
+         do k = 1, km
+            do i = 1, im
+               if (cnvflg(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:
+                     rho = po(i,k)*100. / (rd * to(i,k))
+                     tem = (-omega_u(i,k)) / ((rho * grav))
+                     wu2(i,k) = tem**2
+                     wu2(i,k) = max(wu2(i,k), 0.)
+                  endif
+              endif
+            enddo
+         enddo
+      else
+!     diagnostic method:
+         do k = 2, km1
+            do i = 1, im
+               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-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 + tem2) / ptem1
+                     wu2(i,k) = max(wu2(i,k), 0.)
+                  endif
+               endif
+            enddo
+         enddo
+!       convert to Pa/s for use in closure
+         do k = 1, km
             do i = 1, im
                if (cnvflg(i)) then
                   if(k > kbcon1(i) .and. k < ktcon(i)) then
@@ -1790,10 +1816,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
                endif
             enddo
          enddo
-      endif 
+
+      endif                     !progomega
+
 !
 !  compute updraft velocity average over the whole cumulus
-!
 !> - Calculate the mean updraft velocity within the cloud (wc).
       do i = 1, im
         wc(i) = 0.
@@ -1821,11 +1848,10 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
           val = 1.e-4
           if (wc(i) < val) cnvflg(i)=.false.
         endif
-      enddo
+      enddo      
 c
-
 !> - For progsigma = T, calculate the mean updraft velocity within the cloud (omegac),cast in pressure coordinates.                                                                                                                                  
-      if(progsigma)then                                                                                                                                                            
+      if(progsigma)then
          do i = 1, im
             omegac(i) = 0.
             sumx(i) = 0.
@@ -2912,10 +2938,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
            advfac(i) = min(advfac(i), 1.)
         endif
       enddo
-
+      
 !> - From Bengtsson et al. (2022) \cite Bengtsson_2022 prognostic closure scheme, equation 8, call progsigma_calc() to compute updraft area fraction based on a moisture budget
       if(progsigma)then
-
 !Initial computations, dynamic q-tendency                                                                                                                                               
          if(first_time_step .and. (.not.restart 
      &           .or. sigmab_coldstart))then

physics/CONV/SAMF/samfdeepcnv.meta

@@ -1,7 +1,7 @@
 [ccpp-table-properties]
   name = samfdeepcnv
   type = scheme
-  dependencies = ../../tools/funcphys.f90,../../hooks/machine.F,samfaerosols.F,../progsigma_calc.f90
+  dependencies = ../../tools/funcphys.f90,../../hooks/machine.F,samfaerosols.F,../progsigma_calc.f90,../progomega_calc.f90
 
 ########################################################################
 [ccpp-arg-table]
@@ -321,6 +321,13 @@
   dimensions = ()
   type = logical
   intent = in
+[progomega]
+  standard_name = do_prognostic_updraft_velocity
+  long_name = flag for prognostic omega in cumuls scheme
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
 [cldwrk]
   standard_name = cloud_work_function
   long_name = cloud work function
@@ -453,7 +460,25 @@
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)
   type = real
   kind = kind_phys
-  intent = out
+  intent = inout
+  optional = True
+[omegain]
+  standard_name = prognostic_updraft_velocity_in_convection
+  long_name = convective updraft velocity
+  units = frac
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = True
+[omegaout]
+  standard_name = updraft_velocity_updated_by_physics
+  long_name = convective updraft velocity updated by physics
+  units = frac
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = inout
   optional = True
 [betascu]
   standard_name = tuning_param_for_shallow_cu