Add prognostic cumulus closure description in saSAS

ccpp/data/CCPP_typedefs.F90

@@ -1141,8 +1141,12 @@ subroutine gfs_interstitial_setup_tracers(Interstitial, Model)
     if (Model%ntke > 0) Interstitial%ntkev = Interstitial%nvdiff
 
     if (Model%ntiw > 0) then
-      if (Model%ntclamt > 0) then
+      if (Model%ntclamt > 0 .and. Model%ntsigma <= 0) then
         Interstitial%nn = Model%ntrac - 2
+      elseif (Model%ntclamt <= 0 .and. Model%ntsigma > 0) then
+        Interstitial%nn = Model%ntrac - 2
+      elseif  (Model%ntclamt > 0 .and. Model%ntsigma > 0) then
+        Interstitial%nn = Model%ntrac - 3
       else
         Interstitial%nn = Model%ntrac - 1
       endif
@@ -1162,7 +1166,8 @@ subroutine gfs_interstitial_setup_tracers(Interstitial, Model)
              n /= Model%ntrw  .and. n /= Model%ntsw  .and. n /= Model%ntrnc   .and. &
              n /= Model%ntsnc .and. n /= Model%ntgl  .and. n /= Model%ntgnc   .and. &
              n /= Model%nthl  .and. n /= Model%nthnc .and. n /= Model%ntgv    .and. &
-             n /= Model%nthv  .and. n /= Model%ntccn .and. n /= Model%ntccna )
+             n /= Model%nthv  .and. n /= Model%ntccn .and. n /= Model%ntccna .and.  &
+             n /= Model%ntsigma)
         Interstitial%otsptflag(n) = ltest
         if ( ltest ) then
           tracers = tracers + 1

ccpp/data/GFS_typedefs.F90

@@ -569,6 +569,10 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: pfi_lsan(:,:)   => null()  !< instantaneous 3D flux of ice    nonconvective precipitation (kg m-2 s-1)
     real (kind=kind_phys), pointer :: pfl_lsan(:,:)   => null()  !< instantaneous 3D flux of liquid nonconvective precipitation (kg m-2 s-1)
 
+    !-- prognostic updraft area fraction coupling in convection
+    real (kind=kind_phys), pointer :: tmf (:,:) => null()         !< tmf to be passed from turublence scheme to convection
+    real (kind=kind_phys), pointer :: dqdt_qmicro(:,:) => null()  !< instantanious microphysics tendency to be passed from MP to convection
+
     !--- instantaneous total moisture tendency for smoke coupling:
     real (kind=kind_phys), pointer :: dqdti   (:,:)   => null()  !< rrfs_smoke=true only; instantaneous total moisture tendency (kg/kg/s)
 
@@ -1014,6 +1018,7 @@ module GFS_typedefs
     integer              :: imfshalcnv_gf       = 3 !< flag for scale- & aerosol-aware Grell-Freitas scheme (GSD)
     integer              :: imfshalcnv_ntiedtke = 4 !< flag for new Tiedtke scheme (CAPS)
     logical              :: hwrf_samfdeep           !< flag for HWRF SAMF deepcnv scheme (HWRF)
+    logical              :: progsigma               !< flag for prognostic area fraction in samf ddepcnv scheme (GFS)   
     integer              :: imfdeepcnv      !< flag for mass-flux deep convection scheme
                                             !<     1: July 2010 version of SAS conv scheme
                                             !<           current operational version as of 2016
@@ -1280,6 +1285,7 @@ module GFS_typedefs
     integer              :: ntgv            !< tracer index for graupel particle volume
     integer              :: nthv            !< tracer index for hail particle volume
     integer              :: ntke            !< tracer index for kinetic energy
+    integer              :: ntsigma         !< tracer index for updraft area fraction  
     integer              :: nto             !< tracer index for oxygen ion
     integer              :: nto2            !< tracer index for oxygen
     integer              :: ntwa            !< tracer index for water friendly aerosol
@@ -2690,6 +2696,14 @@ subroutine coupling_create (Coupling, IM, Model)
       Coupling%psurfi_cpl  = clear_val
     endif
 
+    !--prognostic closure - moisture coupling
+    if(Model%progsigma)then
+       allocate(Coupling%dqdt_qmicro (IM,Model%levs))
+       allocate(Coupling%tmf (IM,Model%levs))
+       Coupling%tmf = clear_val
+       Coupling%dqdt_qmicro = clear_val
+    endif
+
     !--- stochastic physics option
     if (Model%do_sppt .or. Model%ca_global) then
       allocate (Coupling%sppt_wts  (IM,Model%levs))
@@ -3153,6 +3167,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
 
     logical              :: hwrf_samfdeep     = .false.               !< flag for HWRF SAMF deepcnv scheme
     logical              :: hwrf_samfshal     = .false.               !< flag for HWRF SAMF shalcnv scheme
+    logical              :: progsigma         = .false.               !< flag for prognostic updraft area fraction closure in saSAS
     logical              :: do_mynnedmf       = .false.               !< flag for MYNN-EDMF
     logical              :: do_mynnsfclay     = .false.               !< flag for MYNN Surface Layer Scheme
     ! DH* TODO - move to MYNN namelist section
@@ -3455,7 +3470,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
                                do_ugwp_v1, do_ugwp_v1_orog_only,  do_ugwp_v1_w_gsldrag,     &
                                var_ric, coef_ric_l, coef_ric_s, hurr_pbl,                   &
                                do_myjsfc, do_myjpbl,                                        &
-                               hwrf_samfdeep, hwrf_samfshal,                                &
+                               hwrf_samfdeep, hwrf_samfshal,progsigma,                      &
                                h2o_phys, pdfcld, shcnvcw, redrag, hybedmf, satmedmf,        &
                                shinhong, do_ysu, dspheat, lheatstrg, lseaspray, cnvcld,     &
                                random_clds, shal_cnv, imfshalcnv, imfdeepcnv, isatmedmf,    &
@@ -4109,6 +4124,12 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%hwrf_samfdeep = hwrf_samfdeep
     Model%hwrf_samfshal = hwrf_samfshal
 
+    if (progsigma .and. imfdeepcnv/=2) then
+       write(*,*) 'Logic error: progsigma requires imfdeepcnv=2'
+       stop
+    end if
+    Model%progsigma = progsigma
+
     if (oz_phys .and. oz_phys_2015) then
        write(*,*) 'Logic error: can only use one ozone physics option (oz_phys or oz_phys_2015), not both. Exiting.'
        stop
@@ -4389,6 +4410,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%ntgv             = get_tracer_index(Model%tracer_names, 'graupel_vol',Model%me, Model%master, Model%debug)
     Model%nthv             = get_tracer_index(Model%tracer_names, 'hail_vol',   Model%me, Model%master, Model%debug)
     Model%ntke             = get_tracer_index(Model%tracer_names, 'sgs_tke',    Model%me, Model%master, Model%debug)
+    Model%ntsigma          = get_tracer_index(Model%tracer_names, 'sigmab',     Model%me, Model%master, Model%debug)
     Model%nqrimef          = get_tracer_index(Model%tracer_names, 'q_rimef',    Model%me, Model%master, Model%debug)
     Model%ntwa             = get_tracer_index(Model%tracer_names, 'liq_aero',   Model%me, Model%master, Model%debug)
     Model%ntia             = get_tracer_index(Model%tracer_names, 'ice_aero',   Model%me, Model%master, Model%debug)
@@ -5987,6 +6009,7 @@ subroutine control_print(Model)
       print *, ' ntgv              : ', Model%ntgv
       print *, ' nthv              : ', Model%nthv
       print *, ' ntke              : ', Model%ntke
+      print *, ' ntsigma           : ', Model%ntsigma
       print *, ' nto               : ', Model%nto
       print *, ' nto2              : ', Model%nto2
       print *, ' ntwa              : ', Model%ntwa
@@ -6228,6 +6251,11 @@ subroutine tbd_create (Tbd, IM, Model)
     allocate (Tbd%hpbl (IM))
     Tbd%hpbl     = clear_val
 
+    if (Model%imfdeepcnv == Model%imfdeepcnv_gf .or. Model%imfdeepcnv == Model%imfdeepcnv_ntiedtke .or. Model%imfdeepcnv == Model%imfdeepcnv_samf .or. Model%imfshalcnv == Model%imfshalcnv_samf) then
+       allocate(Tbd%prevsq(IM, Model%levs))
+       Tbd%prevsq = clear_val
+    endif
+
     if (Model%imfdeepcnv .ge. 0 .or. Model%imfshalcnv .ge. 0) then
        allocate(Tbd%ud_mf(IM, Model%levs))
        Tbd%ud_mf = zero
@@ -6236,12 +6264,11 @@ subroutine tbd_create (Tbd, IM, Model)
     if (Model%imfdeepcnv == Model%imfdeepcnv_gf .or. Model%imfdeepcnv == Model%imfdeepcnv_ntiedtke) then
        allocate(Tbd%forcet(IM, Model%levs))
        allocate(Tbd%forceq(IM, Model%levs))
+       allocate(Tbd%forcet(IM, Model%levs))
        allocate(Tbd%prevst(IM, Model%levs))
-       allocate(Tbd%prevsq(IM, Model%levs))
        Tbd%forcet = clear_val
        Tbd%forceq = clear_val
        Tbd%prevst = clear_val
-       Tbd%prevsq = clear_val
     end if
 
     if (Model%imfdeepcnv == Model%imfdeepcnv_gf) then

ccpp/data/GFS_typedefs.meta

@@ -309,6 +309,14 @@
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)
   type = real
   kind = kind_phys
+[qgrs(:,:,index_of_updraft_area_fraction_in_tracer_concentration_array)]
+  standard_name = prognostic_updraft_area_fraction_in_convection
+  long_name = convective updraft area fraction
+  units = frac
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  active = (index_of_updraft_area_fraction_in_tracer_concentration_array > 0 )
 [qgrs(:,:,index_for_smoke_in_tracer_concentration_array)]
   standard_name = smoke_tracer_concentration
   long_name = concentration of smoke
@@ -559,7 +567,14 @@
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)
   type = real
   kind = kind_phys
-
+[gq0(:,:,index_of_updraft_area_fraction_in_tracer_concentration_array)]
+  standard_name = updraft_area_fraction_updated_by_physics
+  long_name = convective updraft area fraction updated by physics
+  units = frac
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  active = ( index_of_updraft_area_fraction_in_tracer_concentration_array > 0 )
 ########################################################################
 [ccpp-table-properties]
   name = GFS_sfcprop_type
@@ -2330,6 +2345,22 @@
   type = real
   kind = kind_phys
   active = (flag_for_surface_flux_coupling)
+[tmf]
+  standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
+  long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL
+  units = kg kg-1 s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  active = (do_prognostic_updraft_area_fraction)
+[dqdt_qmicro]
+  standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics
+  long_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics
+  units = kg kg-1 s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  active = (do_prognostic_updraft_area_fraction)
 [ca_deep]
   standard_name = cellular_automata_area_fraction_for_deep_convection_from_coupled_process
   long_name = fraction of cellular automata for deep convection
@@ -4509,6 +4540,12 @@
   units = flag
   dimensions = ()
   type = logical
+[progsigma]
+  standard_name = do_prognostic_updraft_area_fraction
+  long_name = do_prognostic_updraft_area_fraction
+  units = flag
+  dimensions = ()
+  type = logical
 [isatmedmf]
   standard_name = choice_of_scale_aware_TKE_moist_EDMF_PBL
   long_name = choice of scale-aware TKE moist EDMF PBL scheme
@@ -5383,6 +5420,12 @@
   units = index
   dimensions = ()
   type = integer
+[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
 [nqrimef]
   standard_name = index_of_mass_weighted_rime_factor_in_tracer_concentration_array
   long_name = tracer index for mass weighted rime factor
@@ -6750,7 +6793,7 @@
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)
   type = real
   kind = kind_phys
-  active = (control_for_deep_convection_scheme == identifier_for_grell_freitas_deep_convection .or. control_for_deep_convection_scheme == identifier_for_new_tiedtke_deep_convection)
+  active = (control_for_deep_convection_scheme == identifier_for_grell_freitas_deep_convection .or. control_for_deep_convection_scheme == identifier_for_new_tiedtke_deep_convection .or. control_for_deep_convection_scheme == identifer_for_scale_aware_mass_flux_deep_convection .or. control_for_shallow_convection_scheme == identifier_for_scale_aware_mass_flux_shallow_convection)
 [cactiv]
   standard_name = counter_for_grell_freitas_convection
   long_name = convective activity memory

ccpp/driver/GFS_restart.F90

@@ -144,6 +144,10 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, &
     if (Model%do_mynnedmf) then
       Restart%num3d = Restart%num3d + 9
     endif
+    !Prognostic area fraction
+    if (Model%progsigma) then
+       Restart%num3d = Restart%num3d + 2
+    endif
 
     if (Model%num_dfi_radar > 0) then
       do itime=1,Model%dfi_radar_max_intervals
@@ -437,6 +441,20 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, &
        num = Model%ntot3d
     endif
 
+    !Prognostic closure
+    if(Model%progsigma)then
+       num = num + 1
+       Restart%name3d(num) = 'sas_3d_qgrs_dsave'
+      do nb = 1,nblks
+        Restart%data(nb,num)%var3p => Tbd(nb)%prevsq(:,:)
+      enddo
+      num = num + 1
+      Restart%name3d(num) = 'sas_3d_dqdt_qmicro'
+      do nb = 1,nblks
+        Restart%data(nb,num)%var3p => Coupling(nb)%dqdt_qmicro(:,:)
+      enddo
+    endif
+
     !--Convection variable used in CB cloud fraction. Presently this
     !--is only needed in sgscloud_radpre for imfdeepcnv == imfdeepcnv_gf.
     if (Model%imfdeepcnv == Model%imfdeepcnv_gf) then
@@ -446,6 +464,7 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, &
         Restart%data(nb,num)%var3p => Tbd(nb)%ud_mf(:,:)
       enddo
     endif
+
     !--- RAP/HRRR-specific variables, 3D
     ! GF
     if (Model%imfdeepcnv == Model%imfdeepcnv_gf) then

stochastic_physics/stochastic_physics_wrapper.F90

@@ -185,11 +185,6 @@ subroutine stochastic_physics_wrapper (GFS_Control, GFS_Data, Atm_block, ierr)
       end if
       ! Consistency check for cellular automata
       if(GFS_Control%do_ca)then
-        ! DH* The current implementation of cellular_automata assumes that all blocksizes are the
-        ! same - abort if this is not the case, otherwise proceed with Atm_block%blksz(1) below
-        if (.not. minval(Atm_block%blksz) == maxblk) then
-           call mpp_error(FATAL, 'Logic errror: cellular_automata not compatible with non-uniform blocksizes')
-        end if
         if(GFS_Control%ca_sgs)then
            allocate(sst         (1:nblks, maxblk))
            allocate(lmsk        (1:nblks, maxblk))