add shallow convection closure updates, add ntsigma in generic files

physics/GFS_DCNV_generic.F90

@@ -18,8 +18,8 @@ end subroutine GFS_DCNV_generic_pre_finalize
     subroutine GFS_DCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, do_cnvgwd, cplchm,  &
                                          gu0, gv0, gt0, gq0, nsamftrac, ntqv,            &
                                          save_u, save_v, save_t, save_q, clw,            &
-                                         ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,   &
-                                         ntgnc, nthl, nthnc, nthv, ntgv,                 &
+                                         ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,&
+                                         ntgl,ntgnc, nthl, nthnc, nthv, ntgv,            &
                                          cscnv, satmedmf, trans_trac, ras, ntrac,        &
                                          dtidx, index_of_process_dcnv, errmsg, errflg)
 
@@ -28,7 +28,7 @@ subroutine GFS_DCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, do_cnvgwd, cplc
       implicit none
 
       integer, intent(in) :: im, levs, nsamftrac, ntqv, index_of_process_dcnv, dtidx(:,:), &
-           ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntrac,ntgnc,nthl,nthnc,nthv,ntgv
+           ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntrac,ntgnc,nthl,nthnc,nthv,ntgv,ntsigma
       logical, intent(in) :: ldiag3d, qdiag3d, do_cnvgwd, cplchm
       real(kind=kind_phys), dimension(:,:),   intent(in)    :: gu0
       real(kind=kind_phys), dimension(:,:),   intent(in)    :: gv0
@@ -74,7 +74,7 @@ subroutine GFS_DCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, do_cnvgwd, cplc
                     n /= ntrw  .and. n /= ntsw  .and. n /= ntrnc   .and. &
                     n /= ntsnc .and. n /= ntgl  .and. n /= ntgnc   .and. &
                     n /= nthl  .and. n /= nthnc .and. n /= nthv    .and. &
-                    n /= ntgv ) then
+                    n /= ntgv .and. n /= ntsigma) then
                   tracers = tracers + 1
                   if(dtidx(100+n,index_of_process_dcnv)>0) then
                      save_q(:,:,n) = clw(:,:,tracers)
@@ -114,7 +114,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, &
       rainc, cldwrk, upd_mf, dwn_mf, det_mf, dtend, dtidx, index_of_process_dcnv, &
       index_of_temperature, index_of_x_wind, index_of_y_wind, ntqv, gq0, save_q,  &
       cnvw, cnvc, cnvw_phy_f3d, cnvc_phy_f3d, flag_for_dcnv_generic_tend,         &
-      ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,                               &
+      ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,                       &
       ntgnc, nthl, nthnc, nthv, ntgv, ntrac,clw,                                  &
       satmedmf, trans_trac, errmsg, errflg)
 
@@ -145,7 +145,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, &
       integer, intent(in) :: dtidx(:,:), index_of_process_dcnv, index_of_temperature, &
            index_of_x_wind, index_of_y_wind, ntqv
       integer, intent(in) :: ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,   &
-                             ntgnc, nthl, nthnc, nthv, ntgv, ntrac
+                             ntgnc, nthl, nthnc, nthv, ntgv, ntsigma, ntrac
       real(kind=kind_phys), dimension(:,:,:), intent(in) :: clw
 
 
@@ -212,7 +212,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, &
                      n /= ntrw  .and. n /= ntsw  .and. n /= ntrnc   .and. &
                      n /= ntsnc .and. n /= ntgl  .and. n /= ntgnc   .and. &
                      n /= nthl  .and. n /= nthnc .and. n /= nthv    .and. &
-                     n /= ntgv ) then
+                     n /= ntgv .and. n /= ntsigma) then
                    tracers = tracers + 1
                    idtend = dtidx(100+n,index_of_process_dcnv)
                    if(idtend>0) then

physics/GFS_DCNV_generic.meta

@@ -190,6 +190,13 @@
   dimensions = ()
   type = integer
   intent = in
+[ntsigma]
+  standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array
+  long_name = tracer index of updraft_area_fraction
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
 [ntrw]
   standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array
   long_name = tracer index for rain water
@@ -670,6 +677,13 @@
   dimensions = ()
   type = integer
   intent = in
+[ntsigma]
+  standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array
+  long_name = tracer index of updraft_area_fraction
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
 [ntrw]
   standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array
   long_name = tracer index for rain water

physics/GFS_MP_generic.meta

@@ -852,7 +852,7 @@
   type = logical
   intent = in
 [dqdt_qmicro]
-  standard_name = instantanious_moisture_tendency_due_to_microphysics
+  standard_name = instantaneous_moisture_tendency_due_to_microphysics
   long_name = moisture tendency due to microphysics
   units = kg kg-1 s-1
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)

physics/GFS_SCNV_generic.F90

@@ -17,14 +17,14 @@ end subroutine GFS_SCNV_generic_pre_finalize
       subroutine GFS_SCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, gu0, gv0, gt0, gq0, &
         save_u, save_v, save_t, save_q, ntqv, nsamftrac, flag_for_scnv_generic_tend,       &
         dtidx, index_of_process_scnv, ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,  &
-        cscnv, satmedmf, trans_trac, ras, ntrac, clw, errmsg, errflg)
+        ntsigma, cscnv, satmedmf, trans_trac, ras, ntrac, clw, errmsg, errflg)
 
         use machine,               only: kind_phys
 
         implicit none
 
         integer, intent(in) :: im, levs, ntqv, nsamftrac, index_of_process_scnv, dtidx(:,:)
-        integer, intent(in) :: ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac
+        integer, intent(in) :: ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac
         logical, intent(in) :: ldiag3d, qdiag3d, flag_for_scnv_generic_tend
         real(kind=kind_phys), dimension(:,:),   intent(in) :: gu0, gv0, gt0
         real(kind=kind_phys), dimension(:,:,:), intent(in) :: gq0
@@ -55,7 +55,7 @@ subroutine GFS_SCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, gu0, gv0, gt0,
                 do n=2,ntrac
                    if ( n /= ntcw  .and. n /= ntiw  .and. n /= ntclamt .and. &
                         n /= ntrw  .and. n /= ntsw  .and. n /= ntrnc   .and. &
-                        n /= ntsnc .and. n /= ntgl  .and. n /= ntgnc) then
+                        n /= ntsnc .and. n /= ntgl  .and. n /= ntgnc .and. n /= ntsigma) then
                       tracers = tracers + 1
                       if(dtidx(100+n,index_of_process_scnv)>0) then
                          save_q(:,:,n) = clw(:,:,tracers)
@@ -97,7 +97,7 @@ subroutine GFS_SCNV_generic_post_run (im, levs, nn, lssav, ldiag3d, qdiag3d, &
         rainc, cnvprcp, cnvprcpb, cnvw_phy_f3d, cnvc_phy_f3d,                      &
         dtend, dtidx, index_of_temperature, index_of_x_wind, index_of_y_wind,      &
         index_of_process_scnv, ntqv, flag_for_scnv_generic_tend,                   &
-        ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,                        &
+        ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,                &
         imfshalcnv, imfshalcnv_sas, imfshalcnv_samf, ntrac,                        &
         cscnv, satmedmf, trans_trac, ras, errmsg, errflg)
 
@@ -106,7 +106,7 @@ subroutine GFS_SCNV_generic_post_run (im, levs, nn, lssav, ldiag3d, qdiag3d, &
       implicit none
 
       integer, intent(in) :: im, levs, nn, ntqv, nsamftrac
-      integer, intent(in) :: ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac
+      integer, intent(in) :: ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac
       logical, intent(in) :: lssav, ldiag3d, qdiag3d, flag_for_scnv_generic_tend
       real(kind=kind_phys),                     intent(in) :: frain
       real(kind=kind_phys), dimension(:,:), intent(in) :: gu0, gv0, gt0
@@ -186,7 +186,7 @@ subroutine GFS_SCNV_generic_post_run (im, levs, nn, lssav, ldiag3d, qdiag3d, &
              do n=2,ntrac
                 if ( n /= ntcw  .and. n /= ntiw  .and. n /= ntclamt .and. &
                      n /= ntrw  .and. n /= ntsw  .and. n /= ntrnc   .and. &
-                     n /= ntsnc .and. n /= ntgl  .and. n /= ntgnc) then
+                     n /= ntsnc .and. n /= ntgl  .and. n /= ntgnc .and. n/= ntsigma) then
                    tracers = tracers + 1
                    idtend = dtidx(100+n,index_of_process_scnv)
                    if(idtend>0) then

physics/GFS_SCNV_generic.meta

@@ -183,6 +183,13 @@
   dimensions = ()
   type = integer
   intent = in
+[ntsigma]
+  standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array
+  long_name = tracer index of updraft_area_fraction
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
 [ntrw]
   standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array
   long_name = tracer index for rain water
@@ -614,6 +621,13 @@
   dimensions = ()
   type = integer
   intent = in
+[ntsigma]
+  standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array
+  long_name = tracer index of updraft_area_fraction
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
 [ntrw]
   standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array
   long_name = tracer index for rain water

physics/progsigma_calc.f90

@@ -1,7 +1,8 @@
 !>\file progsigma
 !! This file contains the subroutine that calculates the prognostic
 !! updraft area fraction that is used for closure computations in 
-!! saSAS deep and shallow convection.
+!! saSAS deep and shallow convection, based on a moisture budget
+!! as described in Bengtsson et al. 2022.
 
 !>\ingroup samfdeepcnv
 !! This subroutine computes a prognostic updraft area fraction
@@ -13,9 +14,8 @@
 !> @{ 
 
       subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
-           del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt,            &
-           qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx,                         &
-           do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep,     &
+           flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,    &
+           delt,qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx,                    &
            ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg)
 !                                                           
 !                                                                                                                                             
@@ -30,12 +30,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
       real,    intent(in)  :: qgrs_dsave(im,km), q(im,km),del(im,km),    &
            qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km),           &
            omega_u(im,km),zeta(im,km),gdx(im)
-      logical, intent(in)  :: flag_init,flag_restart,cnvflg(im)
-      real(kind=kind_phys), intent(in) :: nthresh
-      real(kind=kind_phys), intent(in) :: ca_deep(im)
+      logical, intent(in)  :: flag_init,flag_restart,cnvflg(im),flag_shallow
       real(kind=kind_phys), intent(out):: ca_micro(im)
-      logical, intent(in)  :: do_ca,ca_closure,ca_entr,ca_trigger
-
       real(kind=kind_phys), intent(in) :: sigmain(im,km)
 
 !     intent out
@@ -47,28 +43,29 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
 !     Local variables
       integer              :: i,k,km1
       real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im),   &
-                          mcons(im),zfdqa(im),zform(im,km),              &
+                          mcons(im),fdqa(im),form(im,km),              &
                           qadv(im,km),sigmamax(im)                         
 
 
-      real(kind=kind_phys) :: gcvalmx,ZEPS7,ZZ,ZCVG,mcon,buy2,   &
-                          zfdqb,dtdyn,dxlim,rmulacvg,dp,tem,     &
-                          alpha,DEN
+      real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2,   &
+                          fdqb,dtdyn,dxlim,rmulacvg,dp,tem,     &
+                          alpha,DEN,betascu
       integer :: inbu(im,km)  
 
      !Parameters
      gcvalmx = 0.1
      rmulacvg=10.
-     ZEPS7=1.E-11
+     epsilon=1.E-11
      km1=km-1
      alpha=7000.
+     betascu = 3.0
 
      !Initialization 2D
      do k = 1,km
         do i = 1,im
            sigmaout(i,k)=0.
            inbu(i,k)=0
-           zform(i,k)=0. 
+           form(i,k)=0. 
         enddo
      enddo
 
@@ -80,8 +77,9 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
          termB(i)=0.
          termC(i)=0.
          termD(i)=0.
-         zfdqa(i)=0.
+         fdqa(i)=0.
          mcons(i)=0.
+         ca_micro(i)=0.
       enddo
 
       !Initial computations, place maximum sigmain in sigmab
@@ -94,9 +92,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
              endif
           else
              if(cnvflg(i))then
-                !if(sigmain(i,k)<1.E-5)then
-                !   sigmain(i,k)=0.
-                !endif
                if(sigmain(i,k)>sigmab(i))then
                    sigmab(i)=sigmain(i,k)
                endif
@@ -107,7 +102,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
 
       do i=1,im
          if(sigmab(i) < 1.E-5)then !after advection
-            sigmab(i)=0.                                                                                                             
+            sigmab(i)=0.                                  
          endif
       enddo
 
@@ -180,11 +175,11 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
           do i = 1,im
              if(cnvflg(i))then
                 dp = 1000. * del(i,k)
-                zform(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt)
-                zfdqb=0.5*((zform(i,k)*zdqca(i,k)))
+                form(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt)
+                fdqb=0.5*((form(i,k)*zdqca(i,k)))
                 termC(i)=termC(i)+(float(inbu(i,k))*   &
-                     (zfdqb+zfdqa(i))*hvap*zeta(i,k))
-                zfdqa(i)=zfdqb
+                     (fdqb+fdqa(i))*hvap*zeta(i,k))
+                fdqa(i)=fdqb
              endif
          enddo
       enddo
@@ -193,29 +188,26 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
        do i = 1,im                                                                                                                           
           if(cnvflg(i))then
 
-             DEN=MIN(termC(i)+termB(i),1.E8) !1.E8
-             !DEN=MAX(termC(i)+termB(i),1.E7) !1.E7
-
-             ZCVG=termD(i)*delt
-
+             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)-ZEPS7))   
+                  *MAX(0.0,SIGN(1.0,termC(i)-epsilon))   
 
 
-             ZCVG=MAX(0.0,ZCVG)
+             cvg=MAX(0.0,cvg)
 
-             if(flag_init)then
+             if(flag_init .and. .not. flag_restart)then
                 sigmab(i)=0.03
              else
-                sigmab(i)=(ZZ*(termA(i)+ZCVG))/(DEN+(1.0-ZZ))
+                sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ))
              endif
 
              if(sigmab(i)>0.)then
                 sigmab(i)=MIN(sigmab(i),sigmamax(i))  
                 sigmab(i)=MAX(sigmab(i),0.01)
              endif
-
+             ca_micro(i)=sigmab(i)
           endif!cnvflg
        enddo
 
@@ -226,7 +218,20 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,           &
              endif
           enddo
        enddo
+
+       !Since updraft velocity is much lower in shallow cu region, termC becomes small in shallow cu application, thus the area fraction 
+       !in this regime becomes too large compared with the deep cu region. To address this simply apply a scaling factor for shallow cu 
+       !before computing the massflux to reduce the total strength of the SC MF:
 
+       if(flag_shallow)then
+          do i= 1, im
+             if(cnvflg(i)) then
+                sigmab(i)=sigmab(i)/betascu
+             endif
+          enddo
+       endif
+
+
      end subroutine progsigma_calc
 !> @}                            
 !! @}