WIP: YSU PBL (#1540)

* initial equations and references for YSU PBL

* minor YSU doc changes

* warn for PBL + moisture

* start to YSU

* ignore broken link that actually works

* error on warnings in linkcheck but ignore DOI redirects

* add level to turb visc interfaces to pass to MOST

* rename MYNN pbl params to avoid confusion with ysu pbl params

* add z0 and zref functions to MOST class

* progress on stable PBL height calc

* pbl progress: stable PBL Rib_cr calcs

* stable PBL height calculation done

* spelling fix and other minor cleanup

.github/workflows/docs.yml

@@ -27,7 +27,7 @@ jobs:
       - name: Check links
         run: |
           cd Docs/sphinx_doc
-          make SPHINXOPTS=-v NO_DOXYGEN=TRUE linkcheck
+          make SPHINXOPTS="-v -W --keep-going" NO_DOXYGEN=TRUE linkcheck
 
       - name: Deploy to Webpage
         if: github.event_name == 'push' && github.repository == 'erf-model/ERF' && github.ref == 'refs/heads/development'

Docs/sphinx_doc/Inputs.rst

@@ -825,23 +825,23 @@ List of Parameters
 | **erf.pbl_type**                 | Name of PBL Scheme | "None", "MYNN2.5"   | "None"      |
 |                                  | to be used         |                     |             |
 +----------------------------------+--------------------+---------------------+-------------+
-| **erf.pbl_A1**                   | MYNN Constant A1   | Real                | 1.18        |
+| **erf.pbl_mynn_A1**              | MYNN Constant A1   | Real                | 1.18        |
 +----------------------------------+--------------------+---------------------+-------------+
-| **erf.pbl_A2**                   | MYNN Constant A2   | Real                | 0.665       |
+| **erf.pbl_mynn_A2**              | MYNN Constant A2   | Real                | 0.665       |
 +----------------------------------+--------------------+---------------------+-------------+
-| **erf.pbl_B1**                   | MYNN Constant B1   | Real                | 24.0        |
+| **erf.pbl_mynn_B1**              | MYNN Constant B1   | Real                | 24.0        |
 +----------------------------------+--------------------+---------------------+-------------+
-| **erf.pbl_B2**                   | MYNN Constant B2   | Real                | 15.0        |
+| **erf.pbl_mynn_B2**              | MYNN Constant B2   | Real                | 15.0        |
 +----------------------------------+--------------------+---------------------+-------------+
-| **erf.pbl_C1**                   | MYNN Constant C1   | Real                | 0.137       |
+| **erf.pbl_mynn_C1**              | MYNN Constant C1   | Real                | 0.137       |
 +----------------------------------+--------------------+---------------------+-------------+
-| **erf.pbl_C2**                   | MYNN Constant C1   | Real                | 0.75        |
+| **erf.pbl_mynn_C2**              | MYNN Constant C1   | Real                | 0.75        |
 +----------------------------------+--------------------+---------------------+-------------+
-| **erf.pbl_C3**                   | MYNN Constant C3   | Real                | 0.352       |
+| **erf.pbl_mynn_C3**              | MYNN Constant C3   | Real                | 0.352       |
 +----------------------------------+--------------------+---------------------+-------------+
-| **erf.pbl_C4**                   | MYNN Constant C4   | Real                | 0.0         |
+| **erf.pbl_mynn_C4**              | MYNN Constant C4   | Real                | 0.0         |
 +----------------------------------+--------------------+---------------------+-------------+
-| **erf.pbl_C5**                   | MYNN Constant C5   | Real                | 0.2         |
+| **erf.pbl_mynn_C5**              | MYNN Constant C5   | Real                | 0.2         |
 +----------------------------------+--------------------+---------------------+-------------+
 | **erf.advect_QKE**               | Include advection  | bool                | 1           |
 |                                  | terms in QKE eqn   |                     |             |

Docs/sphinx_doc/conf.py

@@ -163,7 +163,7 @@
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
+# html_static_path = ['_static']
 
 # Add any extra paths that contain custom files (such as robots.txt or
 # .htaccess) here, relative to this directory. These files are copied
@@ -328,9 +328,15 @@
     'https://link.springer.com/article/10.1007/BF00240838',
     'https://onlinelibrary.wiley.com/doi/10.1029/2021MS002904',
     'https://link.springer.com/article/10.1023/B:BOUN.0000020164.04146.98',
-    'https://doi.org/10.1029/RG020i004p00851'
+    'https://doi.org/10.1029/RG020i004p00851',
+    'https://doi.org/10.1002/qj.665'
 ]
 
+linkcheck_allowed_redirects = {
+    r'https://doi\.org/*' : '.*',
+    r'https://dx\.doi\.org/*' : '.*',
+    r'http://dx\.doi\.org/*' : '.*'
+}
 
 # -- Options for Texinfo output -------------------------------------------
 

Docs/sphinx_doc/theory/PBLschemes.rst

@@ -128,4 +128,110 @@ More recent advancements that add significant complexity to the MYNN scheme have
 YSU PBL Model
 -------------
 
-The Yonsei University (YSU) PBL model is another commonly use scheme in WRF. It is not yet supported in ERF, but is under development.
+.. warning::
+
+   Implementation is in progress, this option is not yet supported
+
+The Yonsei University (YSU) PBL model is another commonly use scheme in WRF. It includes nonlocal mixing with  contergradient diffusion within
+the PBL, and a local mixing treatment outside the PBL.
+
+Turbulent diffusion for prognostic variables (:math:`C, u, v, \theta, q_k`), where :math:`q_k` includes all moisture variables and :math:`C`
+any additional scalars (other terms in the equations omitted for brevity):
+
+.. math::
+   \frac{\partial C}{\partial t}
+   = \frac{\partial}{\partial z} \left[
+   K_c \left( \frac{\partial C}{\partial z} - \gamma_c \right)
+   - \overline{\left(w'c' \right)_h} \left( \frac{z}{h} \right)^3
+   \right]
+
+.. note::
+
+   Not applied for vertical velocity?
+
+Where for each variable the turbulent diffusion coefficient :math:`K_c`, countergradient correction :math:`\gamma_c`,
+and entrainment flux at the PBL top :math:`\overline{\left(w'c' \right)_h}` must be diagnosed for each variable.
+The main controlling parameter is the PBL height :math:`h`.
+Notably, a nonlocal model for turbulent diffusion is used for :math:`z \leq h`, but a local model is used for :math:`z \ge h`.
+
+The first step is to determine the PBL height :math:`h`. This is defined as the smallest value of :math:`z` where the bulk
+Richardson number equals the critical value, which is taken to be 0:
+
+.. math::
+
+   {\rm Rib}(z) = \frac{ g \left[ \theta_m(z) - \theta_s\right] }{\theta_{ma} U(z)^2}z
+
+.. math::
+
+   {\rm Rib}(h) = {\rm Rib_{cr}} = 0
+
+where
+
+- :math:`\theta_m` is the moist potential temperature,
+- :math:`\theta_{ma}` is the value at the lowest vertical cell in a column,
+- :math:`U = \sqrt{u^2 + v^2}` is the horizontal wind speed,
+- :math:`\theta_s = \theta_{ma} + \theta_T` is the virtual temperature near the surface,
+- :math:`\theta_T = a\frac{\overline{\left(w'\theta_m' \right)_0}}{w_{s0}}` is the excess virtual temperature near the surface,
+- :math:`a` is a constant taken to be 6.8 per HND06 (matching the :math:`b` constant that appears elsehwere in the YSU model)
+- :math:`\overline{\left(w'\theta_m' \right)_0}` is the surface virtual heat flux (determined from the MOST surface layer model),
+- :math:`w_{s}(z) = \left(u_*^3 + 8 k w_{*b}^3z/h \right)^{1/3}` is a representative velocity scale in the mixed layer, with :math:`w_{s0} = w_s(h/2)` (note this equation matches the WRF implementation and description in H10, but differs from HND06, where :math:`\phi_m` appears in place of the constant 8),
+- :math:`u_*` is the surface frictional velocity scale determined from the MOST surface layer model,
+- :math:`k = 0.4` is the von Karman constant
+- :math:`w_{*b} = \left[ g/\theta_{ma} \overline{\left(w'\theta_m' \right)_0} h \right]^{1/3}` for :math:`\overline{\left(w'\theta_m' \right)_0} > 0`, :math:`w_{*b} = 0` otherwise, is a convective velcoity scale for moist air
+
+In practice, an approximate value of :math:`h` is determined through a two-step process. First, :math:`\theta_T` is set to be zero
+and a provisional value of :math:`h` is estimated. Then this provisional value of :math:`h` is used to compute :math:`\theta_T`,
+which is in turn used to provide an improved estimate of :math:`h`, which is the value used in subsequent calculations.
+
+.. note::
+
+   This two step-process matches the WRF implementation, but this could be extended iteratively to reach convergence.
+
+
+Countergradient corrections are computed as follows:
+
+.. math::
+
+   \gamma_\theta =
+
+.. math::
+   \gamma_u =
+
+.. math::
+   \gamma_v =
+
+.. math::
+   \gamma_{q_k} = \gamma_C = 0
+
+Entrainment fluxes are computed:
+
+.. math::
+   \overline{\left(w'c' \right)_h} =
+
+.. math::
+   \overline{\left(w'c' \right)_h} =
+
+Within the PBL (:math:`z \leq h`),
+
+.. _YSUReferences:
+
+Useful References
+~~~~~~~~~~~~~~~~~
+
+The following references have informed the implementation of the YSU model in ERF:
+
+.. _HP96: https://doi.org/10.1175/1520-0493(1996)124<2322:NBLVDI>2.0.CO;2
+
+- [H10] `Hong, Quarterly Journal of the Royal Meteorological Society, 2010 <https://doi.org/10.1002/qj.665>`_: Most up-to-date YSU model formulation as implemented in WRF, with revisions for stable boundary layers
+
+- [HND06] `Hong, Noh, and Dudhia, Monthly Weather Review, 2006 <https://doi.org/10.1175/MWR3199.1>`_: Initial formulation referred to as the YSU model, adds improved entrainment formulation (relative to NCHR03) to work of TM86 and a few other modifications
+
+- [NCHR03] `Noh, Cheon, Hong, and Raasch, Boundary-Layer Meteorology, 2003 <https://doi.org/10.1023/A:1022146015946>`_: Entrainment effects added to TM86
+
+- [HP96] `Hong and Pan, Monthly Weather Review, 1996 <HP96_>`_: Largely an implementation and evluation of TM86
+
+- [TM86] `Troen and Mahrt, Boundary-Layer Meteorology, 1986 <https://doi.org/10.1007/BF00122760>`_: Initial incorporation of nonlocal counter-graident term in vertical diffusion model
+
+- [WF18] `Wilson and Fovell, Weather and Forecasting, 2018 <https://doi.org/10.1175/WAF-D-17-0109.1>`_: Extension of YSU to handle interplay between radiation and fog, active in WRF with the ``ysu_topdown_pblmix = 1`` option
+
+- The WRF Fortran source code for this `module <https://github.com/wrf-model/WRF/blob/a8eb846859cb39d0acfd1d3297ea9992ce66424a/phys/module_bl_ysu.F>`_ as of Dec. 2023. The ERF implementation supports the same physical models as this WRF implementation, with the exception of the ``ysu_topdown_pblmix = 1`` option from WF18, i.e. the implementation in ERF largely matches the PBL scheme described in H10.

Source/BoundaryConditions/ABLMost.H

@@ -286,6 +286,15 @@ public:
     const amrex::MultiFab*
     get_mac_avg (const int& lev, int comp) { return m_ma.get_average(lev,comp); }
 
+    const amrex::MultiFab*
+    get_t_surf (const int& lev) { return t_surf[lev].get(); }
+
+    amrex::Real
+    get_zref () {return m_ma.get_zref();}
+
+    const amrex::FArrayBox*
+    get_z0 (const int& lev) {return &z_0[lev];}
+
     enum struct FluxCalcType {
         MOENG = 0,      ///< Moeng functional form
         DONELAN,        ///< Donelan functional form

Source/DataStructs/DataStruct.H

@@ -243,6 +243,13 @@ struct SolverChoice {
             turbChoice[lev].init_params(lev,max_level);
         }
 
+        // Warn for PBL models and moisture - these may not yet be compatible
+        for (int lev = 0; lev <= max_level; lev++) {
+            if ((moisture_type != MoistureType::None) && (turbChoice[lev].pbl_type != PBLType::None)) {
+                amrex::Warning("\n*** WARNING: Moisture may not yet be compatible with PBL models, \n    proceed with caution ***");
+            }
+        }
+
         // Which type of refinement
         static std::string coupling_type_string = "OneWay";
         pp.query("coupling_type",coupling_type_string);

Source/DataStructs/TurbStruct.H

@@ -78,15 +78,15 @@ struct TurbChoice {
                     amrex::Print() << "Selected a PBL model and an LES model: " <<
                     "Using PBL for vertical transport, LES for horizontal" << std::endl;
                 }
-                pp.query("pbl_A1", pbl_A1);
-                pp.query("pbl_A2", pbl_A2);
-                pp.query("pbl_B1", pbl_B1);
-                pp.query("pbl_B2", pbl_B2);
-                pp.query("pbl_C1", pbl_C1);
-                pp.query("pbl_C2", pbl_C2);
-                pp.query("pbl_C3", pbl_C3);
-                pp.query("pbl_C4", pbl_C4);
-                pp.query("pbl_C5", pbl_C5);
+                pp.query("pbl_mynn_A1", pbl_mynn_A1);
+                pp.query("pbl_mynn_A2", pbl_mynn_A2);
+                pp.query("pbl_mynn_B1", pbl_mynn_B1);
+                pp.query("pbl_mynn_B2", pbl_mynn_B2);
+                pp.query("pbl_mynn_C1", pbl_mynn_C1);
+                pp.query("pbl_mynn_C2", pbl_mynn_C2);
+                pp.query("pbl_mynn_C3", pbl_mynn_C3);
+                pp.query("pbl_mynn_C4", pbl_mynn_C4);
+                pp.query("pbl_mynn_C5", pbl_mynn_C5);
             }
 
             // Right now, solving the QKE equation is only supported when MYNN PBL is turned on
@@ -163,15 +163,15 @@ struct TurbChoice {
                 } else if (les_type == LESType::Deardorff) {
                     amrex::Error("It is not recommended to use Deardorff LES and a PBL model");
                 }
-                pp.query("pbl_A1", pbl_A1, lev);
-                pp.query("pbl_A2", pbl_A2, lev);
-                pp.query("pbl_B1", pbl_B1, lev);
-                pp.query("pbl_B2", pbl_B2, lev);
-                pp.query("pbl_C1", pbl_C1, lev);
-                pp.query("pbl_C2", pbl_C2, lev);
-                pp.query("pbl_C3", pbl_C3, lev);
-                pp.query("pbl_C4", pbl_C4, lev);
-                pp.query("pbl_C5", pbl_C5, lev);
+                pp.query("pbl_mynn_A1", pbl_mynn_A1, lev);
+                pp.query("pbl_mynn_A2", pbl_mynn_A2, lev);
+                pp.query("pbl_mynn_B1", pbl_mynn_B1, lev);
+                pp.query("pbl_mynn_B2", pbl_mynn_B2, lev);
+                pp.query("pbl_mynn_C1", pbl_mynn_C1, lev);
+                pp.query("pbl_mynn_C2", pbl_mynn_C2, lev);
+                pp.query("pbl_mynn_C3", pbl_mynn_C3, lev);
+                pp.query("pbl_mynn_C4", pbl_mynn_C4, lev);
+                pp.query("pbl_mynn_C5", pbl_mynn_C5, lev);
             }
 
             // Right now, solving the QKE equation is only supported when MYNN PBL is turned on
@@ -237,15 +237,15 @@ struct TurbChoice {
         }
 
         if (pbl_type != PBLType::None) {
-            amrex::Print() << "pbl_A1                      : " << pbl_A1 << std::endl;
-            amrex::Print() << "pbl_A2                      : " << pbl_A2 << std::endl;
-            amrex::Print() << "pbl_B1                      : " << pbl_B1 << std::endl;
-            amrex::Print() << "pbl_B2                      : " << pbl_B2 << std::endl;
-            amrex::Print() << "pbl_C1                      : " << pbl_C1 << std::endl;
-            amrex::Print() << "pbl_C2                      : " << pbl_C2 << std::endl;
-            amrex::Print() << "pbl_C3                      : " << pbl_C3 << std::endl;
-            amrex::Print() << "pbl_C4                      : " << pbl_C4 << std::endl;
-            amrex::Print() << "pbl_C5                      : " << pbl_C5 << std::endl;
+            amrex::Print() << "pbl_mynn_A1                      : " << pbl_mynn_A1 << std::endl;
+            amrex::Print() << "pbl_mynn_A2                      : " << pbl_mynn_A2 << std::endl;
+            amrex::Print() << "pbl_mynn_B1                      : " << pbl_mynn_B1 << std::endl;
+            amrex::Print() << "pbl_mynn_B2                      : " << pbl_mynn_B2 << std::endl;
+            amrex::Print() << "pbl_mynn_C1                      : " << pbl_mynn_C1 << std::endl;
+            amrex::Print() << "pbl_mynn_C2                      : " << pbl_mynn_C2 << std::endl;
+            amrex::Print() << "pbl_mynn_C3                      : " << pbl_mynn_C3 << std::endl;
+            amrex::Print() << "pbl_mynn_C4                      : " << pbl_mynn_C4 << std::endl;
+            amrex::Print() << "pbl_mynn_C5                      : " << pbl_mynn_C5 << std::endl;
         }
     }
 
@@ -277,16 +277,23 @@ struct TurbChoice {
 
     // PBL model
     PBLType pbl_type;
-    // Model coefficients
-    amrex::Real pbl_A1 = 1.18;
-    amrex::Real pbl_A2 = 0.665;
-    amrex::Real pbl_B1 = 24.0;
-    amrex::Real pbl_B2 = 15.0;
-    amrex::Real pbl_C1 = 0.137;
-    amrex::Real pbl_C2 = 0.75;
-    amrex::Real pbl_C3 = 0.352;
-    amrex::Real pbl_C4 = 0.0;
-    amrex::Real pbl_C5 = 0.2;
+    // Model coefficients - MYNN2.5
+    amrex::Real pbl_mynn_A1 = 1.18;
+    amrex::Real pbl_mynn_A2 = 0.665;
+    amrex::Real pbl_mynn_B1 = 24.0;
+    amrex::Real pbl_mynn_B2 = 15.0;
+    amrex::Real pbl_mynn_C1 = 0.137;
+    amrex::Real pbl_mynn_C2 = 0.75;
+    amrex::Real pbl_mynn_C3 = 0.352;
+    amrex::Real pbl_mynn_C4 = 0.0;
+    amrex::Real pbl_mynn_C5 = 0.2;
+
+    // Model coefficients - YSU
+    // TODO: Add parmparse for all of these above
+    amrex::Real pbl_ysu_coriolis_freq = 1.0e-4; // TODO: make this consistent with coriolis forcing (but note WRF always uses 1e-4)
+    bool        pbl_ysu_over_land = false; // TODO: pull from other inputs and make local
+    amrex::Real pbl_ysu_land_Ribcr = 0.25; // Critical Bulk Richardson number of Land for stable conditions
+    amrex::Real pbl_ysu_unst_Ribcr = 0.0;
 
     // QKE stuff - default is not to use it, if MYNN2.5 or YSU PBL is used default is turb transport in Z-direction only
     bool use_QKE = false;

Source/Diffusion/ComputeTurbulentViscosity.cpp

@@ -16,6 +16,7 @@ ComputeTurbulentViscosityPBL (const MultiFab& xvel,
                               const Geometry& geom,
                               const TurbChoice& turbChoice,
                               std::unique_ptr<ABLMost>& most,
+                              int level,
                               const BCRec* bc_ptr,
                               bool /*vert_only*/,
                               const std::unique_ptr<MultiFab>& z_phys_nd);
@@ -412,6 +413,7 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel ,
                                 const std::unique_ptr<MultiFab>& z_phys_nd,
                                 const TurbChoice& turbChoice, const Real const_grav,
                                 std::unique_ptr<ABLMost>& most,
+                                int level,
                                 const BCRec* bc_ptr,
                                 bool vert_only)
 {
@@ -449,6 +451,6 @@ void ComputeTurbulentViscosity (const MultiFab& xvel , const MultiFab& yvel ,
     if (turbChoice.pbl_type != PBLType::None) {
         // NOTE: state_new is passed in for Cons_old (due to ptr swap in advance)
         ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, eddyViscosity,
-                                     geom, turbChoice, most, bc_ptr, vert_only, z_phys_nd);
+                                     geom, turbChoice, most, level, bc_ptr, vert_only, z_phys_nd);
     }
 }

Source/Diffusion/DiffusionSrcForState_N.cpp

@@ -459,7 +459,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
     // Using PBL
     if (l_use_QKE && start_comp <= RhoQKE_comp && end_comp >=RhoQKE_comp) {
         int qty_index = RhoQKE_comp;
-        auto pbl_B1_l = turbChoice.pbl_B1;
+        auto pbl_mynn_B1_l = turbChoice.pbl_mynn_B1;
         bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
         bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) );
         bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
@@ -470,7 +470,7 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain,
         {
             cell_rhs(i, j, k, qty_index) += ComputeQKESourceTerms(i,j,k,u,v,cell_data,cell_prim,
                                                                   mu_turb,cellSizeInv,domain,
-                                                                  pbl_B1_l,tm_arr(i,j,0),
+                                                                  pbl_mynn_B1_l,tm_arr(i,j,0),
                                                                   c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                                                   u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                                                   v_ext_dir_on_zlo, v_ext_dir_on_zhi);