Allow for rebalancing of level 1 on wrfinput pathway

Docs/sphinx_doc/Initialization.rst

@@ -78,6 +78,26 @@ If the sounding is ``ConstantDensity``, then the initial density field is
 uniformly set to 1.0; the potential temperature (and water vapor mixing ratio)
 field(s) are set to the sounding values.
 
+For the ``Ideal``, ``Isentropic`` and ``DryIsentropic`` soundings, the profiles are
+interpolated onto the cell-centered heights of the mesh and each column is then
+rebalanced so that the base state is in discrete hydrostatic equilibrium on that mesh.
+With terrain-fitted coordinates this step matters, because the terrain-following heights
+differ from the nominal levels the 1-D sounding was integrated on.
+
+This is done on **every level**, against that level's own heights, so a refined level is
+hydrostatic on its own mesh rather than inheriting an interpolation of its parent's base
+state.  A refined level whose terrain is read from a text file via
+``erf.terrain_file_name`` re-reads that file at its own resolution, so it genuinely
+resolves topography its parent does not.  A refined region may also cover only part of the
+depth of the domain; see :ref:`subsec:base-state-refined` for why that does not change the
+base state in the cells it does contain, and :ref:`subsec:partial-depth-refinement` for how
+to specify such a region.
+
+The same restrictions apply as for a refined region in any terrain run: the PBL models
+(MYJ, MYNN2.5, MYNN-EDMF, YSU, MRF), the SHOC PBL model, and the column-integral derived
+quantities (``helicity``, ``precipitable``, ``max_reflectivity``, ``mucape``) all need
+entire columns, so they cannot be used on a level whose grids do not span the domain in z.
+
 .. note::
 
    You can optionally replace only the velocity fields (``u``, ``v``, ``w``)
@@ -119,6 +139,9 @@ file **Exec/ERF_Prob.cpp** must still be present for the build.
     We also note that wrfinput files that have later start times than the base level
     can be read and the new level will be initialized once the start time of the finer
     level has been reached.
+    A file may be supplied for each level -- ``erf.nc_init_file_0``,
+    ``erf.nc_init_file_1``, and so on -- but only level 0 is required to have one;
+    see :ref:`sec:nested-wrfinput` below.
 
 * **erf.init_type = Metgrid**
 
@@ -143,6 +166,119 @@ file **Exec/ERF_Prob.cpp** must still be present for the build.
     Optional HSE variables include ``RHO_HSE``, ``T_HSE``, and ``P_HSE``; the base state will be
     calculated if it is not specified.
 
+.. _sec:nested-wrfinput:
+
+Nested Initialization From WRF Input Files
+------------------------------------------
+
+With ``erf.init_type = WRFInput``, a refined level may be initialized either from its
+own wrfinput file or by interpolation from its parent, and it may cover either the full
+depth of the domain or only part of it.  Refinement should be horizontal only, so use
+``amr.ref_ratio_vect`` with 1 in the third slot rather than ``amr.ref_ratio``::
+
+          amr.ref_ratio_vect = 3 3 1
+
+When the finer level is initialized from its own wrfinput file this is enforced: ERF
+aborts if the refinement ratio in z is not 1, because the machinery that reads the file
+onto the finer level assumes the two levels share their vertical grid.
+
+Where the refined region goes
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+If a file is given for the finer level, ERF reads the nest's position and size from that
+file's ``I_PARENT_START``, ``J_PARENT_START`` and ``PARENT_GRID_RATIO`` attributes, and
+by default refines exactly the region the file covers, over the full depth of the domain.
+
+To refine only part of that region, add a refinement indicator.  The box it specifies
+must be contained in the region covered by the file; ERF checks this and aborts with both
+boxes printed if it is not.  The same mechanism is what lets the refined region stop below
+the domain top::
+
+          amr.max_level      = 1
+          amr.ref_ratio_vect = 3 3 1
+
+          erf.init_type      = WRFInput
+          erf.nc_init_file_0 = "wrfinput_d01"
+          erf.nc_init_file_1 = "wrfinput_d02"
+
+          erf.refinement_indicators = box1
+          erf.box1.max_level = 1
+          erf.box1.in_box_lo =  70000.   90000.     0.
+          erf.box1.in_box_hi = 110000.  150000.  6000.
+
+Note that all three components must be given.  Supplying only x and y is the documented
+way to ask for the full depth of the domain (see :ref:`subsec:full-depth-refinement`),
+which is the opposite of what is wanted here.
+
+If no file is given for the finer level, omit ``erf.nc_init_file_1`` and use a refinement
+indicator alone.  The finer level's terrain and state are then interpolated from its
+parent, and its base state is built as described below.  This works whether or not the
+refined region reaches the domain top.
+
+The base state at each level
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The WRF reference state is a closed-form function of height -- piecewise in log-pressure
+across a troposphere, an isothermal layer and an optional stratosphere -- built from the
+six parameters ``T00``, ``P00``, ``TLP``, ``TISO``, ``TLP_STRAT`` and ``P_STRAT``.  ERF
+evaluates that profile at each level's own cell-centered heights and then rebalances each
+column so that the result satisfies :math:`dp_0/dz = -\rho_0 g` discretely on that
+level's mesh.  Every level runs this same construction, which is what makes the base
+states of the levels consistent with one another.
+
+Two properties matter when the refined region does not reach the domain top:
+
+-  The rebalance integrates **upward from the bottom of each column**, so the value in a
+   cell depends only on the cells at or below it.  A refined region that stops partway up
+   therefore gets exactly the base state it would have had if it had been refined to the
+   model top -- the values in the cells it does contain are unchanged, to the last bit.
+
+-  The six parameters are read from the **level 0** file only, and every finer level uses
+   those same values.  The whole base state is built from those six numbers, so levels
+   that disagreed about them could not have base states that agree.  This matters in
+   practice because idealized and hand-built wrfinput files often declare these variables
+   but leave them zero-filled, in which case they fall back to the ERF defaults.  If a
+   nested file carries values that differ from level 0's, ERF prints a warning naming
+   both sets and uses level 0's.
+
+The base states of two levels that each have their own file will still not be identical
+to one another, because each is built on its own terrain and a nested file resolves finer
+topography than its parent.  That difference is set by the difference in surface elevation
+between the files and is carried up the column as an essentially constant hydrostatic
+offset; it does not grow with height.  It is the same difference a full-depth nest already
+has across its lateral faces.
+
+Restrictions
+~~~~~~~~~~~~
+
+-  The refinement ratio in z must be 1 for any level initialized from its own wrfinput
+   file.  Only horizontal refinement has been exercised for the interpolated case as
+   well, so ``amr.ref_ratio_vect`` with 1 in z is the supported configuration either way.
+
+-  The PBL models (MYJ, MYNN2.5, MYNN-EDMF, YSU, MRF), the SHOC PBL model, and the
+   column-integral derived quantities (``helicity``, ``precipitable``,
+   ``max_reflectivity``, ``mucape``) all need entire columns and cannot be used on a
+   level whose grids do not span the domain in z.  Set ``erf.pbl_type = "None"`` for
+   that level and leave those variables out of the plotfile variable list.
+
+-  For a level that has its own wrfinput file and does not reach the domain top, the
+   vertical remap of the **state** is incomplete in the cells near the top of the refined
+   region.  ERF remaps the state from the WRF vertical coordinate onto its own
+   terrain-following grid by searching the column of geopotential heights in the file;
+   both that column and the state being interpolated are stored on the level's own grids,
+   so above the top of the refined region there is nothing to interpolate from and those
+   cells keep the values read from the file.  ERF prints a warning naming the level when
+   this happens.  The base state is not affected.  A level whose state is interpolated
+   from its parent is not affected either.
+
+-  ``erf.box1.in_box_lo`` / ``in_box_hi`` in z are interpreted on the nominal (undeformed)
+   vertical grid, not on the terrain-following heights the wrfinput file supplies.  Check
+   the box ERF reports (``Saving in 'boxes at level'``) and the ``z_phys`` field in the
+   plotfile to confirm the refined region reaches the height you intended.
+
+See :ref:`subsec:partial-depth-refinement` for the mesh-refinement side of this, and
+:ref:`MeshRefinement` for refinement in general.
+
 TKE Initialization
 --------------------
 

Docs/sphinx_doc/Inputs.rst

@@ -2138,9 +2138,12 @@ List of Parameters
 |                                  | these bcs?          |                    | Metgrid;              |
 |                                  |                     |                    | else false            |
 +----------------------------------+---------------------+--------------------+-----------------------+
-| **erf.nc_init_file**             | NetCDF file with    |  String            | NONE                  |
-|                                  | initial mesoscale   |                    |                       |
-|                                  | data                |                    |                       |
+| **erf.nc_init_file_<lev>**       | NetCDF file(s) with |  String or         | NONE                  |
+|                                  | initial mesoscale   |  list of Strings   |                       |
+|                                  | data at level       |                    |                       |
+|                                  | <lev>. Required at  |                    |                       |
+|                                  | level 0; optional   |                    |                       |
+|                                  | at finer levels     |                    |                       |
 +----------------------------------+---------------------+--------------------+-----------------------+
 | **erf.nc_bdy_file**              | NetCDF file with    |  String            | NONE                  |
 |                                  | mesoscale data at   |                    |                       |
@@ -2268,7 +2271,14 @@ Notes
 
 If **erf.init_type = WRFInput** or **erf.init_type = NCFile**,
 the problem is initialized with mesoscale data contained in a NetCDF file,
-provided via ``erf.nc_init_file`` (e.g., "wrfinput_d01").
+provided via ``erf.nc_init_file_0`` (e.g., "wrfinput_d01").
+
+A file may also be given for each refined level -- ``erf.nc_init_file_1`` (e.g.,
+"wrfinput_d02"), ``erf.nc_init_file_2``, and so on.  A file at level 0 is required;
+files at finer levels are optional.  A refined level given no file of its own is
+initialized by interpolation from its parent, and its base state is built from the
+level-0 reference profile evaluated at that level's own heights.  See
+:ref:`sec:nested-wrfinput` for the nesting options and their restrictions.
 
 In addition, if **erf.use_real_bcs = true**, the lateral boundary conditions must be supplied in a NetCDF files
 specified by ``erf.nc_bdy_file`` (e.g., "wrfbdy_d01").  (If **erf.use_real_bcs = false**, no file is read for the

Docs/sphinx_doc/MeshRefinement.rst

@@ -220,6 +220,9 @@ aborts.  Using the two-value form above is the simplest way to guarantee this.
 Similarly, the SHOC PBL model requires that no box be split in the vertical
 direction; see :ref:`subsec:no-vertical-decomposition`.
 
+A refined region that deliberately covers only part of the depth is common in
+terrain runs; see :ref:`subsec:partial-depth-refinement`.
+
 The two-value form applies to statically specified boxes.  For dynamically
 created grids there is a separate, blunter mechanism based on the vertical
 tagging buffer; see :ref:`subsec:dynamic-full-depth`.
@@ -230,6 +233,59 @@ and that the specified indices are snapped outward to the nearest indices aligne
 with the refinement ratio (a message is printed whenever this snapping changes the
 box).
 
+.. _subsec:partial-depth-refinement:
+
+Refining Only Part of the Depth with Terrain
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+A nested LES region rarely needs to be refined all the way to the model top, so a
+refinement box whose vertical extent stops well below ``geometry.prob_hi`` in z is a
+natural thing to want.  This is supported with terrain-fitted coordinates for both of the
+initialization pathways that build a stratified base state:
+
+-  ``erf.init_type = input_sounding``, with the terrain read from a text file via
+   ``erf.terrain_file_name``; and
+
+-  ``erf.init_type = WRFInput``, whether each level has its own wrfinput file or only
+   level 0 does.
+
+For example, with a domain 17500 m deep::
+
+          amr.max_level      = 1
+          amr.ref_ratio_vect = 3 3 1
+
+          erf.refinement_indicators = box1
+          erf.box1.max_level = 1
+          erf.box1.in_box_lo =  70000.   90000.     0.
+          erf.box1.in_box_hi = 110000.  150000.  6000.
+
+Two things to watch, both specific to how the box is specified:
+
+-  All three components must be given.  The two-value form described above fills the
+   vertical extent with the full depth of the domain, which is the opposite of what is
+   wanted here.
+
+-  The z values are located on the nominal (undeformed) vertical grid -- that is, on
+   ``zlevels_stag`` -- and not on the terrain-following heights, so the box may not reach
+   the height you expect.  This matters most with ``erf.init_type = WRFInput``, where the
+   nominal grid is uniform while the heights the wrfinput file supplies are strongly
+   compressed near the surface.  Check the box ERF reports
+   (``Saving in 'boxes at level'``) against the ``z_phys`` field in the plotfile.
+
+With ``erf.init_type = WRFInput``, if a file is given for the finer level then the box
+must be contained in the region that file covers; ERF checks this and aborts with both
+boxes printed if it is not.
+
+The base state is built independently on each level, against that level's own heights,
+and its vertical integration is causal upward -- so a level that stops partway up gets
+exactly the base state it would have had if it had been refined to the model top.  See
+:ref:`subsec:base-state-refined` for that, and :ref:`sec:nested-wrfinput` for the
+wrfinput-specific details.
+
+The restrictions are the same for both pathways: the PBL models and the column-integral
+derived quantities need entire columns and cannot be used on a level whose grids do not
+span the domain in z.  See :ref:`sec:nested-wrfinput`.
+
 Moving Refinement Regions
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 

Docs/sphinx_doc/theory/BaseState.rst

@@ -148,3 +148,48 @@ where :math:`p_0 = 1e5 \, \text{N/m}^2` is the pressure at the base. Hence, we d
     F(p(0)) \equiv p(0) - p_0 + \rho(0)g\frac{\Delta z}{2},
 
 and the Newton-Raphson procedure is the same.
+
+.. _subsec:base-state-refined:
+
+The base state on a refined mesh
+----------------------------------------------------------------------------
+
+With mesh refinement, the base state is constructed **independently on each level**,
+against that level's own cell-centered heights :math:`z_{cc}`, by the same procedure
+described above.  Running the same construction on every level, rather than
+interpolating the coarse answer, is what leaves each level in discrete hydrostatic
+balance on its own mesh.
+
+Before that construction runs, the base state on a refined level is filled by
+conservative interpolation from its parent, and the physical boundary conditions for the
+base state are applied.  That is what gives values to the part of the refined level that
+lies inside the domain but outside the refined grids, and to the ghost cells; the
+per-level construction then overwrites only the refined grids.
+
+Two consequences are worth stating explicitly.
+
+**A refined level need not span the depth of the domain.**  The vertical integration is
+causal upward: the value in a cell depends only on the cells at or below it.  A refined
+region that stops partway up the domain therefore gets exactly the base state it would
+have had if it had been refined all the way to the model top -- the values in the cells it
+does contain are unchanged.  What the integration does require is a valid starting value
+in the lowest cell of each column, which is why the grids must not be decomposed in the
+vertical (see :ref:`subsec:no-vertical-decomposition`); ``amr.refine_grid_layout_z``
+defaults to 0 in ERF, so this holds unless it is overridden.
+
+**The base state does not depend on the coupling type.**  Under two-way coupling ERF
+averages :math:`\det J` down from fine to coarse, so that the volume weighting used when
+averaging the state down telescopes and remains conservative.  The cell-centered heights
+:math:`z_{cc}` are deliberately **not** averaged down.  Were they, the coarse heights
+would change after the coarse base state had already been built against them, leaving that
+base state out of hydrostatic balance with the heights the dynamics subsequently use for
+vertical gradients, Rayleigh damping and the sponge zones.  Not averaging them also keeps
+:math:`z_{cc}` consistent with the nodal heights :math:`z_{nd}` it is derived from, which
+are not averaged down either.
+
+Note that the base states of two levels will not be identical to one another wherever
+their terrain differs -- for instance when a finer level resolves topography that its
+parent does not, either because it re-reads a terrain file at its own resolution or
+because it reads a nested wrfinput file.  Each level is individually hydrostatic; the
+difference between them is set by the difference in surface elevation and is carried up
+the column as an essentially constant offset in pressure.