From 2fc4a64040e36f150af41890270112b1abc3200d Mon Sep 17 00:00:00 2001
From: Lisa Bengtsson <Lisa.Bengtsson@noaa.gov>
Date: Tue, 24 Oct 2023 22:11:54 +0000
Subject: [PATCH 1/4] Updates to bring out tuning parameters in C3 and SAS
 convection schemes

---
 physics/cu_c3_deep.F90     | 26 ++++++++++++++++++--------
 physics/cu_c3_driver.F90   | 16 ++++++++++++----
 physics/cu_c3_driver.meta  | 23 +++++++++++++++++++++++
 physics/cu_c3_sh.F90       | 23 ++++++++++++++---------
 physics/progsigma_calc.f90 | 31 ++++++++++++++++++++-----------
 physics/samfdeepcnv.f      | 19 +++++++++++--------
 physics/samfdeepcnv.meta   | 23 +++++++++++++++++++++++
 physics/samfshalcnv.f      | 19 +++++++++++--------
 physics/samfshalcnv.meta   | 23 +++++++++++++++++++++++
 9 files changed, 155 insertions(+), 48 deletions(-)

diff --git a/physics/cu_c3_deep.F90 b/physics/cu_c3_deep.F90
index 7092840c3..7e907aaba 100644
--- a/physics/cu_c3_deep.F90
+++ b/physics/cu_c3_deep.F90
@@ -97,6 +97,9 @@ subroutine cu_c3_deep_run(        &
               ,tmf           &  ! instantanious tendency from turbulence
               ,qmicro        &  ! instantanious tendency from microphysics
               ,forceqv_spechum & !instantanious tendency from dynamics
+              ,betascu       &  ! Tuning parameter for shallow clouds
+              ,betamcu       &  ! Tuning parameter for mid-level clouds
+              ,betadcu       &  ! Tuning parameter for deep clouds
               ,sigmain       &  ! input area fraction after advection
               ,sigmaout      &  ! updated prognostic area fraction
               ,z1            &  ! terrain
@@ -233,8 +236,8 @@ subroutine cu_c3_deep_run(        &
 
        
        real(kind=kind_phys)                                                            &
-        ,intent (in   )                   ::                           &
-        dtime,ccnclean,fv,r_d
+        ,intent (in   )                   ::                                           &
+        dtime,ccnclean,fv,r_d,betascu,betamcu,betadcu
 
 
 !
@@ -386,13 +389,16 @@ subroutine cu_c3_deep_run(        &
      real(kind=kind_phys), dimension (its:ite) :: pefc
      real(kind=kind_phys) entdo,dp,subin,detdo,entup,                    &
       detup,subdown,entdoj,entupk,detupk,totmas
+     real(kind=kind_phys)                 ::                             &
+          sigmind,sigminm,sigmins
+     parameter(sigmind=0.005,sigmins=0.03,sigminm=0.01)
 
      real(kind=kind_phys), dimension (its:ite) :: lambau,flux_tun,zws,ztexec,zqexec
 !$acc declare create(lambau,flux_tun,zws,ztexec,zqexec)
 
      integer :: jprnt,jmini,start_k22
      logical :: keep_going,flg(its:ite),cnvflg(its:ite)
-     logical :: flag_shallow
+     logical :: flag_shallow,flag_mid
 
 !$acc declare create(flg)
      
@@ -1988,7 +1994,11 @@ subroutine cu_c3_deep_run(        &
 ! equation 8, call progsigma_calc() to compute updraft area fraction based on a moisture budget
                    
       if(progsigma)then
+         flag_mid = .false.
          flag_shallow = .false.
+         if(imid.eq.1)then
+            flag_mid = .true.
+         endif
          do k=kts,ktf
             do i=its,itf
                del(i,k) = delp(i,k)*0.001
@@ -2003,9 +2013,9 @@ subroutine cu_c3_deep_run(        &
             endif
          enddo
          call progsigma_calc(itf,ktf,flag_init,flag_restart,flag_shallow,  &
-              del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,xlv,dtime,           &
-              forceqv_spechum,kbcon,ktop,cnvflg,                           &
-              sigmain,sigmaout,sigmab)        
+              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)        
       endif
 
 !$acc end kernels
@@ -3147,7 +3157,7 @@ subroutine cup_forcing_ens_3d(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2
 !       pcrit,acrit,acritt
      integer, dimension (its:ite)         :: kloc
      real(kind=kind_phys)                                ::                           &
-       a1,a_ave,xff0,xomg,gravinv!,aclim1,aclim2,aclim3,aclim4
+       a1,a_ave,xff0,xomg,gravinv
 
      real(kind=kind_phys), dimension (its:ite) :: ens_adj
 !$acc declare create(kloc,ens_adj)
@@ -5748,7 +5758,7 @@ subroutine calculate_updraft_velocity(its,itf,ktf,ite,kts,kte,ierr,progsigma,
     do k = 2, ktf-1
        do i = 1, itf
          if (ierr(i)==0) then
-           if(k >= kbcon(i) .and. k < ktcon(i))then
+           if(k >= kbcon(i) .and. k < ktcon(i) .and. dbyo(i,k)>0.)then
               gamma   = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
               if(k >= kbcon(i) .and. clw_all(i,k)>0.)then
                buo(i,k) = buo(i,k) - g * qlk(i,k)  
diff --git a/physics/cu_c3_driver.F90 b/physics/cu_c3_driver.F90
index 8592e08f9..0ecb81750 100644
--- a/physics/cu_c3_driver.F90
+++ b/physics/cu_c3_driver.F90
@@ -60,7 +60,8 @@ end subroutine cu_c3_driver_init
       subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
                do_ca,progsigma,cactiv,cactiv_m,g,cp,fv,r_d,xlv,r_v,forcet,      &
                forceqv_spechum,phil,delp,raincv,tmf,qmicro,sigmain,             &
-               qv_spechum,t,cld1d,us,vs,t2di,w,qv2di_spechum,p2di,psuri,        &
+               betascu,betamcu,betadcu,qv_spechum,t,cld1d,us,vs,t2di,w,         &
+               qv2di_spechum,p2di,psuri,                                        &
                hbot,htop,kcnv,xland,hfx2,qfx2,aod_gf,cliw,clcw,ca_deep,rainevap,&
                pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv,                &
                flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend,           &
@@ -97,7 +98,7 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
    logical, intent(in   ) :: flag_init, flag_restart, do_mynnedmf
    logical, intent(in   ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend, &
         do_ca,progsigma
-   real (kind=kind_phys), intent(in) :: g,cp,fv,r_d,xlv,r_v
+   real (kind=kind_phys), intent(in) :: g,cp,fv,r_d,xlv,r_v,betascu,betamcu,betadcu
    logical, intent(in   ) :: ldiag3d
 
    real(kind=kind_phys), intent(inout)                      :: dtend(:,:,:)
@@ -587,7 +588,7 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
       hfx(i)=hfx2(i)*cp*rhoi(i,1)
       qfx(i)=qfx2(i)*xlv*rhoi(i,1)
       dx(i) = sqrt(garea(i))
-     enddo
+     enddo    
 
      do i=its,itf
       do k=kts,kpbli(i)
@@ -669,7 +670,8 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
                          zus,xmbs,kbcons,ktops,k22s,ierrs,ierrcs,                &
 ! Prog closure
                          flag_init, flag_restart,fv,r_d,delp,tmfq,qmicro,        &
-                         forceqv_spechum,sigmain,sigmaout,progsigma,dx,          &
+                         forceqv_spechum,betascu,betamcu,betadcu,sigmain,        &
+                         sigmaout,progsigma,dx,                                  &
 ! output tendencies
                          outts,outqs,outqcs,outus,outvs,cnvwt,prets,cupclws,     &
 ! dimesnional variables
@@ -714,6 +716,9 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
               ,tmfq          &
               ,qmicro        &
               ,forceqv_spechum &
+              ,betascu       &
+              ,betamcu       &
+              ,betadcu       &
               ,sigmain       &
               ,sigmaout      &
               ,ter11         &
@@ -805,6 +810,9 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
               ,tmfq          &
               ,qmicro        &
               ,forceqv_spechum &
+              ,betascu       &
+              ,betamcu       &
+              ,betadcu       &
               ,sigmain       &
               ,sigmaout      &
               ,ter11         &
diff --git a/physics/cu_c3_driver.meta b/physics/cu_c3_driver.meta
index 999b5c2bc..e02116243 100644
--- a/physics/cu_c3_driver.meta
+++ b/physics/cu_c3_driver.meta
@@ -244,6 +244,29 @@
   type = real
   kind = kind_phys
   intent = out
+[betascu]
+  standard_name = tuning_param_for_shallow_cu
+  long_name = tuning param for shallow cu in case prognostic closure is used
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[betamcu]
+  standard_name = tuning_param_for_midlevel_cu
+  long_name = tuning param for midlevel cu in case prognostic closure is used
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[betadcu]
+  standard_name = tuning_param_for_deep_cu
+  long_name = tuning param for deep cu in case prognostic closure is used
+  units = none
+  dimensions = ()
+  type = real
+  intent = in
 [phil]
   standard_name = geopotential
   long_name = layer geopotential
diff --git a/physics/cu_c3_sh.F90 b/physics/cu_c3_sh.F90
index a79e1dfcf..704f2a0fc 100644
--- a/physics/cu_c3_sh.F90
+++ b/physics/cu_c3_sh.F90
@@ -68,7 +68,8 @@ subroutine cu_c3_sh_run (                                            &
                          hfx,qfx,xland,ichoice,tcrit,dtime,         &
                          zuo,xmb_out,kbcon,ktop,k22,ierr,ierrc,     &
                          flag_init, flag_restart,fv,r_d,delp,tmf,qmicro, & 
-                         forceqv_spechum,sigmain,sigmaout,progsigma,dx,  &
+                         forceqv_spechum,betascu,betamcu,betadcu,sigmain,&
+                         sigmaout,progsigma,dx,  &
                          outt,outq,outqc,outu,outv,cnvwt,pre,cupclw,     & ! output tendencies
                          itf,ktf,its,ite, kts,kte,ipr,tropics)  ! dimesnional variables
 !
@@ -131,7 +132,7 @@ subroutine cu_c3_sh_run (                                            &
        
      real(kind=kind_phys)                                                              &
         ,intent (in   )                   ::                           &
-        dtime,tcrit,fv,r_d
+        dtime,tcrit,fv,r_d,betascu,betamcu,betadcu
 !$acc declare sigmaout                                                                                                                                                                                                                      
      real(kind=kind_phys),    dimension (its:,kts:)                              &
         ,intent (out)                     ::                           &
@@ -234,15 +235,18 @@ subroutine cu_c3_sh_run (                                            &
 !$acc       cap_max_increment,lambau,                                       &
 !$acc       kstabi,xland1,kbmax,ktopx)
 
-     logical :: flag_shallow
+     logical :: flag_shallow,flag_mid
      logical, dimension(its:ite) :: cnvflg
      integer                              ::                           &
        kstart,i,k,ki
-     real(kind=kind_phys)                                 ::                           &
+     real(kind=kind_phys)                 ::                           &
       dz,mbdt,zkbmax,                                                  &
       cap_maxs,trash,trash2,frh,el2orc,gravinv
       
-      real(kind=kind_phys) buo_flux,pgeoh,dp,entup,detup,totmas
+     real(kind=kind_phys) buo_flux,pgeoh,dp,entup,detup,totmas
+     real(kind=kind_phys)                 ::                           &
+          sigmind,sigminm,sigmins
+     parameter(sigmind=0.005,sigmins=0.03,sigminm=0.01)
 
      real(kind=kind_phys) xff_shal(3),blqe,xkshal
      character*50 :: ierrc(its:)
@@ -672,13 +676,13 @@ subroutine cu_c3_sh_run (                                            &
               dz=z_cup(i,k)-z_cup(i,k-1)
               ! cloud liquid water
               c1d(i,k)=c1_shal! 0. !.02*up_massdetr(i,k-1)
+              clw_all(i,k)=max(0.,qco(i,k)-trash)
               qrco(i,k)= (qco(i,k)-trash)/(1.+(c0_shal+c1d(i,k))*dz)
               if(qrco(i,k).lt.0.)then  ! hli new test 02/12/19
                  qrco(i,k)=0.
                  !c1d(i,k)=0.
               endif
               pwo(i,k)=c0_shal*dz*qrco(i,k)*zuo(i,k)
-              clw_all(i,k)=qco(i,k)-trash !LB total cloud before rain and detrain
               ! cloud water vapor 
               qco (i,k)= trash+qrco(i,k)
         
@@ -960,6 +964,7 @@ subroutine cu_c3_sh_run (                                            &
 ! equation 8, call progsigma_calc() to compute updraft area fraction based on a moisture budget
       if(progsigma)then
          flag_shallow = .true.
+         flag_mid = .false.
          do k=kts,ktf
             do i=its,itf
                del(i,k) = delp(i,k)*0.001
@@ -974,9 +979,9 @@ subroutine cu_c3_sh_run (                                            &
             endif
          enddo
          call progsigma_calc(itf,ktf,flag_init,flag_restart,flag_shallow,  &
-              del,tmf,qmicro,dbyo,zdqca,omega_u,zeta,xlv,dtime,            &
-              forceqv_spechum,kbcon,ktop,cnvflg,                           &
-              sigmain,sigmaout,sigmab)
+              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)
 
       endif
 
diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90
index c87308602..469df49f6 100644
--- a/physics/progsigma_calc.f90
+++ b/physics/progsigma_calc.f90
@@ -19,10 +19,10 @@ module progsigma
 !! This subroutine computes a prognostic updraft area fracftion
 !! used in the closure computations in the samfshalcnv. scheme
 !!\section gen_progsigma progsigma_calc General Algorithm 
-      subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
-           flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,    &
-           delt,qadv,kbcon1,ktcon,cnvflg,sigmain,sigmaout,           &
-           sigmab)
+      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,          &
+           sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
 !                                                           
 !                                                                                                                                             
       use machine,  only : kind_phys
@@ -32,11 +32,12 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
 
 !     intent in
       integer, intent(in)  :: im,km,kbcon1(im),ktcon(im)
-      real(kind=kind_phys), intent(in)  :: hvap,delt
+      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),    &
            qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km),           &
            omega_u(im,km),zeta(im,km)
-      logical, intent(in)  :: flag_init,flag_restart,cnvflg(im),flag_shallow
+      logical, intent(in)  :: flag_init,flag_restart,cnvflg(im),flag_shallow,flag_mid
       real(kind=kind_phys), intent(in) :: sigmain(im,km)
 
 !     intent out
@@ -53,15 +54,13 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
 
       real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2,   &
                           fdqb,dtdyn,dxlim,rmulacvg,tem,     &
-                          DEN,betascu,betadcu,dp1,invdelt
+                          DEN,dp1,invdelt
 
      !Parameters
       gcvalmx = 0.1
       rmulacvg=10.
       epsilon=1.E-11
       km1=km-1
-      betadcu = 2.0
-      betascu = 8.0
       invdelt = 1./delt
 
      !Initialization 2D
@@ -206,17 +205,27 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
          do i= 1, im
             if(cnvflg(i)) then
                sigmab(i)=sigmab(i)/betascu
-               sigmab(i)=MAX(0.03,sigmab(i))
+               sigmab(i)=MAX(sigmins,sigmab(i))
+            endif
+         enddo
+      elseif(flag_mid)then
+         do i= 1, im
+            if(cnvflg(i)) then
+               sigmab(i)=sigmab(i)/betamcu
+               sigmab(i)=MAX(sigminm,sigmab(i))
             endif
          enddo
       else
          do i= 1, im
             if(cnvflg(i)) then
                sigmab(i)=sigmab(i)/betadcu
-               sigmab(i)=MAX(0.01,sigmab(i))
+               sigmab(i)=MAX(sigmind,sigmab(i))
             endif
          enddo
       endif
+      do i= 1, im
+        sigmab(i) = MIN(0.95,sigmab(i))
+      enddo
 
      end subroutine progsigma_calc
 
diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f
index 8a36fe34c..e8faecf14 100644
--- a/physics/samfdeepcnv.f
+++ b/physics/samfdeepcnv.f
@@ -83,7 +83,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
      &    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, errmsg,errflg)
+     &    rainevap,sigmain,sigmaout,betadcu,betamcu,betascu,            &
+     &    errmsg,errflg)
 !
       use machine , only : kind_phys
       use funcphys , only : fpvs
@@ -100,14 +101,14 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
       real(kind=kind_phys), dimension(:), intent(in) :: fscav
       logical, intent(in)  :: first_time_step,restart,hwrf_samfdeep,    &
      &     progsigma
-      real(kind=kind_phys), intent(in) :: nthresh
+      real(kind=kind_phys), intent(in) :: nthresh,betadcu,betamcu,      &
+     &                                    betascu                       
       real(kind=kind_phys), intent(in) :: ca_deep(:)
       real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:),     &
      &     tmf(:,:,:),q(:,:), prevsq(:,:)
       real(kind=kind_phys), intent(out) :: rainevap(:)
       real(kind=kind_phys), intent(out) :: sigmaout(:,:)
       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
       real(kind=kind_phys), intent(inout) ::   qtr(:,:,:),              &
@@ -213,8 +214,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)
-      real(kind=kind_phys) gravinv,invdelt
-      logical flag_shallow
+      real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins
+      parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01)
+      logical flag_shallow, flag_mid
 c  physical parameters
 !     parameter(grav=grav,asolfac=0.958)
 !     parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp))
@@ -2930,10 +2932,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart,        &
          enddo
 
          flag_shallow = .false.
+         flag_mid = .false.
          call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
-     &        del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
-     &        qadv,kbcon1,ktcon,cnvflg,
-     &        sigmain,sigmaout,sigmab)
+     &        flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,
+     &        delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu,
+     &        sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
       endif
 
 !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer.
diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta
index bed4d655d..d0d39d830 100644
--- a/physics/samfdeepcnv.meta
+++ b/physics/samfdeepcnv.meta
@@ -450,6 +450,29 @@
   type = real
   kind = kind_phys
   intent = out
+[betascu]
+  standard_name = tuning_param_for_shallow_cu
+  long_name = tuning param for shallow cu in case prognostic closure is used
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[betamcu]
+  standard_name = tuning_param_for_midlevel_cu
+  long_name = tuning param for midlevel cu in case prognostic closure is used
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[betadcu]
+  standard_name = tuning_param_for_deep_cu
+  long_name = tuning param for deep cu in case prognostic closure is used
+  units = none
+  dimensions = ()
+  type = real
+  intent = in
 [qlcn]
   standard_name = mass_fraction_of_convective_cloud_liquid_water
   long_name = mass fraction of convective cloud liquid water
diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f
index a7682342f..3869ea6ea 100644
--- a/physics/samfshalcnv.f
+++ b/physics/samfshalcnv.f
@@ -57,7 +57,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
      &     rn,kbot,ktop,kcnv,islimsk,garea,                             &
      &     dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc,                       &
      &     clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal,                & 
-     &     sigmain,sigmaout,errmsg,errflg)
+     &     sigmain,sigmaout,betadcu,betamcu,betascu,errmsg,errflg)
 !
       use machine , only : kind_phys
       use funcphys , only : fpvs
@@ -67,7 +67,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
       integer, intent(in)  :: im, km, itc, ntc, ntk, ntr, ncloud
       integer, intent(in)  :: islimsk(:)
       real(kind=kind_phys), intent(in) :: cliq, cp, cvap,               &
-     &   eps, epsm1, fv, grav, hvap, rd, rv, t0c
+     &   eps, epsm1, fv, grav, hvap, rd, rv, t0c, betascu, betadcu,     &
+     &   betamcu
       real(kind=kind_phys), intent(in) ::  delt
       real(kind=kind_phys), intent(in) :: psp(:), delp(:,:),            &
      &   prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:),            &
@@ -159,8 +160,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
       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)
-      real(kind=kind_phys) gravinv,dxcrtas,invdelt
-      logical flag_shallow
+      real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins,
+     &                     sigminm
+      logical flag_shallow,flag_mid
 c  physical parameters
 !     parameter(g=grav,asolfac=0.89)
 !     parameter(g=grav)
@@ -194,7 +196,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
       parameter(betaw=.03,dxcrtc0=9.e3)
       parameter(h1=0.33333333)
 !  progsigma
-      parameter(dxcrtas=30.e3)
+      parameter(dxcrtas=30.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),
@@ -1974,10 +1976,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap,            &
          enddo
 
          flag_shallow = .true.
+         flag_mid = .false.
          call progsigma_calc(im,km,first_time_step,restart,flag_shallow,
-     &        del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,
-     &        qadv,kbcon1,ktcon,cnvflg,
-     &        sigmain,sigmaout,sigmab)
+     &        flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,
+     &        delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu,
+     &        sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab)
       endif
 
 !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity.
diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta
index c1fffef58..200e33707 100644
--- a/physics/samfshalcnv.meta
+++ b/physics/samfshalcnv.meta
@@ -482,6 +482,29 @@
   type = real
   kind = kind_phys
   intent = out
+[betascu]
+  standard_name = tuning_param_for_shallow_cu
+  long_name = tuning param for shallow cu in case prognostic closure is used
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[betamcu]
+  standard_name = tuning_param_for_midlevel_cu
+  long_name = tuning param for midlevel cu in case prognostic closure is used
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[betadcu]
+  standard_name = tuning_param_for_deep_cu
+  long_name = tuning param for deep cu in case prognostic closure is used
+  units = none
+  dimensions = ()
+  type = real
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP

From 12b400a210854d33e64fd6d211482d9f8ab7add5 Mon Sep 17 00:00:00 2001
From: Lisa Bengtsson <Lisa.Bengtsson@noaa.gov>
Date: Wed, 25 Oct 2023 21:30:26 +0000
Subject: [PATCH 2/4] Ensure prognostic closure is not used at coarse
 resolution

---
 physics/cu_c3_deep.F90    |  2 +-
 physics/cu_c3_driver.F90  | 16 ++++++++++++----
 physics/cu_c3_driver.meta |  7 +++++++
 3 files changed, 20 insertions(+), 5 deletions(-)

diff --git a/physics/cu_c3_deep.F90 b/physics/cu_c3_deep.F90
index 7e907aaba..b7cd5f62d 100644
--- a/physics/cu_c3_deep.F90
+++ b/physics/cu_c3_deep.F90
@@ -5758,7 +5758,7 @@ subroutine calculate_updraft_velocity(its,itf,ktf,ite,kts,kte,ierr,progsigma,
     do k = 2, ktf-1
        do i = 1, itf
          if (ierr(i)==0) then
-           if(k >= kbcon(i) .and. k < ktcon(i) .and. dbyo(i,k)>0.)then
+           if(k >= kbcon(i) .and. k < ktcon(i))then
               gamma   = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
               if(k >= kbcon(i) .and. clw_all(i,k)>0.)then
                buo(i,k) = buo(i,k) - g * qlk(i,k)  
diff --git a/physics/cu_c3_driver.F90 b/physics/cu_c3_driver.F90
index 0ecb81750..5b6be1d6c 100644
--- a/physics/cu_c3_driver.F90
+++ b/physics/cu_c3_driver.F90
@@ -58,7 +58,7 @@ end subroutine cu_c3_driver_init
 !!
 !>\section gen_c3_driver Grell-Freitas Cumulus Scheme Driver General Algorithm
       subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
-               do_ca,progsigma,cactiv,cactiv_m,g,cp,fv,r_d,xlv,r_v,forcet,      &
+               do_ca,progsigma,cnx,cactiv,cactiv_m,g,cp,fv,r_d,xlv,r_v,forcet,  &
                forceqv_spechum,phil,delp,raincv,tmf,qmicro,sigmain,             &
                betascu,betamcu,betadcu,qv_spechum,t,cld1d,us,vs,t2di,w,         &
                qv2di_spechum,p2di,psuri,                                        &
@@ -93,14 +93,14 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
       integer            :: ishallow_g3 ! depend on imfshalcnv
 !-------------------------------------------------------------
    integer      :: its,ite, jts,jte, kts,kte
-   integer, intent(in   ) :: im,km,ntracer
+   integer, intent(in   ) :: im,km,ntracer,cnx
    integer, intent(in   ) :: ichoice_in,ichoicem_in,ichoice_s_in
    logical, intent(in   ) :: flag_init, flag_restart, do_mynnedmf
    logical, intent(in   ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend, &
-        do_ca,progsigma
+        do_ca
    real (kind=kind_phys), intent(in) :: g,cp,fv,r_d,xlv,r_v,betascu,betamcu,betadcu
    logical, intent(in   ) :: ldiag3d
-
+   logical, intent(inout) :: progsigma
    real(kind=kind_phys), intent(inout)                      :: dtend(:,:,:)
 !$acc declare copy(dtend)
    integer, intent(in)                                      :: dtidx(:,:), &
@@ -280,6 +280,14 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
 !$acc end kernels
      endif
 
+
+     if(progsigma)then
+        if(cnx < 384)then
+           progsigma=.false.
+           write(*,*)'Forcing prognostic closure to .false. due to coarse resolution'
+        endif
+     endif
+
      if(ldiag3d) then
        if(flag_for_dcnv_generic_tend) then
          cliw_deep_idx=0
diff --git a/physics/cu_c3_driver.meta b/physics/cu_c3_driver.meta
index e02116243..71a785318 100644
--- a/physics/cu_c3_driver.meta
+++ b/physics/cu_c3_driver.meta
@@ -133,6 +133,13 @@
   units = flag
   dimensions = ()
   type = logical
+  intent = inout
+[cnx]
+  standard_name = number_of_x_points_for_current_cubed_sphere_tile
+  long_name = number of points in x direction for this cubed sphere face
+  units = count
+  dimensions = ()
+  type = integer
   intent = in	
 [cactiv]
   standard_name = counter_for_grell_freitas_convection

From 81563de9686260bc5f1c85fe350d48e56fdf7afc Mon Sep 17 00:00:00 2001
From: Dustin Swales <dustin.swales@noaa.gov>
Date: Thu, 26 Oct 2023 19:17:33 +0000
Subject: [PATCH 3/4] move setting of flag from run to init phase

---
 physics/cu_c3_driver.F90  | 26 ++++++++++++++------------
 physics/cu_c3_driver.meta | 21 ++++++++++++++-------
 2 files changed, 28 insertions(+), 19 deletions(-)

diff --git a/physics/cu_c3_driver.F90 b/physics/cu_c3_driver.F90
index 5b6be1d6c..c911ff5e4 100644
--- a/physics/cu_c3_driver.F90
+++ b/physics/cu_c3_driver.F90
@@ -30,7 +30,8 @@ module cu_c3_driver
 !! \htmlinclude cu_c3_driver_init.html
 !!
       subroutine cu_c3_driver_init(imfshalcnv, imfshalcnv_c3, imfdeepcnv, &
-                          imfdeepcnv_c3,mpirank, mpiroot, errmsg, errflg)
+                          imfdeepcnv_c3,progsigma, cnx, mpirank, mpiroot, &
+                          errmsg, errflg)
 
          implicit none
 
@@ -38,6 +39,8 @@ subroutine cu_c3_driver_init(imfshalcnv, imfshalcnv_c3, imfdeepcnv, &
          integer,                   intent(in) :: imfdeepcnv, imfdeepcnv_c3
          integer,                   intent(in)    :: mpirank
          integer,                   intent(in)    :: mpiroot
+         integer,                   intent(in)    :: cnx
+         logical,                   intent(inout) :: progsigma
          character(len=*),          intent(  out) :: errmsg
          integer,                   intent(  out) :: errflg
 
@@ -45,6 +48,13 @@ subroutine cu_c3_driver_init(imfshalcnv, imfshalcnv_c3, imfdeepcnv, &
          errmsg = ''
          errflg = 0
 
+         if(progsigma)then
+            if(cnx < 384)then
+               progsigma=.false.
+               write(*,*)'Forcing prognostic closure to .false. due to coarse resolution'
+            endif
+         endif
+
       end subroutine cu_c3_driver_init
 
 !
@@ -58,7 +68,7 @@ end subroutine cu_c3_driver_init
 !!
 !>\section gen_c3_driver Grell-Freitas Cumulus Scheme Driver General Algorithm
       subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
-               do_ca,progsigma,cnx,cactiv,cactiv_m,g,cp,fv,r_d,xlv,r_v,forcet,  &
+               do_ca,progsigma,cactiv,cactiv_m,g,cp,fv,r_d,xlv,r_v,forcet,      &
                forceqv_spechum,phil,delp,raincv,tmf,qmicro,sigmain,             &
                betascu,betamcu,betadcu,qv_spechum,t,cld1d,us,vs,t2di,w,         &
                qv2di_spechum,p2di,psuri,                                        &
@@ -93,14 +103,14 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
       integer            :: ishallow_g3 ! depend on imfshalcnv
 !-------------------------------------------------------------
    integer      :: its,ite, jts,jte, kts,kte
-   integer, intent(in   ) :: im,km,ntracer,cnx
+   integer, intent(in   ) :: im,km,ntracer
    integer, intent(in   ) :: ichoice_in,ichoicem_in,ichoice_s_in
    logical, intent(in   ) :: flag_init, flag_restart, do_mynnedmf
    logical, intent(in   ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend, &
         do_ca
    real (kind=kind_phys), intent(in) :: g,cp,fv,r_d,xlv,r_v,betascu,betamcu,betadcu
    logical, intent(in   ) :: ldiag3d
-   logical, intent(inout) :: progsigma
+   logical, intent(in   ) :: progsigma
    real(kind=kind_phys), intent(inout)                      :: dtend(:,:,:)
 !$acc declare copy(dtend)
    integer, intent(in)                                      :: dtidx(:,:), &
@@ -280,14 +290,6 @@ subroutine cu_c3_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
 !$acc end kernels
      endif
 
-
-     if(progsigma)then
-        if(cnx < 384)then
-           progsigma=.false.
-           write(*,*)'Forcing prognostic closure to .false. due to coarse resolution'
-        endif
-     endif
-
      if(ldiag3d) then
        if(flag_for_dcnv_generic_tend) then
          cliw_deep_idx=0
diff --git a/physics/cu_c3_driver.meta b/physics/cu_c3_driver.meta
index 71a785318..801b1e9d7 100644
--- a/physics/cu_c3_driver.meta
+++ b/physics/cu_c3_driver.meta
@@ -49,6 +49,20 @@
   dimensions = ()
   type = integer
   intent = in
+[progsigma]
+  standard_name = do_prognostic_updraft_area_fraction
+  long_name = flag for prognostic sigma in cumuls scheme
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = inout
+[cnx]
+  standard_name = number_of_x_points_for_current_cubed_sphere_tile
+  long_name = number of points in x direction for this cubed sphere face
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP
@@ -133,13 +147,6 @@
   units = flag
   dimensions = ()
   type = logical
-  intent = inout
-[cnx]
-  standard_name = number_of_x_points_for_current_cubed_sphere_tile
-  long_name = number of points in x direction for this cubed sphere face
-  units = count
-  dimensions = ()
-  type = integer
   intent = in	
 [cactiv]
   standard_name = counter_for_grell_freitas_convection

From e861277c1ffe8fdcb1b026240f98077bb7a91473 Mon Sep 17 00:00:00 2001
From: Lisa Bengtsson <Lisa.Bengtsson@noaa.gov>
Date: Thu, 26 Oct 2023 21:10:37 +0000
Subject: [PATCH 4/4] address review comments

---
 physics/cu_c3_sh.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/physics/cu_c3_sh.F90 b/physics/cu_c3_sh.F90
index 704f2a0fc..736292092 100644
--- a/physics/cu_c3_sh.F90
+++ b/physics/cu_c3_sh.F90
@@ -676,7 +676,7 @@ subroutine cu_c3_sh_run (                                            &
               dz=z_cup(i,k)-z_cup(i,k-1)
               ! cloud liquid water
               c1d(i,k)=c1_shal! 0. !.02*up_massdetr(i,k-1)
-              clw_all(i,k)=max(0.,qco(i,k)-trash)
+              clw_all(i,k)=max(0._kind_phys,qco(i,k)-trash)
               qrco(i,k)= (qco(i,k)-trash)/(1.+(c0_shal+c1d(i,k))*dz)
               if(qrco(i,k).lt.0.)then  ! hli new test 02/12/19
                  qrco(i,k)=0.