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

@@ -55,6 +55,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

@@ -776,6 +776,20 @@
         "plot_str": "plot_map.py -ral $NODAL_DATA_DIR/lfric_diagnostics.nc  $PLOT_DIR",
     }) %}
 
+{% 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/common/default_directives.cylc

@@ -10,22 +10,22 @@
         "memory": [6, "GB"],
         "wallclock": 20,
         "build_memory": [1, "GB"],
-        "build_wallclock": 25,
+        "build_wallclock": 45,
         "build_cpus": 6,
         "gen_weights_tasks": 4,
         "gen_weights_memory": [4, "GB"],
 
         "lfric_atm": {
-            "build_wallclock": 45,
+            "build_wallclock": 65,
         },
         "linear_model": {
-            "build_wallclock": 45,
+            "build_wallclock": 65,
         },
         "jules": {
-            "build_wallclock": 35,
+            "build_wallclock": 60,
         },
         "ngarch": {
-            "build_wallclock": 45,
+            "build_wallclock": 65,
         },
         "mesh": {
             "mesh_wallclock": 20,

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 =         48399774
+Inner product checksum theta =         5135493D
+Inner product checksum u =         612F4555
+Inner product checksum mr1 =         4090F4C2
+Inner product checksum mr2 =         35AE1FCE
+Inner product checksum mr3 =         2FA57BD6
+Inner product checksum mr4 =         33F4FD79
+Inner product checksum mr5 =          BEDF8CA
+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 =         483997BC
+Inner product checksum theta =         51354960
+Inner product checksum u =         612F492E
+Inner product checksum mr1 =         4090F9FA
+Inner product checksum mr2 =         35B4A04B
+Inner product checksum mr3 =         2FA5A582
+Inner product checksum mr4 =         33F4FE3E
+Inner product checksum mr5 =          BECFCFF
+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({

rose-stem/templates/runtime/generate_runtime_scripts.cylc

@@ -111,7 +111,7 @@
     {% endif %}
 
         inherit={{local_build_inherit|upper}}
-        execution time limit = PT25M
+        execution time limit = PT40M
         script = """
                  cp -r $CYLC_WORKFLOW_SHARE_DIR/source/lfric_apps $CYLC_TASK_WORK_DIR/apps_local_build_source
                  cp -r $CYLC_WORKFLOW_SHARE_DIR/source/lfric_core $CYLC_TASK_WORK_DIR/core_local_build_source

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

@@ -3718,12 +3718,37 @@ help=Options for Smagorinsky subgrid mixing scheme:
     =_________________________________________________________________
     =Blended scheme + 1D BL above
     = As above, but relaxes to 1D BL above BL
+    =_________________________________________________________________
+    =Blended scheme + cloud-top of for shallow cu
+    = As above, but use cloud-top as blending length-scale only for
+    = shallow cu
     =
-    =  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' ;
+value-titles=3D Smagorinsky, 2D Smagorinsky + 1D BL, Blended scheme + Smag above,
+            =Blended scheme + 1D BL above, Blended scheme with CTH for shallow cu only
+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 by the 1D
+help=If true this caps the magnitude of the blended 1D-3D mixing length
+    =by the 1D 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 cu blending
+help=The maximum permitted cloud top height (in metres) for shallow cumulus
+    =to be diagosed, as then used for method=blending_option=blend_cth_shcu_only
+    =to use the cumulus parcel top height as the length scale in the blending weight
+!kind=default
+sort-key=Panel-A05b
+type=real
 
 [namelist:mixing=mix_factor]
 compulsory=true

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

@@ -31,3 +31,15 @@ def upgrade(self, config, meta_config=None):
         # Add settings
         return config, self.reports
 """
+class vn32_t752(MacroUpgrade):
+    # Upgrade macro for #752 by Adrian Lock
+
+    BEFORE_TAG = "vn3.2"
+    AFTER_TAG = "vn3.2_t752"
+
+    def upgrade(self, config, meta_config=None):
+        # Add settings
+        self.add_setting(
+            config, ["namelist:mixing", "cap_blended_ml"], ".false."
+        )
+        return config, self.reports

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

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

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

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

science/physics_schemes/source/boundary_layer/bdy_expl2.F90

@@ -605,6 +605,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
@@ -720,9 +723,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)
@@ -1935,10 +1935,11 @@ subroutine bdy_expl2 (                                                         &
 !$OMP SCHEDULE(STATIC)                                                         &
 !$OMP DEFAULT(none)                                                            &
 !$OMP private(j,i)                                                             &
-!$OMP SHARED(pdims,l_shallow_cth)
+!$OMP SHARED(pdims,cumulus_cth)
 do j = pdims%j_start, pdims%j_end
   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
 end do
 !$OMP end PARALLEL do
@@ -1950,8 +1951,8 @@ subroutine bdy_expl2 (                                                         &
 !$OMP PARALLEL                                                                 &
 !$OMP DEFAULT(none)                                                            &
 !$OMP private(i,j,k)                                                           &
-!$OMP SHARED(pdims,bl_levels,z_uv,l_shallow_cth,cumulus,ntml,cf_bulk,          &
-!$OMP        shallow_cu_maxtop,cloud_base_found,cloud_top_found)
+!$OMP SHARED(pdims,bl_levels,z_uv,cumulus_cth,cumulus,ntml,cf_bulk,            &
+!$OMP        cloud_base_found,cloud_top_found)
 !$OMP do SCHEDULE(STATIC)
   do j = pdims%j_start, pdims%j_end
     do i = pdims%i_start, pdims%i_end
@@ -1974,7 +1975,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
     end do
@@ -1990,7 +1991,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,                                                        &
@@ -2636,7 +2637,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