From 7d0e742753910ac807e1936748f2e21f63caef52 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Thu, 27 Aug 2026 17:53:29 -0700 Subject: [PATCH 1/2] Allow for rebalancing of level 1 on wrfinput pathway --- Docs/sphinx_doc/Initialization.rst | 116 ++++++ Docs/sphinx_doc/Inputs.rst | 18 +- Docs/sphinx_doc/MeshRefinement.rst | 41 +++ Source/DataStructs/ERF_BaseStateParams.H | 115 ++++++ Source/DataStructs/Make.package | 1 + Source/ERF.H | 21 ++ Source/ERF.cpp | 46 ++- Source/ERF_MakeNewArrays.cpp | 37 +- Source/ERF_MakeNewLevel.cpp | 62 +++- Source/Initialization/ERF_Init1D.cpp | 63 ++-- Source/Initialization/ERF_InitFromMetgrid.cpp | 6 +- .../Initialization/ERF_InitFromWRFInput.cpp | 331 ++++++++++++------ Source/Initialization/ERF_MetgridUtils.H | 79 +---- Source/Utils/ERF_Rebalance.cpp | 28 +- 14 files changed, 741 insertions(+), 223 deletions(-) create mode 100644 Source/DataStructs/ERF_BaseStateParams.H diff --git a/Docs/sphinx_doc/Initialization.rst b/Docs/sphinx_doc/Initialization.rst index c492a1754a..64f10362a1 100644 --- a/Docs/sphinx_doc/Initialization.rst +++ b/Docs/sphinx_doc/Initialization.rst @@ -119,6 +119,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 +146,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-wrfinput` for the mesh-refinement side of this, and +:ref:`MeshRefinement` for refinement in general. + TKE Initialization -------------------- diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 6c74d63e79..bd672745a6 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/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_** | NetCDF file(s) with | String or | NONE | +| | initial mesoscale | list of Strings | | +| | data at level | | | +| | . 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 diff --git a/Docs/sphinx_doc/MeshRefinement.rst b/Docs/sphinx_doc/MeshRefinement.rst index bb34e8236c..4702be9172 100644 --- a/Docs/sphinx_doc/MeshRefinement.rst +++ b/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 when +nesting from WRF input files; see :ref:`subsec:partial-depth-wrfinput`. + 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,44 @@ 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-wrfinput: + +Nesting from WRF Input Files over Only Part of the Depth +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +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 with ``erf.init_type = WRFInput``. This is supported, both when +each level has its own wrfinput file and when 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 interpreted on the nominal (undeformed) vertical grid, not on the + terrain-following heights that a wrfinput file supplies, so the box may not reach the + height you expect. Check the box ERF reports (``Saving in 'boxes at level'``) against + the ``z_phys`` field in the plotfile. + +If a file is given for the finer level, the box must be contained in the region that file +covers; ERF checks this and aborts with both boxes printed if it is not. + +See :ref:`sec:nested-wrfinput` for how the base state is built at each level, why +truncating the refined region does not change the values in the cells it does contain, +and the restrictions that apply (in particular, the PBL models and the column-integral +derived quantities cannot be used on a level whose grids do not span the domain in z). + Moving Refinement Regions ~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/Source/DataStructs/ERF_BaseStateParams.H b/Source/DataStructs/ERF_BaseStateParams.H new file mode 100644 index 0000000000..bf112babd1 --- /dev/null +++ b/Source/DataStructs/ERF_BaseStateParams.H @@ -0,0 +1,115 @@ +/** + * \file ERF_BaseStateParams.H + */ +#ifndef ERF_BASE_STATE_PARAMS_H_ +#define ERF_BASE_STATE_PARAMS_H_ + +#include + +#include +#include + +#include + +/** + * Parameters defining the WRF reference (base) state, along with the layer + * interface heights derived from them. + * + * The base state is piecewise in log-pressure: + * (1) troposphere T = T00 + TLP * ln(p/P00) p > P_iso + * (2) isothermal layer T = TISO P_STRAT < p <= P_iso + * (3) stratosphere T = TISO + TLP_STRAT * ln(p/P_STRAT) p <= P_STRAT + * Each piece inverts to p(z) in closed form under dp/dz = -rho g, so the two + * interface heights are precomputed once and the inversion branches on z. + * + * Because the resulting profile is a pure function of height, it can be + * evaluated on any subset of a column -- it does not require a box that spans + * the domain in the vertical. + * + * NOTE: The defaults are the WRF namelist defaults for base_temp, base_pres, + * base_lapse, iso_temp, base_lapse_strat and base_pres_strat. metgrid + * files do not carry these parameters at all; wrfinput files usually do, + * and are read by read_base_state_params_from_wrfinput. + */ +struct BaseStateParams +{ + amrex::Real T00 = amrex::Real(290.0); + amrex::Real P00 = p_0; + amrex::Real TLP = amrex::Real(50.0); + amrex::Real TISO = amrex::Real(200.0); + amrex::Real TLP_STRAT = amrex::Real(-11.0); + amrex::Real P_STRAT = zero; + + // Derived by set_layer_interfaces() + amrex::Real P_iso = zero; + amrex::Real z_iso = zero; + amrex::Real z_strat = zero; + bool use_strat = false; + + // Whether set_layer_interfaces() has been called, so that consumers can + // catch an un-initialized profile rather than silently using z_iso = 0. + bool is_set = false; + + // + // Compute the layer interfaces and report them. + // + // NOTE: This must be called from outside of any MFIter (and outside of any + // OpenMP parallel region) so the diagnostics are emitted exactly once + // rather than once per box per thread. + // + void set_layer_interfaces () + { + // The troposphere inverts p(z) as sqrt((T00/TLP)^2 - 2*g*z/(TLP*R_d)), which + // folds 1/TLP inside the square root and so discards the sign of TLP. A + // non-positive lapse parameter, or an isothermal temperature at or above the + // surface temperature, would silently select the wrong root instead of failing. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(TLP > zero, + "Base state TLP must be positive"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(TISO < T00, + "Base state TISO must be less than T00"); + + const amrex::Real x_iso = (TISO - T00) / TLP; + P_iso = P00 * std::exp(x_iso); + z_iso = -(R_d/CONST_GRAV) * (T00*x_iso + myhalf*TLP*x_iso*x_iso); + + // The upper stratospheric layer is optional (P_STRAT == 0 or TLP_STRAT == 0 + // disables it) and is only meaningful if it begins above the isothermal layer, + // i.e. if P_STRAT is below the pressure at which the isothermal layer starts. + const bool want_strat = ((P_STRAT > zero) && (TLP_STRAT != zero)); + use_strat = (want_strat && (P_STRAT < P_iso)); + z_strat = (use_strat) ? z_iso + (R_d*TISO/CONST_GRAV)*std::log(P_iso/P_STRAT) + : z_iso; + + // A configured stratospheric layer that lies at or below the isothermal + // transition cannot be represented, so say so rather than dropping it quietly. + if (want_strat && !use_strat) { + amrex::Print() << "WARNING: the base stratospheric layer is being ignored: P_STRAT = " + << P_STRAT << " Pa is not below the pressure at the base of the " + << "isothermal layer, P_iso = " << P_iso << " Pa.\n"; + amrex::Print() << " TLP_STRAT = " << TLP_STRAT << " will have no effect and " + << "the atmosphere above z_iso will be isothermal at TISO = " + << TISO << " K.\n"; + } + + amrex::Print() << "Base state layer interfaces: z_iso = " << z_iso << " m"; + if (use_strat) amrex::Print() << ", z_strat = " << z_strat << " m"; + amrex::Print() << "\n"; + + is_set = true; + } + + // + // True if this holds the same six defining parameters as bsp. Used to warn when + // a nested wrfinput file disagrees with its parent: the levels of an AMR hierarchy + // must share one reference profile, or their base states cannot be consistent + // across the coarse/fine interface. + // + bool same_params_as (const BaseStateParams& bsp) const + { + return (T00 == bsp.T00 ) && (P00 == bsp.P00 ) && + (TLP == bsp.TLP ) && (TISO == bsp.TISO) && + (TLP_STRAT == bsp.TLP_STRAT) && (P_STRAT == bsp.P_STRAT); + } +}; + +#endif diff --git a/Source/DataStructs/Make.package b/Source/DataStructs/Make.package index 85fd922117..ccabafc80e 100644 --- a/Source/DataStructs/Make.package +++ b/Source/DataStructs/Make.package @@ -1,3 +1,4 @@ +CEXE_headers += ERF_BaseStateParams.H CEXE_headers += ERF_DampingStruct.H CEXE_headers += ERF_DataStruct.H CEXE_headers += ERF_DiffStruct.H diff --git a/Source/ERF.H b/Source/ERF.H index 077f78c5ea..bde12c4599 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -41,6 +41,7 @@ #include #include +#include #include #include #include @@ -776,6 +777,13 @@ public: void init_from_wrfinput (int lev, amrex::MultiFab& mf_PSFC); void init_from_metgrid (int lev); void init_from_ncfile (int lev); + + // Rebuild the WRF base state on the grids of base_state_mf at this level, using the + // level-0 reference parameters (wrf_bsp) and this level's z_phys_cc. This is the same + // construction init_from_wrfinput performs at initialization, factored out so that a + // level remade by a regrid gets the same base state rather than one interpolated from + // its parent. Requires z_phys_cc[lev] and wrf_bsp to be valid. + void rebuild_base_state_from_wrfinput (int lev, amrex::MultiFab& base_state_mf); #endif // ERF_USE_NETCDF #ifdef ERF_USE_WINDFARM @@ -931,6 +939,12 @@ private: void initHSE (); void initHSE (int lev); + // Fill base_state[lev] by conservative interpolation from base_state[lev-1], then apply + // the base-state physical bc's. This is how the part of a fine level that lies inside + // the domain but outside the fine grids -- and the fine ghost cells -- get their values, + // so it must run before any construction that overwrites only the fine grids. + void interp_base_state_from_coarse (int lev); + //! Initialize Rayleigh damping profiles at a level void initRayleigh_at_level (const int& lev); @@ -1553,6 +1567,13 @@ private: amrex::Vector ba1d; amrex::Vector ba2d; + // The six parameters defining the WRF analytic reference profile, plus the layer + // interfaces derived from them. These are read from the LEVEL 0 wrfinput file only + // and then shared by every level: the levels of an AMR hierarchy must be built from + // one reference profile, or their base states cannot agree across the coarse/fine + // interface. (A nested wrfinput file that disagrees produces a warning.) + BaseStateParams wrf_bsp; + std::unique_ptr wrf_C1H; std::unique_ptr wrf_C2H; std::unique_ptr wrf_MUB; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 9da8769fdc..c03804172c 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -2057,7 +2057,51 @@ ERF::init_only (int lev, double elapsed_time) } else if (solverChoice.init_type == InitType::WRFInput && nc_init_file[lev].empty()) { - amrex::Abort("This pathway is not quite implemented yet"); + // + // A refined level with no wrfinput file of its own: the level was created by + // ERF's own refinement machinery rather than read from a wrfinput_d0*. Level 0 + // must always have a file (checked in ParameterSanityChecks), so this is lev > 0. + // + // Note the caller has already built the terrain at this level for us -- see the + // ordering in MakeNewLevelFromScratch -- because unlike the has-a-file case there + // is no file here to build z_phys from, and the base state below needs it. + // + AMREX_ALWAYS_ASSERT(lev > 0); + AMREX_ALWAYS_ASSERT(z_phys_cc[lev] != nullptr); + + make_physbcs(lev); + + // Fill the part of this level that lies inside the domain but outside the fine + // grids, and the fine ghost cells, from the coarser level. + interp_base_state_from_coarse(lev); + + // Now build the base state on the fine grids by exactly the construction level 0 + // used: the analytic reference profile from the level-0 parameters, evaluated at + // this level's cell-centered heights, followed by the discrete hydrostatic + // rebalance. Running the same procedure on both levels is what makes their base + // states agree; interpolating the coarse base state instead would leave this level + // not discretely hydrostatic on its own mesh. + rebuild_base_state_from_wrfinput(lev, base_state[lev]); + (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect()); + + // The state itself is interpolated from the coarser level. This must come after the + // base state, since it interpolates perturbational quantities relative to it. + FillCoarsePatch(lev, elapsed_time); + + // The 2D/surface arrays likewise have no file to come from at this level + Interp2DArrays(lev, ba2d[lev], dmap[lev]); + + // PSFC is read from the file at a level that has one; here it can only be + // interpolated from the parent. The surface layer and the LSM both read it. + if (mf_PSFC[lev-1] && mf_PSFC[lev]) { + auto ngv = mf_PSFC[lev]->nGrowVect(); ngv[2] = 0; + InterpFromCoarseLevel(*mf_PSFC[lev], ngv, + IntVect(0,0,0), // do not fill ghost cells outside the domain + *mf_PSFC[lev-1], 0, 0, 1, + geom[lev-1], geom[lev], + refRatio(lev-1), &cell_cons_interp, + domain_bcs_type, BCVars::cons_bc); + } } else if (solverChoice.init_type == InitType::NCFile) { diff --git a/Source/ERF_MakeNewArrays.cpp b/Source/ERF_MakeNewArrays.cpp index b30e3a9c22..452929f398 100644 --- a/Source/ERF_MakeNewArrays.cpp +++ b/Source/ERF_MakeNewArrays.cpp @@ -803,11 +803,21 @@ ERF::init_zphys (int lev, double elapsed_time) } // lev == 0 } else { - // NOTE: If a WRFInput file is NOT provided for a finer level, - // we simply interpolate from the coarse. This is necessary - // since we average_down the terrain (ERF_MakeNewLevel.cpp L351). - // If a WRFInput file IS present, it overwrites the terrain data. - if (lev > 0) { + // If a WRFInput / met_em file is NOT provided for a finer level, we simply + // interpolate the terrain from the coarse level. This is necessary since we + // average_down the terrain (see ERF_MakeNewLevel.cpp). + // + // If a file IS present at this level, the terrain has already been built from + // that file's PH + PHB by init_terrain_from_wrfinput (or its metgrid analogue), + // which runs BEFORE this routine on both the from-scratch and the from-coarse + // paths (init_only precedes init_zphys there because level 0 has to read its + // terrain and its data in the same pass). Interpolating from the coarse level + // here would therefore throw that away and leave the fine level running on + // coarse-derived heights -- so we must not do it. + // + // NOTE: this is the behavior the comment here has always described; the guard + // below is what actually makes it true. + if ( (lev > 0) && nc_init_file[lev].empty() ) { // // First interpolate from coarser level if there is one // NOTE: this interpolater assumes that ALL ghost cells of the coarse MultiFab @@ -962,6 +972,23 @@ ERF::remake_zphys (int lev, std::unique_ptr& temp_zphys_nd) refRatio(lev-1), &node_bilinear_interp, domain_bcs_type, BCVars::cons_bc); + // + // If this level has its own wrfinput / met_em file then its terrain was built + // from that file's PH + PHB, and the coarse interpolation above is only the + // right answer in the parts of the NEW grids that the OLD grids did not cover. + // ParallelCopy writes only where source and destination overlap, so this keeps + // the file-derived terrain wherever we already had it and leaves the + // coarse-interpolated values everywhere else. + // + // Without this, every regrid would quietly demote a nested level back onto + // coarse-interpolated heights -- which is the same thing the guard in + // init_zphys above exists to prevent at initialization. + // + if (!nc_init_file[lev].empty()) { + temp_zphys_nd->ParallelCopy(*z_phys_nd[lev], 0, 0, 1, + z_phys_nd[lev]->nGrowVect(), z_phys_nd[lev]->nGrowVect()); + } + } // lev > 0 std::swap(temp_zphys_nd, z_phys_nd[lev]); diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp index af022e054d..1a2822ae0a 100644 --- a/Source/ERF_MakeNewLevel.cpp +++ b/Source/ERF_MakeNewLevel.cpp @@ -162,10 +162,28 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, // // Note that "time" here is elapsed time, and start_time is the start_time from wrfinput/metgrid files // - init_only(lev, time); - init_zphys(lev, time); - update_terrain_arrays(lev); - make_physbcs(lev); + if (nc_init_file[lev].empty()) { + // + // A refined level with no file of its own has no PH/PHB to build terrain + // from, so its terrain must be interpolated from the coarser level FIRST -- + // the base state init_only builds is a function of the cell-centered heights + // and cannot run before they exist. (init_zphys interpolates from coarse + // exactly when this level has no file, which is the same condition.) + // + init_zphys(lev, time); + update_terrain_arrays(lev); + make_physbcs(lev); + init_only(lev, time); + } else { + // + // A level that does have a file reads its terrain and its data in the same + // pass, so init_only must come first here. + // + init_only(lev, time); + init_zphys(lev, time); + update_terrain_arrays(lev); + make_physbcs(lev); + } } else { // // Note that "time" here is elapsed time, and start_time = 0 when not using wrfinput/metgrid @@ -458,6 +476,20 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, } else { #endif + +#ifdef ERF_USE_NETCDF + // + // A WRFInput level created by refinement rather than read from a file: build its base + // state by the same construction level 0 used, on top of the coarse interpolation done + // above. This must happen BEFORE FillCoarsePatch, which interpolates perturbational + // quantities relative to the base state. + // + if ( (solverChoice.init_type == InitType::WRFInput) && nc_init_file[lev].empty() ) { + rebuild_base_state_from_wrfinput(lev, base_state[lev]); + (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect()); + } +#endif + // // Interpolate the solution data // @@ -698,6 +730,28 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // Impose bc's outside the domain (*physbcs_base[lev])(temp_base_state,0,temp_base_state.nComp(),base_state[lev].nGrowVect()); +#ifdef ERF_USE_NETCDF + // ************************************************************************************************* + // For a WRFInput run, rebuild the base state on the new grids by the SAME construction + // used at initialization -- the analytic reference profile from the level-0 parameters, + // evaluated at this level's heights, followed by the discrete hydrostatic rebalance -- + // rather than leaving it as the conservative interpolation of the parent's base state. + // + // Without this, the first regrid silently replaces a base state that is discretely + // hydrostatic on this level's mesh with one that is not, and the level stops agreeing + // with its parent in the way it did at t = 0. + // + // Note this must come AFTER the FillPatchTwoLevels above (which fills the region of + // this level outside the fine grids, and the ghost cells) and BEFORE the + // FillPatchFineLevel below, which interpolates perturbational quantities relative to + // the new base state. + // ************************************************************************************************* + if (solverChoice.init_type == InitType::WRFInput) { + rebuild_base_state_from_wrfinput(lev, temp_base_state); + (*physbcs_base[lev])(temp_base_state,0,temp_base_state.nComp(),base_state[lev].nGrowVect()); + } +#endif + // ************************************************************************************************* // This will fill the temporary MultiFabs with data from vars_new // NOTE: the momenta here are only used as scratch space, the momenta themselves are not fillpatched diff --git a/Source/Initialization/ERF_Init1D.cpp b/Source/Initialization/ERF_Init1D.cpp index a272422077..437af5fb3a 100644 --- a/Source/Initialization/ERF_Init1D.cpp +++ b/Source/Initialization/ERF_Init1D.cpp @@ -32,8 +32,6 @@ ERF::initHSE (int lev) bool all_boxes_touch_bottom = true; Box domain(geom[lev].Domain()); - int icomp = 0; int ncomp = BaseState::num_comps; - if (lev == 0) { BoxArray ba(base_state[lev].boxArray()); for (int i = 0; i < ba.size(); i++) { @@ -45,26 +43,12 @@ ERF::initHSE (int lev) else { // - // We need to do this interp from coarse level in order to set the values of - // the base state inside the domain but outside of the fine region - // - base_state[lev-1].FillBoundary(geom[lev-1].periodicity()); + // Set the values of the base state inside the domain but outside of the fine region, + // and in the fine ghost cells; the vertical integration below then overwrites only + // the fine grids, and for a box whose klo is in the interior it reads the value this + // leaves in klo-1 as its starting point. // - // NOTE: this interpolater assumes that ALL ghost cells of the coarse MultiFab - // have been pre-filled - this includes ghost cells both inside and outside - // the domain - // - InterpFromCoarseLevel(base_state[lev], base_state[lev].nGrowVect(), - IntVect(0,0,0), // do not fill ghost cells outside the domain - base_state[lev-1], icomp, icomp, ncomp, - geom[lev-1], geom[lev], - refRatio(lev-1), &cell_cons_interp, - domain_bcs_type, BCVars::cons_bc); - - // We need to do this here because the interpolation above may leave corners unfilled - // when the corners need to be filled by, for example, reflection of the fine ghost - // cell outside the fine region but inide the domain. - (*physbcs_base[lev])(base_state[lev],icomp,ncomp,base_state[lev].nGrowVect()); + interp_base_state_from_coarse(lev); } bool is_constant_dz = (solverChoice.mesh_type == MeshType::ConstantDz); @@ -177,6 +161,43 @@ ERF::initHSE (int lev) (*physbcs_base[lev])(base_state[lev],0,base_state[lev].nComp(),base_state[lev].nGrowVect()); } +/** + * Fill base_state[lev] by conservative interpolation from the next coarser level and then + * apply the base-state physical boundary conditions. + * + * This is what gives values to the part of a fine level that lies inside the domain but + * outside the fine grids, and to the fine ghost cells. Any construction that writes only + * the fine grids -- initHSE's vertical integration, or the analytic WRF profile -- must run + * after this, not before. + * + * @param lev Integer specifying the current level; must be > 0 + */ +void +ERF::interp_base_state_from_coarse (int lev) +{ + AMREX_ALWAYS_ASSERT(lev > 0); + + int icomp = 0; int ncomp = BaseState::num_comps; + + base_state[lev-1].FillBoundary(geom[lev-1].periodicity()); + // + // NOTE: this interpolater assumes that ALL ghost cells of the coarse MultiFab + // have been pre-filled - this includes ghost cells both inside and outside + // the domain + // + InterpFromCoarseLevel(base_state[lev], base_state[lev].nGrowVect(), + IntVect(0,0,0), // do not fill ghost cells outside the domain + base_state[lev-1], icomp, icomp, ncomp, + geom[lev-1], geom[lev], + refRatio(lev-1), &cell_cons_interp, + domain_bcs_type, BCVars::cons_bc); + + // We need to do this here because the interpolation above may leave corners unfilled + // when the corners need to be filled by, for example, reflection of the fine ghost + // cell outside the fine region but inside the domain. + (*physbcs_base[lev])(base_state[lev],icomp,ncomp,base_state[lev].nGrowVect()); +} + /** * Initialize density and pressure base state on every active AMR level. */ diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index dca6e0ccc7..80ff09c36e 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -104,7 +104,7 @@ ERF::init_from_metgrid (int lev) AMREX_ALWAYS_ASSERT(metgrid_use_sfc); // Base state parameters from WRF and the layer interfaces derived from them. - MetgridBaseStateParams bsp; + BaseStateParams bsp; bsp.set_layer_interfaces(); // Size the SST and LANDMASK @@ -1368,10 +1368,10 @@ init_base_state_from_metgrid (const bool use_moisture, FArrayBox& z_phys_nd_fab, FArrayBox& z_phys_cc_fab, const FArrayBox& NC_psfc_fab, - const MetgridBaseStateParams& bsp) + const BaseStateParams& bsp) { // Base state parameters and the layer interfaces derived from them. These are - // set once by the caller; see MetgridBaseStateParams for the profile they define. + // set once by the caller; see BaseStateParams for the profile they define. const Real T00 = bsp.T00; const Real P00 = bsp.P00; const Real TLP = bsp.TLP; diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index d2d13582eb..3d18d88f6d 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -130,12 +130,8 @@ init_terrain_from_wrfinput (int lev, * @param[in] mf_PB MultiFab holding WRF data specifying base state pressure (currently unused). * @param[in] mf_ALB MultiFab holding inverse density perturbation data (currently unused). * @param[in] z_phys_cc Cell-centered z-coordinate MultiFab; required, must not be null. - * @param[in] T00 Sea-level base-state temperature. - * @param[in] P00 Sea-level base-state pressure. - * @param[in] TLP Base-state lapse rate. - * @param[in] TISO Isothermal stratosphere temperature. - * @param[in] TLP_STRAT Stratospheric lapse rate. - * @param[in] P_STRAT Pressure at the stratosphere transition. + * @param[in] bsp Reference-state parameters and the layer interfaces derived from them; + * set_layer_interfaces() must already have been called on it. */ void init_base_state_from_wrfinput (const Box& subdomain, @@ -148,12 +144,7 @@ init_base_state_from_wrfinput (const Box& subdomain, MultiFab& mf_PB, MultiFab* mf_ALB, MultiFab* z_phys_cc, - const Real& T00, - const Real& P00, - const Real& TLP, - const Real& TISO, - const Real& TLP_STRAT, - const Real& P_STRAT); + const BaseStateParams& bsp); /** * Read start_time from the first WRF input file. @@ -197,22 +188,20 @@ read_start_time_from_wrfinput (int lev, const std::string& fname) * Read WRF base-state thermodynamic parameters from a NetCDF file. * * @param fname Path to the WRF input file - * @param T00 Sea-level base-state temperature - * @param P00 Sea-level base-state pressure - * @param TLP Base-state lapse rate - * @param TISO Isothermal stratosphere temperature - * @param TLP_STRAT Stratospheric lapse rate - * @param P_STRAT Pressure at the stratosphere transition + * @param bsp Reference-state parameters to fill; any value the file does not supply + * (or supplies as invalid) keeps the value bsp already held. */ void read_base_state_params_from_wrfinput (const std::string& fname, - Real& T00, - Real& P00, - Real& TLP, - Real& TISO, - Real& TLP_STRAT, - Real& P_STRAT) + BaseStateParams& bsp) { + Real& T00 = bsp.T00; + Real& P00 = bsp.P00; + Real& TLP = bsp.TLP; + Real& TISO = bsp.TISO; + Real& TLP_STRAT = bsp.TLP_STRAT; + Real& P_STRAT = bsp.P_STRAT; + if (ParallelDescriptor::IOProcessor()) { auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); @@ -410,16 +399,41 @@ ERF::init_from_wrfinput (int lev, MultiFab& mf_PSFC_lev) MultiFab mf_PH, mf_PB, mf_P; // For geopotential and base state std::unique_ptr mf_ALB; // For density - // Read base state params (used if ALB is not read) - Real T00 = Real(290.0); - Real P00 = p_0; - Real TLP = Real(50.0); - Real TISO = Real(200.0); - Real TLP_STRAT = Real(-11.0); - Real P_STRAT = zero; - read_base_state_params_from_wrfinput(nc_init_file[lev][0], - T00, P00, TLP, TISO, - TLP_STRAT, P_STRAT); + // ************************************************************************** + // Read the six parameters that define the analytic reference profile. + // + // These are read from the LEVEL 0 file only, and every finer level then uses the + // level-0 values. The whole base state is built from these six numbers, so two + // levels that disagree about them cannot have base states that agree across the + // coarse/fine interface -- and a nested wrfinput file very often does disagree, + // because idealized and hand-built files routinely declare T00/P00/TLP but leave + // them zero-filled, in which case read_base_state_params_from_wrfinput falls back + // to the ERF defaults. Previously each level read its own file and could silently + // end up on a different profile than its parent. + // + // If the file at this level does carry its own values and they differ from level + // zero's, say so rather than picking one silently. + // ************************************************************************** + if (lev == 0) { + read_base_state_params_from_wrfinput(nc_init_file[lev][0], wrf_bsp); + wrf_bsp.set_layer_interfaces(); + } else { + BaseStateParams lev_bsp; + read_base_state_params_from_wrfinput(nc_init_file[lev][0], lev_bsp); + if (!lev_bsp.same_params_as(wrf_bsp)) { + Print() << "WARNING: the base state parameters in " << nc_init_file[lev][0] + << " differ from those at level 0:\n"; + Print() << " level " << lev << " file: (T00, P00, TLP, TISO, TLP_STRAT, P_STRAT) = (" + << lev_bsp.T00 << ", " << lev_bsp.P00 << ", " << lev_bsp.TLP << ", " + << lev_bsp.TISO << ", " << lev_bsp.TLP_STRAT << ", " << lev_bsp.P_STRAT << ")\n"; + Print() << " level 0 : (T00, P00, TLP, TISO, TLP_STRAT, P_STRAT) = (" + << wrf_bsp.T00 << ", " << wrf_bsp.P00 << ", " << wrf_bsp.TLP << ", " + << wrf_bsp.TISO << ", " << wrf_bsp.TLP_STRAT << ", " << wrf_bsp.P_STRAT << ")\n"; + Print() << " Using the level 0 values at every level, so that the base states " + "of the two levels are built from one reference profile.\n"; + } + } + AMREX_ALWAYS_ASSERT(wrf_bsp.is_set); // Temporary MFs for derived quantities auto& ba = lev_new[Vars::cons].boxArray(); @@ -1328,6 +1342,39 @@ ERF::init_from_wrfinput (int lev, MultiFab& mf_PSFC_lev) int klo = geom[lev].Domain().smallEnd(2); int khi = geom[lev].Domain().bigEnd(2); + // + // KNOWN LIMITATION for a level whose grids do not reach the top of the domain + // (a nest refined over only the lower part of the column). + // + // The remap below searches the wrfinput column (mf_PH + mf_PHB) for the interval + // bracketing each ERF cell-centered height, and interpolates cons_tmp between the + // two bracketing source cells. All three of those MultiFabs live on this level's + // grids with no ghost cells in z, so above the top of a grid there is no source + // geopotential AND no source state to interpolate from. The cells near the top of + // such a grid therefore keep the values read straight from the file rather than the + // remapped ones. (Before the box bound was introduced, the search instead ran to + // the top of the domain and read past the end of those FABs.) + // + // Fixing this properly means reading PH, PHB and the state over a box grown upward + // in z beyond the grid, remapping, and then discarding the extra planes. + // + { + int max_grid_khi = -1; + for (int i = 0; i < lev_new[Vars::cons].boxArray().size(); ++i) { + max_grid_khi = std::max(max_grid_khi, lev_new[Vars::cons].boxArray()[i].bigEnd(2)); + } + if (max_grid_khi < khi) { + Print() << "WARNING: the grids at level " << lev << " reach only k = " << max_grid_khi + << " of a domain whose top is k = " << khi << ".\n"; + Print() << " The vertical remap of the wrfinput state onto the ERF grid " + "cannot be completed in the cells near the top of those grids, because " + "neither the source geopotential nor the source state is available " + "above them; those cells retain the values read from the file.\n"; + Print() << " The BASE STATE is not affected -- it is an analytic function " + "of height and is built (and rebalanced) correctly over the whole patch.\n"; + } + } + MultiFab cons_tmp(lev_new[Vars::cons].boxArray(), lev_new[Vars::cons].DistributionMap(), ncons, 0); MultiFab xvel_tmp(lev_new[Vars::xvel].boxArray(), lev_new[Vars::xvel].DistributionMap(), 1 , 0); MultiFab yvel_tmp(lev_new[Vars::yvel].boxArray(), lev_new[Vars::yvel].DistributionMap(), 1 , 0); @@ -1341,6 +1388,21 @@ ERF::init_from_wrfinput (int lev, MultiFab& mf_PSFC_lev) const Box& bxx = mfi.tilebox(IntVect(1,0,0)); const Box& bxy = mfi.tilebox(IntVect(0,1,0)); + // + // mf_PH and mf_PHB are defined on the zvel BoxArray with ZERO ghost cells in + // z (see their define() calls above), so within this box they only hold the + // nodes [vbx.smallEnd(2), vbx.bigEnd(2)+1]. The bracket searches below read + // ph_arr(...,lk+1), so lk must stop at vbx.bigEnd(2). Bounding those searches + // by the *domain* top instead would walk off the end of the FAB for any box + // that does not reach the top of the domain -- which is precisely the case on + // a refined level whose patch covers only the lower part of the domain. + // + // For a box that does span the domain in z this is exactly the old bound, so + // full-height runs are unaffected. + // + const Box& vbx = mfi.validbox(); + const int khi_src = amrex::min(khi, vbx.bigEnd(2)); + const Array4< Real>& cons_arr = lev_new[Vars::cons].array(mfi); const Array4< Real>& xvel_arr = lev_new[Vars::xvel].array(mfi); const Array4< Real>& yvel_arr = lev_new[Vars::yvel].array(mfi); @@ -1365,7 +1427,7 @@ ERF::init_from_wrfinput (int lev, MultiFab& mf_PSFC_lev) bool found = false; int kend = kstart; - for (int lk(kstart+1); lk<=khi; ++lk) { + for (int lk(kstart+1); lk<=khi_src; ++lk) { Real z_hi_src = Real(0.5) * ( ph_arr(i,j,lk ) + phb_arr(i,j,lk ) + ph_arr(i,j,lk+1) + phb_arr(i,j,lk+1)) / CONST_GRAV; @@ -1409,7 +1471,7 @@ ERF::init_from_wrfinput (int lev, MultiFab& mf_PSFC_lev) bool found = false; int kend = kstart; - for (int lk(kstart+1); lk<=khi; ++lk) { + for (int lk(kstart+1); lk<=khi_src; ++lk) { Real z_hi_src = Real(0.25) * ( ph_arr(ii ,j,lk ) + phb_arr(ii ,j,lk ) + ph_arr(ii ,j,lk+1) + phb_arr(ii ,j,lk+1) + @@ -1454,7 +1516,7 @@ ERF::init_from_wrfinput (int lev, MultiFab& mf_PSFC_lev) bool found = false; int kend = kstart; - for (int lk(kstart+1); lk<=khi; ++lk) { + for (int lk(kstart+1); lk<=khi_src; ++lk) { Real z_hi_src = Real(0.25) * ( ph_arr(i,jj ,lk ) + phb_arr(i,jj ,lk ) + ph_arr(i,jj ,lk+1) + phb_arr(i,jj ,lk+1) + @@ -1520,23 +1582,12 @@ ERF::init_from_wrfinput (int lev, MultiFab& mf_PSFC_lev) // ************************************************************************** // Initialize the base state // ************************************************************************** - MultiFab r_hse (base_state[lev], make_alias, BaseState::r0_comp, 1); - MultiFab p_hse (base_state[lev], make_alias, BaseState::p0_comp, 1); - MultiFab pi_hse(base_state[lev], make_alias, BaseState::pi0_comp, 1); - MultiFab th_hse(base_state[lev], make_alias, BaseState::th0_comp, 1); - MultiFab qv_hse(base_state[lev], make_alias, BaseState::qv0_comp, 1); - - init_base_state_from_wrfinput(boxes_at_level[lev][0], l_rdOcp, - p_hse, pi_hse, th_hse, qv_hse, r_hse, - mf_PB, mf_ALB.get(), z_phys_cc[lev].get(), - T00, P00, TLP, TISO, TLP_STRAT, P_STRAT); + rebuild_base_state_from_wrfinput(lev, base_state[lev]); - // FillBoundary to populate the internal ghost cells (no averaging in above call) - r_hse.FillBoundary(geom[lev].periodicity()); - p_hse.FillBoundary(geom[lev].periodicity()); - pi_hse.FillBoundary(geom[lev].periodicity()); - th_hse.FillBoundary(geom[lev].periodicity()); - qv_hse.FillBoundary(geom[lev].periodicity()); + // rho0 is consumed below by read_and_convert_from_wrfbdy (via scale_bdy_normal_by_rho0), + // including one cell outside the domain -- which is why the routine above fills the + // base state's ghost cells as well as its valid region. + MultiFab r_hse (base_state[lev], make_alias, BaseState::r0_comp, 1); // ******************************************************************************************* // Initialize the bdy data @@ -1770,12 +1821,8 @@ ERF::init_from_wrfinput (int lev, MultiFab& mf_PSFC_lev) * @param mf_PB MultiFab holding WRF data specifying base state pressure; currently unused * @param mf_ALB MultiFab holding inverse density perturbation data; currently unused * @param z_phys_cc Cell-centered z-coordinate MultiFab; required, must not be null - * @param T00 Sea-level base-state temperature - * @param P00 Sea-level base-state pressure - * @param TLP Base-state lapse rate - * @param TISO Isothermal stratosphere temperature - * @param TLP_STRAT Stratospheric lapse rate - * @param P_STRAT Pressure at the stratosphere transition + * @param bsp Reference-state parameters and the layer interfaces derived from + * them; set_layer_interfaces() must already have been called */ void init_base_state_from_wrfinput (const Box& subdomain, @@ -1788,15 +1835,9 @@ init_base_state_from_wrfinput (const Box& subdomain, MultiFab& /*mf_PB*/, MultiFab* /*mf_ALB*/, MultiFab* z_phys_cc, - const Real& T00, - const Real& P00, - const Real& TLP, - const Real& TISO, - const Real& TLP_STRAT, - const Real& P_STRAT) + const BaseStateParams& bsp) { const auto& dom_lo = lbound(subdomain); - const auto& dom_hi = ubound(subdomain); // The analytic inversion below is a function of the true cell-centered // height, so z_phys_cc is required here -- it is dereferenced unconditionally @@ -1808,35 +1849,25 @@ init_base_state_from_wrfinput (const Box& subdomain, // (1) troposphere T = T00 + TLP * ln(p/P00) p > P_iso // (2) isothermal layer T = TISO P_STRAT < p <= P_iso // (3) stratosphere T = TISO + TLP_STRAT * ln(p/P_STRAT) p <= P_STRAT - // Each piece inverts to p(z) in closed form under dp/dz = -rho g; here we - // precompute the two interface heights so the inversion can branch on z. + // Each piece inverts to p(z) in closed form under dp/dz = -rho g. The two + // interface heights were precomputed once by bsp.set_layer_interfaces(), which + // is also where the parameters were validated and reported; here we only need + // to copy them into locals that the device kernel can capture. // ************************************************************************** - const Real x_iso = (TISO - T00) / TLP; - const Real P_iso = P00 * std::exp(x_iso); - const Real z_iso = -(R_d/CONST_GRAV) * (T00*x_iso + myhalf*TLP*x_iso*x_iso); - - // The upper stratospheric layer is optional (P_STRAT == 0 or TLP_STRAT == 0 - // disables it) and is only meaningful if it begins above the isothermal layer, - // i.e. if P_STRAT is below the pressure at which the isothermal layer starts. - const bool want_strat = (P_STRAT > zero) && (TLP_STRAT != zero); - const bool use_strat = want_strat && (P_STRAT < P_iso); - const Real z_strat = (use_strat) ? z_iso + (R_d*TISO/CONST_GRAV)*std::log(P_iso/P_STRAT) - : z_iso; - - // A configured stratospheric layer that lies at or below the isothermal - // transition cannot be represented, so say so rather than dropping it quietly. - if (want_strat && !use_strat) { - Print() << "WARNING: the WRF stratospheric layer is being ignored: P_STRAT = " - << P_STRAT << " Pa is not below the pressure at the base of the " - << "isothermal layer, P_iso = " << P_iso << " Pa.\n"; - Print() << " TLP_STRAT = " << TLP_STRAT << " will have no effect and " - << "the atmosphere above z_iso will be isothermal at TISO = " - << TISO << " K.\n"; - } - - Print() << "WRF base state layer interfaces: z_iso = " << z_iso << " m"; - if (use_strat) Print() << ", z_strat = " << z_strat << " m"; - Print() << "\n"; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(bsp.is_set, + "init_base_state_from_wrfinput requires that " + "set_layer_interfaces() has been called on bsp"); + + const Real T00 = bsp.T00; + const Real P00 = bsp.P00; + const Real TLP = bsp.TLP; + const Real TISO = bsp.TISO; + const Real TLP_STRAT = bsp.TLP_STRAT; + const Real P_STRAT = bsp.P_STRAT; + const Real P_iso = bsp.P_iso; + const Real z_iso = bsp.z_iso; + const Real z_strat = bsp.z_strat; + const bool use_strat = bsp.use_strat; #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -1913,7 +1944,6 @@ init_base_state_from_wrfinput (const Box& subdomain, // discretely satisfy dp0/dz = -rho0 g on the ERF grid // ************************************************************************** int k_dom_lo = dom_lo.z; - int k_dom_hi = dom_hi.z; // The vertical integration below is seeded with the analytic profile in the // lowest cell of each column, z_cc(i,j,klo) and p_hse(i,j,klo), not with @@ -1934,7 +1964,31 @@ init_base_state_from_wrfinput (const Box& subdomain, int klo = bx.smallEnd(2); int khi = bx.bigEnd(2); - AMREX_ALWAYS_ASSERT((klo == k_dom_lo) && (khi == k_dom_hi)); + // + // The integration below is bottom-up: the value in cell k depends only on + // cells at or below k. A box that stops below the top of the region is + // therefore perfectly well defined -- it produces exactly the values the + // full-height column would have had in the cells it does contain -- so we + // deliberately do NOT require khi to reach the top. That is what lets a + // refined level cover only the lower part of the domain, which is the normal + // way to nest an LES patch inside a mesoscale parent. Cells above such a + // patch simply do not exist at this level; they are covered by the parent. + // + // klo is a different matter: each column is seeded from the analytic profile + // in its lowest cell and then marched upward, so a box whose klo is in the + // interior would restart the discrete integration there and put a kink in + // p_0 at the box boundary. (With amr.refine_grid_layout_z = 0, which is the + // ERF default set in main.cpp, grids are never chopped in z and this holds.) + // + // NOTE: TileNoZ() above guarantees that klo/khi are the *box's* z extent + // rather than a tile's, so each column is integrated exactly once. + // + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(klo == k_dom_lo, + "init_base_state_from_wrfinput requires boxes that " + "reach the bottom of the region being initialized: " + "each column is seeded there and integrated upward. " + "Set erf.max_grid_size_z large enough that the grids " + "are not decomposed in z."); bx.makeSlab(2,klo); const Array4& r_hse_arr = r_hse.array(mfi); @@ -1999,6 +2053,52 @@ init_base_state_from_wrfinput (const Box& subdomain, } // mfi } +/** + * Build the WRF base state on the grids of base_state_mf at this level. + * + * This is the construction init_from_wrfinput performs at initialization, factored out so + * that a level remade by a regrid can run exactly the same procedure rather than inheriting + * a base state interpolated from its parent. Every level uses the SAME reference + * parameters (wrf_bsp, read from the level-0 file), evaluated at this level's own + * cell-centered heights -- that is what makes the base states of the levels consistent + * with one another. + * + * @param lev Level to build the base state at + * @param base_state_mf Base state MultiFab to fill (base_state[lev], or a temporary + * standing in for it during a regrid) + */ +void +ERF::rebuild_base_state_from_wrfinput (int lev, MultiFab& base_state_mf) +{ + AMREX_ALWAYS_ASSERT(z_phys_cc[lev] != nullptr); + + MultiFab r_hse (base_state_mf, make_alias, BaseState::r0_comp, 1); + MultiFab p_hse (base_state_mf, make_alias, BaseState::p0_comp, 1); + MultiFab pi_hse(base_state_mf, make_alias, BaseState::pi0_comp, 1); + MultiFab th_hse(base_state_mf, make_alias, BaseState::th0_comp, 1); + MultiFab qv_hse(base_state_mf, make_alias, BaseState::qv0_comp, 1); + + // mf_PB / mf_ALB are unused by the routine (the base state is built entirely from the + // six reference parameters and z_phys_cc), so there is nothing to hand it here. + // + // We pass the level's Domain rather than boxes_at_level[lev][0] because the only thing + // the routine takes from this box is the k index at which each column starts, and the + // two agree there; the Domain is also always available, whereas boxes_at_level[lev] is + // only populated when the refinement was specified as an explicit box. + MultiFab dummy_PB; + init_base_state_from_wrfinput(geom[lev].Domain(), solverChoice.rdOcp, + p_hse, pi_hse, th_hse, qv_hse, r_hse, + dummy_PB, nullptr, z_phys_cc[lev].get(), + wrf_bsp); + + // FillBoundary to populate the internal ghost cells (no averaging in above call) + r_hse.FillBoundary(geom[lev].periodicity()); + p_hse.FillBoundary(geom[lev].periodicity()); + pi_hse.FillBoundary(geom[lev].periodicity()); + th_hse.FillBoundary(geom[lev].periodicity()); + qv_hse.FillBoundary(geom[lev].periodicity()); +} + /** * Helper function for verifying the top boundary is valid and computing the bottom boundary. * @@ -2212,6 +2312,20 @@ init_terrain_from_wrfinput (int /*lev*/, int klo = z_face_box.smallEnd()[2]; int khi = z_face_box.bigEnd()[2]; + // + // z_top is the height of the top of the *domain*, so it is only the correct + // value for the node at the top of the domain. khi above is the top node of + // `subdomain`, which for a refined level is boxes_at_level[lev][0] -- and that + // box need not reach the domain top: a nested LES patch is routinely refined + // over only the lower part of the column. Forcing z_top at khi in that case + // would stretch the fine column all the way to the model top and produce a + // terrain mapping that has nothing to do with the parent's. So only pin the + // top node to z_top when this really is the top of the domain; otherwise keep + // the height averaged from the file, and extrapolate above it from that. + // + const int khi_dom = surroundingNodes(geom.Domain(),2).bigEnd(2); + const bool box_at_top = (khi == khi_dom); + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { int ii = std::max(std::min(i,ihi),ilo); @@ -2238,16 +2352,23 @@ init_terrain_from_wrfinput (int /*lev*/, nc_ph_arr (ii,jm,khi-1) + nc_ph_arr (im,jm,khi-1) + nc_phb_arr(ii,jj,khi-1) + nc_phb_arr(im,jj,khi-1) + nc_phb_arr(ii,jm,khi-1) + nc_phb_arr(im,jm,khi-1) ) / CONST_GRAV; - // Extrapolate linearly above the top of the domain -- note that - // z_phys_nd has more than one ghost node in the vertical, so this - // must depend on k rather than filling every ghost node with one value - z_arr(i, j, k) = z_top + static_cast(k-khi) * (z_top - z_khim1); + Real z_khi = Real(0.25) * ( nc_ph_arr (ii,jj,khi ) + nc_ph_arr (im,jj,khi ) + + nc_ph_arr (ii,jm,khi ) + nc_ph_arr (im,jm,khi ) + + nc_phb_arr(ii,jj,khi ) + nc_phb_arr(im,jj,khi ) + + nc_phb_arr(ii,jm,khi ) + nc_phb_arr(im,jm,khi ) ) / CONST_GRAV; + // Extrapolate linearly above the top node -- note that z_phys_nd has + // more than one ghost node in the vertical, so this must depend on k + // rather than filling every ghost node with one value. The reference + // is z_top only if this box reaches the top of the domain; for a patch + // that stops lower down it is that patch's own top height. + Real z_ref = (box_at_top) ? z_top : z_khi; + z_arr(i, j, k) = z_ref + static_cast(k-khi) * (z_ref - z_khim1); } else if (k == khi) { - z_arr(i, j, k) = Real(0.25) * ( nc_ph_arr (ii,jj,k) + nc_ph_arr (im,jj,k) + - nc_ph_arr (ii,jm,k) + nc_ph_arr (im,jm,k) + - nc_phb_arr(ii,jj,k) + nc_phb_arr(im,jj,k) + - nc_phb_arr(ii,jm,k) + nc_phb_arr(im,jm,k) ) / CONST_GRAV; - z_arr(i, j, k) = z_top; + Real z_khi = Real(0.25) * ( nc_ph_arr (ii,jj,k) + nc_ph_arr (im,jj,k) + + nc_ph_arr (ii,jm,k) + nc_ph_arr (im,jm,k) + + nc_phb_arr(ii,jj,k) + nc_phb_arr(im,jj,k) + + nc_phb_arr(ii,jm,k) + nc_phb_arr(im,jm,k) ) / CONST_GRAV; + z_arr(i, j, k) = (box_at_top) ? z_top : z_khi; } else { // Note: wrfinput geopotentials ph, phb are only staggered in the vertical, i.e., // they have dims (bottom_top_stag, south_north, west_east). On k==klo, we diff --git a/Source/Initialization/ERF_MetgridUtils.H b/Source/Initialization/ERF_MetgridUtils.H index 2570b67522..e4be8e5613 100644 --- a/Source/Initialization/ERF_MetgridUtils.H +++ b/Source/Initialization/ERF_MetgridUtils.H @@ -12,6 +12,7 @@ #include #include #include +#include /** * Read metgrid NetCDF data and populate ERF initialization FABs. @@ -116,82 +117,6 @@ init_msfs_from_metgrid (const bool metgrid_debug_msf, amrex::FArrayBox& NC_MSFV_fab, amrex::FArrayBox& NC_MSFM_fab); -/** - * Parameters defining the WRF reference (base) state, along with the layer - * interface heights derived from them. - * - * The base state is piecewise in log-pressure: - * (1) troposphere T = T00 + TLP * ln(p/P00) p > P_iso - * (2) isothermal layer T = TISO P_STRAT < p <= P_iso - * (3) stratosphere T = TISO + TLP_STRAT * ln(p/P_STRAT) p <= P_STRAT - * Each piece inverts to p(z) in closed form under dp/dz = -rho g, so the two - * interface heights are precomputed once and the inversion branches on z. - * - * NOTE: metgrid files do not carry these parameters, so we default to the WRF - * namelist defaults for base_temp, base_pres, base_lapse, iso_temp, - * base_lapse_strat and base_pres_strat, respectively. - */ -struct MetgridBaseStateParams -{ - amrex::Real T00 = amrex::Real(290.0); - amrex::Real P00 = p_0; - amrex::Real TLP = amrex::Real(50.0); - amrex::Real TISO = amrex::Real(200.0); - amrex::Real TLP_STRAT = amrex::Real(-11.0); - amrex::Real P_STRAT = zero; - - // Derived by set_layer_interfaces() - amrex::Real P_iso = zero; - amrex::Real z_iso = zero; - amrex::Real z_strat = zero; - bool use_strat = false; - - // - // Compute the layer interfaces and report them. - // - // NOTE: This must be called from outside of any MFIter (and outside of any - // OpenMP parallel region) so the diagnostics are emitted exactly once - // rather than once per box per thread. - // - void set_layer_interfaces () - { - // The troposphere inverts p(z) as sqrt((T00/TLP)^2 - 2*g*z/(TLP*R_d)), which - // folds 1/TLP inside the square root and so discards the sign of TLP. A - // non-positive lapse parameter, or an isothermal temperature at or above the - // surface temperature, would silently select the wrong root instead of failing. - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(TLP > zero, - "Base state TLP must be positive"); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(TISO < T00, - "Base state TISO must be less than T00"); - - const amrex::Real x_iso = (TISO - T00) / TLP; - P_iso = P00 * std::exp(x_iso); - z_iso = -(R_d/CONST_GRAV) * (T00*x_iso + myhalf*TLP*x_iso*x_iso); - - // The upper stratospheric layer is optional (P_STRAT == 0 or TLP_STRAT == 0 - // disables it) and is only meaningful if it begins above the isothermal layer, - // i.e. if P_STRAT is below the pressure at which the isothermal layer starts. - const bool want_strat = ((P_STRAT > zero) && (TLP_STRAT != zero)); - use_strat = (want_strat && (P_STRAT < P_iso)); - z_strat = (use_strat) ? z_iso + (R_d*TISO/CONST_GRAV)*std::log(P_iso/P_STRAT) - : z_iso; - - // A configured stratospheric layer that lies at or below the isothermal - // transition cannot be represented, so say so rather than dropping it quietly. - if (want_strat && !use_strat) { - amrex::Print() << "WARNING: the base stratospheric layer is being ignored: P_STRAT = " - << P_STRAT << " Pa is not below the pressure at the base of the " - << "isothermal layer, P_iso = " << P_iso << " Pa.\n"; - amrex::Print() << " TLP_STRAT = " << TLP_STRAT << " will have no effect and " - << "the atmosphere above z_iso will be isothermal at TISO = " - << TISO << " K.\n"; - } - - amrex::Print() << "Base state layer interfaces: z_iso = " << z_iso << " m"; - if (use_strat) amrex::Print() << ", z_strat = " << z_strat << " m"; - amrex::Print() << "\n"; - } -}; /** * Initialize hydrostatic base-state quantities from metgrid data. @@ -212,7 +137,7 @@ init_base_state_from_metgrid (const bool use_moisture, amrex::FArrayBox& z_phys_nd_fab, amrex::FArrayBox& z_phys_cc_fab, const amrex::FArrayBox& NC_psfc_fab, - const MetgridBaseStateParams& bsp); + const BaseStateParams& bsp); /** * Interpolate a value with a Lagrange polynomial stencil. diff --git a/Source/Utils/ERF_Rebalance.cpp b/Source/Utils/ERF_Rebalance.cpp index 51c9fa5734..fa7d76128f 100644 --- a/Source/Utils/ERF_Rebalance.cpp +++ b/Source/Utils/ERF_Rebalance.cpp @@ -35,13 +35,35 @@ rebalance_columns (MultiFab& rho, // int ncomp = cons.nComp(); int k_dom_lo = geom.Domain().smallEnd(2); - int k_dom_hi = geom.Domain().bigEnd(2); for (MFIter mfi(rho,TileNoZ()); mfi.isValid(); ++mfi) { Box bx = mfi.tilebox(); int klo = bx.smallEnd(2); int khi = bx.bigEnd(2); - AMREX_ALWAYS_ASSERT((klo == k_dom_lo) && (khi == k_dom_hi)); + + // + // This is a bottom-up integration: the value in cell k depends only on cells at + // or below k. A box that stops below the top of the domain is therefore perfectly + // well defined -- it produces exactly the values the full-height column would have + // had in the cells it does contain -- so we deliberately do NOT require + // khi == geom.Domain().bigEnd(2) here. That matters for a refined level whose + // patch covers only the lower part of the domain, which is the normal way to nest + // an LES region inside a mesoscale parent. + // + // What the integration does require is a valid starting value in its lowest cell. + // The use_sfc seeding below marches up from p_0 at z = 0, so it is only meaningful + // for a box that reaches the ground; a box whose klo is in the interior has no + // surface to start from. (With amr.refine_grid_layout_z = 0, which is the ERF + // default set in main.cpp, grids are never chopped in z and this always holds.) + // + // NOTE: TileNoZ() above guarantees that klo/khi are the *box's* z extent rather + // than a tile's, so each column is integrated exactly once. + // + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!use_sfc || (klo == k_dom_lo), + "rebalance_columns with use_sfc requires boxes that " + "reach the bottom of the domain: the integration is " + "seeded from p_0 at the surface. Set erf.max_grid_size_z " + "large enough that the grids are not decomposed in z."); bx.makeSlab(2,klo); const Array4< Real>& rho_arr = rho.array(mfi); @@ -53,7 +75,7 @@ rebalance_columns (MultiFab& rho, ParallelFor(bx, [=,RdoCp_d=RdoCp] AMREX_GPU_DEVICE (int i, int j, int /*k*/) noexcept { - // integrate from surface to domain top + // Integrate upward from the bottom of this box to its top Real dz, F, C; Real rho_tot_hi, rho_tot_lo; Real z_lo, z_hi; From 78db2671336a18e64c536adb1deec14bb8a9271b Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 28 Aug 2026 07:28:34 -0700 Subject: [PATCH 2/2] Enable multilevel base state rebalancing --- Docs/sphinx_doc/Initialization.rst | 22 ++++++++++- Docs/sphinx_doc/MeshRefinement.rst | 55 ++++++++++++++++++---------- Docs/sphinx_doc/theory/BaseState.rst | 45 +++++++++++++++++++++++ Source/ERF.cpp | 15 +++++++- Source/ERF_MakeNewLevel.cpp | 30 +++++++++++++-- 5 files changed, 140 insertions(+), 27 deletions(-) diff --git a/Docs/sphinx_doc/Initialization.rst b/Docs/sphinx_doc/Initialization.rst index 64f10362a1..3e4fc4064b 100644 --- a/Docs/sphinx_doc/Initialization.rst +++ b/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``) @@ -256,7 +276,7 @@ Restrictions 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-wrfinput` for the mesh-refinement side of this, and +See :ref:`subsec:partial-depth-refinement` for the mesh-refinement side of this, and :ref:`MeshRefinement` for refinement in general. TKE Initialization diff --git a/Docs/sphinx_doc/MeshRefinement.rst b/Docs/sphinx_doc/MeshRefinement.rst index 4702be9172..277e3f690a 100644 --- a/Docs/sphinx_doc/MeshRefinement.rst +++ b/Docs/sphinx_doc/MeshRefinement.rst @@ -220,8 +220,8 @@ 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 when -nesting from WRF input files; see :ref:`subsec:partial-depth-wrfinput`. +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 @@ -233,16 +233,23 @@ 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-wrfinput: +.. _subsec:partial-depth-refinement: -Nesting from WRF Input Files over Only Part of the Depth -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +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 with ``erf.init_type = WRFInput``. This is supported, both when -each level has its own wrfinput file and when only level 0 does. For example, with a -domain 17500 m deep:: +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 @@ -258,18 +265,26 @@ Two things to watch, both specific to how the box is specified: vertical extent with the full depth of the domain, which is the opposite of what is wanted here. -- The z values are interpreted on the nominal (undeformed) vertical grid, not on the - terrain-following heights that a wrfinput file supplies, so the box may not reach the - height you expect. Check the box ERF reports (``Saving in 'boxes at level'``) against - the ``z_phys`` field in the plotfile. - -If a file is given for the finer level, the box must be contained in the region that file -covers; ERF checks this and aborts with both boxes printed if it is not. - -See :ref:`sec:nested-wrfinput` for how the base state is built at each level, why -truncating the refined region does not change the values in the cells it does contain, -and the restrictions that apply (in particular, the PBL models and the column-integral -derived quantities cannot be used on a level whose grids do not span the domain in z). +- 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 ~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/Docs/sphinx_doc/theory/BaseState.rst b/Docs/sphinx_doc/theory/BaseState.rst index 428e0990d0..2eddd4a4c1 100644 --- a/Docs/sphinx_doc/theory/BaseState.rst +++ b/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. diff --git a/Source/ERF.cpp b/Source/ERF.cpp index c03804172c..256965face 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -662,13 +662,24 @@ ERF::InitData_post () setPlotVariables2D("plot2d_vars_2", plot2d_var_names_2); // - // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one and if two way coupling + // Make sure that detJ is the average of the data on a finer level if there is one and if two way coupling // if (SolverChoice::mesh_type != MeshType::ConstantDz) { if (solverChoice.coupling_type == CouplingType::TwoWay) { + // NOTE: z_phys_cc is deliberately NOT averaged down. Every level's base state is + // built to be in discrete hydrostatic balance against that level's own + // cell-centered heights, so replacing the coarse heights with the average of + // the fine ones -- after the base states have already been built -- leaves the + // coarse base state out of balance with the heights the dycore then uses for + // vertical gradients, Rayleigh damping and the sponge zones. z_phys_cc is also + // derived from z_phys_nd, which is not averaged down either, so averaging only + // the cell-centered heights made the two disagree inside the refined region. + // + // detJ IS still averaged down: AverageDownTo weights (rho S) by detJ_cc before + // averaging and divides by it afterwards, so the coarse detJ must be the average + // of the fine detJ for that average-down to telescope and stay conservative. for (int crse_lev = finest_level-1; crse_lev >= 0; crse_lev--) { average_down( *detJ_cc[crse_lev+1], *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev)); - average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev)); } } for (int crse_lev = finest_level-1; crse_lev >= 0; crse_lev--) { diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp index 1a2822ae0a..237de08c47 100644 --- a/Source/ERF_MakeNewLevel.cpp +++ b/Source/ERF_MakeNewLevel.cpp @@ -367,13 +367,24 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, update_terrain_arrays(lev); // - // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one + // Make sure that detJ is the average of the data on a finer level if there is one // *and* if there is two-way coupling // if ( (SolverChoice::mesh_type != MeshType::ConstantDz) && (solverChoice.coupling_type == CouplingType::TwoWay) ) { + // NOTE: z_phys_cc is deliberately NOT averaged down. Every level's base state is + // built to be in discrete hydrostatic balance against that level's own + // cell-centered heights, so replacing the coarse heights with the average of + // the fine ones -- after the base states have already been built -- leaves the + // coarse base state out of balance with the heights the dycore then uses for + // vertical gradients, Rayleigh damping and the sponge zones. z_phys_cc is also + // derived from z_phys_nd, which is not averaged down either, so averaging only + // the cell-centered heights made the two disagree inside the refined region. + // + // detJ IS still averaged down: AverageDownTo weights (rho S) by detJ_cc before + // averaging and divides by it afterwards, so the coarse detJ must be the average + // of the fine detJ for that average-down to telescope and stay conservative. for (int crse_lev = lev-1; crse_lev >= 0; crse_lev--) { average_down( *detJ_cc[crse_lev+1], *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev)); - average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev)); } } @@ -682,14 +693,25 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp update_terrain_arrays(lev); // ******************************************************************************************** - // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one + // Make sure that detJ is the average of the data on a finer level if there is one // Note that this shouldn't be necessary because the fine grid is created by interpolation // from the coarse ... but just in case ... // ******************************************************************************************** if ( (SolverChoice::mesh_type != MeshType::ConstantDz) && (solverChoice.coupling_type == CouplingType::TwoWay) ) { + // NOTE: z_phys_cc is deliberately NOT averaged down. Every level's base state is + // built to be in discrete hydrostatic balance against that level's own + // cell-centered heights, so replacing the coarse heights with the average of + // the fine ones -- after the base states have already been built -- leaves the + // coarse base state out of balance with the heights the dycore then uses for + // vertical gradients, Rayleigh damping and the sponge zones. z_phys_cc is also + // derived from z_phys_nd, which is not averaged down either, so averaging only + // the cell-centered heights made the two disagree inside the refined region. + // + // detJ IS still averaged down: AverageDownTo weights (rho S) by detJ_cc before + // averaging and divides by it afterwards, so the coarse detJ must be the average + // of the fine detJ for that average-down to telescope and stay conservative. for (int crse_lev = lev-1; crse_lev >= 0; crse_lev--) { average_down( *detJ_cc[crse_lev+1], *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev)); - average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev)); } }