diff --git a/.codespell-ignore-words b/.codespell-ignore-words index 20008eda4..909e4a13e 100644 --- a/.codespell-ignore-words +++ b/.codespell-ignore-words @@ -1,4 +1,7 @@ -latice +Blocs +BOUN +inout parms pres -ptd +rime +wth diff --git a/Docs/sphinx_doc/Applications_Requirements.rst b/Docs/sphinx_doc/Applications_Requirements.rst index eb599267d..5ceff1b4c 100644 --- a/Docs/sphinx_doc/Applications_Requirements.rst +++ b/Docs/sphinx_doc/Applications_Requirements.rst @@ -232,7 +232,7 @@ These standards should be articulated at a sufficient level of detail to guide c Implementation Details for Specific Code Components =================================================== -This section provides detailed information on core elements of the ERF code that enable its required functionality. Details of both the formulation of methods and their implementation and useability within the code will be added to the existing higher-level descriptions as the project progresses. +This section provides detailed information on core elements of the ERF code that enable its required functionality. Details of both the formulation of methods and their implementation and usability within the code will be added to the existing higher-level descriptions as the project progresses. 1. Governing Equation Set ------------------------- diff --git a/Docs/sphinx_doc/BoundaryConditions.rst b/Docs/sphinx_doc/BoundaryConditions.rst index 1bfdea0ae..5df9d59ae 100644 --- a/Docs/sphinx_doc/BoundaryConditions.rst +++ b/Docs/sphinx_doc/BoundaryConditions.rst @@ -81,7 +81,7 @@ As an example, xlo.theta = 300. xlo.scalar = 2. -sets the boundary condtion type at the low x face to be an inflow with xlo.type = “Inflow”. +sets the boundary condition type at the low x face to be an inflow with xlo.type = “Inflow”. xlo.velocity = 1. 0. 0. sets all three componentns the inflow velocity, xlo.density = 1. sets the inflow density, @@ -133,7 +133,7 @@ Couette regression test example, in which we specify We also note that in the case of a "slipwall" boundary condition in a simulation with non-zero viscosity specified, the "foextrap" boundary condition enforces zero strain at the wall. -The keywork "MOST" is an ERF-specific boundary type; see :ref:`sec:MOST` for more information. +The keyword "MOST" is an ERF-specific boundary type; see :ref:`sec:MOST` for more information. It is important to note that external Dirichlet boundary data should be specified as the value on the face of the cell bounding the domain, even for cell-centered @@ -192,7 +192,7 @@ and :math:`n` is the minimum number of grid points from a lateral boundary. Sponge zone domain BCs ---------------------- -ERF provides the capability to apply sponge zones at the boundaries to prevent spurious reflections that otherwise occur at the domain boundaries if standard extrapolation boundary condition is used. The sponge zone is implemented as a source term in the governing equations, which are active in a volumteric region at the boundaries that is specified by the user in the inputs file. Currently the target condition to which the sponge zones should be forced towards is to be specifed by the user in the inputs file. +ERF provides the capability to apply sponge zones at the boundaries to prevent spurious reflections that otherwise occur at the domain boundaries if standard extrapolation boundary condition is used. The sponge zone is implemented as a source term in the governing equations, which are active in a volumteric region at the boundaries that is specified by the user in the inputs file. Currently the target condition to which the sponge zones should be forced towards is to be specified by the user in the inputs file. .. math:: diff --git a/Docs/sphinx_doc/Discretizations.rst b/Docs/sphinx_doc/Discretizations.rst index 39325c4a4..04ff52071 100644 --- a/Docs/sphinx_doc/Discretizations.rst +++ b/Docs/sphinx_doc/Discretizations.rst @@ -474,12 +474,12 @@ The goal is to compute eddy viscosity at the *cell centers* and interpolated the .. math:: \begin{matrix} - S_{12} = & \frac{1}{4}\left\lbrack S_{12i,j - \frac{1}{2}} + S_{12i,j + \frac{1}{2}} + S_{12i + 1,j - \frac{1}{2}} + S_{12i + 1,j + \frac{1}{2}} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrouding the cell}\end{smallmatrix} \\ - S_{21} = & \frac{1}{4}\left\lbrack S_{21i - \frac{1}{2},j} + S_{21i + \frac{1}{2},j} + S_{21i - \frac{1}{2},j + 1} + S_{21i + \frac{1}{2},j + 1} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrouding the cell}\end{smallmatrix} \\ - S_{13} = & \frac{1}{4}\left\lbrack S_{13i,k - \frac{1}{2}} + S_{13i,k + \frac{1}{2}} + S_{13i + 1,k - \frac{1}{2}} + S_{13i + 1,k + \frac{1}{2}} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrouding the cell}\end{smallmatrix} \\ - S_{31} = & \frac{1}{4}\left\lbrack S_{31i - \frac{1}{2},k} + S_{31i + \frac{1}{2},k} + S_{31i - \frac{1}{2},k + 1} + S_{31i + \frac{1}{2},k + 1} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrouding the cell}\end{smallmatrix} \\ - S_{23} = & \frac{1}{4}\left\lbrack S_{23j,k - \frac{1}{2}} + S_{23j,k + \frac{1}{2}} + S_{23j + 1,k - \frac{1}{2}} + S_{23j + 1,k + \frac{1}{2}} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrouding the cell}\end{smallmatrix} \\ - S_{32} = & \frac{1}{4}\left\lbrack S_{32j - \frac{1}{2},k} + S_{32j + \frac{1}{2},k} + S_{32j - \frac{1}{2},k + 1} + S_{32j + \frac{1}{2},k + 1} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrouding the cell}\end{smallmatrix} + S_{12} = & \frac{1}{4}\left\lbrack S_{12i,j - \frac{1}{2}} + S_{12i,j + \frac{1}{2}} + S_{12i + 1,j - \frac{1}{2}} + S_{12i + 1,j + \frac{1}{2}} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrounding the cell}\end{smallmatrix} \\ + S_{21} = & \frac{1}{4}\left\lbrack S_{21i - \frac{1}{2},j} + S_{21i + \frac{1}{2},j} + S_{21i - \frac{1}{2},j + 1} + S_{21i + \frac{1}{2},j + 1} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrounding the cell}\end{smallmatrix} \\ + S_{13} = & \frac{1}{4}\left\lbrack S_{13i,k - \frac{1}{2}} + S_{13i,k + \frac{1}{2}} + S_{13i + 1,k - \frac{1}{2}} + S_{13i + 1,k + \frac{1}{2}} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrounding the cell}\end{smallmatrix} \\ + S_{31} = & \frac{1}{4}\left\lbrack S_{31i - \frac{1}{2},k} + S_{31i + \frac{1}{2},k} + S_{31i - \frac{1}{2},k + 1} + S_{31i + \frac{1}{2},k + 1} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrounding the cell}\end{smallmatrix} \\ + S_{23} = & \frac{1}{4}\left\lbrack S_{23j,k - \frac{1}{2}} + S_{23j,k + \frac{1}{2}} + S_{23j + 1,k - \frac{1}{2}} + S_{23j + 1,k + \frac{1}{2}} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrounding the cell}\end{smallmatrix} \\ + S_{32} = & \frac{1}{4}\left\lbrack S_{32j - \frac{1}{2},k} + S_{32j + \frac{1}{2},k} + S_{32j - \frac{1}{2},k + 1} + S_{32j + \frac{1}{2},k + 1} \right\rbrack = \begin{smallmatrix} \text{Average of the 4 edges} \\ \text{surrounding the cell}\end{smallmatrix} \end{matrix} Note that: diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index e84316607..230b1fcaf 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -457,7 +457,7 @@ Notes * | If **erf.fixed_mri_dt_ratio** is specified but is not an even positive integer * | If **erf.fixed_dt** and **erf.fast_fixed_dt** are specified and the ratio of **fixed_dt** to **fast_fixed_dt** is not an even positive integer - * | If **erf.fixed_dt** and **erf.fast_fixed_dt** and **erf.fixed_mri_dt_ratio** are all specified but are inconsitent + * | If **erf.fixed_dt** and **erf.fast_fixed_dt** and **erf.fixed_mri_dt_ratio** are all specified but are inconsistent * | Once the slow timestep is set and the inputs are allowed per the above criteria, the fast timestep is computed in one of several ways: @@ -568,7 +568,7 @@ List of Parameters +-------------------------------+------------------+----------------+----------------+ | **erf.profile_int** | Interval (number)| Integer | -1 | | | of steps between | | | -| | ouputs | | | +| | outputs | | | +-------------------------------+------------------+----------------+----------------+ | **erf.interp_profiles_to_cc** | Interpolate all | Boolean | true | | | outputs to cell | | | diff --git a/Docs/sphinx_doc/MOST.rst b/Docs/sphinx_doc/MOST.rst index cbd33578d..58312106f 100644 --- a/Docs/sphinx_doc/MOST.rst +++ b/Docs/sphinx_doc/MOST.rst @@ -9,7 +9,7 @@ MOST Boundaries Monin-Obukhov similarity theory (MOST) is used to describe the atmospheric surface layer (ASL), the lowest part of the atmospheric boundary layer. The implementation of MOST in ERF follows that in `AMR-Wind `_, which is based on the surface layer profiles presented in `P. van der Laan, et al., Wind Energy, 2017 `_ and `D. Etling, "Modeling the vertical ABL structure", 1999 `_. -MOST theory assumes that the ASL is in a steady state and horizontally homogenous, and kinematic fluxes due to turbulent transport (:math:`\overline{u^{'}w^{'}}`, :math:`\overline{v^{'}w^{'}}`, and :math:`\overline{\theta^{'}w^{'}}`) are constant with height. +MOST theory assumes that the ASL is in a steady state and horizontally homogeneous, and kinematic fluxes due to turbulent transport (:math:`\overline{u^{'}w^{'}}`, :math:`\overline{v^{'}w^{'}}`, and :math:`\overline{\theta^{'}w^{'}}`) are constant with height. :math:`\Phi_m` and :math:`\Phi_h` are the nondimensional wind shear and temperature gradient, respectively, which are assumed to follow universal similarity laws based on dimensional arguments. With these assumptions, the MOST theory can be written as: @@ -120,7 +120,7 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in of the velocity terms from the form of the equations presented in Moeng to match the form implemented in AMR-Wind. -#. These local flux values are used to populate values in the ghost cells that will lead to appropiate fluxes, assuming the fluxes are computed from the turbulent transport coefficients (in the vertical direction, if applicable) :math:`K_{m,v}` and :math:`K_{\theta,v}` as follows: +#. These local flux values are used to populate values in the ghost cells that will lead to appropriate fluxes, assuming the fluxes are computed from the turbulent transport coefficients (in the vertical direction, if applicable) :math:`K_{m,v}` and :math:`K_{\theta,v}` as follows: .. math:: @@ -143,7 +143,7 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in (\rho \theta)_{z} = \frac{(\rho \theta)_{i,j,1} - (\rho \theta)_{i,j,0}}{\Delta z} (\rho \theta)_{i,j,-n} = (\rho \theta)_{i,j,0} - (\rho \theta)_{z} n \Delta z . - Finally, it must be noted that complex terrain will modify the surface normal and tangent vectors. Consequently, the MOST implentation with terrain will require local vector rotations. While the ERF dycore accounts for + Finally, it must be noted that complex terrain will modify the surface normal and tangent vectors. Consequently, the MOST implementation with terrain will require local vector rotations. While the ERF dycore accounts for terrain metrics when computing fluxes (e.g. for advection, diffusion, etc.), the impact of terrain metrics on MOST is still a work in progress. Therefore, running with terrain (``erf.use_terrain = true``) and with MOST (``zlo.type = "Most"``) should be cautioned. diff --git a/Docs/sphinx_doc/Visualization.rst b/Docs/sphinx_doc/Visualization.rst index 6b2842fe2..8fb46bc9c 100644 --- a/Docs/sphinx_doc/Visualization.rst +++ b/Docs/sphinx_doc/Visualization.rst @@ -30,11 +30,11 @@ To open a plotfile #. If you have run the ERF executable with terrain, then the mapped grid information will be stored as nodal data. Choose the "point data" called "nu", then click on "Warp by Vector" - which can be found via Filters-->Alphabetical. This wil then plot data onto the mapped grid + which can be found via Filters-->Alphabetical. This will then plot data onto the mapped grid locations. #. Under the "Cell Arrays" field, select a variable (e.g., "x_velocity") and click - "Apply". Note that the default number of refinement levels loaded and vizualized is 1. + "Apply". Note that the default number of refinement levels loaded and visualized is 1. Change to the required number of AMR level before clicking "Apply". #. For "Representation" select "Surface". diff --git a/Docs/sphinx_doc/coc.rst b/Docs/sphinx_doc/coc.rst index 749d1bc5b..ce3d7e2ba 100644 --- a/Docs/sphinx_doc/coc.rst +++ b/Docs/sphinx_doc/coc.rst @@ -6,7 +6,7 @@ Code of Conduct Our Pledge ---------- -In the interest of fostering an open and welcoming environment, we as contributors and maintainers pledge to making participation in our project and our community a harassment-free experience for everyone, regardless of age, body size, disability, ethnicity, sex characteristics, gender identity and expression, level of experience, education, socio-economic status, nationality, personal appearance, race, religion, or sexual identity and orientation. +In the interest of fostering an open and welcoming environment, we as contributors and maintainers pledge to making participation in our project and our community a harassment-free experience for everyone, regardless of age, body size, disability, ethnicity, sex characteristics, gender identity and expression, level of experience, education, socioeconomic status, nationality, personal appearance, race, religion, or sexual identity and orientation. Our Standards ------------- diff --git a/Docs/sphinx_doc/testing.rst b/Docs/sphinx_doc/testing.rst index 28198ee2e..e19fc9d3a 100644 --- a/Docs/sphinx_doc/testing.rst +++ b/Docs/sphinx_doc/testing.rst @@ -3,7 +3,7 @@ Testing and Verification ------------------------ -Testing and verfication of ERF can be performed using CTest, which is included in the CMake build system. If one builds ERF with CMake, the testing suite, and the verification suite, can be enabled during the CMake configure step. +Testing and verification of ERF can be performed using CTest, which is included in the CMake build system. If one builds ERF with CMake, the testing suite, and the verification suite, can be enabled during the CMake configure step. An example ``cmake`` configure command performed in the ``Build`` directory in ERF is shown below with options relevant to the testing suite: diff --git a/Docs/sphinx_doc/theory/DNSvsLES.rst b/Docs/sphinx_doc/theory/DNSvsLES.rst index 9273a9bad..1b5f09470 100644 --- a/Docs/sphinx_doc/theory/DNSvsLES.rst +++ b/Docs/sphinx_doc/theory/DNSvsLES.rst @@ -73,7 +73,7 @@ respectively (Martin et al., Theoret. Comput. Fluid Dynamics (2000)). compressible turbulence and scalar transport", PoF (1991); Martin et al., Subgrid-scale models for compressible large-eddy simulations", Theoret. Comput. Fluid Dynamics (2000). -When substituted back into the filtered equtions, the gradient transport LES models take exactly the same form as the +When substituted back into the filtered equations, the gradient transport LES models take exactly the same form as the molecular transport terms, but with the constant molecular transport coefficients replaced by turbulent equivalents (e.g. :math:`\mu` becomes the turbulent viscosity, :math:`\mu_{t}`). Therefore, when the code is run in LES mode, the :ref:`equation set` remains the same, diff --git a/Docs/sphinx_doc/theory/PBLschemes.rst b/Docs/sphinx_doc/theory/PBLschemes.rst index f5cbe28ae..100ffd2c8 100644 --- a/Docs/sphinx_doc/theory/PBLschemes.rst +++ b/Docs/sphinx_doc/theory/PBLschemes.rst @@ -63,7 +63,7 @@ coefficients are then computed: K_{m,v} = l q S_m, K_{q,v} = 3 l q S_m, K_{\theta, v} = l q S_\theta where :math:`S_m` and :math:`S_\theta` are stability parameters thaat -account for bouyancy effects. These +account for buoyancy effects. These coefficients are then applied in evaluating the vertical component of turbulent fluxes in a similar manner as is described for the :ref:`Smagorinsky LES model`. Computation of the stability parameters @@ -172,12 +172,12 @@ where - :math:`U = \sqrt{u^2 + v^2}` is the horizontal wind speed, - :math:`\theta_s = \theta_{ma} + \theta_T` is the virtual temperature near the surface, - :math:`\theta_T = a\frac{\overline{\left(w'\theta_m' \right)_0}}{w_{s0}}` is the excess virtual temperature near the surface, -- :math:`a` is a constant taken to be 6.8 per HND06 (matching the :math:`b` constant that appears elsehwere in the YSU model) +- :math:`a` is a constant taken to be 6.8 per HND06 (matching the :math:`b` constant that appears elsewhere in the YSU model) - :math:`\overline{\left(w'\theta_m' \right)_0}` is the surface virtual heat flux (determined from the MOST surface layer model), - :math:`w_{s}(z) = \left(u_*^3 + 8 k w_{*b}^3z/h \right)^{1/3}` is a representative velocity scale in the mixed layer, with :math:`w_{s0} = w_s(h/2)` (note this equation matches the WRF implementation and description in H10, but differs from HND06, where :math:`\phi_m` appears in place of the constant 8), - :math:`u_*` is the surface frictional velocity scale determined from the MOST surface layer model, - :math:`k = 0.4` is the von Karman constant -- :math:`w_{*b} = \left[ g/\theta_{ma} \overline{\left(w'\theta_m' \right)_0} h \right]^{1/3}` for :math:`\overline{\left(w'\theta_m' \right)_0} > 0`, :math:`w_{*b} = 0` otherwise, is a convective velcoity scale for moist air +- :math:`w_{*b} = \left[ g/\theta_{ma} \overline{\left(w'\theta_m' \right)_0} h \right]^{1/3}` for :math:`\overline{\left(w'\theta_m' \right)_0} > 0`, :math:`w_{*b} = 0` otherwise, is a convective velocity scale for moist air In practice, an approximate value of :math:`h` is determined through a two-step process. First, :math:`\theta_T` is set to be zero and a provisional value of :math:`h` is estimated. Then this provisional value of :math:`h` is used to compute :math:`\theta_T`, @@ -228,7 +228,7 @@ The following references have informed the implementation of the YSU model in ER - [NCHR03] `Noh, Cheon, Hong, and Raasch, Boundary-Layer Meteorology, 2003 `_: Entrainment effects added to TM86 -- [HP96] `Hong and Pan, Monthly Weather Review, 1996 `_: Largely an implementation and evluation of TM86 +- [HP96] `Hong and Pan, Monthly Weather Review, 1996 `_: Largely an implementation and evaluation of TM86 - [TM86] `Troen and Mahrt, Boundary-Layer Meteorology, 1986 `_: Initial incorporation of nonlocal counter-graident term in vertical diffusion model diff --git a/Docs/sphinx_doc/theory/WindFarmModels.rst b/Docs/sphinx_doc/theory/WindFarmModels.rst index ed709e697..ee679529d 100644 --- a/Docs/sphinx_doc/theory/WindFarmModels.rst +++ b/Docs/sphinx_doc/theory/WindFarmModels.rst @@ -157,7 +157,7 @@ The following are the inputs required for simulations with Fitch, EWP models. erf.windfarm_spec_table = "wind-turbine_1WT.tbl" 1. ``erf.windfarm_type`` has to be one of the supported models - ``Fitch``, ``EWP``. -2. ``erf.windfarm_loc_type`` is a variable to specify how the wind turbine locations in the wind farm is specified. If using the latitude and longitude of the turbine location, this has to be ``lat_lon`` or if using x and y co-ordinates to specify the turbine locations, this input is ``x_y``. +2. ``erf.windfarm_loc_type`` is a variable to specify how the wind turbine locations in the wind farm is specified. If using the latitude and longitude of the turbine location, this has to be ``lat_lon`` or if using x and y coordinates to specify the turbine locations, this input is ``x_y``. - If using ``lat_lon`` format, ``erf.latitude_lo`` and ``erf.longitude_lo`` specifies the latitude and longitude of the lower bottom corner of the domain box. ie. if the domain box is specified as @@ -185,7 +185,7 @@ The following are the inputs required for simulations with Fitch, EWP models. 36.0171171171 -99.0168 1 35.7828828829 -98.9705333333 1 - - For the x-y format, an example is as below. Each line specifies the x and y co-ordinates of the wind turbine location in metres + - For the x-y format, an example is as below. Each line specifies the x and y coordinates of the wind turbine location in metres .. code-block:: cpp diff --git a/Exec/DevTests/Bomex/GNUmakefile b/Exec/DevTests/Bomex/GNUmakefile index f9b76e850..ef670d8f7 100644 --- a/Exec/DevTests/Bomex/GNUmakefile +++ b/Exec/DevTests/Bomex/GNUmakefile @@ -4,7 +4,7 @@ PRECISION = DOUBLE # Profiling PROFILE = FALSE -TINY_PROFILE = TRUE +TINY_PROFILE = FALSE COMM_PROFILE = FALSE TRACE_PROFILE = FALSE MEM_PROFILE = FALSE diff --git a/Exec/DevTests/Bomex/input_Kessler b/Exec/DevTests/Bomex/input_Kessler index 8e5af436f..cc018cbdf 100644 --- a/Exec/DevTests/Bomex/input_Kessler +++ b/Exec/DevTests/Bomex/input_Kessler @@ -1,5 +1,9 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -stop_time = 43200 +#stop_time = 21600 # 6 hours +#stop_time = 3600 # 1 hour +#stop_time = 7200 # 2 hours + +stop_time = 4800 # 80 minutes amrex.fpe_trap_invalid = 1 @@ -25,26 +29,27 @@ zhi.type = "SlipWall" #zhi.theta_grad = 0.00365 # TIME STEP CONTROL -erf.fixed_dt = 0.5 # fixed time step depending on grid resolution +erf.no_substepping = 1 +erf.fixed_dt = 0.075 # fixed time step depending on grid resolution # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass -erf.v = 1 # verbosity in ERF.cpp +erf.v = 0 # verbosity in ERF.cpp amr.v = 1 # verbosity in Amr.cpp -erf.data_log = "surf" "mean" "flux" "subgrid" -erf.profile_int = 120 +erf.data_log = "surf" "mean" "flux" # "subgrid" +erf.profile_int = 800 # (every minute with dt = 0.075) # REFINEMENT / REGRIDDING amr.max_level = 0 # maximum level number allowed # CHECKPOINT FILES -erf.check_file = chk # root name of checkpoint file -erf.check_int = 600 # number of timesteps between checkpoints +erf.check_file = chk # root name of checkpoint file +erf.check_int = 8000 # number of timesteps between checkpoints # PLOTFILES -erf.plot_file_1 = plt # prefix of plotfile name -erf.plot_int_1 = 120 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta x_velocity y_velocity z_velocity pressure temp theta qt qp qv qc qi +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 8000 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta x_velocity y_velocity z_velocity pressure temp theta qt qp qv qc qsat # SOLVER CHOICE erf.alpha_T = 0.0 @@ -58,7 +63,8 @@ erf.dryscal_vert_adv_type = WENOZ5 erf.moistscal_horiz_adv_type = WENOZ5 erf.moistscal_vert_adv_type = WENOZ5 -erf.moisture_model = "Kessler" +erf.moisture_model = "Kessler_NoRain" +erf.buoyancy_type = 1 erf.molec_diff_type = "None" @@ -70,8 +76,8 @@ erf.Cs = 0.17 #erf.sigma_k = 1.0 #erf.Ce = 0.1 -erf.Pr_t = 0.333333 -erf.Sc_t = 0.333333 +erf.Pr_t = 0.33333333333333 +erf.Sc_t = 0.33333333333333 erf.init_type = "input_sounding" erf.init_sounding_ideal = true @@ -80,7 +86,7 @@ erf.add_custom_rhotheta_forcing = true erf.add_custom_moisture_forcing = true erf.add_custom_geostrophic_profile = true erf.add_custom_w_subsidence = true -erf.custom_forcing_uses_primitive_vars = false +erf.custom_forcing_uses_primitive_vars = true # Higher values of perturbations lead to instability # Instability seems to be coming from BC diff --git a/Exec/DevTests/Bomex/prob.cpp b/Exec/DevTests/Bomex/prob.cpp index 2908afba5..e6f2d286f 100644 --- a/Exec/DevTests/Bomex/prob.cpp +++ b/Exec/DevTests/Bomex/prob.cpp @@ -66,7 +66,7 @@ Problem::init_custom_pert ( const Box& xbx, const Box& ybx, const Box& zbx, - Array4 const& /*state*/, + Array4 const& state, Array4 const& state_pert, Array4 const& x_vel_pert, Array4 const& y_vel_pert, @@ -83,6 +83,8 @@ Problem::init_custom_pert ( { const bool use_moisture = (sc.moisture_type != MoistureType::None); + const Real rdOcp = sc.rdOcp; + ParallelForRNG(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k, const RandomEngine& engine) noexcept { // Geometry @@ -99,11 +101,32 @@ Problem::init_custom_pert ( const Real zc = 0.5 * (prob_lo[2] + prob_hi[2]); const Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); - - // Add temperature perturbations + // + // Add temperature perturbations -- we want to keep pressure constant + // so these effectively end up as density perturbations + // if ((z <= parms.pert_ref_height) && (parms.T_0_Pert_Mag != 0.0)) { + + Real rhotheta = state(i,j,k,RhoTheta_comp); + Real rho = state(i,j,k,Rho_comp); + Real qv = state(i,j,k,RhoQ1_comp) / rho; + Real Told = getTgivenRandRTh(rho,rhotheta,qv); + Real P = getPgivenRTh(rhotheta,qv); + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 - state_pert(i, j, k, RhoTheta_comp) = (rand_double*2.0 - 1.0)*parms.T_0_Pert_Mag; //*state_pert(i, j, k, Rho_comp); + Real Tpert = (rand_double*2.0 - 1.0)*parms.T_0_Pert_Mag; + Real Tnew = Told + Tpert; + + Real theta_new = getThgivenPandT(Tnew,P,rdOcp); + Real rhonew = getRhogivenThetaPress(theta_new,P,rdOcp,qv); + state_pert(i, j, k, Rho_comp) = rhonew - rho; + + // Note we do not perturb this + state_pert(i, j, k, RhoTheta_comp) = 0.0; + + // Instead of perturbing (rho theta) we perturb T and hold (rho theta) fixed, + // which ends up being stored as a perturbation in rho + // state_pert(i, j, k, RhoTheta_comp) = (rand_double*2.0 - 1.0)*parms.T_0_Pert_Mag; } // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain @@ -119,9 +142,17 @@ Problem::init_custom_pert ( if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; state_pert(i, j, k, RhoQ2_comp) = 0.0; - if ((z <= parms.pert_ref_height) && (parms.qv_0_Pert_Mag != 0.0)) { + if ((z <= parms.pert_ref_height) && (parms.qv_0_Pert_Mag != 0.0)) + { + Real rhoold = state(i,j,k,Rho_comp); + Real rhonew = rhoold + state_pert(i,j,k,Rho_comp); + + Real qvold = state(i,j,k,RhoQ1_comp) / rhoold; + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 - state_pert(i, j, k, RhoQ1_comp) = (rand_double*2.0 - 1.0)*parms.qv_0_Pert_Mag; //*state_pert(i, j, k, Rho_comp); + Real qvnew = qvold + (rand_double*2.0 - 1.0)*parms.qv_0_Pert_Mag; + + state_pert(i, j, k, RhoQ1_comp) = rhonew * qvnew - rhoold * qvold; } } }); @@ -312,8 +343,8 @@ Problem::update_w_subsidence (const Real& /*time*/, Real z_0 = 0.0; // (z_phys_cc) ? zlevels[0] : prob_lo[2] + 0.5 * dx[2]; Real slope1 = parms.wbar_sub_max / (parms.wbar_cutoff_max - z_0); Real slope2 = -parms.wbar_sub_max / (parms.wbar_cutoff_min - parms.wbar_cutoff_max); - wbar[0] = 0.0; - for (int k = 1; k <= khi; k++) { + // wbar[0] = 0.0; + for (int k = 0; k <= khi; k++) { const Real z_cc = (z_phys_cc) ? zlevels[k] : prob_lo[2] + (k+0.5)* dx[2]; if (z_cc <= parms.wbar_cutoff_max) { wbar[k] = slope1 * (z_cc - z_0); diff --git a/Exec/DevTests/MultiBlock/inputs_amrex b/Exec/DevTests/MultiBlock/inputs_amrex index 60f3228e9..453573c24 100644 --- a/Exec/DevTests/MultiBlock/inputs_amrex +++ b/Exec/DevTests/MultiBlock/inputs_amrex @@ -13,7 +13,7 @@ amr.v = 1 # verbosity in Amr.cpp # REFINEMENT / REGRIDDING amr.max_level = 0 # maximum level number allowed -# TRAP OVER/UNDEFLOW +# TRAP OVER/UNDERFLOW amrex.fpe_trap_invalid = 1 #amrex.fpe_trap_zero = 1 #amrex.fpe_trap_overflow = 1 diff --git a/Exec/RegTests/DensityCurrent/README b/Exec/RegTests/DensityCurrent/README index 97be5cd2c..f03f8eec2 100644 --- a/Exec/RegTests/DensityCurrent/README +++ b/Exec/RegTests/DensityCurrent/README @@ -7,7 +7,7 @@ using a number of different codes can be found in the Straka 1993 paper. See: > Straka, J.M., Wilhelmson, R.B., Wicker, L.J., Anderson, J.R. and Droegemeier, K.K. (1993), Numerical solutions of a non-linear density current: A benchmark -solution and comparisons. Int. J. Numer. Meth. Fluids, 17: 1-22. +solution and comparisons. Int. J. Numerical Methods Fluids, 17: 1-22. https://doi.org/10.1002/fld.1650170103 Compare with `WRF/test/em_grav2d_x` diff --git a/Source/Advection/AdvectionSrcForMom.cpp b/Source/Advection/AdvectionSrcForMom.cpp index 6c2601b7a..39589d4db 100644 --- a/Source/Advection/AdvectionSrcForMom.cpp +++ b/Source/Advection/AdvectionSrcForMom.cpp @@ -9,7 +9,7 @@ using namespace amrex; /** * Function for computing the advective tendency for the momentum equations * This routine has explicit expressions for all cases (terrain or not) when - * the horizontal and vertial spatial orders are <= 2, and calls more specialized + * the horizontal and vertical spatial orders are <= 2, and calls more specialized * functions when either (or both) spatial order(s) is greater than 2. * * @param[in] bxx box over which the x-momentum is updated diff --git a/Source/Advection/AdvectionSrcForState.cpp b/Source/Advection/AdvectionSrcForState.cpp index 59a01fe78..a7379af6e 100644 --- a/Source/Advection/AdvectionSrcForState.cpp +++ b/Source/Advection/AdvectionSrcForState.cpp @@ -8,7 +8,7 @@ using namespace amrex; /** * Function for computing the advective tendency for the update equations for rho and (rho theta) * This routine has explicit expressions for all cases (terrain or not) when - * the horizontal and vertial spatial orders are <= 2, and calls more specialized + * the horizontal and vertical spatial orders are <= 2, and calls more specialized * functions when either (or both) spatial order(s) is greater than 2. * * @param[in] bx box over which the scalars are updated @@ -96,7 +96,7 @@ AdvectionSrcForRho (const Box& bx, /** * Function for computing the advective tendency for the update equations for all scalars other than rho and (rho theta) * This routine has explicit expressions for all cases (terrain or not) when - * the horizontal and vertial spatial orders are <= 2, and calls more specialized + * the horizontal and vertical spatial orders are <= 2, and calls more specialized * functions when either (or both) spatial order(s) is greater than 2. * * @param[in] bx box over which the scalars are updated if no external boundary conditions @@ -105,7 +105,7 @@ AdvectionSrcForRho (const Box& bx, * @param[in] avg_xmom x-component of time-averaged momentum defined in this routine * @param[in] avg_ymom y-component of time-averaged momentum defined in this routine * @param[in] avg_zmom z-component of time-averaged momentum defined in this routine - * @param[in] cell_prim primtive form of scalar variales, here only potential temperature theta + * @param[in] cell_prim primitive form of scalar variables, here only potential temperature theta * @param[out] advectionSrc tendency for the scalar update equation * @param[in] detJ Jacobian of the metric transformation (= 1 if use_terrain is false) * @param[in] cellSizeInv inverse of the mesh spacing diff --git a/Source/BoundaryConditions/BoundaryConditions_realbdy.cpp b/Source/BoundaryConditions/BoundaryConditions_realbdy.cpp index c4eb07528..b548bda72 100644 --- a/Source/BoundaryConditions/BoundaryConditions_realbdy.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_realbdy.cpp @@ -63,7 +63,7 @@ ERF::fill_from_realbdy (const Vector& mfs, const auto& dom_lo = lbound(domain); const auto& dom_hi = ubound(domain); - // Offset only applys to cons (we may fill a subset of these vars) + // Offset only applies to cons (we may fill a subset of these vars) int offset = (var_idx == Vars::cons) ? icomp_cons : 0; // Loop over each component diff --git a/Source/BoundaryConditions/ERF_FillPatcher.cpp b/Source/BoundaryConditions/ERF_FillPatcher.cpp index 4be4106dd..0b80bcb26 100644 --- a/Source/BoundaryConditions/ERF_FillPatcher.cpp +++ b/Source/BoundaryConditions/ERF_FillPatcher.cpp @@ -9,10 +9,10 @@ using namespace amrex; * * @param[in] fba BoxArray of data to be filled at fine level * @param[in] fdm DistributionMapping of data to be filled at fine level - * @param[in] fgeom container of geometry infomation at fine level + * @param[in] fgeom container of geometry information at fine level * @param[in] cba BoxArray of data to be filled at coarse level * @param[in] cdm DistributionMapping of data to be filled at coarse level - * @param[in] cgeom container of geometry infomation at coarse level + * @param[in] cgeom container of geometry information at coarse level * @param[in] nghost number of ghost cells to be filled * @param[in] ncomp number of components to be filled * @param[in] interp interpolation operator to be used @@ -46,10 +46,10 @@ ERFFillPatcher::ERFFillPatcher (BoxArray const& fba, DistributionMapping const& * * @param[in] fba BoxArray of data to be filled at fine level * @param[in] fdm DistributionMapping of data to be filled at fine level - * @param[in] fgeom container of geometry infomation at fine level + * @param[in] fgeom container of geometry information at fine level * @param[in] cba BoxArray of data to be filled at coarse level * @param[in] cdm DistributionMapping of data to be filled at coarse level - * @param[in] cgeom container of geometry infomation at coarse level + * @param[in] cgeom container of geometry information at coarse level * @param[in] nghost number of ghost cells to be filled * @param[in] ncomp number of components to be filled * @param[in] interp interpolation operator to be used @@ -65,7 +65,7 @@ void ERFFillPatcher::Define (BoxArray const& fba, DistributionMapping const& fdm AMREX_ALWAYS_ASSERT(nghost_set <= 0); AMREX_ALWAYS_ASSERT(nghost <= nghost_set); - // Set data memebers + // Set data members m_fba = fba; m_cba = cba; m_fdm = fdm; m_cdm = cdm; m_fgeom = fgeom; m_cgeom = cgeom; diff --git a/Source/BoundaryConditions/MOSTAverage.H b/Source/BoundaryConditions/MOSTAverage.H index b78be715a..1965b85fc 100644 --- a/Source/BoundaryConditions/MOSTAverage.H +++ b/Source/BoundaryConditions/MOSTAverage.H @@ -163,7 +163,7 @@ protected: //-------------------------------------------- const amrex::Vector m_geom; // Geometry at each level amrex::Vector> m_fields; // Ptr to fields to be averaged - amrex::Vector m_z_phys_nd; // Ptr to terrain heigh coords + amrex::Vector m_z_phys_nd; // Ptr to terrain height coords // General vars for multiple or all policies //-------------------------------------------- diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index b15bd2410..684d7ac90 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -1004,7 +1004,7 @@ struct moeng_flux // Note: The surface mean shear stress is decomposed into tau_xz by // multiplying the modeled shear stress (rho*ustar^2) with // a factor of umean/wsp_mean for directionality; this factor - // modifies the demoninator from what is in Moeng 1984. + // modifies the denominator from what is in Moeng 1984. amrex::Real wsp = sqrt(velx*velx+vely*vely); amrex::Real num1 = wsp * umean; amrex::Real num2 = wsp_mean * (velx-umean); @@ -1076,7 +1076,7 @@ struct moeng_flux // Note: The surface mean shear stress is decomposed into tau_yz by // multiplying the modeled shear stress (rho*ustar^2) with // a factor of vmean/wsp_mean for directionality; this factor - // modifies the demoninator from what is in Moeng 1984. + // modifies the denominator from what is in Moeng 1984. amrex::Real wsp = sqrt(velx*velx+vely*vely); amrex::Real num1 = wsp * vmean; amrex::Real num2 = wsp_mean * (vely-vmean); diff --git a/Source/DataStructs/TurbPertStruct.H b/Source/DataStructs/TurbPertStruct.H index 5d5021248..40224515c 100644 --- a/Source/DataStructs/TurbPertStruct.H +++ b/Source/DataStructs/TurbPertStruct.H @@ -83,7 +83,7 @@ struct TurbulentPerturbation { // Perturbation update frequency check. // This function will trigger calc_tpi_meanMag() and calc_tpi_amp(), when - // the local elasped time is greater than the perturbation frequency + // the local elapsed time is greater than the perturbation frequency void calc_tpi_update (const int lev, const amrex::Real dt, amrex::MultiFab& mf_xvel, @@ -117,7 +117,7 @@ struct TurbulentPerturbation { #define USE_VOLUME_AVERAGE // Perturbation box mean velocity magnitude calculation - // This is pulled into the strucutre to also utilize during runtime + // This is pulled into the structure to also utilize during runtime void calc_tpi_meanMag (const int boxIdx, const int lev, amrex::MultiFab& mf_xvel, @@ -240,7 +240,7 @@ struct TurbulentPerturbation { //#define ECKERT_FORMULATION #define BULK_RICHARDSON_FORMULATION - // Perturbation amplitude calcluation + // Perturbation amplitude calculation void calc_tpi_amp (const int boxIdx) { amrex::Real g = CONST_GRAV; diff --git a/Source/DataStructs/TurbStruct.H b/Source/DataStructs/TurbStruct.H index 31a16e391..006246508 100644 --- a/Source/DataStructs/TurbStruct.H +++ b/Source/DataStructs/TurbStruct.H @@ -104,8 +104,8 @@ struct TurbChoice { pp.query("Sc_t", Sc_t); // Compute relevant forms of diffusion parameters - Pr_t_inv = 1.0 / Pr_t; - Sc_t_inv = 1.0 / Sc_t; + Pr_t_inv = amrex::Real(1.0) / Pr_t; + Sc_t_inv = amrex::Real(1.0) / Sc_t; pp.query("Ce" , Ce); pp.query("Ce_wall" , Ce_wall); @@ -189,8 +189,8 @@ struct TurbChoice { pp.query("Sc_t",Sc_t, lev); // Compute relevant forms of diffusion parameters - Pr_t_inv = 1.0 / Pr_t; - Sc_t_inv = 1.0 / Sc_t; + Pr_t_inv = amrex::Real(1.0) / Pr_t; + Sc_t_inv = amrex::Real(1.0) / Sc_t; pp.query("Ce" , Ce, lev); pp.query("Ce_wall" , Ce_wall, lev); @@ -259,8 +259,8 @@ struct TurbChoice { // Smagorinsky CI coefficient amrex::Real CI = 0.0; // Smagorinsky Turbulent Prandtl Number - amrex::Real Pr_t = 1. / 3.; - amrex::Real Pr_t_inv = 3.0; + amrex::Real Pr_t = amrex::Real(1.0) / amrex::Real(3.0); + amrex::Real Pr_t_inv = amrex::Real(3.0); // Smagorinsky Turbulent Schmidt Number amrex::Real Sc_t = 1.0; amrex::Real Sc_t_inv = 1.0; diff --git a/Source/Diffusion/ComputeStrain_N.cpp b/Source/Diffusion/ComputeStrain_N.cpp index 2906c1247..ccf153ecc 100644 --- a/Source/Diffusion/ComputeStrain_N.cpp +++ b/Source/Diffusion/ComputeStrain_N.cpp @@ -32,7 +32,7 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, const BCRec* bc_ptr, const GpuArray& dxInv, const Array4& mf_m, const Array4& mf_u, const Array4& mf_v) { - // Conver domain to each index type to test if we are on dirichlet boundary + // Convert domain to each index type to test if we are on dirichlet boundary Box domain_xy = convert(domain, tbxxy.ixType()); Box domain_xz = convert(domain, tbxxz.ixType()); Box domain_yz = convert(domain, tbxyz.ixType()); diff --git a/Source/Diffusion/ComputeStrain_T.cpp b/Source/Diffusion/ComputeStrain_T.cpp index bc22b56a4..e8b5bc2d4 100644 --- a/Source/Diffusion/ComputeStrain_T.cpp +++ b/Source/Diffusion/ComputeStrain_T.cpp @@ -42,7 +42,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, const Array4& mf_u, const Array4& mf_v) { - // Conver domain to each index type to test if we are on dirichlet boundary + // Convert domain to each index type to test if we are on dirichlet boundary Box domain_xy = convert(domain, tbxxy.ixType()); Box domain_xz = convert(domain, tbxxz.ixType()); Box domain_yz = convert(domain, tbxyz.ixType()); diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index cf681c12a..cfe8a9cee 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -168,13 +168,13 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv; } - Real E = cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp); - Real strat = l_abs_g * dtheta_dz * l_inv_theta0; // stratification + Real E = cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp); + Real stratification = l_abs_g * dtheta_dz * l_inv_theta0; // stratification Real length; - if (strat <= eps) { + if (stratification <= eps) { length = DeltaMsf; } else { - length = 0.76 * std::sqrt(E / strat); + length = 0.76 * std::sqrt(E / stratification); // mixing length should be _reduced_ for stable stratification length = amrex::min(length, DeltaMsf); // following WRF, make sure the mixing length isn't too small @@ -207,7 +207,7 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, } } - // Extrapolate Kturb in x/y, fill remaining elements (relevent to lev==0) + // Extrapolate Kturb in x/y, fill remaining elements (relevant to lev==0) //*********************************************************************************** int ngc(1); Real inv_Pr_t = turbChoice.Pr_t_inv; @@ -336,9 +336,13 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22, { int indx = n; int indx_v = indx + offset; + mu_turb(i,j,k,indx) = mu_turb(i,j,k,EddyDiff::Mom_h) * fac_ptr[indx-1]; + // NOTE: Theta_v has already been set for Deardorff - if (!(indx_v == EddyDiff::Theta_v && use_KE)) mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); + if (!(indx_v == EddyDiff::Theta_v && use_KE)) { + mu_turb(i,j,k,indx_v) = mu_turb(i,j,k,indx); + } }); break; } diff --git a/Source/Diffusion/PBLModels.H b/Source/Diffusion/PBLModels.H index 7deaa92a9..94e0baf00 100644 --- a/Source/Diffusion/PBLModels.H +++ b/Source/Diffusion/PBLModels.H @@ -114,7 +114,7 @@ ComputeQKESourceTerms (int i, int j, int k, // Production (We store mu_turb, which is 0.5*K_turb) source_term += 4.0*K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz); - // Bouyancy + // Buoyancy source_term -= 2.0*(CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz; // Dissipation diff --git a/Source/EOS.H b/Source/EOS.H index c2e8c6828..133f45598 100644 --- a/Source/EOS.H +++ b/Source/EOS.H @@ -11,7 +11,7 @@ * * @params[in] pressure * @params[in] temperature - * @params[in] rd0cp ratio of R_d to c_p + * @params[in] rdOcp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenPandT(const amrex::Real T, const amrex::Real P, const amrex::Real rdOcp) @@ -41,7 +41,7 @@ amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, * * @params[in] rho density * @params[in] T temperature - * @params[in] rd0cp ratio of R_d to c_p + * @params[in] rdOcp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const amrex::Real rdOcp, const amrex::Real qv=0.0) @@ -69,7 +69,7 @@ amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv = 0.) * * @params[in] theta potential temperature * @params[in] p pressure - * @params[in] rd0cp ratio of R_d to c_p + * @params[in] rdOcp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, const amrex::Real rdOcp, const amrex::Real qv=0.0) @@ -84,7 +84,7 @@ amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, * * @params[in] theta potential temperature * @params[in] p pressure - * @params[in] rd0cp ratio of R_d to c_p + * @params[in] rdOcp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getdPdRgivenConstantTheta(const amrex::Real rho, const amrex::Real theta, const amrex::Real qv=0.0) @@ -97,7 +97,7 @@ amrex::Real getdPdRgivenConstantTheta(const amrex::Real rho, const amrex::Real t /** * Function to return the Exner function pi given pressure * @params[in] p pressure - * @params[in] rd0cp ratio of R_d to c_p + * @params[in] rdOcp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getExnergivenP(const amrex::Real P, const amrex::Real rdOcp) @@ -110,7 +110,7 @@ amrex::Real getExnergivenP(const amrex::Real P, const amrex::Real rdOcp) * Function to return the Exner function pi given densith times potential temperature * * @params[in] rhotheta density times potential temperature - * @params[in] rd0cp ratio of R_d to c_p + * @params[in] rdOcp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=0.0 ) @@ -131,7 +131,7 @@ amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=0.0) { // diagnostic relation for the full pressure // see https://erf.readthedocs.io/en/latest/theory/NavierStokesEquations.html - // For cases with mositure theta = theta_m / (1 + R_v/R_d*qv) + // For cases with moisture, theta = theta_m / (1 + R_v/R_d*qv) return std::pow(p*std::pow(p_0, Gamma-1), iGamma) * iR_d / (1.0 + R_v/R_d*qv) ; } diff --git a/Source/ERF.H b/Source/ERF.H index 4ca77a03a..e7c5cd9a2 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -153,6 +153,8 @@ public: void write_1D_profiles (amrex::Real time); void write_1D_profiles_stag (amrex::Real time); + amrex::Real cloud_fraction (amrex::Real time); + // Fill the physical boundary conditions for cell-centered velocity (diagnostic only) void FillBdyCCVels (amrex::Vector& mf_cc_vel); @@ -241,7 +243,6 @@ public: #ifdef ERF_USE_WW3_COUPLING void read_waves (int lev); - void send_waves (int lev); #endif // Interface for advancing the data at one level by one "slow" timestep diff --git a/Source/ERF.cpp b/Source/ERF.cpp index fa9eb9ca3..3428cb417 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -807,13 +807,6 @@ ERF::InitData () if (is_it_time_for_action(istep[0], t_new[0], dt[0], sum_interval, sum_per)) { sum_integrated_quantities(t_new[0]); - if (cc_profiles) { - // all variables cell-centered - write_1D_profiles(t_new[0]); - } else { - // some variables staggered - write_1D_profiles_stag(t_new[0]); - } } // We only write the file at level 0 for now @@ -1010,6 +1003,16 @@ ERF::InitData () setRecordDataInfo(i,datalogname[i]); } + if (restart_chkfile.empty() && profile_int > 0) { + if (cc_profiles) { + // all variables cell-centered + write_1D_profiles(t_new[0]); + } else { + // some variables staggered + write_1D_profiles_stag(t_new[0]); + } + } + if (pp.contains("sample_point_log") && pp.contains("sample_point")) { int lev = 0; diff --git a/Source/ERF_Constants.H b/Source/ERF_Constants.H index 0fa6852ae..d224c5192 100644 --- a/Source/ERF_Constants.H +++ b/Source/ERF_Constants.H @@ -82,8 +82,8 @@ constexpr amrex::Real crain = b_rain / 4.0; constexpr amrex::Real csnow = b_snow / 4.0; constexpr amrex::Real cgrau = b_grau / 4.0; -constexpr amrex::Real latvap = 2.5e6; // Latent heat of vaporization (J/kg) -constexpr amrex::Real latice = 3.337e5; // latent heat of fusion (J/kg) +constexpr amrex::Real lat_vap = 2.5e6; // Latent heat of vaporization (J/kg) +constexpr amrex::Real lat_ice = 3.337e5; // latent heat of fusion (J/kg) constexpr amrex::Real Rd_on_Rv = R_d/R_v; constexpr amrex::Real tmelt = 273.15; // melting temp. constexpr amrex::Real h2otrip = 273.15; // Triple point temperature of water (K) diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index e216a44f3..8aef07062 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -53,7 +53,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, //******************************************************************************************** // This allocates all kinds of things, including but not limited to: solution arrays, - // terrain arrays and metrics,a nd base state. + // terrain arrays, metric terms and base state. // ******************************************************************************************* init_stuff(lev, ba, dm, lev_new, lev_old); @@ -278,7 +278,7 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, //******************************************************************************************** // This allocates all kinds of things, including but not limited to: solution arrays, - // terrain arrays and metrics,a nd base state. + // terrain arrays, metric terms and base state. // ******************************************************************************************* init_stuff(lev, ba, dm, vars_new[lev], vars_old[lev]); diff --git a/Source/IO/ERF_Write1DProfiles.cpp b/Source/IO/ERF_Write1DProfiles.cpp index d8c4ec0de..62f016099 100644 --- a/Source/IO/ERF_Write1DProfiles.cpp +++ b/Source/IO/ERF_Write1DProfiles.cpp @@ -16,14 +16,11 @@ ERF::write_1D_profiles (Real time) { BL_PROFILE("ERF::write_1D_profiles()"); - if (verbose <= 0) - return; - int datwidth = 14; int datprecision = 6; int timeprecision = 13; // e.g., 1-yr LES: 31,536,000 s with dt ~ 0.01 ==> min prec = 10 - if (verbose > 0 && NumDataLogs() > 1) + if (NumDataLogs() > 1) { // Define the 1d arrays we will need Gpu::HostVector h_avg_u, h_avg_v, h_avg_w; @@ -46,7 +43,7 @@ ERF::write_1D_profiles (Real time) h_avg_p, h_avg_pu, h_avg_pv, h_avg_pw); } - if (NumDataLogs() > 3) { + if (NumDataLogs() > 3 && time > 0.) { derive_stress_profiles(h_avg_tau11, h_avg_tau12, h_avg_tau13, h_avg_tau22, h_avg_tau23, h_avg_tau33, h_avg_sgshfx, @@ -133,7 +130,7 @@ ERF::write_1D_profiles (Real time) } // if good } // NumDataLogs - if (NumDataLogs() > 3) { + if (NumDataLogs() > 3 && time > 0.) { std::ostream& data_log3 = DataLog(3); if (data_log3.good()) { // Write the average stresses @@ -152,9 +149,9 @@ ERF::write_1D_profiles (Real time) << std::endl; } // loop over z } // if good - } // NumDataLogs + } // if (NumDataLogs() > 3) } // if IOProcessor - } // if verbose + } // if (NumDataLogs() > 1) } /** @@ -266,8 +263,16 @@ void ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVecto } fab_arr(i, j, k, 2) = ksgs; Real kturb = 0.0; +#if 1 if (l_use_Turb) kturb = eta_arr(i,j,k,EddyDiff::Mom_h); fab_arr(i, j, k, 3) = kturb; +#else + // Here we hijack the "Kturb" variable name to print out the resolved kinetic energy + Real upert = u_cc_arr(i,j,k) - h_avg_u[k]; + Real vpert = v_cc_arr(i,j,k) - h_avg_v[k]; + Real wpert = w_cc_arr(i,j,k) - h_avg_w[k]; + fab_arr(i, j, k, 3) = 0.5 * (upert*upert + vpert*vpert + wpert*wpert); +#endif fab_arr(i, j, k, 4) = u_cc_arr(i,j,k) * u_cc_arr(i,j,k); // u*u fab_arr(i, j, k, 5) = u_cc_arr(i,j,k) * v_cc_arr(i,j,k); // u*v fab_arr(i, j, k, 6) = u_cc_arr(i,j,k) * w_cc_arr(i,j,k); // u*w @@ -400,7 +405,18 @@ void ERF::derive_diag_profiles(Gpu::HostVector& h_avg_u , Gpu::HostVecto h_avg_wqv[k] /= area_z; h_avg_wqc[k] /= area_z; h_avg_wqr[k] /= area_z; h_avg_qi[k] /= area_z; h_avg_qs[k] /= area_z; h_avg_qg[k] /= area_z; } + +#if 0 + // Here we print the integrated total kinetic energy as computed in the 1D profile above + Real sum = 0.; + Real dz = geom[0].ProbHi(2) / static_cast(h_avg_u_size); + for (int k = 0; k < h_avg_u_size; ++k) { + sum += h_avg_kturb[k] * h_avg_rho[k] * dz; + } + amrex::Print() << "ITKE " << time << " " << sum << " using " << h_avg_u_size << " " << dz << std::endl; +#endif } + void ERF::derive_stress_profiles (Gpu::HostVector& h_avg_tau11, Gpu::HostVector& h_avg_tau12, Gpu::HostVector& h_avg_tau13, Gpu::HostVector& h_avg_tau22, diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp index 78a61f9ec..0370afa54 100644 --- a/Source/IO/ERF_WriteScalarProfiles.cpp +++ b/Source/IO/ERF_WriteScalarProfiles.cpp @@ -126,8 +126,8 @@ ERF::sum_integrated_quantities (Real time) // The first data log only holds scalars if (NumDataLogs() > 0) { - int nd = 0; - std::ostream& data_log1 = DataLog(nd); + int n_d = 0; + std::ostream& data_log1 = DataLog(n_d); if (data_log1.good()) { if (time == 0.0) { data_log1 << std::setw(datwidth) << " time"; @@ -169,6 +169,68 @@ ERF::sum_integrated_quantities (Real time) } } +Real +ERF::cloud_fraction (Real time) +{ + BL_PROFILE("ERF::cloud_fraction()"); + + Real sum = 0.0; + + int lev = 0; + // This holds all of qc + MultiFab qc(vars_new[lev][Vars::cons],make_alias,RhoQ2_comp,1); + + int direction = 2; // z-direction + Box const& domain = geom[lev].Domain(); + + auto const& qc_arr = qc.const_arrays(); + + // qc_2d is an BaseFab holding the max value over the column + auto qc_2d = ReduceToPlane(direction, domain, qc, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) -> int + { + if (qc_arr[box_no](i,j,k) > 0) { + return 1; + } else { + return 0; + } + }); + + auto* p = qc_2d.dataPtr(); + + Long numpts = qc_2d.numPts(); + + AMREX_ASSERT(numpts < Long(std::numeric_limits::max)); + +#if 1 + ParallelDescriptor::ReduceIntMax(p,static_cast(numpts)); + + // Sum over component 0 + Long num_cloudy = qc_2d.template sum(0); + +#else + // + // We need this if we allow domain decomposition in the vertical + // but for now we leave it commented out + // + Long num_cloudy = Reduce::Sum(numpts, + [=] AMREX_GPU_DEVICE (Long i) -> Long { + if (p[i] == 1) { + return 1; + } else { + return 0; + } + }); + ParallelDescriptor::ReduceLongSum(num_cloudy); +#endif + + Real num_total = qc_2d.box().d_numPts(); + + Real cloud_frac = num_cloudy / num_total; + + return cloud_frac; +} + /** * Utility function for sampling MultiFab data at a specified cell index. * diff --git a/Source/IO/NCWpsFile.H b/Source/IO/NCWpsFile.H index be8fc9577..268968a04 100644 --- a/Source/IO/NCWpsFile.H +++ b/Source/IO/NCWpsFile.H @@ -42,7 +42,7 @@ enum class NC_Data_Dims_Type { // NDArray is the datatype designed to hold any data, including scalars, multidimensional // arrays, that read from the NetCDF file. // -// The data read from NetCDF file are stored in a continuous memory, and the data layout is descriped +// The data read from NetCDF file are stored in a continuous memory, and the data layout is described // by using a vector (shape). AMRex Box can be constructed using the data shape information, and MultiFab // data array can be setup using the data that stored in the NDArray. // diff --git a/Source/Initialization/ERF_init_TurbPert.cpp b/Source/Initialization/ERF_init_TurbPert.cpp index beb2c924d..28aacd70e 100644 --- a/Source/Initialization/ERF_init_TurbPert.cpp +++ b/Source/Initialization/ERF_init_TurbPert.cpp @@ -22,7 +22,7 @@ ERF::turbPert_constants(const int lev) void ERF::turbPert_update (const int lev, const Real local_dt) { - // Grabing data from velocity field + // Grabbing data from velocity field auto& lev_new = vars_new[lev]; // Accessing local data diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 71dcc1748..b950d9c49 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -574,16 +574,16 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, * @param y_vel_fab FArrayBox holding the y-velocity data to initialize * @param z_vel_fab FArrayBox holding the z-velocity data to initialize * @param z_phys_nd_fab FArrayBox holding nodal z coordinate data for terrain - * @param NC_hgt_fab Vector of FArrayBox obects holding metgrid data for terrain height - * @param NC_ght_fab Vector of FArrayBox objects holding metgrid data for height of cell centers - * @param NC_xvel_fab Vector of FArrayBox obects holding metgrid data for x-velocity - * @param NC_yvel_fab Vector of FArrayBox obects holding metgrid data for y-velocity - * @param NC_zvel_fab Vector of FArrayBox obects holding metgrid data for z-velocity - * @param NC_temp_fab Vector of FArrayBox obects holding metgrid data for temperature - * @param NC_rhum_fab Vector of FArrayBox obects holding metgrid data for relative humidity - * @param NC_pres_fab Vector of FArrayBox obects holding metgrid data for pressure - * @param theta_fab Vector of FArrayBox obects holding potential temperature calculated from temperature and pressure - * @param mxrat_fab Vector of FArrayBox obects holding vapor mixing ratio calculated from relative humidity + * @param NC_hgt_fab Vector of FArrayBox objects holding metgrid data for terrain height + * @param NC_ght_fab Vector of FArrayBox objects holding metgrid data for height of cell centers + * @param NC_xvel_fab Vector of FArrayBox objects holding metgrid data for x-velocity + * @param NC_yvel_fab Vector of FArrayBox objects holding metgrid data for y-velocity + * @param NC_zvel_fab Vector of FArrayBox objects holding metgrid data for z-velocity + * @param NC_temp_fab Vector of FArrayBox objects holding metgrid data for temperature + * @param NC_rhum_fab Vector of FArrayBox objects holding metgrid data for relative humidity + * @param NC_pres_fab Vector of FArrayBox objects holding metgrid data for pressure + * @param theta_fab Vector of FArrayBox objects holding potential temperature calculated from temperature and pressure + * @param mxrat_fab Vector of FArrayBox objects holding vapor mixing ratio calculated from relative humidity * @param fabs_for_bcs Vector of Vector of FArrayBox objects holding MetGridBdyVars at each met_em time. */ void diff --git a/Source/Microphysics/EulerianMicrophysics.H b/Source/Microphysics/EulerianMicrophysics.H index ba5f72735..de056e20d 100644 --- a/Source/Microphysics/EulerianMicrophysics.H +++ b/Source/Microphysics/EulerianMicrophysics.H @@ -1,5 +1,5 @@ /*! @file EulerianMicrophysics.H - * \brief Containes the Eulerian microphysics class + * \brief Contains the Eulerian microphysics class */ #ifndef EULERIANMICROPHYSICS_H diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index f00d8759b..2c67ffdd9 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -98,7 +98,7 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); //------- Autoconversion/accretion - Real qcc, autor, accrr; + Real qcc, auto_r, accrr; Real dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; Real pressure = pres_array(i,j,k); @@ -116,7 +116,7 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2)); //Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); - // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature + // If water vapor content exceeds saturation value, then vapor condenses to water and latent heat is released, increasing temperature if (qv_array(i,j,k) > qsat) { dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac)); } @@ -142,14 +142,14 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) if (qc_array(i,j,k) > 0.0) { qcc = qc_array(i,j,k); - autor = 0.0; + auto_r = 0.0; if (qcc > qcw0) { - autor = alphaelq; + auto_r = alphaelq; } accrr = 0.0; accrr = 2.2 * std::pow(qp_array(i,j,k) , 0.875); - dq_clwater_to_rain = dtn *(accrr*qcc + autor*(qcc - qcw0)); + dq_clwater_to_rain = dtn *(accrr*qcc + auto_r*(qcc - qcw0)); // If the amount of change is more than the amount of qc present, then dq = qc dq_clwater_to_rain = std::min(dq_clwater_to_rain, qc_array(i,j,k)); @@ -212,7 +212,7 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) //Real fac = qsat*4093.0*L_v/(Cp_d*std::pow(tabs_array(i,j,k)-36.0,2)); Real fac = qsat*L_v*L_v/(Cp_d*R_v*tabs_array(i,j,k)*tabs_array(i,j,k)); - // If water vapor content exceeds saturation value, then vapor condenses to waterm and latent heat is released, increasing temperature + // If water vapor content exceeds saturation value, then vapor condenses to water and latent heat is released, increasing temperature if (qv_array(i,j,k) > qsat){ dq_vapor_to_clwater = std::min(qv_array(i,j,k), (qv_array(i,j,k)-qsat)/(1.0 + fac)); } diff --git a/Source/Microphysics/LagrangianMicrophysics.H b/Source/Microphysics/LagrangianMicrophysics.H index a1b6a39c8..9daafd242 100644 --- a/Source/Microphysics/LagrangianMicrophysics.H +++ b/Source/Microphysics/LagrangianMicrophysics.H @@ -1,5 +1,5 @@ /*! @file LagrangianMicrophysics.H - * \brief Containes the Lagrangian microphysics class + * \brief Contains the Lagrangian microphysics class */ #ifndef LAGRANGIANMICROPHYSICS_H diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index dd3bc87b4..fc8d480f0 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -1,5 +1,5 @@ /*! @file Microphysics.H - * \brief Containes the base class for microphysics + * \brief Contains the base class for microphysics */ #ifndef MICROPHYSICS_H diff --git a/Source/Microphysics/SAM/Precip.cpp b/Source/Microphysics/SAM/Precip.cpp index 6582c5290..a97b3b86e 100644 --- a/Source/Microphysics/SAM/Precip.cpp +++ b/Source/Microphysics/SAM/Precip.cpp @@ -81,7 +81,7 @@ SAM::Precip (const SolverChoice& sc) Real dqc, dqca, dqi, dqia, dqp; Real dqpr, dqps, dqpg; - Real autor, autos; + Real auto_r, autos; Real accrcr, accrcs, accris, accrcg, accrig; // Work to be done for autoc/accr or evap @@ -114,9 +114,9 @@ SAM::Precip (const SolverChoice& sc) accrig = 0.0; if (qcc > qcw0) { - autor = alphaelq; + auto_r = alphaelq; } else { - autor = 0.0; + auto_r = 0.0; } if (qii > qci0) { @@ -140,7 +140,7 @@ SAM::Precip (const SolverChoice& sc) } // Autoconversion & accretion (sink for cloud comps) - dqca = dtn * autor * (qcc-qcw0); + dqca = dtn * auto_r * (qcc-qcw0); dprc = dtn * accrcr * qcc * std::pow(qpr, powr1); dpsc = dtn * accrcs * qcc * std::pow(qps, pows1); dpgc = dtn * accrcg * qcc * std::pow(qpg, powg1); diff --git a/Source/Particles/ERFPCUtils.cpp b/Source/Particles/ERFPCUtils.cpp index 0f83345ae..8f9841fb0 100644 --- a/Source/Particles/ERFPCUtils.cpp +++ b/Source/Particles/ERFPCUtils.cpp @@ -23,15 +23,15 @@ void ERFPC::massDensity ( MultiFab& a_mf, a_mf.setVal(0.0); ParticleToMesh( *this, a_mf, a_lev, - [=] AMREX_GPU_DEVICE ( const ERFPC::ParticleTileType::ConstParticleTileDataType& ptd, + [=] AMREX_GPU_DEVICE ( const ERFPC::ParticleTileType::ConstParticleTileDataType& ptr, int i, Array4 const& rho) { - auto p = ptd.m_aos[i]; + auto p = ptr.m_aos[i]; ParticleInterpolator::Linear interp(p, plo, dxi); interp.ParticleToMesh ( p, rho, 0, a_comp, 1, [=] AMREX_GPU_DEVICE ( const ERFPC::ParticleType&, int) { - auto mass = ptd.m_rdata[ERFParticlesRealIdxSoA::mass][i]; + auto mass = ptr.m_rdata[ERFParticlesRealIdxSoA::mass][i]; return mass*inv_cell_volume; }); }); diff --git a/Source/Prob/init_density_hse_dry.H b/Source/Prob/init_density_hse_dry.H index 5bebcf6aa..c4dc5cdcb 100644 --- a/Source/Prob/init_density_hse_dry.H +++ b/Source/Prob/init_density_hse_dry.H @@ -2,7 +2,7 @@ * Initialize hydrostatically balanced density * * Calls init_isentropic_hse_terrain() or init_isentropic_hse() for cases with - * and without terrain. Hydrostatic equlibrium (HSE) is satisfied discretely. + * and without terrain. Hydrostatic equilibrium (HSE) is satisfied discretely. * Note that these routines presume that qv==0 when evaluating the EOS. Both * density and pressure in HSE are calculated but only the density is used at * this point. diff --git a/Source/Radiation/Aero_rad_props.H b/Source/Radiation/Aero_rad_props.H index 18d3d21fe..acedb4900 100644 --- a/Source/Radiation/Aero_rad_props.H +++ b/Source/Radiation/Aero_rad_props.H @@ -12,7 +12,7 @@ #include "Mam4_constituents.H" // -// AerRadProps class convers the aerosol mass to bulk optical properties for short wave +// AerRadProps class converts the aerosol mass to bulk optical properties for short wave // and long wave radiation computations // class AerRadProps { diff --git a/Source/Radiation/Linear_interpolate.H b/Source/Radiation/Linear_interpolate.H index 421b9d2ee..cfa644dd4 100644 --- a/Source/Radiation/Linear_interpolate.H +++ b/Source/Radiation/Linear_interpolate.H @@ -26,7 +26,7 @@ class LinInterp // Initialize a variable of type(interp_type) with weights for linear interpolation. // this variable can then be used in calls to lininterp1d and lininterp2d. - // yin is a 1d array of length nin of locations to interpolate from - this array must + // yin is a 1d array of length n_in of locations to interpolate from - this array must // be monotonic but can be increasing or decreasing // yout is a 1d array of length nout of locations to interpolate to, this array need // not be ordered @@ -38,7 +38,7 @@ class LinInterp // cyclic mapping - these default to 0 and 360. YAKL_INLINE static void init (const real1d& yin, - const int& nin, + const int& n_in, const real1d& yout, const int& nout, const InterpMethod& extrap_method, @@ -71,7 +71,7 @@ class LinInterp interp_wgts.wgts = real1d("wgts", nout); interp_wgts.wgtn = real1d("wgtn", nout); - parallel_for(SimpleBounds<1>(nin-1), YAKL_LAMBDA (int j) + parallel_for(SimpleBounds<1>(n_in-1), YAKL_LAMBDA (int j) { if(yin(j) > yin(j+1)) { printf("init: inputs are not monotonic!\n"); @@ -103,9 +103,9 @@ class LinInterp interp_wgts.wgts(j) = 0.; interp_wgts.wgtn(j) = 0.; } - else if (yout(j) > yin(nin)) { - interp_wgts.jjm(j) = nin; - interp_wgts.jjp(j) = nin; + else if (yout(j) > yin(n_in)) { + interp_wgts.jjm(j) = n_in; + interp_wgts.jjp(j) = n_in; interp_wgts.wgts(j) = 0.; interp_wgts.wgtn(j) = 0.; } @@ -120,9 +120,9 @@ class LinInterp interp_wgts.wgts(j) = 0.; interp_wgts.wgtn(j) = 0.; } - else if (yout(j) < yin(nin)) { - interp_wgts.jjm(j) = nin; - interp_wgts.jjp(j) = nin; + else if (yout(j) < yin(n_in)) { + interp_wgts.jjm(j) = n_in; + interp_wgts.jjp(j) = n_in; interp_wgts.wgts(j) = 0.; interp_wgts.wgtn(j) = 0.; } @@ -143,9 +143,9 @@ class LinInterp interp_wgts.wgts(j) = 1.; interp_wgts.wgtn(j) = 0.; } - else if (yout(j) > yin(nin)) { - interp_wgts.jjm(j) = nin; - interp_wgts.jjp(j) = nin; + else if (yout(j) > yin(n_in)) { + interp_wgts.jjm(j) = n_in; + interp_wgts.jjp(j) = n_in; interp_wgts.wgts(j) = 1.; interp_wgts.wgtn(j) = 0.; } @@ -160,9 +160,9 @@ class LinInterp interp_wgts.wgts(j) = 1.; interp_wgts.wgtn(j) = 0.; } - else if (yout(j) <= yin(nin)) { - interp_wgts.jjm(j) = nin; - interp_wgts.jjp(j) = nin; + else if (yout(j) <= yin(n_in)) { + interp_wgts.jjm(j) = n_in; + interp_wgts.jjp(j) = n_in; interp_wgts.wgts(j) = 1.; interp_wgts.wgtn(j) = 0.; } @@ -174,28 +174,28 @@ class LinInterp // For values which extend beyond boundaries, set weights // for circular boundaries // - dyinwrap = yin(1) + (cmax-cmin) - yin(nin); - avgdyin = abs(yin(nin)-yin(1))/(nin-1.); + dyinwrap = yin(1) + (cmax-cmin) - yin(n_in); + avgdyin = abs(yin(n_in)-yin(1))/(n_in-1.); auto ratio = dyinwrap/avgdyin; if (ratio < 0.9 || ratio > 1.1) { - printf("ratio is too large, dyinwrap= %13.6e, avgdyin= %13.6e, yin(1) = %13.6e, yin(nin)= %13.6e\n", - dyinwrap, avgdyin, yin(1), yin(nin)); + printf("ratio is too large, dyinwrap= %13.6e, avgdyin= %13.6e, yin(1) = %13.6e, yin(n_in)= %13.6e\n", + dyinwrap, avgdyin, yin(1), yin(n_in)); } if(increasing) { parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j) { if (yout(j) <= yin(1)) { - interp_wgts.jjm(j) = nin; + interp_wgts.jjm(j) = n_in; interp_wgts.jjp(j) = 1; interp_wgts.wgts(j) = (yin(1)-yout(j))/dyinwrap; - interp_wgts.wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap; + interp_wgts.wgtn(j) = (yout(j)+(cmax-cmin) - yin(n_in))/dyinwrap; } - else if (yout(j) > yin(nin)) { - interp_wgts.jjm(j) = nin; + else if (yout(j) > yin(n_in)) { + interp_wgts.jjm(j) = n_in; interp_wgts.jjp(j) = 1; interp_wgts.wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap; - interp_wgts.wgtn(j) = (yout(j)-yin(nin))/dyinwrap; + interp_wgts.wgtn(j) = (yout(j)-yin(n_in))/dyinwrap; } }); } @@ -203,16 +203,16 @@ class LinInterp parallel_for(SimpleBounds<1>(nout), YAKL_LAMBDA (int j) { if (yout(j) > yin(1)) { - interp_wgts.jjm(j) = nin; + interp_wgts.jjm(j) = n_in; interp_wgts.jjp(j) = 1; interp_wgts.wgts(j) = (yin(1)-yout(j))/dyinwrap; - interp_wgts.wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap; + interp_wgts.wgtn(j) = (yout(j)+(cmax-cmin) - yin(n_in))/dyinwrap; } - else if (yout(j) <= yin(nin)) { - interp_wgts.jjm(j) = nin; + else if (yout(j) <= yin(n_in)) { + interp_wgts.jjm(j) = n_in; interp_wgts.jjp(j) = 1; interp_wgts.wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap; - interp_wgts.wgtn(j) = (yout(j)+(cmax-cmin)-yin(nin))/dyinwrap; + interp_wgts.wgtn(j) = (yout(j)+(cmax-cmin)-yin(n_in))/dyinwrap; } }); } @@ -222,7 +222,7 @@ class LinInterp // Loop though output indices finding input indices and weights // if(increasing) { - parallel_for(SimpleBounds<2>(nout, nin-1), YAKL_LAMBDA (int j, int jj) + parallel_for(SimpleBounds<2>(nout, n_in-1), YAKL_LAMBDA (int j, int jj) { if (yout(j) > yin(jj) && yout(j) <= yin(jj+1)) { interp_wgts.jjm(j) = jj; @@ -234,7 +234,7 @@ class LinInterp }); } else { - parallel_for(SimpleBounds<2>(nout, nin-1), YAKL_LAMBDA (int j, int jj) + parallel_for(SimpleBounds<2>(nout, n_in-1), YAKL_LAMBDA (int j, int jj) { if (yout(j) <= yin(jj) && yout(j) > yin(jj+1)) { interp_wgts.jjm(j) = jj; diff --git a/Source/Radiation/Mam4_aero.H b/Source/Radiation/Mam4_aero.H index fb6e55c7d..fa4fd1254 100644 --- a/Source/Radiation/Mam4_aero.H +++ b/Source/Radiation/Mam4_aero.H @@ -43,13 +43,13 @@ class Mam4_aer { real1d crefwlwr; // complex refractive index for water infrared real1d crefwlwi; - // These are defined as module level variables to aviod allocation-deallocation in a loop + // These are defined as module level variables to avoid allocation-deallocation in a loop real3d dgnumdry_m; // number mode dry diameter for all modes real3d dgnumwet_m; // number mode wet diameter for all modes real3d qaerwat_m; // aerosol water (g/g) for all modes real3d wetdens_m; - // aerosol constitutents + // aerosol constituents MamConstituents mam_consti; public: // initialization @@ -108,7 +108,7 @@ class Mam4_aer { int nsw_bands = water_ref_file.getDimSize( "sw_band" ); if (nswbands != nsw_bands || nlwbands != nlw_bands) - amrex::Print() << "ERROR - file and bandwith values do not match: " + amrex::Print() << "ERROR - file and bandwidth values do not match: " << "\n nswbands: " << nswbands << "; nsw_bands: " << nsw_bands << "\n nlwbands: " << nlwbands << "; nlw_bands: " << nlw_bands << "\n"; @@ -553,7 +553,7 @@ class Mam4_aer { } // modal_aero_sw // - // calcualtes aerosol lw radiative properties + // calculates aerosol lw radiative properties // void modal_aero_lw (int list_idx, real dt, const real2d& pdeldry, const real2d& pmid, const real2d& temperature, const real2d& qt, diff --git a/Source/Radiation/Mam4_constituents.H b/Source/Radiation/Mam4_constituents.H index 901edf925..a1c420d53 100644 --- a/Source/Radiation/Mam4_constituents.H +++ b/Source/Radiation/Mam4_constituents.H @@ -207,7 +207,7 @@ class MamConstituents { aerosollist.resize(N_DIAG); ma_list.resize(N_DIAG); - // Set the list_id fields which distinquish the climate and diagnostic lists + // Set the list_id fields which distinguish the climate and diagnostic lists for(auto i = 0; i < N_DIAG; ++i) { aerosollist[i].list_id = ""; gaslist[i].list_id = ""; diff --git a/Source/Radiation/Modal_aero_wateruptake.H b/Source/Radiation/Modal_aero_wateruptake.H index 2900550a9..22ae2cfba 100644 --- a/Source/Radiation/Modal_aero_wateruptake.H +++ b/Source/Radiation/Modal_aero_wateruptake.H @@ -88,7 +88,7 @@ class ModalAeroWateruptake { qaerwat = qaerwat_m; wetdens = wetdens_m; - // retreive aerosol properties + // retrieve aerosol properties for (auto m = 1; m <= nmodes; ++m) { yakl::memset(dryvolmr, 0.); @@ -261,7 +261,7 @@ class ModalAeroWateruptake { } // modes } - // calculates equlibrium radius r of haze droplets as function of + // calculates equilibrium radius r of haze droplets as function of // dry particle mass and relative humidity s using kohler solution // given in pruppacher and klett (eqn 6-35) diff --git a/Source/Radiation/Optics.cpp b/Source/Radiation/Optics.cpp index c74b9f2c4..6fbc25ea0 100644 --- a/Source/Radiation/Optics.cpp +++ b/Source/Radiation/Optics.cpp @@ -152,7 +152,7 @@ void Optics::get_cloud_optics_sw (int ncol, int nlev, int nbnd, } // Copy to output arrays, converting to optical depth, single scattering - // albedo, and assymmetry parameter from the products that the CAM routines + // albedo, and asymmetry parameter from the products that the CAM routines // return. Make sure we do not try to divide by zero... parallel_for(SimpleBounds<3>(nbnd, ncol, nlev), YAKL_LAMBDA (int iband, int icol, int ilev) { @@ -401,7 +401,7 @@ void Optics::set_aerosol_optics_sw (int icall, int ncol, int nlev, int nswbands, ssa_out(icol,ilev,iband) = 1.; } - // Extract assymmetry parameter from the product-defined fields + // Extract asymmetry parameter from the product-defined fields if (tau_w(icol,ilev,iband) > 0) { asm_out(icol,ilev,iband) = tau_w_g(icol,ilev,iband) / tau_w(icol,ilev,iband); } else { diff --git a/Source/Radiation/Phys_prop.H b/Source/Radiation/Phys_prop.H index 54c30ba34..ce07505a1 100644 --- a/Source/Radiation/Phys_prop.H +++ b/Source/Radiation/Phys_prop.H @@ -74,7 +74,7 @@ class PhysProp { real dgnum; // geometric dry mean diameter of the number distribution for aerosol mode real dgnumlo; // lower limit of dgnum real dgnumhi; // upper limit of dgnum - real rhcrystal; // crystalization relative humidity for mode + real rhcrystal; // crystallization relative humidity for mode real rhdeliques; // deliquescence relative humidity for mode }; diff --git a/Source/Radiation/Radiation.H b/Source/Radiation/Radiation.H index f05504a74..6e4539b08 100644 --- a/Source/Radiation/Radiation.H +++ b/Source/Radiation/Radiation.H @@ -1,10 +1,10 @@ /* * RTE-RRTMGP radiation model interface to ERF - * The orginal code is developed by RobertPincus, and the code is open source available at: + * The original code is developed by RobertPincus, and the code is open source available at: * https://github.com/earth-system-radiation/rte-rrtmgp * Please reference to the following paper, * https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001621 - * NOTE: we use the C++ version of RTE-RRTMGP, which is the implemention of the original Fortran + * NOTE: we use the C++ version of RTE-RRTMGP, which is the implementation of the original Fortran * code using C++ YAKL for CUDA, HiP and SYCL application by E3SM ECP team, the C++ version * of the rte-rrtmgp code is located at: * https://github.com/E3SM-Project/rte-rrtmgp @@ -143,7 +143,7 @@ class Radiation { // rrtmgp Rrtmgp radiation; - // optics radation properties + // optics radiation properties Optics optics; // aerosol optics properties @@ -185,7 +185,7 @@ class Radiation { // Disabled when value is less than 0. real fixed_total_solar_irradiance = -1.; - // The RRTMGP warnings are printed when the state variables need to be limted, + // The RRTMGP warnings are printed when the state variables need to be limited, // such as when the temperature drops too low. This is not normally an issue, // but in aquaplanet and RCE configurations these situations occur much more // frequently, so this flag was added to be able to disable those messages. diff --git a/Source/Radiation/Radiation.cpp b/Source/Radiation/Radiation.cpp index e0b9a62a1..cefd52354 100644 --- a/Source/Radiation/Radiation.cpp +++ b/Source/Radiation/Radiation.cpp @@ -1,6 +1,6 @@ /* * RTE-RRTMGP radiation model interface to ERF - * The orginal code is developed by RobertPincus, and the code is open source available at: + * The original code is developed by RobertPincus, and the code is open source available at: * https://github.com/earth-system-radiation/rte-rrtmgp * Please reference to the following paper, * https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001621 @@ -256,8 +256,8 @@ void Radiation::initialize (const MultiFab& cons_in, ncol, nlev, nrh, top_lev, aero_names, zi, pmid, pint, tmid, qt, geom_radius); - amrex::Print() << "LW coefficents file: " << rrtmgp_coefficients_file_lw - << "\nSW coefficents file: " << rrtmgp_coefficients_file_sw + amrex::Print() << "LW coefficients file: " << rrtmgp_coefficients_file_lw + << "\nSW coefficients file: " << rrtmgp_coefficients_file_sw << "\nFrequency (timesteps) of Shortwave Radiation calc: " << dt << "\nFrequency (timesteps) of Longwave Radiation calc: " << dt << "\nDo aerosol radiative calculations: " << do_aerosol_rad << std::endl; @@ -310,7 +310,7 @@ void Radiation::run () // Needed for shortwave aerosol; //int nday, nnight; // Number of daylight columns - int1d day_indices("day_indices", ncol), night_indices("night_indices", ncol); // Indicies of daylight coumns + int1d day_indices("day_indices", ncol), night_indices("night_indices", ncol); // Indices of daylight coumns // Flag to carry (QRS,QRL)*dp across time steps. // TODO: what does this mean? @@ -576,9 +576,9 @@ void Radiation::radiation_driver_sw (int ncol, const real3d& gas_vmr, // and earth-sun distance real2d solar_irradiance_by_gpt("solar_irradiance_by_gpt",ncol,nswgpts); - // Gathered indicies of day and night columns + // Gathered indices of day and night columns // chunk_column_index = day_indices(daylight_column_index) - int1d day_indices("day_indices",ncol), night_indices("night_indices", ncol); // Indicies of daylight coumns + int1d day_indices("day_indices",ncol), night_indices("night_indices", ncol); // Indices of daylight coumns real1d coszrs_day("coszrs_day", ncol); real2d albedo_dir_day("albedo_dir_day", nswbands, ncol), albedo_dif_day("albedo_dif_day", nswbands, ncol); @@ -791,7 +791,7 @@ void Radiation::set_daynight_indices (const real1d& coszrs, const int1d& day_ind { // Loop over columns and identify daytime columns as those where the cosine // solar zenith angle exceeds zero. Note that we wrap the setting of - // day_indices in an if-then to make sure we are not accesing day_indices out + // day_indices in an if-then to make sure we are not accessing day_indices out // of bounds, and stopping with an informative error message if we do for some reason. int1d iday("iday", 1); int1d inight("inight",1); diff --git a/Source/Radiation/Rrtmgp.H b/Source/Radiation/Rrtmgp.H index 31a37a7f6..45bc233a5 100644 --- a/Source/Radiation/Rrtmgp.H +++ b/Source/Radiation/Rrtmgp.H @@ -1,6 +1,6 @@ /* * RTE-RRTMGP radiation model interface to ERF - * The orginal code is developed by RobertPincus, and the code is open source available at: + * The orgiinal code is developed by RobertPincus, and the code is open source available at: * https://github.com/earth-system-radiation/rte-rrtmgp * For details of the radiation algorithm, please reference to the following paper, * https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001621 diff --git a/Source/Radiation/m2005_effradius.H b/Source/Radiation/m2005_effradius.H index b33fd37c9..b66142d75 100644 --- a/Source/Radiation/m2005_effradius.H +++ b/Source/Radiation/m2005_effradius.H @@ -68,7 +68,7 @@ void m2005_effradius (const real2d& ql, const real2d& nl, const real2d& qi, res = 500.; } - // from Hugh Morrision: rhos/917 accouts for assumptions about + // from Hugh Morrision: rhos/917 accounts for assumptions about // ice density in the Mitchell optics. des(i,j) = res*rhos/917.*2.; diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index f8c6baa59..668710512 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/Source/TimeIntegration/ERF_TimeStep.cpp @@ -87,6 +87,12 @@ ERF::timeStep (int lev, Real time, int /*iteration*/) timeStep(lev+1, strt_time_for_fine, i); } } + + if (verbose && lev == 0) { + if (solverChoice.moisture_type != MoistureType::None) { + amrex::Print() << "Cloud fraction " << time << " " << cloud_fraction(time) << std::endl; + } + } } /** diff --git a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp index 118c33908..d15af909a 100644 --- a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp +++ b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp @@ -28,7 +28,7 @@ void make_fast_coeffs (int /*level*/, MultiFab& fast_coeffs, Vector& S_stage_data, // S_bar = S^n, S^* or S^** const MultiFab& S_stage_prim, - const MultiFab& pi_stage, // Exner function evaluted at least stage + const MultiFab& pi_stage, // Exner function evaluated at least stage const amrex::Geometry geom, bool l_use_moisture, bool l_use_terrain, diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index b985380d0..9e3360319 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -245,14 +245,14 @@ void erf_slow_rhs_post (int level, int finest_level, const GpuArray scomp_slow = { 2,0,0,0}; const GpuArray ncomp_slow = {nsv,0,0,0}; - { - BL_PROFILE("rhs_post_7"); + // ************************************************************************** + // Note that here we do copy only the "slow" variables, not (rho) or (rho theta) + // ************************************************************************** ParallelFor(tbx, ncomp_slow[IntVars::cons], [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) { const int n = scomp_slow[IntVars::cons] + nn; cur_cons(i,j,k,n) = new_cons(i,j,k,n); }); - } // end profile // We have projected the velocities stored in S_data but we will use // the velocities stored in S_scratch to update the scalars, so @@ -354,7 +354,7 @@ void erf_slow_rhs_post (int level, int finest_level, if (l_use_terrain) { DiffusionSrcForState_T(tbx, domain, start_comp, num_comp, exp_most, u, v, - cur_cons, cur_prim, cell_rhs, + new_cons, cur_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, z_nd, ax_arr, ay_arr, az_arr, detJ_arr, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, @@ -362,7 +362,7 @@ void erf_slow_rhs_post (int level, int finest_level, mu_turb, dc, tc, tm_arr, grav_gpu, bc_ptr_d, use_most); } else { DiffusionSrcForState_N(tbx, domain, start_comp, num_comp, exp_most, u, v, - cur_cons, cur_prim, cell_rhs, + new_cons, cur_prim, cell_rhs, diffflux_x, diffflux_y, diffflux_z, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, diss, diff --git a/Source/Utils/Water_vapor_saturation.H b/Source/Utils/Water_vapor_saturation.H index f2b4d9220..13cd8381b 100644 --- a/Source/Utils/Water_vapor_saturation.H +++ b/Source/Utils/Water_vapor_saturation.H @@ -75,8 +75,8 @@ public: AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static void no_ip_hltalt (const amrex::Real& t, amrex::Real hltalt) { - hltalt = latvap; - // Account for change of latvap with t above freezing where + hltalt = lat_vap; + // Account for change of lat_vap with t above freezing where // constant slope is given by -2369 j/(kg c) = cpv - cw if (t >= tmelt) hltalt = hltalt - 2369.0*(t-tmelt); } @@ -110,7 +110,7 @@ public: weight = 1.0; } - hltalt = hltalt + weight*latice; + hltalt = hltalt + weight*lat_ice; } } @@ -220,7 +220,7 @@ public: SatMethods::wv_sat_qsat_ice(t, p, es, qs); // For pure ice, just add latent heats. - hltalt = latvap + latice; + hltalt = lat_vap + lat_ice; enthalpy = tq_enthalpy(t, qs, hltalt);