Improved blending of 1D and 3D turbulence schemes

documentation/source/science_guide/turbulence_schemes/bl_scheme_doc.rst

@@ -1986,7 +1986,7 @@ length, :math:`S` is the wind shear, :math:`f_\chi(Ri)` is the stability
 function and :math:`\chi` represents conserved heat and moisture
 variables, or momentum. Both schemes use the same stability function,
 and both schemes can use the full 3D shear for :math:`S`. Therefore the
-only difference is in the mixing length, which is calculated as
+only difference is in the mixing length, which is initially calculated as
 
 .. math:: :label: eq-lblend
 
@@ -2119,18 +2119,27 @@ boundary layer depth and type <sec_types>`)
 but alternatively (``blending_option``\ :math:`=`\ 4) this can be
 restricted to strictly shallow cumulus clouds, defined as contiguously
 cloudy levels (cloud fraction greater than SC_CFTOL) with cloud top
-height below input parameter ``shallow_cu_maxtop``. Note that the
+height below the input parameter ``shallow_cu_maxtop``. Note that the
 diagnosis of shallow cumulus from the diagnosis parcel ascent (that was
-used to identify a cumulus regime) was found frequently to indicate deep
+used to identify a cumulus regime) has frequently been found to indicate deep
 convection even when the resolved clouds were shallow because the
 diagnosis parcel, being undilute, would penetrate to the tropopause.
-However, having decided the regime is shallow convection, we do still
-set :math:`z_{\mathrm{turb}}` to the diagnosis parcel top height because, for
-current km-scale configurations (without a cumulus convection
-parametrization), it was found that the resulting stronger parametrized
-vertical mixing was beneficial for the development of the convection,
-and that without this a widespread stratiform cloud layer could develop
-instead.
+Under ``blending_option``\ :math:`=`\ 4, then, if such a shallow cumulus
+regime is diagnosed, :math:`z_{\mathrm{turb}}` is set to the resolved cloud
+top height.
+
+One less desirable feature of the combination of the forms of the blended
+mixing length :eq:`eq-lblend` and :math:`W_{1D}` :eq:`eq-tanh` functions
+is that the resulting mixing length does not simply increase monotonically
+towards :math:`l_{\mathrm{bl}}` as :math:`z_{\mathrm{turb}}` reduces (and
+the turbulence becomes entirely unresolved).   To illustrate this, consider
+the case of :math:`\Delta x= z_{\mathrm{turb}} = 1` km, which gives
+:math:`W_{1D} \approx 0.9` (see :numref:`Figure %s <fig-blend>`\ a).
+Ignoring the reduction in length scales towards the surface, we find
+:math:`l_{\mathrm{bl}}=150` m and :math:`l_{\mathrm{smag}}=200` m such that
+:math:`l_{\mathrm{blend}} = 155` m.  A simple way to avoid this is to cap
+the blended length scale explicitly by :math:`l_{\mathrm{bl}}` via the
+namelist switch ``cap_blended_ml``.
 
 Above the boundary layer top, `Boutle et al. (2014)`_ aimed
 for any free atmospheric mixing to be done by the 3D Smagorinsky scheme.

interfaces/physics_schemes_interface/source/support/um_physics_init_mod.f90

@@ -194,14 +194,17 @@ module um_physics_init_mod
                                         heavy_rain_evap_fac_in =>            &
                                                 heavy_rain_evap_fac
 
-  use mixing_config_mod,         only : smagorinsky,                 &
-                                        mixing_method => method,     &
-                                        method_3d_smag,              &
-                                        method_2d_smag,              &
-                                        method_blend_smag_fa,        &
-                                        method_blend_1dbl_fa,        &
-                                        mix_factor_in => mix_factor, &
-                                        leonard_term
+  use mixing_config_mod, only : smagorinsky,                               &
+                                mixing_method => method,                   &
+                                method_3d_smag,                            &
+                                method_2d_smag,                            &
+                                method_blend_smag_fa,                      &
+                                method_blend_1dbl_fa,                      &
+                                method_blend_cth_shcu_only,                &
+                                cap_blended_ml_in => cap_blended_ml,       &
+                                shallow_cu_maxtop_in => shallow_cu_maxtop, &
+                                mix_factor_in => mix_factor,               &
+                                leonard_term
 
   use radiation_config_mod,      only : topography, topography_horizon
 
@@ -370,6 +373,7 @@ subroutine um_physics_init()
          l_noice_in_turb, l_use_var_fixes,                                 &
          i_interp_local_cf_dbdz, tke_diag_fac, a_ent_2, dec_thres_cloud,   &
          dec_thres_cu, near_neut_z_on_l, blend_gridindep_fa,               &
+         blend_cth_shcu_only, shallow_cu_maxtop, cap_blended_ml,           &
          specified_fluxes_tstar, buoy_integ_low, num_sweeps_bflux,         &
          l_use_sml_dsc_fixes, l_converge_ga, improved_tke_diag
     use cloud_inputs_mod, only: i_cld_vn, forced_cu, i_rhcpt, i_cld_area,  &
@@ -1600,6 +1604,9 @@ subroutine um_physics_init()
       mix_factor = real(mix_factor_in, r_um)
       turb_startlev_vert  = 2
       turb_endlev_vert    = bl_levels
+      if (mixing_method /= method_3d_smag) then
+        cap_blended_ml = cap_blended_ml_in
+      end if
 
       ! Options which are bespoke to the choice of scheme
       select case ( mixing_method )
@@ -1622,6 +1629,11 @@ subroutine um_physics_init()
         l_subfilter_horiz = .true.
         l_subfilter_vert  = .true.
         blending_option   = blend_gridindep_fa
+      case( method_blend_cth_shcu_only )
+        l_subfilter_horiz = .true.
+        l_subfilter_vert  = .true.
+        blending_option   = blend_cth_shcu_only
+        shallow_cu_maxtop = real(shallow_cu_maxtop_in, r_bl)
       end select
 
     else ! not Smagorinsky

rose-stem/app/lfric_atm/opt/README

@@ -56,6 +56,8 @@ Science settings
 - ral3 - regional configuration
 - ral3_ens - extra config for regional ensembles (including stochastic physics)
 - ral3_scm - single column regional physics (includes convection scheme)
+- ral_blendshcu - development version of the turbulence blending for shallow
+                  cumulus regimes
 - ral_mixmol_opts - extra options for regional models when using the MoL
                     transport scheme and mixed solver configuration
 - ral_opts - extra options for regional models, including LBC blending weights

rose-stem/app/lfric_atm/opt/rose-app-blendshcu.conf

@@ -0,0 +1,5 @@
+[namelist:mixing]
+cap_blended_ml=.true.
+method='blend_cth_shcu_only'
+shallow_cu_maxtop=3000.0
+smagorinsky=.true.

rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc

@@ -788,6 +788,20 @@
         "xios_min_buffer_size": 2097152,
     }) %}
 
+{% elif task_ns.conf_name == "ral_dev-seuk_MG" %}
+
+    {% do task_dict.update({
+        "opt_confs": ["l70_40km","ral3","blendshcu","ral_opts","hh_solve"],
+        "resolution": "seuk_MG",
+        "DT": 60,
+        "tsteps": 120,
+        "crun": 2,
+        "crun_compare": false,
+        "mpi_parts": 16,
+        "kgo_checks": ["checksum"],
+        "plot_str": "plot_map.py -ral $NODAL_DATA_DIR/lfric_diagnostics.nc  $PLOT_DIR",
+    }) %}
+
 {% elif task_ns.conf_name == "ral3_mixmol-seuk_MG" %}
 
     {% do task_dict.update({

rose-stem/site/meto/groups/groups_lfric_atm.cylc

@@ -37,6 +37,7 @@
     "lfric_atm_ral_azspice_developer": [
         "lfric_atm_ral3-seuk_MG_azspice_gnu_fast-debug-32bit",
         "lfric_atm_ral3_ens-seuk_MG_azspice_gnu_fast-debug-32bit",
+        "lfric_atm_ral_dev-seuk_MG_azspice_gnu_fast-debug-32bit",
         "lfric_atm_ral3_mixmol-seuk_MG_azspice_gnu_fast-debug-32bit",
     ],
     "lfric_atm_scm_azspice_developer": [
@@ -140,6 +141,7 @@
     "lfric_atm_ral_ex1a_developer": [
         "lfric_atm_ral3-seuk_MG_ex1a_cce_fast-debug-32bit",
         "lfric_atm_ral3_ens-seuk_MG_ex1a_cce_fast-debug-32bit",
+        "lfric_atm_ral_dev-seuk_MG_ex1a_cce_fast-debug-32bit",
         "lfric_atm_ral3_mixmol-seuk_MG_ex1a_cce_fast-debug-32bit",
     ],
     "lfric_atm_ral_ex1a_weekly": [

rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_ral_dev-seuk_MG_azspice_gnu_fast-debug-32bit.txt

@@ -0,0 +1,9 @@
+Inner product checksum rho =         4839977E
+Inner product checksum theta =         5135493C
+Inner product checksum u =         612F4556
+Inner product checksum mr1 =         4090F41B
+Inner product checksum mr2 =         35AC7EDC
+Inner product checksum mr3 =         2FA559EC
+Inner product checksum mr4 =         33F4FDAE
+Inner product checksum mr5 =          BEDF860
+Inner product checksum mr6 =                0

rose-stem/site/meto/kgos/lfric_atm/ex1a/checksum_lfric_atm_ral_dev-seuk_MG_ex1a_cce_fast-debug-32bit.txt

@@ -0,0 +1,9 @@
+Inner product checksum rho =         483997BF
+Inner product checksum theta =         5135495F
+Inner product checksum u =         612F491D
+Inner product checksum mr1 =         4090F962
+Inner product checksum mr2 =         35B6E523
+Inner product checksum mr3 =         2FA5A395
+Inner product checksum mr4 =         33F4FE66
+Inner product checksum mr5 =          BECFC40
+Inner product checksum mr6 =                0

rose-stem/site/meto/lfric_atm/tasks_lfric_atm_ex1a.cylc

@@ -64,6 +64,12 @@
         "memory": [24, "GB"],
     }) %}
 
+{% elif task_ns.conf_name == "ral_dev-seuk_MG" %}
+
+    {% do task_dict.update({
+        "memory": [24, "GB"],
+    }) %}
+
 {% elif task_ns.conf_name == "ral3_mixmol-seuk_MG" %}
 
     {% do task_dict.update({

science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf

@@ -3768,7 +3768,7 @@ type=real
 
 [namelist:mixing=method]
 compulsory=true
-description=Smagorinsky subgrid mixing scheme option
+description=Smagorinsky (Smag) subgrid mixing scheme option
 !enumeration=true
 help=Options for Smagorinsky subgrid mixing scheme:
     =_________________________________________________________________
@@ -3786,12 +3786,42 @@ help=Options for Smagorinsky subgrid mixing scheme:
     =_________________________________________________________________
     =Blended scheme + 1D BL above
     = As above, but relaxes to 1D BL above BL
+    =_________________________________________________________________
+    =Blended scheme only using cloud top in shallow clouds
+    = As above, but use cloud-top height as blending length-scale only
+    = for shallow cumulus clouds
     =
-    =  WARNING: This option currently does not conserve energy and as
-    =  a result may be unstable.
+    =  WARNING: This option will only conserve if
+    =           conservative_diffusion is true
 sort-key=Panel-A05
-value-titles=3D Smagorinsky, 2D Smagorinsky + 1D BL, Blended scheme + Smag above, Blended scheme + 1D BL above
-values='3d_smag', '2d_smag', 'blend_smag_fa','blend_1dbl_fa'
+trigger=namelist:mixing=shallow_cu_maxtop: 'blend_cth_shcu_only' ;
+       =namelist:mixing=cap_blended_ml: 'blend_smag_fa', 'blend_1dbl_fa', 'blend_cth_shcu_only' ;
+value-titles=3D Smag, 2D Smag + 1D BL, Blended scheme + Smag above,
+            =Blended scheme + 1D BL above, Blended scheme using cloud top only in shallow clouds
+values='3d_smag', '2d_smag', 'blend_smag_fa','blend_1dbl_fa','blend_cth_shcu_only'
+
+[namelist:mixing=cap_blended_ml]
+compulsory=true
+description=Cap the blended mixing length at the 1D
+help=In the the blended subgrid mixing scheme (namelist:mixing=method) the physical
+    =mixing length from the one-dimensional scheme is combined with the grid-dependent
+    =Smagorinsky mixing length.  This option caps this blended mixing length to be at
+    =most the one dimensional mixing length
+!kind=default
+sort-key=Panel-A05a
+type=logical
+
+[namelist:mixing=shallow_cu_maxtop]
+compulsory=true
+description=Max. permitted cloud top height for shallow cloud blending
+help=The maximum permitted cloud top height (in metres) for shallow cumulus clouds
+    =to be diagnosed.  This is used in the blending weight length-scale calculation
+    =when "Blended scheme using cloud top only in shallow clouds" is selected for
+    =the blended subgrid mixing scheme (namelist:mixing=method).
+!kind=default
+range=0.0:
+sort-key=Panel-A05b
+type=real
 
 [namelist:mixing=mix_factor]
 compulsory=true

science/gungho/rose-meta/lfric-gungho/versions.py

@@ -28,8 +28,6 @@ def upgrade(self, config, meta_config=None):
         # Add settings
         return config, self.reports
 """
-
-
 class vn32_t634(MacroUpgrade):
     """Upgrade macro for ticket #634 by Ian Boutle."""
 
@@ -61,3 +59,18 @@ def upgrade(self, config, meta_config=None):
         self.add_setting(config, ["namelist:mixing", "fullstress"], ".false.")
 
         return config, self.reports
+
+
+class vn32_t752(MacroUpgrade):
+    # Upgrade macro for #752 by Adrian Lock
+
+    BEFORE_TAG = "vn3.2_t479"
+    AFTER_TAG = "vn3.2_t752"
+
+    def upgrade(self, config, meta_config=None):
+        # Add settings
+        self.add_setting(
+            config, ["namelist:mixing", "cap_blended_ml"], ".false.")
+        self.add_setting(
+            config, ["namelist:mixing", "shallow_cu_maxtop"], "3000.0")
+        return config, self.reports

science/gungho/unit-test/kernel/diffusion/diffusion_momentum_flux_kernel_mod_test.pf

@@ -128,7 +128,9 @@ contains
                               leonard_kl    = 4.0_r_def,         &
                               conservative_diffusion = .true.,   &
                               density_weighted = .true.,         &
-                              max_diff_factor=1.0_r_def )
+                              max_diff_factor  = 1.0_r_def,      &
+                              cap_blended_ml   = .false.,        &
+                              shallow_cu_maxtop= 0.0_r_def )
 
     ! Variables for testing momentum_smagorinsky
     ! Get sizes of everything

science/gungho/unit-test/kernel/diffusion/divergence_momentum_flux_kernel_mod_test.pf

@@ -115,7 +115,9 @@ contains
                               leonard_kl    = 4.0_r_def,         &
                               conservative_diffusion = .true.,   &
                               density_weighted = .true.,         &
-                              max_diff_factor=1.0_r_def )
+                              max_diff_factor  = 1.0_r_def,      &
+                              cap_blended_ml   = .false.,        &
+                              shallow_cu_maxtop= 0.0_r_def )
 
     ! Variables for testing momentum_smagorinsky
     ! Get sizes of everything

science/gungho/unit-test/kernel/diffusion/leonard_term_kl_kernel_mod_test.pf

@@ -55,8 +55,10 @@ contains
                               leonard_kl    = 4.0_r_def,         &
                               conservative_diffusion = .false.,  &
                               density_weighted = .false.,        &
-                              max_diff_factor = 1.0_r_def,       &
-                              fullstress    = .false. )
+                              max_diff_factor  = 1.0_r_def,      &
+                              fullstress       = .false.,        &
+                              cap_blended_ml   = .false.,        &
+                              shallow_cu_maxtop= 0.0_r_def )
 
     call feign_finite_element_config(      &
         cellshape=cellshape_quadrilateral, &

science/gungho/unit-test/kernel/diffusion/tracer_smagorinsky_diff_kernel_mod_test.pf

@@ -56,8 +56,10 @@ contains
                               leonard_kl    = 1.0_r_def,         &
                               conservative_diffusion = .false.,  &
                               density_weighted = .false.,        &
-                              max_diff_factor = 1.0_r_def,       &
-                              fullstress    = .false. )
+                              max_diff_factor  = 1.0_r_def,      &
+                              fullstress       = .false.,        &
+                              cap_blended_ml   = .false.,        &
+                              shallow_cu_maxtop= 0.0_r_def )
 
     call feign_finite_element_config(      &
         cellshape=cellshape_quadrilateral, &

science/physics_schemes/source/boundary_layer/bdy_expl2.F90

@@ -602,6 +602,9 @@ subroutine bdy_expl2 (                                                         &
            2:bl_levels),                                                       &
                                 ! Mixing length for heat (m),
                                 ! held on theta and rho levels, resp.
+   cumulus_cth(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end),           &
+                                ! cf_bulk cloud top height used in blended BL
+                                ! for shallow convection
    fm_3d(tdims%i_start:tdims%i_end,tdims%j_start:tdims%j_end,bl_levels),       &
                                 ! stability function for momentum transport
                                 ! level 1 value is dummy
@@ -717,9 +720,6 @@ subroutine bdy_expl2 (                                                         &
                                ! Flag for having reached
                                ! the top of the turbulently mixed
                                ! layer.
- l_shallow_cth(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end),           &
-                               ! Flag to indicate shallow convection based on
-                               ! cf_bulk cloud top height
  cloud_base_found(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end),        &
                                ! Flag for having reached cloud base
  cloud_top_found(pdims%i_start:pdims%i_end,pdims%j_start:pdims%j_end)
@@ -1829,20 +1829,21 @@ subroutine bdy_expl2 (                                                         &
 !$OMP SCHEDULE(STATIC)                                                         &
 !$OMP DEFAULT(none)                                                            &
 !$OMP private(i)                                                               &
-!$OMP SHARED(pdims,l_shallow_cth)
+!$OMP SHARED(pdims,cumulus_cth)
 do i = pdims%i_start, pdims%i_end
-  l_shallow_cth(i,j) = .false.
+  cumulus_cth(i,j) = zero ! initialise to zero for case of no cloud
+                          ! - will then blend using BL depth
 end do
 !$OMP end PARALLEL do
 
 if (blending_option == blend_cth_shcu_only) then
-  ! only going to use the parcel top as the length scale in blending
+  ! going to use cloud top height as the length scale in blending only
   ! if the convection is shallow, where we define shallow convection here
   ! as the resolved cloud top (from cf_bluk) being below shallow_cu_maxtop
 !$OMP PARALLEL                                                                 &
 !$OMP DEFAULT(none)                                                            &
 !$OMP private(i,k)                                                             &
-!$OMP SHARED(pdims,bl_levels,z_uv,l_shallow_cth,cumulus,ntml,cf_bulk,          &
+!$OMP SHARED(pdims,bl_levels,z_uv,cumulus_cth,cumulus,ntml,cf_bulk,            &
 !$OMP        shallow_cu_maxtop,cloud_base_found,cloud_top_found)
 !$OMP do SCHEDULE(STATIC)
   do i = pdims%i_start, pdims%i_end
@@ -1863,7 +1864,7 @@ subroutine bdy_expl2 (                                                         &
             cf_bulk(i,j,k+1) < sc_cftol ) then
         ! got to cloud-top
         cloud_top_found(i,j) = .true.
-        l_shallow_cth(i,j)   = z_uv(i,j,k+1) < shallow_cu_maxtop
+        cumulus_cth(i,j)     = z_uv(i,j,k+1)
       end if
     end do
 !$OMP end do
@@ -1878,7 +1879,7 @@ subroutine bdy_expl2 (                                                         &
    bl_levels,k_log_layr,BL_diag,                                               &
 ! in fields
    sigma_h,flandg,dvdzm,ri,rho_wet_tq,z_uv,z_tq,z0m_eff_gb,zhnl,zhpar,zhsc,    &
-   zdsc_base,ntpar,ntml_nl,ntdsc,nbdsc,l_shallow_cth,rmlmax2,rneutml_sq,       &
+   zdsc_base,ntpar,ntml_nl,ntdsc,nbdsc,cumulus_cth,rmlmax2,rneutml_sq,         &
    delta_smag,                                                                 &
 ! in/out fields
    cumulus,weight_1dbl,                                                        &
@@ -2475,7 +2476,7 @@ subroutine bdy_expl2 (                                                         &
       bl_type_6(i,j) = one
       ! Label this shallow regime as 2.0_r_bl here, to be able to identify it
       ! in diagnostics_bl, but the "cumulus" stash output will still be 1.0
-      if (l_shallow_cth(i,j)) bl_type_6(i,j) = 2.0_r_bl
+      if (cumulus_cth(i,j) < shallow_cu_maxtop) bl_type_6(i,j) = 2.0_r_bl
     end if
   end if
 end do