ERF¶

ERF solves the compressible Navier-Stokes on a Arakawa C-grid for large-scale weather modeling.

ERF is built on AMReX, an adaptive mesh refinement software framework, which provides the underlying software infrastructure and performance portability. Visit the AMReX documentation and AMRex tutorials for more information.

ERF is designed to run on machines from laptops to multicore CPU and hybrid CPU/GPU systems.

In addition to this documentation, there is API documentation for ERF generated by Doxygen.

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.

Our Standards¶

Examples of behavior that contributes to creating a positive environment include:

Using welcoming and inclusive language
Being respectful of differing viewpoints and experiences
Gracefully accepting constructive criticism
Focusing on what is best for the community
Showing empathy towards other community members

Examples of unacceptable behavior by participants include:

The use of sexualized language or imagery and unwelcome sexual attention or advances
Trolling, insulting/derogatory comments, and personal or political attacks
Public or private harassment
Publishing others’ private information, such as a physical or electronic address, without explicit permission
Other conduct which could reasonably be considered inappropriate in a professional setting

Our Responsibilities¶

Project maintainers are responsible for clarifying the standards of acceptable behavior and are expected to take appropriate and fair corrective action in response to any instances of unacceptable behavior.

Project maintainers have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, issues, and other contributions that are not aligned to this Code of Conduct, or to ban temporarily or permanently any contributor for other behaviors that they deem inappropriate, threatening, offensive, or harmful.

Scope¶

This Code of Conduct applies both within project spaces and in public spaces when an individual is representing the project or its community. Examples of representing a project or community include using an official project e-mail address, posting via an official social media account, or acting as an appointed representative at an online or offline event. Representation of a project may be further defined and clarified by project maintainers.

Enforcement¶

Instances of abusive, harassing, or otherwise unacceptable behavior may be reported by contacting the project team. All complaints will be reviewed and investigated and will result in a response that is deemed necessary and appropriate to the circumstances. The project team is obligated to maintain confidentiality with regard to the reporter of an incident. Further details of specific enforcement policies may be posted separately.

Project maintainers who do not follow or enforce the Code of Conduct in good faith may face temporary or permanent repercussions as determined by other members of the project’s leadership.

Attribution¶

This Code of Conduct is adapted from the Contributor Covenant, version 1.4, available at https://www.contributor-covenant.org/version/1/4/code-of-conduct/

For answers to common questions about this code of conduct, see https://www.contributor-covenant.org/faq

Doxygen Documentation¶

For details about the ERF code structure and APIs, see the Doxygen documentation

Getting Started¶

Downloading the code¶

First, make sure that git is installed on your machine.

Then download the ERF repository by typing:

git clone https://github.com/erf-model/ERF.git

Or, to automatically include the necessary submodules when downloading ERF, type:

git clone --recursive https://github.com/erf-model/ERF.git

Git Submodules¶

When using the submodule to build, it is ideal to properly update and match what is in the repository. Depending on Git version, different commands and options to ensure these match. An example workflow is to run git pull to get the latest commit on your current branch, and then run git submodule update to explicitly update the submodule. This should work for all versions of git which support submodules.

The following example demonstrates a shorter form that combines both commands and requires Git 2.14 or newer:

# Replaces your git pull to use both the updated code and the updated submodule
git pull --recurse-submodules

The following example demonstrates setting defaults in the config file for the repository and requires Git 2.15 or newer:

# Set this once
git config submodule.recurse true
# When configured as true, this will use both the updated code and the updated submodule
git pull
# This option will override any configuration and not update the submodule
git pull --no-recurse-submodules

These example also apply to git checkout. For more details see Git Tools Submodules: https://git-scm.com/book/en/v2/Git-Tools-Submodules

Building¶

The ERF code is dependent on AMReX, and uses the radiation model (RTE-RRTMGP) which is based on YAKL C++ implementation for heterogeneous computing infrastructure (which are all available as submodules in the ERF repo). ERF can be built using either GNU Make or CMake, however, if radiation model is activated, only CMake build system is supported.

Minimum Requirements¶

ERF requires a C++ compiler that supports the C++17 standard and a C compiler that supports the C99 standard. Building with GPU support may be done with CUDA, HIP, or SYCL. For CUDA, ERF requires versions >= 11.0. For HIP and SYCL, only the latest compilers are supported. Prerequisites for building with GNU Make include Python (>= 2.7, including 3) and standard tools available in any Unix-like environments (e.g., Perl and sed). For building with CMake, the minimal requirement is version 3.18.

Note

While ERF is designed to work with SYCL, we do not make any guarantees that it will build and run on your Intel platform.

Note

ERF was successfully compiled with the Intel compiler suite (e.g., icx version 2024.1.0). However, for older versions, it may be necessary to use reduced compiler optimization (-O1) to avoid an internal compiler error. For example, ERF was successfully compiled with icpc version 19.1.2.254, with -O2 (CMAKE_BUILD_TYPE = RelWithDebInfo) but TimeIntegration/ERF_advance_dycore.cpp had to be manually compiled with -O1. Your mileage may vary.

Paradigm¶

ERF uses the paradigm that different executables are built in different subdirectories within the Exec directory. When using gmake (see below), the user/developer should build in the directory of the selected problem. When using cmake (see below), separate executables are built for all of the problem directories listed in Exec/CMakeLists.txt. The problem directories within Exec are sorted into 1) science-relevant setups, such as ABL for modeling the atmospheric boundary layer or DensityCurrent for running the standard density current test case, etc, 2) regression tests in Exec/RegTests that are used for testing specific known aspects of the code functionality, such as boundary conditions or Rayleigh damping, and 3) tests for features under development in Exec/DevTests, such as moving terrain. There is a README in each problem directory that describes the purpose/role of that problem.

GNU Make¶

The GNU Make system is best for use on large computing facility machines and production runs. With the GNU Make implementation, the build system will inspect the machine and use known compiler optimizations explicit to that machine if possible. These explicit settings are kept up-to-date by the AMReX project.

Using the GNU Make build system involves first setting environment variables for the directories of the dependencies of ERF (AMReX, RTE-RRTMGP, and YAKL); note, RTE-RRTMGP, and YAKL are only required if running with radiation. All dependencies are provided as git submodules in ERF and can be populated by using git submodule init; git submodule update in the ERF repo, or before cloning by using git clone --recursive <erf_repo>. Although submodules of these projects are provided, they can be placed externally as long as the <REPO_HOME> environment variables for each dependency is set correctly. An example of setting the <REPO_HOME> environment variables in the user’s .bashrc is shown below:

export ERF_HOME=${HOME}/ERF
export AMREX_HOME=${ERF_HOME}/Submodules/AMReX

The GNU Make system is set up to use the path to AMReX submodule by default, so it is not necessary to set these paths explicitly, unless it is desired to do so. It is also possible to use an external version of AMReX, downloaded by running

git clone https://github.com/amrex-codes/amrex.git

in which case the AMREX_HOME environment variable must point to the location where AMReX has been downloaded, which will take precedence over the default path to the submodule. If using bash shell,

export AMREX_HOME=/path/to/external/amrex

or if using tcsh,

setenv AMREX_HOME /path/to/external/amrex

cd to the desired build directory, e.g. ERF/Exec/RegTests/IsentropicVortex/

Edit the GNUmakefile; options include

Option name	Description	Possible values	Default value
COMP	Compiler (gnu or intel)	gnu / intel	None
USE_MPI	Whether to enable MPI	TRUE / FALSE	FALSE
USE_OMP	Whether to enable OpenMP	TRUE / FALSE	FALSE
USE_CUDA	Whether to enable CUDA	TRUE / FALSE	FALSE
USE_HIP	Whether to enable HIP	TRUE / FALSE	FALSE
USE_SYCL	Whether to enable SYCL	TRUE / FALSE	FALSE
USE_NETCDF	Whether to enable NETCDF	TRUE / FALSE	FALSE
USE_HDF5	Whether to enable HDF5	TRUE / FALSE	FALSE
USE_PARTICLES	Whether to enable particles	TRUE / FALSE	FALSE
USE_WARM_NO_PRECIP	Whether to use warm moisture	TRUE / FALSE	FALSE
USE_MULTIBLOCK	Whether to enable multiblock	TRUE / FALSE	FALSE
DEBUG	Whether to use DEBUG mode	TRUE / FALSE	FALSE
PROFILE	Include profiling info	TRUE / FALSE	FALSE
TINY_PROFILE	Include tiny profiling info	TRUE / FALSE	FALSE
COMM_PROFILE	Include comm profiling info	TRUE / FALSE	FALSE
TRACE_PROFILE	Include trace profiling info	TRUE / FALSE	FALSE

Note

Do not set both USE_OMP and USE_CUDA to true.

Information on using other compilers can be found in the AMReX documentation at https://amrex-codes.github.io/amrex/docs_html/BuildingAMReX.html .

Make the executable by typing
```
make
```
The name of the resulting executable (generated by the GNUmake system) encodes several of the build characteristics, including dimensionality of the problem, compiler name, and whether MPI and/or OpenMP were linked with the executable. Thus, several different build configurations may coexist simultaneously in a problem folder. For example, the default build in ERF/Exec/RegTests/IsentropicVortex will look like ERF3d.gnu.MPI.ex, indicating that this is a 3-d version of the code, made with COMP=gnu, and USE_MPI=TRUE.

Job info¶

The build information can be accessed by typing

./ERF*ex --describe

in the directory where the executable has been built.

CMake¶

CMake is often preferred by developers of ERF; CMake allows for building as well as easy testing and verification of ERF through the use of CTest which is included in CMake.

Compiling with CMake involves an additional configure step before using the make command and it is expected that the user has cloned the ERF repo with the --recursive option or performed git submodule init; git submodule update in the ERF repo to populate its submodules.

ERF provides example scripts for CMake configuration in the /path/to/ERF/Build directory. Once the CMake configure step is done, the make command will build the executable.

An example CMake configure command to build ERF with MPI is listed below:

cmake -DCMAKE_BUILD_TYPE:STRING=Release \
      -DERF_ENABLE_MPI:BOOL=ON \
      -DCMAKE_CXX_COMPILER:STRING=mpicxx \
      -DCMAKE_C_COMPILER:STRING=mpicc \
      -DCMAKE_Fortran_COMPILER:STRING=mpifort \
      .. && make

Typically, a user will create a build directory in the project directory and execute the configuration from said directory (cmake <options> ..) before building. Note that CMake is able to generate makefiles for the Ninja build system as well which will allow for faster building of the executable(s).

Analogous to GNU Make, the list of cmake directives is as follows:

Option name

Description

Possible values

Default value

CMAKE_BUILD_TYPE

Whether to use DEBUG

Release / Debug

Release

ERF_ENABLE_MPI

Whether to enable MPI

TRUE / FALSE

FALSE

ERF_ENABLE_OPENMP

Whether to enable OpenMP

TRUE / FALSE

FALSE

ERF_ENABLE_CUDA

Whether to enable CUDA

TRUE / FALSE

FALSE

ERF_ENABLE_HIP

Whether to enable HIP

TRUE / FALSE

FALSE

ERF_ENABLE_SYCL

Whether to enable SYCL

TRUE / FALSE

FALSE

ERF_ENABLE_NETCDF

Whether to enable NETCDF

TRUE / FALSE

FALSE

ERF_ENABLE_HDF5

Whether to enable HDF5

TRUE / FALSE

FALSE

ERF_ENABLE_PARTICLES

Whether to enable particles

TRUE / FALSE

FALSE

ERF_ENABLE_WARM_NO_PRECIP

Whether to use warm moisture

TRUE / FALSE

FALSE

ERF_ENABLE_MULTIBLOCK

Whether to enable multiblock

TRUE / FALSE

FALSE

ERF_ENABLE_RADIATION

Whether to enable radiation

TRUE / FALSE

FALSE

ERF_ENABLE_TESTS

Whether to enable tests

TRUE / FALSE

FALSE

ERF_ENABLE_FCOMPARE

Whether to enable fcompare

TRUE / FALSE

FALSE

Mac with CMake¶

Tested with macOS 12.7 (Monterey) using cmake (3.27.8), open-mpi (5.0.0), and pkg-config (0.29.2) installed with the homebrew package manager. HDF5 and NetCDF will be compiled from source. The instructions below should be version agnostic.

HDF5 (tested with v1.14.3)

Download latest source package from hdfgroup.org
Extract source code tar xzf hdf5-<version>.tar.gz
Create build directory cd hdf5-<version> && mkdir build && cd build
Configure for your system ../configure --prefix=/usr/local --enable-parallel
Build make -j8 and sudo make install

NetCDF (tested with v4.9.2)

Download latest source package from ucar.edu
(Optional) install Zstd compression library brew install zstd
Create build directory cd netcdf-c-4.9.2 && mkdir build && cd build
Configure for your system ../configure --enable-parallel CC=mpicc CXX=mpicxx LDFLAGS="-L/opt/homebrew/Cellar/zstd/1.5.5/lib" CPPFLAGS="-I/opt/homebrew/Cellar/zstd/1.5.5/include" (omit the LDFLAGS and CPPFLAGS if you do not have Zstd installed) – note that you may encounter cmake errors if you do not have pkg-config installed
Build make -j8 and sudo make install

ERF (tested with commit 40e64ed35ebc080ad61d08aea828330dfbdbc162)

Get latest source code git clone --recursive git@github.com:erf-model/ERF.git
Create build directory cd ERF && mkdir MyBuild && cd MyBuild
Configure with cmake and build

cmake -DCMAKE_INSTALL_PREFIX:PATH=./install \
   -DCMAKE_CXX_COMPILER:STRING=mpicxx \
   -DCMAKE_C_COMPILER:STRING=mpicc \
   -DCMAKE_Fortran_COMPILER:STRING=mpifort \
   -DCMAKE_BUILD_TYPE:STRING=RelWithDebInfo \
   -DERF_DIM:STRING=3 \
   -DERF_ENABLE_MPI:BOOL=ON \
   -DERF_ENABLE_TESTS:BOOL=ON \
   -DERF_ENABLE_FCOMPARE:BOOL=ON \
   -DERF_ENABLE_DOCUMENTATION:BOOL=OFF \
   -DERF_ENABLE_NETCDF:BOOL=ON \
   -DERF_ENABLE_HDF5:BOOL=ON \
   -DCMAKE_EXPORT_COMPILE_COMMANDS:BOOL=ON \
   .. && make -j8

Perlmutter (NERSC)¶

Recall the GNU Make system is best for use on large computing facility machines and production runs. With the GNU Make implementation, the build system will inspect the machine and use known compiler optimizations explicit to that machine if possible. These explicit settings are kept up-to-date by the AMReX project.

For Perlmutter at NERSC, look at the general instructions for building ERF using GNU Make, and then you can initialize your environment by loading these modules:

module load PrgEnv-gnu
module load cudatoolkit

Then build ERF as, for example (specify your own path to the AMReX submodule in ERF/Submodules/AMReX):

make -j 4 COMP=gnu USE_MPI=TRUE USE_OMP=FALSE USE_CUDA=TRUE AMREX_HOME=/global/u2/d/dwillcox/dev-erf/ERF/Submodules/AMReX

Finally, you can prepare your SLURM job script, using the following as a guide:

#!/bin/bash

## specify your allocation (with the _g) and that you want GPU nodes
#SBATCH -A m4106_g
#SBATCH -C gpu

## the job will be named "ERF" in the queue and will save stdout to erf_[job ID].out
#SBATCH -J ERF
#SBATCH -o erf_%j.out

## set the max walltime
#SBATCH -t 10

## specify the number of nodes you want
#SBATCH -N 2

## we use the same number of MPI ranks per node as GPUs per node
#SBATCH --ntasks-per-node=4
#SBATCH --gpus-per-node=4
#SBATCH --gpu-bind=none

# pin to closest NIC to GPU
export MPICH_OFI_NIC_POLICY=GPU

# use GPU-aware MPI
#GPU_AWARE_MPI=""
GPU_AWARE_MPI="amrex.use_gpu_aware_mpi=1"

# the -n argument is (--ntasks-per-node) * (-N) = (number of MPI ranks per node) * (number of nodes)
# set ordering of CUDA visible devices inverse to local task IDs for optimal GPU-aware MPI
srun -n 8 --cpus-per-task=32 --cpu-bind=cores bash -c "
  export CUDA_VISIBLE_DEVICES=\$((3-SLURM_LOCALID));
  ./ERF3d.gnu.MPI.CUDA.ex inputs_wrf_baseline max_step=100 ${GPU_AWARE_MPI}" \
> test.out

To submit your job script, do sbatch [your job script] and you can check its status by doing squeue -u [your username].

Running¶

The input file specified on the command line is a free-format text file, one entry per row, that specifies input data processed by the AMReX ParmParse module.

This file needs to be specified along with the executable as an argv option, for example:

mpirun -np 64 ./ERF3d.xxx.yyy.ex inputs

Also, any entry that can be specified in the inputs file can also be specified on the command line; values specified on the command line override values in the inputs file, e.g.:

mpirun -np 64 ./ERF3d.gnu.DEBUG.MPI.ex inputs amr.restart=chk0030 erf.use_gravity=true

See Inputs for details on run-time options that can be specified. If running on a Mac and getting errors like SIGILL Invalid, privileged, or ill-formed instruction, see the note on that page about runtime error-checking options.

# To be added later #.. include:: tutorials.rst

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.

An example cmake configure command performed in the Build directory in ERF is shown below with options relevant to the testing suite:

cmake -DCMAKE_INSTALL_PREFIX:PATH=./install \
      -DCMAKE_BUILD_TYPE:STRING=Release \
      -DERF_ENABLE_MPI:BOOL=ON \
      -DCMAKE_CXX_COMPILER:STRING=mpicxx \
      -DCMAKE_C_COMPILER:STRING=mpicc \
      -DCMAKE_Fortran_COMPILER:STRING=mpifort \
      -DERF_ENABLE_FCOMPARE:BOOL=ON \
      -DERF_ENABLE_TESTS:BOOL=ON \
      -DERF_USE_CPP:BOOL=ON \
      ..

While performing a cmake -LAH .. command will give descriptions of every option for the CMake project. Descriptions of particular options regarding the testing suite are listed below:

ERF_ENABLE_FCOMPARE – builds the fcompare utility from AMReX as well as the executable(s), to allow for testing differences between plot files

ERF_ENABLE_TESTS – enables the base level regression test suite that will check whether each test will run its executable to completion successfully

Building the Tests¶

Once the user has performed the CMake configure step, the make command will build every executable required for each test. In this step, it is highly beneficial for the user to use the -j option for make to build source files in parallel.

Running the Tests¶

Once the test executables are built, CTest also creates working directories for each test within the Build directory where plot files will be output, etc. This directory is analogous to the source location of the tests in Tests/test_files.

To run the test suite, run ctest in the Build directory. CTest will run the tests and report their exit status. Useful options for CTest are -VV which runs in a verbose mode where the output of each test can be seen. -R where a regex string can be used to run specific sets of tests. -j where CTest will bin pack and run tests in parallel based on how many processes each test is specified to use and fit them into the amount of cores available on the machine. -L where the subset of tests containing a particular label will be run. Output for the last set of tests run is available in the Build directory in Tests/Temporary/LastTest.log.

Adding Tests¶

Developers are encouraged to add tests to ERF and in this section we describe how the tests are organized in the CTest framework. The locations (relative to the ERF code base) of the tests are in Tests. To add a test, first create a problem directory with a name in Exec/RegTests/<prob_name> (for problems to be used as regression tests) or Exec/DevTests/<prob_name> (for problems testing features under development), depending on which type of test is being added. Prepare a suitable input file. As an example, the TaylorGreenVortex problem with input file Exec/RegTests/TaylorGreenVortex/inputs_ex solves a simple advection-diffusion problem. The corresponding regression tests are driven by the input files Tests/test_files/TaylorGreenAdvecting/TaylorGreenAdvecting.i and Tests/test_files/TaylorGreenAdvectingDiffusing/TaylorGreenAdvectingDiffusing.i.

Any file in the test directory will be copied during CMake configure to the test’s working directory. The input files meant for regression test run only until a few time steps. The reference solution that the regression test will refer to should be placed in Tests/ERFGoldFiles/<test_name>. Next, edit the Exec/CMakeLists.txt and Tests/CTestList.cmake files, add the problem and the corresponding tests to the list. Note that there are different categories of tests and if your test falls outside of these categories, a new function to add the test will need to be created. After these steps, your test will be automatically added to the test suite database when doing the CMake configure with the testing suite enabled.

Inputs¶

The ERF executable reads run-time information from an inputs file which you name on the command line. This section describes the inputs which can be specified either in the inputs file or on the command line. A value specified on the command line will override a value specified in the inputs file.

Problem Geometry¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
geometry.prob_lo	physical location of low corner of the domain	Real	must be set
geometry.prob_hi	physical location of high corner of the domain	Real	must be set
geometry.is_periodic	is the domain periodic in this direction	0 if false, 1 if true	0 0 0

Examples of Usage¶

geometry.prob_lo = 0 0 0 defines the low corner of the domain at (0,0,0) in physical space.
geometry.prob_hi = 1.e8 2.e8 2.e8 defines the high corner of the domain at (1.e8,2.e8,2.e8) in physical space.
geometry.is_periodic = 0 1 0 says the domain is periodic in the y-direction only.

Domain Boundary Conditions¶

List of Parameters¶

Parameter	Definition	Default
xlo.type	boundary type of xlo face	must be set if not periodic
xhi.type	boundary type of xhi face	must be set if not periodic
ylo.type	boundary type of ylo face	must be set if not periodic
yhi.type	boundary type of yhi face	must be set if not periodic
zlo.type	boundary type of zlo face	must be set if not periodic
zhi.type	boundary type of zhi face	must be set if not periodic

Resolution¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
amr.n_cell	number of cells in each direction at the coarsest level	Integer > 0	must be set
amr.max_level	number of levels of refinement above the coarsest level	Integer >= 0	must be set
amr.ref_ratio	ratio of coarse to fine grid spacing between subsequent levels	2 / 3 / 4 (one per level)	2 for all levels
amr.ref_ratio_vect	ratio of coarse to fine grid spacing between subsequent levels	3 integers (one per dir) 2 / 3 / 4	2 for all directions
amr.regrid_int	how often to regrid	Integer > 0 (if negative, no regridding)	-1
amr.regrid_on_restart	should we regrid immediately after restarting	0 or 1	0
amr.iterate_grids	do we iterate on the grids?	true, false	true

Note: if amr.max_level = 0 then you do not need to set amr.ref_ratio or amr.regrid_int.

Examples of Usage¶

amr.n_cell = 32 64 64

would define the domain to have 32 cells in the x-direction, 64 cells in the y-direction, and 64 cells in the z-direction at the coarsest level.
amr.max_level = 2

would allow a maximum of 2 refined levels in addition to the coarse level. Note that these additional levels will only be created only if the tagging criteria are such that cells are flagged as needing refinement. The number of refined levels in a calculation must be \(\leq\) amr.max_level, but can change in time and need not always be equal to amr.max_level.
amr.ref_ratio = 2 3

would set factor 2 refinement between levels 0 and 1, and factor 3 refinement between levels 1 and 2. Note that you must have at least amr.max_level values of amr.ref_ratio (Additional values may appear in that line and they will be ignored).
amr.ref_ratio_vect = 2 4 3

would set factor {2 in x-dir, 4 in y-dir, 3 in z-dir} refinement between all adjacent levels. Note that you must specify 3 values, one for each coordinate direction.
amr.regrid_int = 2 2

tells the code to regrid every 2 steps. Thus in this example, new level-1 grids will be created every 2 level-0 time steps, and new level-2 grids will be created every 2 level-1 time steps.

Grid Stretching¶

This automatically activates erf.use_terrain. By default, the problem-specific terrain is initialized to be flat at an elevation of z=0. These inputs are used to automatically generate the staggered z levels used to calculate the grid metric transformation. Alternatively, arbitrary z levels may be specified with the erf.terrain_z_levels parameter, which should vary from 0 (at the surface) to the domain height.

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.grid_stretching_ratio	scaling factor applied to delta z at each level	Real > 1	0 (no grid stretching)
erf.initial_dz	vertical grid spacing for the cell above the bottom surface	Real > 0	must be set if grid stretching ratio is set
erf.terrain_z_levels	nominal staggered z levels	List of Real	NONE

Notes¶

If both erf.terrain_z_levels and erf.grid_stretching_ratio are specified, the simple grid stretching will be ignored.
The number of input erf.terrain_z_levels must be equal amr.n_cell in the z direction + 1.

Examples of Usage¶

erf.grid_stretching_ratio = 1.025
erf.initial_dz = 5.0

the first cell center would be at z=2.5

Regridding¶

Overview¶

The user defines how to tag individual cells at a given level for refinement. This list of tagged cells is sent to a grid generation routine, which uses the Berger-Rigoutsos algorithm to create rectangular grids that contain the tagged cells.

See Mesh Refinement for more details on how to specify regions for refinement.

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
amr.regrid_file	name of file from which to read the grids	text	no file
amr.grid_eff	grid efficiency at coarse level at which grids are created	Real > 0, < 1	0.7
amr.n_error_buf	radius of additional tagging around already tagged cells	Integer >= 0	1
amr.max_grid_size	maximum size of a grid in any direction	Integer > 0	32
amr.max_grid_size	maximum size	Integer	32
amr.blocking_factor	grid size must be a multiple of this	Integer > 0	2
amr.refine_grid_layout	refine grids more if # of processors \(>\) # of grids	0 if false, 1 if true	1

Notes¶

amr.n_error_buf, amr.max_grid_size and amr.blocking_factor can be read in as a single value which is assigned to every level, or as multiple values, one for each level
amr.max_grid_size at every level must be even
amr.blocking_factor at every level must be a power of 2
the domain size amr.n_cell must be a multiple of amr.blocking_factor at level 0
amr.max_grid_size must be a multiple of amr.blocking_factor at every level

Examples of Usage¶

amr.regrid_file = fixed_grids

In this case the list of grids at each fine level are contained in the file fixed_grids, which will be read during the gridding procedure. These grids must not violate the amr.max_grid_size criterion. The rest of the gridding procedure described below will not occur if amr.regrid_file is set.
amr.grid_eff = 0.9

During the grid creation process, at least 90% of the cells in each grid at the level at which the grid creation occurs must be tagged cells. Note that this is applied at the coarsened level at which the grids are actually made, and before amr.max_grid_size is imposed.
amr.max_grid_size = 64

The final grids will be no longer than 64 cells on a side at every level.
amr.max_grid_size = 64 32 16

The final grids will be no longer than 64 cells on a side at level 0, 32 cells on a side at level 1, and 16 cells on a side at level 2.
amr.blocking_factor = 32

The dimensions of all the final grids will be multiples of 32 at all levels.
amr.blocking_factor = 32 16 8

The dimensions of all the final grids will be multiples of 32 at level 0, multiples of 16 at level 1, and multiples of 8 at level 2.

Gridding and Load Balancing¶

The ERF gridding and load balancing strategy is based on that in AMReX. See the Gridding section of the AMReX documentation for details.

Simulation Time¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
max_step	maximum number of level 0 time steps	Integer >= 0	-1
start_time	starting simulation time	Real >= 0	0.0
stop_time	final simulation time	Real >= 0	Very Large

Notes¶

To control the number of time steps, you can limit by the maximum number of level-0 time steps (max_step), or the final simulation time (stop_time), or both. The code will stop at whichever criterion comes first. Note that if the code reaches stop_time then the final time step will be shortened so as to end exactly at stop_time, not pass it.

Examples of Usage¶

max_step = 1000
stop_time = 1.0

will end the calculation when either the simulation time reaches 1.0 or the number of level-0 steps taken equals 1000, whichever comes first.

Time Step¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.no_substepping	Should we turn off substepping in time?	int (0 or 1)	0
erf.cfl	CFL number for hydro	Real > 0 and <= 1	0.8
erf.fixed_dt	set level 0 dt as this value regardless of cfl or other settings	Real > 0	unused if not set
erf.fixed_fast_dt	set fast dt as this value	Real > 0	only relevant if use_native_mri is true
erf.fixed_mri_dt_ratio	set fast dt as slow dt / this ratio	even int > 0	only relevant if no_substepping is 0
erf.init_shrink	factor by which to shrink the initial dt	Real > 0 and <= 1	1.0
erf.change_max	factor by which dt can grow in subsequent steps	Real >= 1	1.1

Notes¶

The time step controls work somewhat differently depending on whether one is using acoustic substepping in time; this is determined by the value of no_substepping.
If erf.no_substepping = 1 there is only one time step to be calculated, and fixed_fast_dt and fixed_mri_dt_ratio are not used.
- If erf.fixed_dt is also specified, the timestep will be set to fixed_dt.
- If erf.fixed_dt is not specified, the timestep will be computed using the CFL condition for compressible flow. If erf.cfl is specified, that CFL value will be used. If not, the default value will be used.
If erf.no_substepping = 0 we must determine both the slow and fast timesteps.
- If erf.fixed_dt is specified, the slow timestep will be set to fixed_dt.
- If erf.fixed_dt is not set, the slow timestep will be computed using the CFL condition for incompressible flow. If erf.cfl is specified, that CFL value will be used. If not, the default value will be used.
- There are several consistency checks before the fast timestep is computed. Specifically, if any of the following are true the code will abort while reading the inputs.
  - 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
- 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:
  - If erf.fixed_fast_dt is specified, the fast timestep will be set to fixed_fast_dt.
  - If erf.fixed_mri_dt_ratio is specified and erf.fixed_fast_dt is not specified, the fast timestep will be the slow timestep divided by fixed_mri_dt_ratio.
  - If neither erf.fixed_mri_dt_ratio nor erf.fixed_fast_dt is specified, then the fast timestep will be computed using the CFL condition for compressible flow, then adjusted (reduced if necessary) as above so that the ratio of slow timestep to fine timestep is an even integer. If erf.cfl is specified, that CFL value will be used. If not, the default value will be used.

Examples of Usage of Additional Parameters¶

erf.init_shrink = 0.01

sets the initial slow time step to 1% of what it would be otherwise. Note that if erf.init_shrink \(\neq 1\) and fixed_dt is specified, then the first time step will in fact be erf.init_shrink * erf.fixed_dt.
erf.change_max = 1.1

allows the slow time step to increase by no more than 10% in this case. Note that the time step can shrink by any factor; this only controls the extent to which it can grow.

Restart Capability¶

See Checkpoint / Restart for how to control the checkpoint/restart capability.

PlotFiles¶

See Plotfiles for how to control the types and frequency of plotfile generation.

Screen Output¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
amr.v	verbosity of Amr.cpp	0 or 1	0
erf.v	verbosity of ERF.cpp	0 or 1	0
erf.sum_interval	if \(> 0,\) how often (in level-0 time steps) to compute and print integral quantities	Integer	-1

Examples of Usage¶

erf.sum_interval = 2

if erf.sum_interval \(> 0\) then the code computes and prints certain integral quantities, such as total mass, momentum and energy in the domain every erf.sum_interval level-0 steps. In this example the code will print these quantities every two coarse time steps. The print statements have the form

TIME= 1.91717746 MASS= 1.792410279e+34

for example. If this line is commented out then it will not compute and print these quantities.

Diagnostic Outputs¶

If erf.v is set then one or more additional output files may be requested. These include (1) a surface time history file, (2) a history of mean profiles, (3) a history of vertical flux profiles (i.e., variances and covariances), and (4) a history of modeled subgrid stresses. The number of specified output filenames dictates the level of output. E.g., specifying 3 filenames will give outputs (1), (2), and (3). Data files are only written if erf.profile_int > 0. This output functionality has not been implemented for terrain.

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.datalog	Output filename(s)	Up to four strings	NONE
erf.profile_int	Interval (number) of steps between ouputs	Integer	-1
erf.interp_profiles_to_cc	Interpolate all outputs to cell centers	Boolean	true

By default, all profiles are planar-averaged quantities \(\langle\cdot\rangle\) interpolated to cell centers. Setting erf.interp_profiles_to_cc = false will keep vertically staggered quantities on z faces (quantities already at cell centers or on x/y faces will remain at those locations). Note that all output quantities–whether cell-centered or face-centered–will be output on the staggered grid. The user should discard the highest z level (corresponding to the z-dir amr.n_cell + 1) for cell-centered quantities. Staggered quantities are indicated below.

The requested output files have the following columns:

Surface time history
1. Time (s)
2. Friction velocity, \(u_*\) (m/s)
3. Surface-layer potential temperature scale, \(\theta_*\) (K)
4. Obukhov length, \(L\) (m)
Mean flow profiles
1. Time (s)
2. Height (m)
3. X-velocity, \(\langle u \rangle\) (m/s)
4. Y-velocity, \(\langle v \rangle\) (m/s)
5. Z-velocity, \(\langle w \rangle\) (m/s) – staggered
6. Dry air density, \(\langle \rho \rangle\) (kg/m3)
7. Total (moist) potential temperature, \(\langle \theta \rangle\) (K)
8. Turbulent kinetic energy (TKE), \(\langle k \rangle\) (m2/s2) for the subgrid model
Vertical flux profiles
1. Time (s)
2. Height (m)
3. X-velocity variance, \(\langle u^\prime u^\prime \rangle\) (m2/s2)
4. X,Y-velocity covariance, \(\langle u^\prime v^\prime \rangle\) (m2/s2)
5. X,Z-velocity covariance, \(\langle u^\prime w^\prime \rangle\) (m2/s2) – staggered
6. Y-velocity variance, \(\langle v^\prime v^\prime \rangle\) (m2/s2)
7. Y,Z-velocity covariance, \(\langle v^\prime w^\prime \rangle\) (m2/s2) – staggered
8. Z-velocity variance, \(\langle w^\prime w^\prime \rangle\) (m2/s2) – staggered
9. X-direction heat flux, \(\langle u^\prime \theta^\prime \rangle\) (K m/s)
10. Y-direction heat flux, \(\langle v^\prime \theta^\prime \rangle\) (K m/s)
11. Z-direction heat flux, \(\langle w^\prime \theta^\prime \rangle\) (K m/s) – staggered
12. Temperature variance, \(\langle \theta^\prime \theta^\prime \rangle\) (K m/s)
13. X-direction turbulent transport of TKE, \(\langle u_i^\prime u_i^\prime u^\prime \rangle\) (m3/s3) – Note: \(u_i u_i = uu + vv + ww\)
14. Y-direction turbulent transport of TKE, \(\langle u_i^\prime u_i^\prime v^\prime \rangle\) (m3/s3)
15. Z-direction turbulent transport of TKE, \(\langle u_i^\prime u_i^\prime w^\prime \rangle\) (m3/s3) – staggered
16. X-direction pressure transport of TKE, \(\langle p^\prime u^\prime \rangle\) (m3/s3)
17. Y-direction pressure transport of TKE, \(\langle p^\prime v^\prime \rangle\) (m3/s3)
18. Z-direction pressure transport of TKE, \(\langle p^\prime w^\prime \rangle\) (m3/s3) – staggered
Modeled subgrid-scale (SGS) profiles
1. SGS stress tensor component, \(\tau_{11}\) (m2/s2)
2. SGS stress tensor component, \(\tau_{12}\) (m2/s2)
3. SGS stress tensor component, \(\tau_{13}\) (m2/s2) – staggered
4. SGS stress tensor component, \(\tau_{22}\) (m2/s2)
5. SGS stress tensor component, \(\tau_{23}\) (m2/s2) – staggered
6. SGS stress tensor component, \(\tau_{33}\) (m2/s2)
7. SGS heat flux, \(\tau_{\theta w}\) (K m/s)
8. SGS turbulence dissipation, \(\epsilon\) (m2/s3)

Advection Schemes¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.dycore_horiz_adv_type	Horizontal advection type for dycore vars	see below	Upwind_3rd
erf.dycore_vert_adv_type	Vertical advection type for dycore vars	see below	Upwind_3rd
erf.dryscal_horiz_adv_type	Horizontal advection type for dry scalars	see below	Upwind_3rd
erf.dryscal_vert_adv_type	Vertical advection type for dry scalars	see below	Upwind_3rd
erf.moistscal_horiz_adv_type	Horizontal advection type for moist scalars	see below	WENO3
erf.moistscal_vert_adv_type	Vertical advection type for moist scalars	see below	WENO3
erf.use_efficient_advection	Use efficient advection scheme for scalars	true/false	false

The allowed advection types for the dycore variables are “Centered_2nd”, “Upwind_3rd”, “Blended_3rd4th”, “Centered_4th”, “Upwind_5th”, “Blended_5th6th”, and “Centered_6th”.

The allowed advection types for the dry and moist scalars are “Centered_2nd”, “Upwind_3rd”, “Blended_3rd4th”, “Centered_4th”, “Upwind_5th”, “Blended_5th6th”, “Centered_6th” and in addition, “WENO3”, “WENOZ3”, “WENOMZQ3”, “WENO5”, and “WENOZ5.”

Note: if using WENO schemes, the horizontal and vertical advection types must be set to the same string.

The efficient advection schemes for dry and moist scalars exploit the substages of the time advancing RK3 scheme by using lower order schemes in the first two substages and the solver’s choice of scheme in the final stage. Based on CPU-only runtimes on Perlmutter for the scalar advection routine, the approximate computational savings for the scalar advection schemes are as follows when using efficient advection option: roughly 30% for Centered_4th and Centered_6th, 35% for Upwind_5th, roughly 45% for WENO5 and WENOZ5, and roughly 60% for Upwind_3rd, WENO3, WENOZ3, and WENOMZQ3.

Diffusive Physics¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.alpha_T	Diffusion coeff. for temperature	Real	0.0
erf.alpha_C	Diffusion coeff. for scalar	Real	0.0
erf.rho0_trans	Reference density to compute const. rho*Alpha	Real	1.0
erf.les_type	Using an LES model, and if so, which type?	“None”, “Smagorinsky”, “Deardorff”	“None”
erf.molec_diff_type	Using molecular viscosity and diffusivity?	“None”, “Constant”, or “ConstantAlpha”	“None”
erf.dynamicViscosity	Viscous coeff. if DNS	Real	0.0
erf.Cs	Constant Smagorinsky coeff.	Real	0.0
erf.Pr_t	Turbulent Prandtl Number	Real	1.0
erf.Sc_t	Turbulent Schmidt Number	Real	1.0
erf.use_NumDiff	Use 6th order numerical diffusion	“true”, “false”	“false”
erf.NumDiffCoeff	Coefficient for 6th order numerical diffusion	Real [0.0, 1.0]	0.0

Note: in the equations for the evolution of momentum, potential temperature and advected scalars, the diffusion coefficients are written as \(\mu\), \(\rho \alpha_T\) and \(\rho \alpha_C\), respectively.

If we set erf.molec_diff_type to Constant, then

erf.dynamicViscosity is used as the value of \(\mu\) in the momentum equation, and
erf.alpha_T is multiplied by erf.rho0_trans to form the coefficient for potential temperature, and
erf.alpha_C is multiplied by erf.rho0_trans to form the coefficient for an advected scalar.

If we set erf.molec_diff_type to ConstantAlpha, then

the dynamic viscosity in the momentum equation is assumed to have the form \(\mu = \rho \alpha_M\) where \(\alpha_M\) is a momentum diffusivity constant with units of kinematic viscosity, calculated as erf.dynamicViscosity divided by erf.rho0_trans; this diffusivity is multiplied by the instantaneous local density \(\rho\) to form the coefficient in the momentum equation; and
erf.alpha_T is multiplied by the instantaneous local density \(\rho\) to form the coefficient for potential temperature, and
erf.alpha_C is multiplied by the instantaneous local density \(\rho\) to form the coefficient for an advected scalar.

PBL Scheme¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.pbl_type	Name of PBL Scheme to be used	“None”, “MYNN2.5”	“None”
erf.pbl_mynn_A1	MYNN Constant A1	Real	1.18
erf.pbl_mynn_A2	MYNN Constant A2	Real	0.665
erf.pbl_mynn_B1	MYNN Constant B1	Real	24.0
erf.pbl_mynn_B2	MYNN Constant B2	Real	15.0
erf.pbl_mynn_C1	MYNN Constant C1	Real	0.137
erf.pbl_mynn_C2	MYNN Constant C1	Real	0.75
erf.pbl_mynn_C3	MYNN Constant C3	Real	0.352
erf.pbl_mynn_C4	MYNN Constant C4	Real	0.0
erf.pbl_mynn_C5	MYNN Constant C5	Real	0.2
erf.advect_QKE	Include advection terms in QKE eqn	bool	1
erf.diffuse_QKE_3D	Include horizontal turb. diffusion terms in QKE eqn.	bool	0

Note that the MYNN2.5 scheme must be used in conjunction with a MOST boundary condition at the surface (Zlo) boundary.

If the PBL scheme is activated, it determines the turbulent diffusivity in the vertical direction. If an LES model is also specified, it determines only the horizontal turbulent diffusivity.

Right now, the QKE equation is solved if and only if the MYNN2.5 PBL model is selected. In that transport equation, it is optional to advect QKE, and to apply LES diffusive transport for QKE in the horizontal directions (the vertical component is always computed as part of the PBL scheme).

Forcing Terms¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.abl_driver_type	Type of external forcing term	None, PressureGradient GeostrophicWind	None
erf.abl_pressure_grad	Pressure gradient forcing term (only if abl.driver_type = PressureGradient)	3 Reals	(0.,0.,0.)
erf.abl_geo_wind	Geostrophic forcing term (only if abl.driver_type = GeostrophicWind)	3 Reals	(0.,0.,0.)
erf.use_gravity	Include gravity in momentum update? If true, there is buoyancy	true / false	false
erf.use_coriolis	Include Coriolis forcing	true / false	false
erf.use_rayleigh_damping	Include explicit Rayleigh damping	true / false	false

In addition, custom forcings or tendencies may be defined on a problem-specific basis. This affords additional flexibility in defining the RHS source term as a function of time and/or height. Implementation entails modifying problem source code inside the Exec directory and overriding the update_*_sources() function(s).

Parameter	Definition	Acceptable Values	Default
erf.custom_forcing_uses_primitive_vars	User-defined source terms set the tendency of primitive variables instead of conserved quantities (rho*prim_var)	true or false	false
erf.add_custom_rhotheta_forcing	Apply the user-defined temperature source term	true or false	false

Initialization¶

ERF can be initialized in different ways. These are listed below:

Custom initialization:
Several problems under Exec are initialized in a custom manner. The state and velocity components are specific to the problem. These problems are meant for demonstration and do not include any terrain or map scale factors.
Initialization using a NetCDF file:
Problems in ERF can be initialized using a NetCDF file containing the mesoscale data. The state and velocity components of the ERF domain are ingested from the mesoscale data. This is a more realistic problem with real atmospheric data used for initialization. The typical filename used for initialization is wrfinput_d01, which is the outcome of running ideal.exe or real.exe of the WPS/WRF system. These problems are run with both terrain and map scale factors.
Initialization using an input_sounding file:
Problems in ERF can be initialized using an input_sounding file containing the vertical profile. This file has the same format as used by ideal.exe executable in WRF. Using this option for initialization, running ideal.exe and reading from the resulting wrfinput_d01 file are not needed. This option is used for initializing ERF domain to a horizontally homogeneous mesoscale state and does not include terrain or map scale factors.

In addition, there is a run-time option to project the initial velocity field to make it divergence-free. To take advantage of this option, the code must be built with USE_POISSON_SOLVE = TRUE in the GNUmakefile if using gmake, or with -DERF_ENABLE_POISSON_SOLVE:BOOL=ON in the cmake.sh file if using cmake.

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.init_type	Initialization type	“custom”, “ideal”, “real”, “input_sounding”	“custom”
erf.input_sounding_file	Path to WRF-style input sounding file	String	“input_sounding”
erf.init_sounding_ideal	Perform initialization like WRF’s ideal.exe	true or false	false
erf.use_real_bcs	If init_type is real or metgrid, do we want to use these bcs?	true or false	true if if init_type is real or metgrid; else false
erf.nc_init_file	NetCDF file with initial mesoscale data	String	NONE
erf.nc_bdy_file	NetCDF file with mesoscale data at lateral boundaries	String	NONE
erf.project_initial_velocity	project initial velocity?	Integer	1

Notes¶

If erf.init_type = ideal, the problem is initialized with mesoscale data contained in a NetCDF file, provided via erf.nc_init_file. The mesoscale data are horizontally homogeneous, i.e., there is variation only in the vertical direction.

If erf.init_type = real, the problem is initialized with mesoscale data contained in a NetCDF file, provided via erf.nc_init_file. The mesoscale data are realistic with variation in all three directions. In addition, the lateral boundary conditions must be supplied in a NetCDF files specified by erf.nc_bdy_file = wrfbdy_d01

If erf.init_type = input_sounding, a WRF-style input sounding is read from erf.input_sounding_file. This text file includes any set of levels that goes at least up to the model top height. The first line includes the surface pressure [hPa], potential temperature [K], and water vapor mixing ratio [g/kg]. Each subsequent line has five input values: height [m above sea level], dry potential temperature [K], vapor mixing ratio [g/kg], x-direction wind component [m/s], and y-direction wind component [m/s]. Please pay attention to the units of pressure and mixing ratio. If erf.init_sounding_ideal = true, then moist and dry conditions throughout the air column are determined by integrating the hydrostatic equation from the surface.

If erf.init_type = custom or erf.init_type = input_sounding, erf.nc_init_file and erf.nc_bdy_file do not need to be set.

Setting erf.project_initial_velocity = 1 will have no effect if the code is not built with ERF_USE_POISSON_SOLVE defined.

Map Scale Factors¶

Map scale factors are always present in the evolution equations, but the values default to one unless specified in the initialization when erf.init_type = real.

There is an option to test the map scale factors by setting erf.test_mapfactor = true; this arbitrarily sets the map factors to 0.5 in order to test the implementation.

Terrain¶

ERF allows the use to specify whether terrain-fitted coordinates should be used by setting erf.use_terrain (default false). If terrain-fitted coordinates are chosen, they are defined to be static (default) or moving by setting erf.terrain_type. If using terrain, the user also has the option to specify one of three methods for defining how the terrain-fitted coordinates given the topography:

Basic Terrain Following (BTF):
The influence of the terrain decreases linearly with height.
Smoothed Terrain Following (STF):
Small-scale terrain structures are progressively smoothed out of the coordinate system as height increases.
Sullivan Terrain Following (name TBD):
The influence of the terrain decreases with the cube of height.

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.use_terrain	use terrain-fitted coordinates?	true / false	false
erf.terrain_type	static or moving?	Static / Moving	Static
erf.terrain_smoothing	specify terrain following	0, 1, 2	0

Examples of Usage¶

erf.terrain_smoothing = 0
BTF is used when generating the terrain following coordinate.
erf.terrain_smoothing = 1
STF is used when generating the terrain following coordinate.
erf.terrain_smoothing = 2
Sullivan TF is used when generating the terrain following coordinate.

Moisture¶

ERF has several different moisture models. The models that are currently implemented are Eulerian models; however, ERF has the capability for Lagrangian models when compiled with particles.

The following run-time options control how the full moisture model is used.

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.moisture_model	Name of moisture model	“SAM”, “Kessler”, “FastEddy”	“Null”
erf.do_cloud	use basic moisture model	true / false	true
erf.do_precip	include precipitation in treatment of moisture	true / false	true

Runtime Error Checking¶

Through AMReX functionality, ERF supports options to raise errors when NaNs, division by zero, and overflow errors are detected. These checks are activated at runtime using the input parameters below.

Note

When running on Macs using the Apple-Clang compilers with optimization (DEBUG = FALSE in the GNUmakefile), these checks may lead to false positives due to optimizations performed by the compiler and the flags should be turned off. It is still possible to run with these error checks with Apple-Clang debug builds.

List of Parameters¶

Parameter	Definition	Acceptable Values
amrex.fpe_trap_invalid	Raise errors for NaNs	0 / 1
amrex.fpe_trap_zero	Raise errors for divide by zero	0 / 1
amrex.fpe_trap_overflow	Raise errors for overflow	0 / 1

Best Practices¶

Please note this section is a work in progress and aims to offer guidelines but, in general, your mileage may vary.

Large-Eddy Simulations¶

Advection Scheme
- A WRF-like configuration is generally robust, but may under-resolve turbulence features.
```
erf.dycore_horiz_adv_type  = "Upwind_5th"
erf.dycore_vert_adv_type   = "Upwind_3rd"
erf.dryscal_horiz_adv_type = "Upwind_5th"
erf.dryscal_vert_adv_type  = "Upwind_3rd"
```
- Centered difference schemes will generally give some non-physical numerical noise, clearly visible in the free atmosphere, but may also better resolve turbulence features. With Centered_2nd, the simulation may remain numerically stable but without any upwinding or numerical diffusion, these results should be carefully interpreted.
- For higher-order central differencing alone (i.e., without any added upwinding), at least 5% numerical diffusion should be included to stabilize the solution; this was tested with Centered_6th. Note that this does not necessarily kill the numerical noise and is only for numerical stability. These options are identical to WRF’s diff_6th_opt (default: off) and diff_6th_factor (default: 12%) options.
```
erf.use_NumDiff  = true
erf.NumDiffCoeff = 0.05
```
Time Integration
- Split timestepping offers some computational cost savings but still does not allow you to run with an incompressible time-step size.
- The acoustic CFL should be less than 0.5, with 4–6 fast timesteps (substeps) according to WRF best practices.
```
erf.fixed_dt           = 0.06  # slow timestep

# These are equivalent and result in a fixed fast timestep size
#   if dx=10, speed of sound ~ 350 m/s
erf.fixed_fast_dt      = 0.01  # ==> CFL~0.35
#erf.fixed_mri_dt_ratio = 6

# Alternatively, let ERF chose the fast timestep
#erf.cfl                = 0.5
```
- We note that ERF LESs with up to 10 fast timesteps have successfully been run but your mileage may vary.

Single-Column Model¶

Currently, ERF does not have the ability to run a true single-column model (SCM). The grid size in the lateral directions must have a minimum number of cells. This will give comparable results, e.g.:
```
geometry.prob_extent = 25  25  400
amr.n_cell           =  4   4   64
geometry.is_periodic =  1   1    0
```
An SCM was successfully run with third-order advection in the horizontal and vertical.

Prognostic Equations (Dry)¶

The following partial differential equations governing dry compressible flow are solved in ERF for mass, momentum, potential temperature, and scalars:

\[ \begin{align}\begin{aligned}\frac{\partial \rho}{\partial t} &= - \nabla \cdot (\rho \mathbf{u}),\\\frac{\partial (\rho \mathbf{u})}{\partial t} &= - \nabla \cdot (\rho \mathbf{u} \mathbf{u}) - \nabla p^\prime + \delta_{i,3}\mathbf{B} - \nabla \cdot \tau + \mathbf{F},\\\frac{\partial (\rho \theta)}{\partial t} &= - \nabla \cdot (\rho \mathbf{u} \theta) + \nabla \cdot ( \rho \alpha_{T}\ \nabla \theta) + F_{\rho \theta},\\\frac{\partial (\rho C)}{\partial t} &= - \nabla \cdot (\rho \mathbf{u} C) + \nabla \cdot (\rho \alpha_{C}\ \nabla C)\end{aligned}\end{align} \]

where

\(\tau\) is the viscous stress tensor,

\[\tau_{ij} = -2\mu \sigma_{ij},\]

with \(\sigma_{ij} = S_{ij} -D_{ij}\) being the deviatoric part of the strain rate, and

\[S_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right), \hspace{24pt} D_{ij} = \frac{1}{3} S_{kk} \delta_{ij} = \frac{1}{3} (\nabla \cdot \mathbf{u}) \delta_{ij},\]

\(\mathbf{F}\) and \(F_{\rho \theta}\) are the forcing terms described in Physical Forcings,
\(\mathbf{g} = (0,0,-g)\) is the gravity vector,
the potential temperature \(\theta\) is defined from temperature \(T\) and pressure \(p\) as

\[\theta = T \left( \frac{p_0}{p} \right)^{R_d / c_p}.\]

pressure and density are defined as perturbations from a hydrostatically stratified background state, i.e.

\[p = \overline{p}(z) + p^\prime \hspace{24pt} \rho = \overline{\rho}(z) + \rho^\prime\]

with

\[\frac{d \overline{p}}{d z} = - \overline{\rho} g\]

Assumptions¶

The assumptions involved in deriving these equations from first principles are:

Continuum behavior
Ideal gas behavior (\(p = \rho R_d T\)) with constant specific heats (\(c_p,c_v\))
Constant mixture molecular weight (therefore constant \(R_d\))
Viscous heating is negligible
No chemical reactions, second order diffusive processes or radiative heat transfer
Newtonian viscous stress with no bulk viscosity contribution (i.e., \(\kappa S_{kk} \delta_{ij}\))
Depending on the simulation mode, the transport coefficients \(\mu\), \(\rho\alpha_C\), and \(\rho\alpha_T\) may correspond to the molecular transport coefficients, turbulent transport coefficients computed from an LES or PBL model, or a combination. See the sections on DNS vs. LES modes and PBL schemes for more details.

Diagnostic Relationships¶

In order to close the above prognostic equations, a relationship between the pressure and the other state variables must be specified. This is obtained by re-expressing the ideal gas equation of state in terms of \(\theta\):

\[p = \left( \frac{\rho R_d \theta}{p_0^{R_d / c_p}} \right)^\gamma = p_0 \left( \frac{\rho R_d \theta}{p_0} \right)^\gamma\]

Nomenclature¶

Here \(\rho, T, \theta\), and \(p\) are the density, temperature, potential temperature and pressure, respectively; these variables are all defined at cell centers. \(C\) is an advected quantity, i.e., a tracer, also defined at cell centers. \(\mathbf{u}\) and \((\rho \mathbf{u})\) are the velocity and momentum, respectively, and are defined on faces.

\(R_d\) and \(c_p\) are the gas constant and specific heat capacity for dry air respectively, and \(\gamma = c_p / (c_p - R_d)\) . \(p_0\) is a reference value for pressure.

Prognostic Equations (Moist)¶

Model 1: Warm Moisture with no Precipitation¶

With this model, which is analogous to that in FASTEddy, we consider a mixture of dry air \(\rho_d\) and nonprecipitating water vapor \(\rho_v\), assumed to be a perfect ideal gas with constant heat capacities \(C_{vd}\), \(C_{vv}\), \(C_{pd}\), \(C_{pv}\), and (non-precipitating) cloud water \(\rho_c\).

Neglecting the volume occupied by all water not in vapor form, we have

\[p = p_d + p_v = \rho_d R_d T + \rho_v R_v T\]

where \(p_d\) and \(p_v\) are the partial pressures of dry air and water vapor, respectively, and \(R_d\) and \(R_v\) are the gas constants for dry air and water vapor, respectively.

We define the mixing ratio of each moist component, \(q_s\), as the mass density of species \(s\) relative to the density of dry air, i.e. \(q_s = \frac{\rho_s}{\rho_d}\).

Governing Equations¶

The governing equations for this model are

\[ \begin{align}\begin{aligned}\frac{\partial \rho_d}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u})\\\frac{\partial (\rho_d \mathbf{u})}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \mathbf{u}) - \frac{1}{1 + q_v + q_c} ( \nabla p^\prime + \delta_{i,3}\mathbf{B} ) - \nabla \cdot \tau + \mathbf{F}\\\frac{\partial (\rho_d \theta_d)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \theta_d) + \nabla \cdot ( \rho_d \alpha_{T}\ \nabla \theta_d) + \frac{\theta_d L_v}{T_d C_{pd}} f_{cond}\\\frac{\partial (\rho_d q_v)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} q_vi) + \nabla \cdot (\rho_d \alpha \nabla q_v) - f_{cond}\\\frac{\partial (\rho_d q_c)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} q_c) + \nabla \cdot (\rho_d \alpha \nabla q_c) + f_{cond}\end{aligned}\end{align} \]

Here \(L_v\) is the latent heat of vaporization, \(\theta_d\) is the (dry) potential temperature \(\mathbf{B}\) is the buoyancy force, which is defined in Buoyancy.

The pressure perturbation is computed as

\[p^\prime = p_0 \left( \frac{R_d \rho_d \theta_m}{p_0} \right)^\gamma - p_0\]

where \(\gamma = C_{pd} / C_{vd}\) and

\[\theta_m = \theta_d (1 + \frac{R_v}{R_d} q_v)\]

is the moist potential temperature. We note that this is the only place \(\theta_m\) is used; we evolve \(\theta_d\) above.

Model 2: Full Moisture Including Precipitation¶

With this model, in addition to dry air \(\rho_d\) and nonprecipitating water vapor \(\rho_v\), assumed to be a perfect ideal gas with constant heat capacities \(C_{vd}\), \(C_{vv}\), \(C_{pd}\), \(C_{pv}\), we include non-precipitating condensates \(\rho_c + \rho_i\), and precipitating condensates \(\rho_p = \rho_{rain} + \rho_{snow} + \rho_{graupel}\). Here \(\rho_c\) is the density of cloud water and \(\rho_i\) is the density of cloud ice, and we define the sum of all non-precipitating moist quantities to be \(\rho_T = \rho_v + \rho_c + \rho_i\). All condensates are treated as incompressible; cloud water and ice have constant heat capacities \(C_p\) and \(C_i\), respectively.

Neglecting the volume occupied by all water not in vapor form, we have

\[p = p_d + p_v = \rho_d R_d T + \rho_v R_v T\]

where \(p_d\) and \(p_v\) are the partial pressures of dry air and water vapor, respectively, and \(R_d\) and \(R_v\) are the gas constants for dry air and water vapor, respectively.

We define the mixing ratio of each moist component, \(q_s\), as the mass density of species \(s\) relative to the density of dry air, i.e. \(q_s = \frac{\rho_s}{\rho_d}\).

We define the total potential temperature

\[\theta = \frac{\sum_s \rho_s \theta_s}{\sum_s \rho_s} \approx (\theta_d + q_v \theta_v + q_i \theta_i + q_c \theta_c).\]

and write the EOS as

\[T = \theta (\frac{p}{p_0})^\frac{R^\star}{C_p^\star}\]

or

\[p = p_0 (\frac{\Pi}{C_p^\star})^{\frac{C_p^\star}{R^\star}}\]

where \(p_0\) is the reference pressure. and

\[\Pi = C_p^\star (\frac{p}{\alpha p_0})^\frac{R^\star}{C_p^\star}\]

with \(\alpha = \frac{R^\star}{p}(\frac{p}{p_0})^\frac{R^\star}{c_p^\star} \theta\)

here, \(R^\star = R_{d} + q_v R_{v} + q_i R_{i} + q_p R_{p}\), and \(C_p^\star = C_{pd} + q_v C_{pv} + q_i C_{pi} + q_p C_{pp}\).

\(R_d\), \(R_v\), \(R_i\), and \(R_p\) are the gas constants for dry air, water vapor, cloud ice, precipitating condensates, respectively. \(C_{pd}\), \(C_{pv}\), \(C_{pi}\), and \(C_{pp}\) are the specific heats for dry air, water vapor, cloud ice, and precipitating condensates, respectively.

Governing Equations¶

We assume that all species have same average speed, Then the governing equations become

\[ \begin{align}\begin{aligned}\frac{\partial \rho_d}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} + \mathbf{F}_\rho)\\\frac{\partial (\rho_d \mathbf{u})}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \mathbf{u} + \mathbf{F}_u) - \frac{1}{1 + q_T + q_p} \nabla p^\prime - \nabla \cdot \tau + \mathbf{F} + \delta_{i,3}\mathbf{B}\\\frac{\partial (\rho_d \theta)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \theta + F_{\theta}) + \nabla \cdot ( \rho_d \alpha_{T}\ \nabla \theta) + F_Q\\\frac{\partial (\rho_d q_T)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} q_T +F_{q_{T}}) - Q\\\frac{\partial (\rho_d q_p)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} q_p + F_{q_{p}}) + Q\end{aligned}\end{align} \]

In this set of equations, the subgrid turbulent parameterization effects are included with fluxes \(F_\rho\), \(F_u\), \(F_C\), \(F_{\theta}\), \(F_{q_{T}}\), \(F_{q_{r}}\). \(\mathbf{F}\) stands for the external force, and \(Q\) and \(F_Q\) represent the mass and energy transformation of water vapor to/from water through condensation/evaporation, which is determined by the microphysics parameterization processes. \(\mathbf{B}\) is the buoyancy force, which is defined in Buoyancy.

Initialization¶

To initialize a background (base) state for a simulation with a hydrostatic atmosphere, the hydrostatic equation balancing the pressure gradient and gravity must be satisfied. This section details the procedure for the initialization of the background state. The procedure is similar for the cases with and without moisture, the only difference being that the density for the cases with moisture has to be the total density \(\rho = \rho_d(1 + q_t)\), where \(\rho_d\) is the dry density, and \(q_t\) is the total mass mixing ratio – water vapor and liquid water, instead of the dry density \(\rho_d\) for cases without moisture.

Computation of the dry density¶

We express the total pressure \(p\) as:

\[p = \rho_d R_d T_v.\]

By definition, we have:

\[T_v = \theta_m\left(\frac{p}{p_0}\right)^{R_d/C_p},\]

\[T = \theta_d\left(\frac{p}{p_0}\right)^{R_d/C_p}.\]

This gives:

\[p = \rho_d R_d \theta_m\left(\frac{p}{p_0}\right)^{R_d/C_p},\]

which, on rearranging, yields:

\[p = p_0\left(\frac{\rho_d R_d \theta_m}{p_0}\right)^\gamma.\]

To obtain \(\theta_m\), consider the density of dry air:

\[\rho_d = \frac{p - p_v}{R_d T}.\]

Substituting for \(\rho_d\) from the above equation, we get:

\[\frac{p}{T_v} = \frac{p - p_v}{T},\]

which implies:

\[\frac{T_v}{T} = \frac{p}{p - p_v} = \frac{1}{1-\left(\cfrac{p_v}{p}\right)}.\]

We also have:

\[q_v = \frac{\rho_v}{\rho_d} = \frac{p_v}{R_v T}\frac{R_d T}{p-p_v} = \frac{r p_v}{p - p_v},\]

where \(p_v\) is the partial pressure of water vapor, \(r = R_d/R_v \approx 0.622\), and \(q_v\) is the vapor mass mixing ratio—the ratio of the mass of vapor to the mass of dry air. Rearranging and using \(q_v \ll r\), we get:

\[\frac{p_v}{p} = \frac{1}{1 + \left(\cfrac{r}{q_v}\right)} \approx \frac{q_v}{r},\]

which, on substitution in the equation for \(\frac{T_v}{T}\), gives:

\[\frac{T_v}{T} = \frac{1}{1 - \left(\cfrac{q_v}{r}\right)}.\]

As \(q_v \ll 1\), a binomial expansion, ignoring higher-order terms, gives:

\[T_v = T\left(1 + \frac{R_v}{R_d}q_v\right).\]

Hence, the density of dry air is given by:

\[\rho_d = \frac{p}{R_d T_v} = \frac{p}{R_d T\left(1 + \cfrac{R_v}{R_d}q_v\right)}.\]

Initialization with a second-order integration of the hydrostatic equation¶

We have the hydrostatic equation given by

\[\frac{\partial p}{\partial z} = -\rho g,\]

where \(\rho = \rho_d(1 + q_t)\), \(\rho_d\) is the dry density, and \(q_t\) is the total mass mixing ratio – water vapor and liquid water. Using an average value of \(\rho\) for the integration from \(z = z(k-1)\) to \(z(k)\), we get

\[p(k) = p(k-1) - \frac{(\rho(k-1) + \rho(k))}{2} g\Delta z.\]

The density at a point is a function of the pressure, potential temperature, and relative humidity. The latter two quantities are computed using user-specified profiles, and hence, for simplicity, we write \(\rho(k) = f(p(k))\). Hence

\[p(k) = p(k-1) - \frac{\rho(k-1)}{2}g\Delta z - \frac{f(p(k))}{2}g\Delta z.\]

Now, we define

\[F(p(k)) \equiv p(k) - p(k-1) + \frac{\rho(k-1)}{2}g\Delta z + \frac{f(p(k))}{2}g\Delta z = 0.\]

This is a non-linear equation in \(p(k)\). Consider a Newton-Raphson iteration (where \(n\) denotes the iteration number) procedure

\[F(p+\delta p) \approx F(p) + \delta p \frac{\partial F}{\partial p} = 0,\]

which implies

\[\delta p = -\frac{F}{F'},\]

with the gradient being evaluated as

\[F' = \frac{F(p+\epsilon) - F(p)}{\epsilon},\]

and the iteration update is given by

\[p^{n+1} = p^n + \delta p.\]

For the first cell (\(k=0\)), which is at a height of \(z = \frac{\Delta z}{2}\) from the base, we have

\[p(0) = p_0 - \rho(0)g\frac{\Delta z}{2},\]

where \(p_0 = 1e5 \, \text{N/m}^2\) is the pressure at the base. Hence, we define

\[F(p(0)) \equiv p(0) - p_0 + \rho(0)g\frac{\Delta z}{2},\]

and the Newton-Raphson procedure is the same.

Buoyancy¶

ERF has three options for how to define the buoyancy force. Even in the absence of moisture these expressions are not equivalent.

Density of the mixture¶

The total density in a cell containing air, water vapor, liquid water and precipitates is given by

\[\rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V},\]

where \(m_a\) is the mass of dry air, \(m_v\) is the mass of water vapor, \(m_c\) is the mass of liquid water, and \(m_p\) is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component

\[q_i \equiv \frac{m_i}{m_a}.\]

Using this we can write

\[\rho = m_a\frac{(1 + q_v + q_c + q_p)}{V} = \rho_d(1 + q_v + q_c + q_p),\]

where \(\rho_d \equiv \cfrac{m_a}{V}\) is the density of dry air.

Type 1¶

One version of the buoyancy force is expressed simply as

\[\mathbf{B} = \rho^\prime \mathbf{g}\]

\[\rho^\prime = \rho_{total} - \rho_0\]

where the total density \(\rho_{total} = \rho_d(1 + q_v + q_c + q_p)\) is the sum of dry and moist components and \(\rho_0\) is the total density for the background state. For eg., a usual scenario is that of a background state that contains only air and vapor and no cloud water or precipitates. For such a state, the total background density \(\rho_0 = \rho_{d_0}(1 + q_{v_0})\), where \(\rho_{d_0}\) and \(q_{v_0}\) are the background dry density and vapor mixing ratio respectively. As a check, we observe that \(\rho^\prime_0 = 0\), which means that the background state is not buoyant.

Type 2¶

The second option for the buoyancy force is

\[\mathbf{B} = -\rho_0 \mathbf{g} ( 0.61 q_v^\prime - q_c^\prime - q_i^\prime - q_p^\prime + \frac{T^\prime}{\bar{T}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) )\]

To derive this expression, we define \(T_v = T (1 + 0.61 q_v − q_c − q_i - q_p)\), then we can write

\[p = \rho (R_d q_d + R_v q_v) T = \rho R_d T (1 + 0.61 q_v − q_c − q_i - q_p ) = \rho R_d T_v\]

Starting from \(p = \rho R_d T_v\) and neglecting \(\frac{p^\prime}{\bar{p}}\), we now write

\[\frac{\rho^\prime}{\overline{\rho}} = -\frac{T_v^\prime}{\overline{T_v}}\]

and define

\[T_v^\prime = T_v - \overline{T_v} \approx \overline{T} [ 0.61 q_v^\prime - (q_c^\prime + q_i^\prime + q_p^\prime)] + (T - \overline{T}) [1+ 0.61 \bar{q_v} - \bar{q_c} - \bar{q_i} - \bar{q_p} ] .\]

where we have retained only first order terms in perturbational quantities.

Then

\[\mathbf{B} = \rho^\prime \mathbf{g} = -\overline{\rho} \frac{\overline{T}}{\overline{T_v}} \mathbf{g} [ 0.61 q_v^\prime - q_c^\prime - q_i^\prime - q_p^\prime ) + \frac{T^\prime}{\overline{T_v}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ]\]

where the overbar represents a horizontal average of the current state and the perturbation is defined relative to that average.

Again keeping only the first order terms in the mass mixing ratios, we can simplify this to

\[\mathbf{B} = \rho^\prime \mathbf{g} = -\rho_0 \mathbf{g} [ 0.61 q_v^\prime - q_c^\prime + q_i^\prime + q_p^\prime + \frac{T^\prime}{\overline{T}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ]\]

We note that this reduces to Type 3 if the horizontal averages of the moisture terms are all zero.

Type 3¶

The third formulation of the buoyancy term assumes that the horizontal averages of the moisture quantities are negligible, which removes the need to compute horizontal averages of these quantities. This reduces the Type 2 expression to the following:

\[\mathbf{B} = \rho^\prime \mathbf{g} \approx -\rho_0 \mathbf{g} ( \frac{T^\prime}{\overline{T}} + 0.61 q_v - q_c - q_i - q_p)\]

We note that this version of the buoyancy force matches that given in Marat F. Khairoutdinov and David A. Randall’s paper (J. Atm Sciences, 607, 1983) if we neglect \(\frac{p^\prime}{\bar{p_0}}\).

Type 4¶

This expression for buoyancy is from khairoutdinov2003cloud and bryan2002benchmark.

\[\begin{equation} \mathbf{B} = \rho'\mathbf{g} \approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg), \end{equation}\]

The derivation follows. The total density is given by \(\rho = \rho_d(1 + q_v + q_c + q_p)\), which can be written as

\[\rho = \frac{p (1 + q_v + q_c + q_p)}{R_dT\Bigg(1 + \cfrac{R_v}{R_d}q_v\Bigg)}\]

This can be written using binomial expansion as

\[\begin{split}\begin{align*} \rho &= \frac{p}{R_dT} (1 + q_v + q_c + q_p)\Bigg(1 + \frac{R_v}{R_d}q_v\Bigg)^{-1} \\ &= \frac{p}{R_dT} (1 + q_v + q_c + q_p)\Bigg(1 - \frac{R_v}{R_d}q_v + O(q_v^2)\Bigg) \\ &= \frac{p}{R_dT}\Bigg(1 + q_v + q_c + q_p - \frac{R_v}{R_d}q_v + \text{H.O.T. such as } O(q_v^2) + O(q_vq_c)\Bigg) \\ &\approx \frac{p}{R_dT}\Bigg(1 + q_v + q_c + q_p - \frac{R_v}{R_d}q_v\Bigg) \end{align*}\end{split}\]

Taking log on both sides, we get

\[\log{\rho} = \log{p} - \log{R_d} - \log{T} + \log(1 - 0.61 q_v + q_c + q_p)\]

Taking derivative gives

\[\frac{\rho'}{\rho} = \frac{p'}{p} - \frac{T'}{T} + \frac{(-0.61 q_v' + q_c' + q_p')}{(1 - 0.61 q_v + q_c + q_p)}\]

Using \(- 0.61 q_v + q_c + q_p \ll 1\), we have

\[\frac{\rho'}{\rho} = \frac{p'}{p} - \frac{T'}{T} + (-0.61 q_v' + q_c' + q_p')\]

Since the background values of cloud water and precipitate mass mixing ratios – \(q_c\) and \(q_p\) are zero, we have \(q_c' = q_c\) and \(q_p' = q_p\). Hence, we have

\[\begin{equation} \rho'\approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg), \end{equation}\]

Microphysics model¶

Kessler Microphysics model¶

The Kessler microphysics model is a simple version of cloud microphysics which has precipitation only in the form of rain. Hence \(q_p = q_r\). Governing equations for the microphysical quantities for Kessler microphysics from gabervsek2012dry are

\[\frac{\partial q_v}{\partial t} = -C_c + E_c + E_r\]

\[\frac{\partial q_c}{\partial t} = C_c - E_c - (A_c + K_c)\]

\[\frac{\partial q_p}{\partial t} = \frac{1}{\overline{\rho}}\frac{\partial}{\partial z}(\overline{\rho}Vq_p) + (A_c + K_c) - E_r\]

\[\frac{\partial q_t}{\partial t} = \frac{\partial q_v}{\partial t} + \frac{\partial q_c}{\partial t} = E_r - (A_c + K_c)\]

where \(C_c\) is the rate of condensation of water vapor to cloud water, \(E_c\) is the rate of evaporation of cloud water to water vapor, \(A_c\) is the autoconversion of cloud water to rain, \(K_c\) is the accretion of cloud water to rain drops, \(E_r\) is the evaporation of rain to water vapor and \(F_r\) is the sedimentation of rain. The parametrization used is given in klemp1978simulation, and is given below. Note that in all the equations, \(p\) is specified in millibars and \(\overline{\rho}\) is specified in g cm \(^{-3}\). The parametrization of the source terms are given below.

Rate of condensation of water vapor/evaporation of cloud water¶

From klemp1978simulation, we have the following expressions.

If the air is not saturated, i.e. \(q_v > q_{vs}\)

\[C_c = \frac{q_v - q_{vs}}{1 + \cfrac{q_{vs}^*4093L}{C_p(T-36)^2}}\]

If the air is not saturated, i.e. \(q_v < q_{vs}\), then cloud water evaporates to water vapor at the rate

\[E_c = \frac{q_{vs} - q_v}{1 + \cfrac{q_{vs}^*4093L}{C_p(T-36)^2}}\]

Rate of autoconversion of cloud water into rain¶

The rate of autoconversion of cloud water into rain is given by

\[A_c = k_1(q_c - a)\]

where \(k_1 = 0.001\) s^-1, \(a = 0.001\) kg kg^-1.

Rate of accretion of cloud water into rain water drops¶

The rate of accretion of cloud water into rain water drops is given by

\[K_c = k_2q_cq_r^{0.875}\]

where \(k_2= 2.2\) s^-1.

The rate of evaporation of rain into water vapor¶

The rate of evaporation of rain into water vapor is given by

\[E_r = \cfrac{1}{\overline{\rho}}\cfrac{(1- q_v/q_s)C(\overline{\rho}q_r)^{0.525}}{5.4\times10^5 + 2.55\times10^6/(\overline{p}q_s)},\]

where the ventilation factor \(C = 1.6 + 124.9(\overline{\rho}q_r)^{0.2046}\).

Terminal fall velocity of rain¶

The terminal fall velocity of rain is given by

\[V = 3634(\overline{\rho}q_r)^{0.1346}\Bigg(\cfrac{\overline{\rho}}{\rho_0}\Bigg)^{-\frac{1}{2}}~\text{[cm/s]}\]

Single Moment Microphysics Model¶

The conversion rates among the moist hydrometeors are parameterized assuming that

\[\frac{\partial N_{m}}{\partial D} = n_{m}\left(D_{m}\right) = N_{0m} exp \left(-\lambda _{m} D_{m}\right)\]

where \(N_{0m}\) is the intercept parameter, \(D_{m}\) is the diameters, and

\[\lambda_{m} = (\frac{\pi \rho_{m} N_{0m}}{q_{m}\rho})^{0.25}\]

where \(\rho_{m}\) is the density of moist hydrometeors. Assuming that the particle terminal velocity

\[v_{m} \left( D_{m},p \right) = a_{m}D_{m}^{b_{m}}\left(\frac{\rho_{0}}{\rho}\right)^{0.5}\]

The total production rates including the contribution from aggregation, accretion, sublimation, melting, bergeron process, freezing and autoconversion are listed below without derivation. For details, please refer to Yuh-Lang Lin et al (J. Climate Appl. Meteor, 22, 1065, 1983) and Marat F. Khairoutdinov and David A. Randall’s (J. Atm Sciences, 607, 1983). The implementation of microphysics model in ERF is similar to the that in the SAM code (http://rossby.msrc.sunysb.edu/~marat/SAM.html)

Accretion¶

There are several different type of accretional growth mechanisms that need to be included; these describe the interaction of water vapor and cloud water with rain water.

The accretion of cloud water forms in either the dry or wet growth rate can be written as:

\[Q_{gacw} = \frac{\pi E_{GW}n_{0G}q_{c}\Gamma(3.5)}{4\lambda_{G}^{3.5}}(\frac{4g\rho G}{3C_{D}\rho})^{0.5}\]

The accretion of raindrops by accretion of cloud water is

\[Q_{racw} = \frac{\pi E_{RW}n_{0R}\alpha q_{c}\Gamma(3+b)}{4\lambda_{R}^{3+b}}(\frac{\rho_{0}}{\rho})^{1/2}\]

The bergeron Process¶

The cloud water transform to snow by deposition and rimming can be written as

\[Q_{sfw} = N_{150}\left(\alpha_{1}m_{150}^{\alpha_{2}}+\pi E_{iw}\rho q_{c}R_{150}^{2}U_{150}\right)\]

Autoconversion¶

The collision and coalescence of cloud water to from raindrops is parameterized as following

\[Q_{raut} = \rho\left(q_{c}-q_{c0}\right)^{2}\left[1.2 \times 10^{-4}+{1.569 \times 10^{-12}N_{1}/[D_{0}(q_{c}-q_{c0})]}\right]^{-1}\]

Evaporation¶

The evaporation rate of rain is

\[Q_{revp} = 2\pi(S-1)n_{0R}[0.78\lambda_{R}^{-2}+0.31S_{c}^{1/3}\Gamma[(b+5)/2]a^{1/2}\mu^{-1/2}(\frac{\rho_{0}}{\rho})^{1/4}\lambda_{R}^{(b+5)/2}](\frac{1}{\rho})(\frac{L_{v}^{2}}{K_{0}R_{w}T^{2}}+\frac{1}{\rho r_{s}\psi})^{-1}\]

Simulation Modes: DNS and LES¶

DNS¶

When running in Direct Numerical Simulation (DNS) mode, it is assumed that the mesh is fine enough to resolve the Kolmogorov scales of turbulence. Therefore, in DNS mode the previously described equations are solved directly with no additional models being required. The transport coefficients correspond to the molecular transport coefficients, with one of these assumptions:

constant transport coefficients \(\mu\), \(\rho\alpha_C\), and \(\rho\alpha_T\), or
constant \(\nu = \mu / \rho\), \(\alpha_T\) and constant \(\alpha_C\),
each of which is then multiplied by \(\rho\) in the construction of the diffusive terms.

Note that scaling arguments from simple kinetic theory indicate that the quantities \(\mu\), \(\rho\alpha_C\), and \(\rho\alpha_T\) are all independent of density and pressure but weakly dependent on temperature (\(T^{1/2}\) scaling), justifying use of constant parameters over the range of conditions of interest.

LES¶

When running in Large Eddy Simulation (LES) mode, it is assumed that the Kolmogorov scales are not resolved. As a result, the numerical discretization acts as a filter on the governing equations, resulting in the following set of filtered equations:

\[ \begin{align}\begin{aligned}\frac{\partial \overline{\rho}}{\partial t} &= - \nabla \cdot (\overline{\rho} \mathbf{\tilde{u}}),\\\frac{\partial (\overline{\rho} \mathbf{\tilde{u}})}{\partial t} &= - \nabla \cdot (\overline{\rho} \mathbf{\tilde{u}} \mathbf{\tilde{u}}) - \nabla \overline{p} + \overline{\rho} \mathbf{g} - \nabla \cdot \overline{\tau} + \mathbf{\overline{F}} &- \nabla \cdot (\overline{\rho} \mathbf{\widetilde{u u}} - \overline{\rho}\mathbf{\tilde{u}\tilde{u}} ) ,\\\frac{\partial (\overline{\rho} \tilde{\theta})}{\partial t} &= - \nabla \cdot (\overline{\rho} \mathbf{\tilde{u}} \tilde{\theta}) + \nabla \cdot \left( \overline{\rho} \alpha_{T} \nabla \tilde{\theta} \right) &- \nabla \cdot (\overline{\rho} {\widetilde{\mathbf{u} \theta}} - \overline{\rho}\mathbf{\tilde{u}}\tilde{\theta} ) ,\\\frac{\partial (\overline{\rho} \tilde{C})}{\partial t} &= - \nabla \cdot (\overline{\rho} \mathbf{\tilde{u}} \tilde{C}) + \nabla \cdot \left( \overline{\rho} \alpha_{C} \nabla \tilde{C} \right) &- \nabla \cdot (\overline{\rho} \widetilde{\mathbf{u} C} - \overline{\rho}\mathbf{\tilde{u}}\tilde{C} ) ,\end{aligned}\end{align} \]

where overbars indicate filtering and tildes indicate density-weighted (Favre) filtering (e.g., \(\tilde{\theta} = \overline{\rho \theta} / \overline{\rho}\)). When the code is run in LES mode, all variables correspond to their appropriate filtered version.

In the above equations, the final term in each of the momentum, potential temperature, and scalar equations is unclosed due to containing a filtered nonlinear function of the state quantities. These terms represent the effect of turbulent transport at unresolved scales. LES models attempt to account for these terms by invoking a gradient transport hypothesis, which assumes that turbulent transport acts similarly to molecular transport in that quantities are transported down their resolved gradients:

\[ \begin{align}\begin{aligned}\overline{\rho} {\widetilde{\mathbf{u} \theta}} - \overline{\rho}\mathbf{\tilde{u}}\tilde{\theta} &= -\frac{\mu_t}{Pr_t} \nabla \tilde{\theta}\\\overline{\rho} \widetilde{\mathbf{u} C} - \overline{\rho}\mathbf{\tilde{u}}\tilde{C} &= -\frac{\mu_t}{Sc_t} \nabla \tilde{C}\\\overline{\rho} \mathbf{\widetilde{u u}} - \overline{\rho}\mathbf{\tilde{u}\tilde{u}} &= \tau^{sfs}\end{aligned}\end{align} \]

\[ \begin{align}\begin{aligned}\tau^{sfs}_{ij} - \frac{\delta_{ij}}{3} \tau^{sfs}_{kk} = -2 \mu_t \tilde{\sigma}_{ij}\\\tau^{sfs}_{kk} = 2 \mu_t \frac{C_I}{C_s^2} (2 \tilde{S}_{ij} \tilde{S}_{ij})^{1/2}.\end{aligned}\end{align} \]

The model coefficients \(C_s, C_I, Pr_t, Sc_t\) have nominal values of 0.16, 0.09, 0.7, and 0.7, respectively (Martin et al., Theoret. Comput. Fluid Dynamics (2000)).

Note

Presently, we assume \(C_I =0\). This term is similar to the bulk viscosity term for molecular transport and should be added if the bulk viscosity term is added. It is believed to be small for low-Mach number flows, but there is some discussion in the literature about this topic. See Moin et al., “A dynamic subgrid-scale model for 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 molecular transport terms, but with the constant molecular transport coefficients replaced by turbulent equivalents (e.g. \(\mu\) becomes the turbulent viscosity, \(\mu_{t}\)). Therefore, when the code is run in LES mode, the equation set remains the same, but all variables are interpereted as the appropriate filtered version and the turbulent transport coefficients are used. Molecular transport is omitted by default in the present LES implementation because the molecular transport coefficients are insignificant compared to turbulent transport for most LES grids. However, for fine LES grids, molecular transport and LES models may both be activated, with the effective transport coefficient being the sum of the molecular and turbulent coefficients.

It should also be noted that filtering affects the computation of pressure from density and potential temperature, but the nonlinearity in the equation of state is weak for \(\gamma = 1.4\), so the subfilter contribution is neglected:

\[\overline{p} = \overline{ \left( \frac{\rho R_d \theta}{p_0^{R_d / c_p}} \right)^\gamma} \approx \left( \frac{\overline{\rho} R_d \tilde{\theta}}{p_0^{R_d / c_p}} \right)^\gamma.\]

ERF offers two LES options: Smagorinsky and Deardorff models, which differ in how \(\mu_{t}\) is computed.

Smagorinsky Model¶

\[\mu_{t} = (C_s \Delta)^2 (\sqrt{2 \tilde{S} \tilde{S}}) \overline{\rho}\]

\(C_s\) is the Smagorinsky constant and \(\Delta\) is the cube root of cell volume, the representative mesh spacing.

\[\tau_{ij} = -2\mu_{t} \tilde{\sigma}_{ij} = -K \tilde{\sigma}_{ij}\]

where \(K = 2\mu_{t}\)

In the Smagorinsky model, modeling of \(\mu_{t}\) does not account for the turbulent kinetic energy (TKE) corresponding to unresolved scales and no extra equation for TKE is solved.

Deardorff Model¶

Unlike the Smagorinsky model, the Deardorff model accounts for the contribution of TKE in modeling \(\mu_{t}\) and a prognostic equation for TKE is solved. The turbulent viscosity is computed as:

\[\mu_t = \overline{\rho} C_k \ell (e^{sfs})^{1/2},\]

where the mixing length \(\ell = \Delta\) for unstable (or neutral) stratification, otherwise the mixing length is reduced as per Deardorff 1980:

\[\ell = \frac{0.76 (e^{sfs})^{1/2}}{\left(\frac{g}{\theta_0}\frac{\partial\theta}{\partial z}\right)^{1/2}}.\]

The potential temperature gradient in the denominator dictates the stratification, which is scaled by a reference virtual potential temperature \(\theta_0\). The mixing length is set to a maximum of \(\Delta\) to prevent unbounded behavior under weakly stable conditions.

Then the equation solved to determine \(e^{sfs}\), the subfilter contribution to TKE, is:

\[\frac{\partial \overline{\rho} e^{sfs}}{\partial t} = - \nabla \cdot (\overline{\rho} \mathbf{\tilde{u}} \tilde{e}^{sfs}) - \tau_{ij} \frac{\partial \tilde{u}_i}{\partial x_j} + \frac{g}{\theta_0} \tau_{\theta w} + \nabla \cdot \left( \frac{\mu_t}{\sigma_k} \nabla e^{sfs} \right) - \overline{\rho} C_\epsilon \frac{(e^{sfs})^{3/2}}{\overline{\Delta}}.\]

where \(\sigma_k\) is a constant model coefficient representing the ratio of turbulent viscosity to turbulent diffusivity of TKE that should be order unity (e.g., Moeng 1984 uses TKE diffusivity of \(2 \mu_t\)), we have used the downgradient diffusion assumption

\[\frac{\partial\left\langle \left( u_{n}^{'}\rho e + u_{n}^{'}p^{'} \right) \right\rangle}{\partial x_{n}} = -\nabla \cdot \left( \frac{\mu_t}{\sigma_k} \nabla e^{sfs} \right),\]

the eddy diffusivity of heat is

\[K_H = \left(1 + \frac{2\ell}{\Delta}\right) \mu_t,\]

and the SFS heat flux is

\[\tau_{\theta i} = -K_H \frac{\partial\theta}{\partial x_i}.\]

The RHS terms of the TKE transport equation correspond to advection, shear production, buoyant production, diffusion, and dissipation.

PBL Schemes¶

Planetary Boundary Layer (PBL) schemes are used to model unresolved transport in the vertical direction within the planetary boundary layer when mesh resolutions are too coarse to resolve any of the turbulent eddies responsible for this transport (~1 km grid resolution or larger). The PBL scheme is used to provide closure for vertical turbulent fluxes (i.e., \(\widetilde{w'\phi'} = \widetilde{w\phi} - \widetilde{w}\widetilde{\phi}\), for any quantity \(\phi\)). PBL schemes may be used in conjunction with an LES model that specifies horizontal turbulent transport, in which case the vertical component of the LES model is ignored.

Right now, the only PBL scheme supported in ERF is the Mellor-Yamada-Nakanishi-Niino Level 2.5 model, largely matching the original forumulation proposed by Nakanishi and Niino in a series of papers from 2001 to 2009.

MYNN Level 2.5 PBL Model¶

In this model, the vertical turbulent diffusivities are computed in a local manner based on a gradient diffusion approach with coefficients computed based on a transported turbulent kinetic energy value. The implementation and description here largely follows Nakanishi and Niino, Journal of Meteorological Society of Japan, 2009, but has also been influenced by the full series of papers that led to the development of this model and by a few documents published since then, as listed in the Useful References section below.

The prognostic equation for \(q^2 = \widetilde{u_i u_i} - \widetilde{u}_i\widetilde{u}_i\) is

\[\frac{\partial \bar{\rho} q^2}{\partial t} + \left[ \frac{\partial \bar{\rho} \widetilde{u}_i q^2}{\partial x_i} \right] = \frac{\partial}{\partial z} \left(K_{q,v} \frac{\partial q^2}{\partial z} \right) + 2\bar{\rho} \left(-\widetilde{u'w'} \frac{\partial \widetilde{u}}{\partial z} - \widetilde{v'w'}\frac{\partial \widetilde{v}}{\partial z} + \beta g \widetilde{w'\theta'} - \frac{q^3}{B_1 l} \right)\]

where \(B_1\) is a model parameter, \(\beta\) is the thermal expansion coefficient and l is a lengthscale. The vertical turbulent transport 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 \(S_m\) and \(S_\theta\) are stability parameters thaat account for bouyancy effects. These coefficients are then applied in evaluating the vertical component of turbulent fluxes in a similar manner as is described for the Smagorinsky LES model. Computation of the stability parameters and lengthscale depend on the Obukhov length and surface heat flux, which are obtained from the MOST Boundaries. Further detail on these computations can be found in the cited works. Several model coefficients are required, with default values in ERF taken from the work of Nakanishi and Niino.

Useful References¶

The following references have informed the implementation of the MYNN PBL model in ERF:

Mellor, Journal of the Atmospheric Sciences, 1973: Introduces a PBL model based on \(q^2\)
Mellor and Yamada, Journal of the Atmospheric Sciences, 1974: Introduces PBL Model Hierarchy (Levels 1-4)
Mellor and Yamada, Reviews of Geophysics and Space Physics, 1982: Introduces Level 2.5 Model
Nakanishi, Boundary-Layer Meteorology, 2001: Fits new model coefficients and proposes new diagnostic equation for the length scale
Nakanishi and Niino, Boundary-Layer Meteorology, 2004: Extends the MYNN PBL modeling framework for moist conditions
Nakanishi and Niino, Boundary-Layer Meteorology, 2006: Numerical stability improvements for the MYNN PBL modeling framework
Nakanishi and Niino, Journal of the Meteorological Society of Japan, 2009: Summary of MYNN model development, re-evaluation of coefficients, and additional demonstration cases
Skamarock et al., A Description of the Advanced Research WRF Model Version 4, 2021: Description of the models implemented in WRF
Olson et al., A Description of the MYNN-EDMF Scheme and the Coupling to Other Components in WRF–ARW, 2019: Description of more recent advancements upon the MYNN model
Juliano et al., Monthly Weather Review, 2022: Description of a 3D generalization Mellor-Yamada PBL models

Discussions with Branko Kosovic (NCAR) and Joseph B. Olson (NOAA) have also played a major role in informing the implementation of MYNN PBL models in ERF.

MYNN-EDMF Level 2.5 PBL Model¶

More recent advancements that add significant complexity to the MYNN scheme have been incorporated into WRF, as described in Olson et al. 2019. These advancements are not included in ERF, but may be in the future.

YSU PBL Model¶

Warning

Implementation is in progress, this option is not yet supported

The Yonsei University (YSU) PBL model is another commonly use scheme in WRF. It includes nonlocal mixing with contergradient diffusion within the PBL, and a local mixing treatment outside the PBL.

Turbulent diffusion for prognostic variables (\(C, u, v, \theta, q_k\)), where \(q_k\) includes all moisture variables and \(C\) any additional scalars (other terms in the equations omitted for brevity):

\[\frac{\partial C}{\partial t} = \frac{\partial}{\partial z} \left[ K_c \left( \frac{\partial C}{\partial z} - \gamma_c \right) - \overline{\left(w'c' \right)_h} \left( \frac{z}{h} \right)^3 \right]\]

Note

Not applied for vertical velocity?

Where for each variable the turbulent diffusion coefficient \(K_c\), countergradient correction \(\gamma_c\), and entrainment flux at the PBL top \(\overline{\left(w'c' \right)_h}\) must be diagnosed for each variable. The main controlling parameter is the PBL height \(h\). Notably, a nonlocal model for turbulent diffusion is used for \(z \leq h\), but a local model is used for \(z \ge h\).

The first step is to determine the PBL height \(h\). This is defined as the smallest value of \(z\) where the bulk Richardson number equals the critical value, which is taken to be 0:

\[{\rm Rib}(z) = \frac{ g \left[ \theta_m(z) - \theta_s\right] }{\theta_{ma} U(z)^2}z\]

\[{\rm Rib}(h) = {\rm Rib_{cr}} = 0\]

where

\(\theta_m\) is the moist potential temperature,
\(\theta_{ma}\) is the value at the lowest vertical cell in a column,
\(U = \sqrt{u^2 + v^2}\) is the horizontal wind speed,
\(\theta_s = \theta_{ma} + \theta_T\) is the virtual temperature near the surface,
\(\theta_T = a\frac{\overline{\left(w'\theta_m' \right)_0}}{w_{s0}}\) is the excess virtual temperature near the surface,
\(a\) is a constant taken to be 6.8 per HND06 (matching the \(b\) constant that appears elsehwere in the YSU model)
\(\overline{\left(w'\theta_m' \right)_0}\) is the surface virtual heat flux (determined from the MOST surface layer model),
\(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 \(w_{s0} = w_s(h/2)\) (note this equation matches the WRF implementation and description in H10, but differs from HND06, where \(\phi_m\) appears in place of the constant 8),
\(u_*\) is the surface frictional velocity scale determined from the MOST surface layer model,
\(k = 0.4\) is the von Karman constant
\(w_{*b} = \left[ g/\theta_{ma} \overline{\left(w'\theta_m' \right)_0} h \right]^{1/3}\) for \(\overline{\left(w'\theta_m' \right)_0} > 0\), \(w_{*b} = 0\) otherwise, is a convective velcoity scale for moist air

In practice, an approximate value of \(h\) is determined through a two-step process. First, \(\theta_T\) is set to be zero and a provisional value of \(h\) is estimated. Then this provisional value of \(h\) is used to compute \(\theta_T\), which is in turn used to provide an improved estimate of \(h\), which is the value used in subsequent calculations.

Note

This two step-process matches the WRF implementation, but this could be extended iteratively to reach convergence.

Countergradient corrections are computed as follows:

\[\gamma_\theta =\]

\[\gamma_u =\]

\[\gamma_v =\]

\[\gamma_{q_k} = \gamma_C = 0\]

Entrainment fluxes are computed:

\[\overline{\left(w'c' \right)_h} =\]

Within the PBL (\(z \leq h\)),

Useful References¶

The following references have informed the implementation of the YSU model in ERF:

[H10] Hong, Quarterly Journal of the Royal Meteorological Society, 2010: Most up-to-date YSU model formulation as implemented in WRF, with revisions for stable boundary layers
[HND06] Hong, Noh, and Dudhia, Monthly Weather Review, 2006: Initial formulation referred to as the YSU model, adds improved entrainment formulation (relative to NCHR03) to work of TM86 and a few other modifications
[NCHR03] Noh, Cheon, Hong, and Raasch, Boundary-Layer Meteorology, 2003: Entrainment effects added to TM86
[HP96] Hong and Pan, Monthly Weather Review, 1996: Largely an implementation and evluation of TM86
[TM86] Troen and Mahrt, Boundary-Layer Meteorology, 1986: Initial incorporation of nonlocal counter-graident term in vertical diffusion model
[WF18] Wilson and Fovell, Weather and Forecasting, 2018: Extension of YSU to handle interplay between radiation and fog, active in WRF with the ysu_topdown_pblmix = 1 option
The WRF Fortran source code for this module as of Dec. 2023. The ERF implementation supports the same physical models as this WRF implementation, with the exception of the ysu_topdown_pblmix = 1 option from WF18, i.e. the implementation in ERF largely matches the PBL scheme described in H10.

Physical Forcings¶

ERF includes the following forcing terms as options:

Buoyancy¶

If

use_gravity == true

then buoyancy is included in the momentum equations in the form

\[(0, 0, -\rho^\prime g)\]

Coriolis Forcing¶

If

use_coriolis == true

then Coriolis forcing is included in the momentum equations, i.e. :

\[\mathbf{F} = (C_f \; (\rho v \sin{\phi} - \rho w \cos{\phi}), -C_f \; \rho u \sin{\phi}, C_f \; \rho u \cos{\phi})\]

where \(C_f = 4 \pi / P_{rot}\) is the Coriolis factor with \(P_{rot}\) the rotational period (measured in seconds), and \(\phi\) the latitude.

There is no dependence on the radial distance from the center of the earth, thus the curvature of the earth is neglected.

Rayleigh Damping¶

Rayleigh damping can be imposed on any or all of \(u, v, w, T\) and is controlled by setting

rayleigh_damp_U = true
rayleigh_damp_V = true
rayleigh_damp_W = true
rayleigh_damp_T = true

in the inputs file. When one or more of those is true, explicit Rayleigh damping is included in the energy and/or momentum equations as described in Section 4.4.3 of the WRF Model Version 4 documentation (p40), i.e. :

\[\mathbf{F} = - \tau(z) \rho \; (u - \overline{u}, v - \overline{v}, 0)\]

and

\[F_{\rho \theta} = - \tau(z) \rho (\theta - \overline{\theta})\]

where \((\overline{u}, \overline{v}, 0)\) is the reference state velocity, typically defined as the initial horizontally homogeneous fields in idealized simulations, and \(\overline{\theta}\) is the reference state potential temperature. As in the WRF model, the reference state vertical velocity is assumed to be zero.

Problem-Specific Forcing¶

There are two ways to specify background conditions to drive the simulation:

Pressure Gradient¶

If

abl_driver_type == "PressureGradient"

then

\[\mathbf{F} = (\nabla p_{x,ext}, \nabla p_{y,ext}, \nabla p_{z,ext})\]

where \((\nabla p_{x,ext}, \nabla p_{y,ext}, \nabla p_{z,ext})\) are user-specified through erf.abl_pressure_grad.

Geostrophic Forcing¶

If

abl_driver_type == "GeostrophicWind"

then geostrophic forcing is included in the forcing terms, i.e.

\[\mathbf{F} = (-C_f \; v_{geo}, C_f \; u_{geo}, 0)\]

where \(C_f = 4 \pi / P_{rot}\) is the Coriolis factor with \(P_{rot}\) the rotational period (measured in seconds), and the geostrophic wind \((u_{geo}, v_{geo}, 0)\) is user-specified through erf.abl_geo_wind. Note that if geostrophic forcing is enabled, Coriolis forcing must also be included.

Particles¶

ERF has the option to include Lagrangian particles in addition to the mesh-based solution. Currently there are two particle types available in ERF: tracer_particles and hydro_particles. The particle functionality is very simple and meant for demonstration. The particles are initialized one per mesh cell in a vertical plane at \(i = 3\) for tracer particles and a horizontal plane at \(k = 23\) for hydro particles. The tracer particles are advected by the velocity field; the hydro particles fall with a velocity determined by gravity minus drag.

However, the AMReX particle data structure is very general and particles may take on a number of different roles in future.

To enable the use of particles, one must set

USE_PARTICLES = TRUE

in the GNUmakefile if using gmake, or add

-DERF_ENABLE_PARTICLES:BOOL=ON \

to the cmake command if using cmake. (See, e.g., Build/cmake_with_particles.sh)

One must also set

erf.use_tracer_particles = 1

and / or

erf.use_hydro_particles = 1

in the inputs file or on the command line at runtime.

Caveat: the particle information is currently output when using the AMReX-native plotfile format, but not when using netcdf. Writing particles into the netcdf files is a WIP.

To see an example of using the particle functionality, build the executable using gmake in Exec/DevTests/ParticlesOverWoA.

To visualize the number of particles per cell as a mesh-based variable, add tracer_particle_count (if you have set erf.use_tracer_particles) and hydro_particle_count (if you have set erf.use_tracer_particles) to the line in the inputs file that begins

erf.plot_vars_1 =

Wind farm models¶

Introduction¶

ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently only the Fitch model (Fitch et al. 2012) is supported.

Fitch model¶

The Fitch model for wind farms introduced in Fitch et al. 2012 models the effect of wind farms (See Fig. 1) as source terms in the governing equations for the horizontal components of momentum (i.e., \(x\) and \(y\) momentum) and the turbulent kinetic energy (TKE). The wind turbine is discretized only in the vertical (ie. z direction). At a given cell \((i,j,k)\), the source terms in the governing equations are

\[\begin{split}\frac{\partial u_{ijk}}{\partial t} &= \frac{u_{ijk}}{|V|_{ijk}}\frac{\partial |V|_{ijk}}{\partial t} \\ \frac{\partial v_{ijk}}{\partial t} &= \frac{v_{ijk}}{|V|_{ijk}}\frac{\partial |V|_{ijk}}{\partial t} \\ \frac{\partial \text{TKE}_{ijk}}{\partial t} &= \frac{0.5N^{ij}C_{TKE}(|V|_{ijk})|V|_{ijk}^3A_{ijk}}{z_{k+1}-z_k}\end{split}\]

where

\[\frac{\partial |V|_{ijk}}{\partial t} = \frac{0.5N^{ij}C_T(|V|_{ijk})|V|_{ijk}^2A_{ijk}}{z_{k+1}-z_k}\]

where u and v are horizontal components of velocity, |V| is the velocity magnitude, \(N^{ij}\) is the number of turbines in cell \((i,j)\), \(C_T\) is the coefficient of thrust of the turbines, \(C_{TKE}\) is the fraction of energy converted to turbulent kinetic energy – both of which are functions of the velocity magnitude, and \(A_{ijk}\) is the area intersected by the swept area of the turbine between levels \(z=z_k\) and \(z= z_{k+1}\), and is explained in the next section.

Intersected area \(A_{ijk}\)¶

Consider \(A_k^{k+1}\) – the area intersected by the swept area of the wind turbine between \(z=z_k\) and \(z = z_{k+1}\). We have (see Figs. 2 and 3 below)

\[A_k = \frac{\pi R^2}{2} - A_{ks}\]

where \(A_{ks}\) is the area of the segment of the circle as shown in Fig. 3. We have from geometry, \(d_k = \min(|z_k - z_c|,R)\) is the perpendicular distance of the center of the turbine to \(z = z_k\), where \(z_c\) is the height of the center of the turbine from the ground. The area of the segment is

\[A_{ks} = R^2\theta_k - d_k\sqrt{R^2 - d_k^2}\]

where \(\theta_k = \cos^{-1}\left(\frac{d_k}{R}\right)\).

Hence, we have the intersected area \(A_{ijk}\equiv A_k^{k+1}\) as

\[\begin{split}A_k^{k+1} = \begin{cases} |A_k - A_{k+1}| & \text{if } (z_k - z_c)(z_{k+1}-z_c) > 0 \\ |A_k + A_{k+1}| & \text{if } (z_k - z_c)(z_{k+1}-z_c) \le 0 \\ \end{cases}\end{split}\]

An example of the Fitch model is in Exec/Fitch

_images/WindFarm_Fitch.png — Horizontal view of the wind farm in the Fitch model – shows a wind farm in cell (i,j) with 5 wind turbines. The turbines are discretized only in the vertical direction.¶

_images/WindTurbine_Fitch.png — The vertical discretization of the wind turbine in the Fitch model.¶

_images/FitchModel_A_ijk.png — The area terminology in the Fitch model. The circle represents the area swept by the turbine blades.¶

Constants and Units¶

Units¶

The following units are used in ERF:

Name

Units

Description

\(t\)

\(s\)

time

\(\rho\)

\(kg/m^3\)

density

\(\mathbf{u}\)

\(m/s\)

velocity

\(p\)

\(Pa\)

pressure

\(T\)

\(K\)

temperature

\(\theta\)

\(K\)

potential temperature

Constants¶

Name

Value

Description

\(R_d\)

287.0 J/kg/K

gas constant for dry air

\(c_p\)

1004.5 J/kg/K

specific heat capacity

\(p_0\)

1.0e5 Pa

reference pressure

\(g\)

9.81 m/s/s

gravity

\(\gamma\)

1.4

\(c_p /(c_p-R_d)\)

These constants are defined in the file Source/ERF_Constants.H

Arakawa C-Grid¶

Variables are located on Arakawa C-grid as pictured in the images below.

Map Factors (Dry)¶

ERF supports the use of isotropic projections only, e.g. Lambert conformal, polar stereographic and Mercator. Like WRF, ERF implements the projections using map factors, which are defined as the ratio of the distance in computational space to the corresponding distance on the earth’s surface.

With map factors, the dry governing equations have the following form

\[ \begin{align}\begin{aligned}\frac{\partial \rho}{\partial t} &= - m^2 \nabla_H \cdot (\rho \mathbf{u_H} m^{-1}) - m \frac{\partial (\rho w m^{-1})}{\partial z}\\\frac{\partial (\rho u)}{\partial t} &= - m^2 \nabla_H \cdot (\rho \mathbf{u_H} u m^{-1}) - m \frac{\partial (\rho w u m^{-1})}{\partial z} - m \frac{\partial p^\prime}{\partial x} + \nabla \cdot \tau + \mathbf{F},\\\frac{\partial (\rho v)}{\partial t} &= - m^2 \nabla_H \cdot (\rho \mathbf{u_H} v m^{-1}) - m \frac{\partial (\rho w v m^{-1})}{\partial z} - m \frac{\partial p^\prime}{\partial y} + \nabla \cdot \tau + \mathbf{F},\\\frac{\partial (\rho w) }{\partial t} &= - m^2 \nabla_H \cdot (\rho \mathbf{u_H} w m^{-1}) - m \frac{\partial (\rho w w m^{-1})}{\partial z} - \frac{\partial p^\prime}{\partial z} + \rho^\prime \mathbf{g} + \nabla \cdot \tau + F_z,\\\frac{\partial (\rho \theta)}{\partial t} &= - m^2 \nabla_H \cdot (\rho \mathbf{u_H} \theta m^{-1}) - m \frac{\partial (\rho w \theta m^{-1})}{\partial z} + \nabla \cdot ( \rho \alpha_{T}\ m \nabla \theta) + F_{\rho \theta},\\\frac{\partial (\rho C)}{\partial t} &= - m^2 \nabla_H \cdot (\rho \mathbf{u_H} C m^{-1}) - m \frac{\partial (\rho w C m^{-1})}{\partial z} + \nabla \cdot (\rho \alpha_{C}\ m \nabla C)\end{aligned}\end{align} \]

where \(u_H\) is the velocity in the horizontal (lateral) only, \(\nabla_H\) is the horizontal (lateral) gradient, \(m\) is the map factor at the appropriate spatial location and \(m^{-1} = 1 / m\).

The viscous stress tensor \(\tau\) is modified via the strain rates \(S_{ij}\):

\[ \begin{align}\begin{aligned}S_{11} &= m^2*\left[ \partial_x (um^{-1}) - \frac{h_\xi}{h_\zeta}\partial_z (um^{-1}) \right]\\S_{22} &= m^2*\left[ \partial_y (vm^{-1}) - \frac{h_\eta}{h_\zeta}\partial_z (vm^{-1})) \right]\\S_{33} &= \frac{1}{h_\zeta}\partial_z w\\S_{12} &= \frac{m^2}{2} * \left[ \partial_y (m^{-1}u) + \partial_x (vm^{-1}) - \frac{h_\eta}{h_\zeta} \partial_z (um^{-1}) - \frac{h_\xi}{h_\zeta}\partial_z (vm^{-1}) \right]\\S_{13} &= \frac{1}{2} * \left[ \frac{1}{h_\zeta}\partial_z u + m * \left( \partial_x w - \frac{h_\xi}{h_\zeta} \partial_z w \right) \right]\\S_{23} &= \frac{1}{2} * \left[ \frac{1}{h_\zeta}\partial_z v + m * \left( \partial_y w - \frac{h_\eta}{h_\zeta} \partial_z w \right) \right]\end{aligned}\end{align} \]

When LES models are used, the cell volume is used to compute the eddy viscosities. The cell volume must be adjusted using the map factors:

\[cellVol = \frac{1}{\Delta x * m_x * \Delta y * m_y * \Delta z}\]

Time Advance¶

To advance the solution in time, ERF uses a 3rd order Runge-Kutta method with acoustic sub-stepping in each Runge-Kutta stage, following the approach of Klemp, Skamarock and Dudhia (2007)

Specifically, for

\[\frac{d \mathbf{S}}{dt} = f(\mathbf{S})\]

where \(\mathbf{S}\) is the solution vector, we solve

\[ \begin{align}\begin{aligned}\mathbf{S}^{*} &=& \mathbf{S}^n + \frac{1}{3} \Delta t f(\mathbf{S}^n)\\\mathbf{S}^{**} &=& \mathbf{S}^n + \frac{1}{2} \Delta t f(\mathbf{S}^{*})\\\mathbf{S}^{n+1} &=& \mathbf{S}^n + \Delta t f(\mathbf{S}^{**})\end{aligned}\end{align} \]

Acoustic Sub-stepping¶

We sub-step the acoustic modes within each Runge-Kutta stage.

We first recall the equations in the following form, here defining \(\mathbf{R}\) for each equation to include all additional terms that contribute to the time evolution.

\[ \begin{align}\begin{aligned}\frac{\partial \rho}{\partial t} &=& - \nabla \cdot (\rho \mathbf{u}) = R_\rho,\\\frac{\partial (\rho \mathbf{u})}{\partial t} &=& - \nabla \cdot (\rho \mathbf{u} \mathbf{u} + p^\prime I) + {\mathbf F}_\mathbf{u} = \mathbf{R}_\mathbf{u}\\\frac{\partial (\Theta)}{\partial t} &=& - \nabla \cdot (\mathbf{u} \Theta) + \rho \alpha_{T}\ \nabla^2 \theta = R_{\Theta},\\\frac{\partial (\rho C)}{\partial t} &=& - \nabla \cdot (\rho \mathbf{u} C) + \rho \alpha_{C}\ \nabla^2 C = R_{\rho C},\end{aligned}\end{align} \]

where we have defined \(\mathbf{U} = (U,V,W) = \rho \mathbf{u} = (\rho u, \rho v, \rho w)\) , \(\Theta = \rho \theta\) and \(\mathbf{F}_\mathbf{U} = (F_U, F_V, F_W) = \rho^\prime \mathbf{g} + \nabla \cdot \tau + \mathbf{F}\)

Using the relation \(\nabla p = \gamma R_d \pi \nabla \Theta,\) where the Exner function \(\pi = (p/p_0)^\kappa\) with \(\kappa = R_d / c_p,\) we can re-write the unapproximated momentum equations

\[ \begin{align}\begin{aligned}\frac{\partial U}{\partial t} &=& - \nabla \cdot (\mathbf{u} U) - \gamma R_d \pi \frac{\partial \Theta^\prime}{\partial x} + F_u\\\frac{\partial V}{\partial t} &=& - \nabla \cdot (\mathbf{u} V) - \gamma R_d \pi \frac{\partial \Theta^\prime}{\partial y} + F_v\\\frac{\partial W}{\partial t} &=& - \nabla \cdot (\mathbf{u} W) - \gamma R_d \pi \frac{\partial \Theta^\prime}{\partial z} - g (\overline{\rho} \frac{\pi^\prime}{\overline{\pi}} - \rho^\prime) + F_w\end{aligned}\end{align} \]

We then define new perturbational quantities, e.g., \(\mathbf{U}^{\prime \prime} = \mathbf{U} - \mathbf{U}^t\) where \(\mathbf{U}^t\) is the momentum at the most recent time reached by the “large” time step, which could be \(t^{n}\) or one of the intermediate Runge-Kutta steps. Then the acoustic substepping evolves the equations in the form

\[ \begin{align}\begin{aligned}U^{\prime \prime, \tau + \delta \tau} - U^{\prime \prime, \tau} &= \delta \tau \left( -\gamma R_d \pi^t \frac{\partial \Theta^{\prime \prime, \tau}}{\partial x} + R^t_U \right)\\V^{\prime \prime, \tau + \delta \tau} - V^{\prime \prime, \tau} &= \delta \tau \left( -\gamma R_d \pi^t \frac{\partial \Theta^{\prime \prime, \tau}}{\partial y} + R^t_V \right)\end{aligned}\end{align} \]

\[\begin{split}W^{\prime \prime, \tau + \delta \tau} - W^{\prime \prime, \tau} = \delta \tau \biggl( &-\gamma R_d \pi^t \frac{\partial (\beta_1 \Theta^{\prime \prime, \tau} + \beta_2 \Theta^{\prime \prime, \tau + \delta \tau} ) }{\partial z} \\ & - g \overline{\rho} \frac{R_d}{c_v} \frac{\pi^t}{\overline{\pi}} \frac{ (\beta_1 \Theta^{\prime \prime, \tau} + \beta_2 \Theta^{\prime \prime, \tau + \delta \tau} )}{\Theta^t} \\ & + g (\beta_1 \rho^{\prime \prime, \tau} + \beta_2 \rho^{\prime \prime, \tau + \delta \tau } ) \\ & + R^t_W \biggr)\end{split}\]

\[\Theta^{\prime \prime, \tau + \delta \tau} - \Theta^{\prime \prime, \tau} = \delta \tau \left( -\frac{\partial (U^{\prime \prime, \tau + \delta \tau} \theta^t)}{\partial x} -\frac{\partial (V^{\prime \prime, \tau + \delta \tau} \theta^t)}{\partial y} -\frac{\partial \left(( \beta_1 W^{\prime \prime, \tau} + \beta_2 W^{\prime \prime, \tau + \delta \tau} ) \theta^t\right)}{\partial z} + R^t_{\Theta} \right)\]

\[\rho^{\prime \prime, \tau + \delta \tau} - \rho^{\prime \prime, \tau} = \delta \tau \left( - \frac{\partial U^{\prime \prime, \tau + \delta \tau }}{\partial x} - \frac{\partial V^{\prime \prime, \tau + \delta \tau }}{\partial y} - \frac{\partial (\beta_1 W^{\prime \prime, \tau} + \beta_2 W^{\prime \prime, \tau + \delta \tau})}{\partial z} + R^t_{\rho} \right)\]

where \(\beta_1 = 0.5 (1 - \beta_s)\) and \(\beta_2 = 0.5 (1 + \beta_s)\) with \(\beta_s = 0.1\). \(\beta_s\) is the acoustic step off-centering coefficient and 0.1 is the typical WRF value. This off-centering is intended to provide damping of both horizontally and vertically propagating sound waves by biasing the time average toward the future time step.

To solve the coupled system, we first evolve the equations for \(U^{\prime \prime, \tau + \delta \tau}\) and \(V^{\prime \prime, \tau + \delta \tau}\) explicitly using \(\Theta^{\prime \prime, \tau}\) which is already known. We then solve a tridiagonal system for \(W^{\prime \prime, \tau + \delta \tau}\), and once \(W^{\prime \prime, \tau + \delta \tau}\) is known, we update \(\rho^{\prime \prime, \tau + \delta \tau}\) and \(\Theta^{\prime \prime, \tau + \delta \tau}.\)

In addition to the acoustic off-centering, divergence damping is also applied to control horizontally propagating sound waves.

\[p^{\prime\prime,\tau*} = p^{\prime\prime,\tau} + \beta_d \left( p^{\prime\prime,\tau} + p^{\prime\prime,\tau-\delta\tau} \right)\]

where \(\tau*\) is the forward projected value used in RHS of the acoustic substepping equations for horizontal momentum. According to Skamarock et al, This is equivalent to including a horizontal diffusion term in the continuity equation. A typical damping coefficient of \(\beta_d = 0.1\) is used, as in WRF.

Discretizations¶

Last update: 2021-11-24

NOTE: For the sake of simplicity, the discretized equations mention time levels \(n\) and \(n+1\). They should be treated as initial and final states of each RK3 stage.

Staggered Grids¶

The staggered grids indicating where different variables are located.

XY Plane¶

YZ Plane¶

Mass Conservation¶

\[\begin{split}\begin{align} \rho_{i,j,k}^{n + 1} = \rho_{i,j,k}^{n} - \Delta t & \left\{ \frac{1}{\Delta x} \left\lbrack \left( \rho u \right)_{i + 1,j,k}^{n} - \left( \rho u \right)_{i.j,k}^{n} \right\rbrack \right. \\ & \hspace{-5pt} + \frac{1}{\Delta y} \left\lbrack \left( \rho v \right)_{i,j + 1,k}^{n} - \left( \rho v \right)_{i.j,k}^{n} \right\rbrack \\ & \left. + \frac{1}{\Delta z} \left\lbrack \left( \rho w \right)_{i,j,k + 1}^{n} - \left( \rho w \right)_{i,j,k}^{n} \right\rbrack \right\} \end{align}\end{split}\]

Contributions from different directions¶

Advection Contribution to DNS/LES¶

U Momentum¶

\[\begin{split}\begin{align} \left( \rho u \right)_{i,j,k}^{n + 1} & = \left( \rho u \right)_{i,j,k}^{n} - \\ \Delta t & \left\{ \frac{1}{2\Delta x}\ \left\lbrack \left( \left( \rho u \right)_{i + 1,j,k}^{n} + \left( \rho u \right)_{i,j,k}^{n} \right)u_{i + \frac{1}{2},j,k}^{n} - \left( \left( \rho u \right)_{i,j,k}^{n} + \left( \rho u \right)_{i - 1,j,k}^{n} \right)u_{i - \frac{1}{2},j,k}^{n} \right\rbrack \right. \\ & + \frac{1}{2\Delta y} \left\lbrack \left( \left( \rho v \right)_{i,j + 1,k}^{n} + \left( \rho v \right)_{i - 1,j + 1,k}^{n} \right)u_{i,j + \frac{1}{2},k}^{n} - \left( \left( \rho v \right)_{i,j,k}^{n} + \left( \rho v \right)_{i - 1,j,k}^{n} \right)u_{i,j - \frac{1}{2},k}^{n} \right\rbrack \\ & + \left. \frac{1}{2\Delta z} \left\lbrack \left( \left( \rho w \right)_{i,j,k + 1}^{n} + \left( \rho w \right)_{i - 1,j,k + 1}^{n} \right)u_{i,j,k + \frac{1}{2}}^{n} - \left( \left( \rho w \right)_{i,j,k}^{n} + \left( \rho w \right)_{i - 1,j,k}^{n} \right)u_{i,j,k - \frac{1}{2}}^{n} \right\rbrack \right\} \\ & - \frac{\Delta t}{\Delta x}\left\lbrack p_{i,\ j,\ k}^{n} - p_{i - 1,\ j,\ k}^{n} \right\rbrack \\ \end{align}\end{split}\]

Contributions from different directions¶

V Momentum¶

\[\begin{split}\begin{align} \left( \rho v \right)_{i,j,k}^{n + 1} & = \left( \rho v \right)_{i,j,k}^{n} - \\ \Delta t & \left\{ \frac{1}{2\Delta x}\ \left\lbrack \left( \left( \rho u \right)_{i + 1,j,k}^{n} + \left( \rho u \right)_{i + 1,j - 1,k}^{n} \right)v_{i + \frac{1}{2},j,k}^{n} - \left( \left( \rho u \right)_{i,j,k}^{n} + \left( \rho u \right)_{i,j - 1,k}^{n} \right)v_{i - \frac{1}{2},j,k}^{n} \right\rbrack \right. \\ & + \frac{1}{2\Delta y} \left\lbrack \left( \left( \rho v \right)_{i,j + 1,k}^{n} + \left( \rho v \right)_{i,j,k}^{n} \right)v_{i,j + \frac{1}{2},k}^{n} - \left( \left( \rho v \right)_{i,j,k}^{n} + \left( \rho v \right)_{i,j - 1,k}^{n} \right)v_{i,j - \frac{1}{2},k}^{n} \right\rbrack \\ & + \left. \ \frac{1}{2\Delta z} \left\lbrack \left( \left( \rho w \right)_{i,j,k + 1}^{n} + \left( \rho w \right)_{i,j - 1,k + 1}^{n} \right)v_{i,j,k + \frac{1}{2}}^{n} - \left( \left( \rho w \right)_{i,j,k}^{n} + \left( \rho w \right)_{i,j - 1,k}^{n} \right)v_{i,j,k - \frac{1}{2}}^{n} \right\rbrack \right\} \\ & - \frac{\Delta t}{\Delta y}\left\lbrack p_{i,j,\ k}^{n} - p_{i,\ j - 1,\ k}^{n} \right\rbrack \\ \end{align}\end{split}\]

Contributions from different directions¶

W Momentum¶

\[\begin{split}\begin{align} \left( \rho w \right)_{i,j,k}^{n + 1} & = \left( \rho w \right)_{i,j,k}^{n} - \\ \Delta t & \left\{ \frac{1}{2\Delta x}\ \left\lbrack \left( \left( \rho u \right)_{i + 1,j,k}^{n} + \left( \rho u \right)_{i + 1,j,k - 1}^{n} \right)w_{i + \frac{1}{2},j,k}^{n} - \left( \left( \rho u \right)_{i,j,k}^{n} + \left( \rho u \right)_{i,j,k - 1}^{n} \right)w_{i - \frac{1}{2},j,k}^{n} \right\rbrack \right. \\ & + \frac{1}{2\Delta y} \left\lbrack \left( \left( \rho v \right)_{i,j + 1,k}^{n} + \left( \rho v \right)_{i,j + 1,k - 1}^{n} \right)w_{i,j + \frac{1}{2},k}^{n} - \left( \left( \rho v \right)_{i,j,k}^{n} + \left( \rho v \right)_{i,j,k - 1}^{n} \right)w_{i,j - \frac{1}{2},k}^{n} \right\rbrack \\ & + \left. \ \frac{1}{2\Delta z} \left\lbrack \left( \left( \rho w \right)_{i,j,k + 1}^{n} + \left( \rho w \right)_{i,j,k}^{n} \right)w_{i,j,k + \frac{1}{2}}^{n} - \left( \left( \rho w \right)_{i,j,k}^{n} + \left( \rho w \right)_{i,j,k - 1}^{n} \right)w_{i,j,k - \frac{1}{2}}^{n} \right\rbrack \right\} \\ & - \frac{\Delta t}{\Delta z}\left\lbrack p_{i,\ j,\ k}^{n} - p_{i,\ j,\ \ k - 1}^{n} \right\rbrack\ + \ \Delta t \left\lbrack \rho_{i,j,k - \ \frac{1}{2}}^{n} \right\rbrack g \\ \end{align}\end{split}\]

Contributions from different directions¶

Potential Temperature Advection¶

\[\begin{split}\begin{align} \left( \rho \theta \right)_{i,j,k}^{n + 1} = \left( \rho \theta \right)_{i,j,k}^{n} - \Delta t & \left\{ \frac{1}{\Delta x}\ \left\lbrack \left( \rho u \right)_{i + 1,j,k}^{n} \theta_{i + \frac{1}{2},j,k}^{n} - \left( \rho u \right)_{i,j,k}^{n} \theta_{i - \frac{1}{2},j,k}^{n} \right\rbrack \right. \\ & + \frac{1}{\Delta y} \left\lbrack \left( \rho v \right)_{i,j + 1,k}^{n} \theta_{i,j + \frac{1}{2},k}^{n} - \left( \rho v \right)_{i,j,k}^{n} \theta_{i,j - \frac{1}{2},k}^{n} \right\rbrack \\ & + \left. \frac{1}{\Delta z} \left\lbrack \left( \rho w \right)_{i,j,k + 1}^{n} \theta_{i,j,k + \frac{1}{2}}^{n} - \left( \rho w \right)_{i,j,k}^{n} \theta_{i,j,k - \frac{1}{2}}^{n} \right\rbrack \right\} \\ \end{align}\end{split}\]

Contributions from different directions¶

Scalar Advection¶

\[\begin{split}\begin{align} \left( \rho C \right)_{i,j,k}^{n + 1} = \left( \rho C \right)_{i,j,k}^{n} - \Delta t & \left\{ \frac{1}{\Delta x}\ \left\lbrack \left( \rho u \right)_{i + 1,j,k}^{n} C_{i + \frac{1}{2},j,k}^{n} - \left( \rho u \right)_{i,j,k}^{n} C_{i - \frac{1}{2},j,k}^{n} \right\rbrack \right. \\ & + \frac{1}{\Delta y} \left\lbrack \left( \rho v \right)_{i,j + 1,k}^{n} C_{i,j + \frac{1}{2},k}^{n} - \left( \rho v \right)_{i,j,k}^{n} C_{i,j - \frac{1}{2},k}^{n} \right\rbrack \\ & + \left. \frac{1}{\Delta z} \left\lbrack \left( \rho w \right)_{i,j,k + 1}^{n} C_{i,j,k + \frac{1}{2}}^{n} - \left( \rho w \right)_{i,j,k}^{n} C_{i,j,k - \frac{1}{2}}^{n} \right\rbrack \right\} \\ \end{align}\end{split}\]

Contributions from different directions¶

Diagnostic Variables¶

\[p_{i, j, k}^n = (\rho_{i, j, k}^n R_d \theta_{i, j, k}^n / p_0^{R_d / c_p} )^\gamma\]

\[T_{i, j, k}^n = \frac{p_{i, j, k}^n}{ \rho_{i, j, k}^n R_d}\]

Here \(\rho_{i, j, k}^n, T_{i, j, k}^n, \theta_{i, j, k}^n\), and \(p_{i, j, k}^n\) are the density, temperature, potential temperature and pressure, respectively; these variables are all defined at cell centers of cell indexed by \((i, j, k)\) and at time level \(n\).

\(R_d\) and \(c_p\) are the gas constant and specific heat capacity for dry air respectively, and \(\gamma = c_p / (c_p - R_d)\) . \(p_0\) is a reference value for pressure.

Differencing of Different Orders¶

Interpolation Methods¶

The midpoint values given above \(q_{m - \frac{1}{2}}\), where \(q\) may be \([u, v, w, \theta, C]\), \(m = i, j, k\), \(U_d = [u, v, w]\) for \([x, y, z]\) directions respectively, the following interpolation schemes are available:

\[\begin{split}\begin{array}{lll} \left. q_{m - \frac{1}{2}} \right|^{2nd} = \frac{1}{2}\left( q_{m} + q_{m - 1} \right) & & \\ \left. q_{m - \frac{1}{2}} \right|^{4th} = \frac{7}{12}\left( q_{m} + q_{m - 1} \right) & \hspace{-5pt} - \frac{1}{12}\left( q_{m + 1} + q_{m - 2} \right) & \\ \left. q_{m - \frac{1}{2}} \right|^{6th} = \frac{37}{60}\left( q_{m} + q_{m - 1} \right) & \hspace{-5pt} - \frac{2}{15}\left( q_{m + 1} + q_{m - 2} \right) & \hspace{-5pt} + \frac{1}{60}\left( q_{m + 2} + q_{m - 3} \right) \\ & & \\ \left. q_{m - \frac{1}{2}} \right|^{3rd} = \left. q_{m - \frac{1}{2}} \right|^{4th} & \hspace{-5pt} + \beta_{\text{up}}\frac{U_{d}}{\left| U_{d} \right|}\frac{1}{12}\left\lbrack \left( q_{m + 1} - q_{m - 2} \right) \right.\ & \hspace{-5pt} - 3\left. \ \left( q_{m} - q_{m - 1} \right) \right\rbrack \\ & & \\ \left. q_{m - \frac{1}{2}} \right|^{5th} = \left. q_{m - \frac{1}{2}} \right|^{6th} & \hspace{-5pt} - \beta_{\text{up}}\frac{U_{d}}{\left| U_{d} \right|}\frac{1}{60}\left\lbrack \left( q_{m + 2} - q_{m - 1} \right) \right.\ & \hspace{-5pt} - 5\left( q_{m + 1} - q_{m - 2} \right) + 10\left. \left( q_{m} - q_{m - 1} \right) \right\rbrack \end{array}\end{split}\]

An extra blending factor (\(\beta_{\text{up}}\)) has been added to allow control over the amount of upwinding added to the scheme. This hybrid scheme has been demonstrated by Sauer and Muñoz-Esparza (2020, JAMES).

Refs:

Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Liu, Z., Berner, J., … Huang, X. -yu. (2019). A Description of the Advanced Research WRF Model Version 4 (No. NCAR/TN-556+STR). doi:10.5065/1dfh-6p97

Sauer, J. A., & Muñoz-Esparza, D. (2020). The FastEddy® resident-GPU accelerated large-eddy simulation framework: Model formulation, dynamical-core validation and performance benchmarks. Journal of Advances in Modeling Earth Systems, 12, e2020MS002100. doi:10.1029/2020MS002100

WENO Methods¶

Additionally, weighted essentially non-oscillatory (WENO) schemes are available for \(3rd\) and \(5th\) order interpolation. The general formulation is as follows:

\[\begin{split}\begin{array}{ll} q_{m + \frac{1}{2}} = \sum_{n=1}^{N} w_n q_{m + \frac{1}{2}}^{(n)} & \\ w_{n} = \frac{\hat{w}_{n}}{\sum_{l=1}^{N} \hat{w}_{l}} & \hat{w}_{l} = \frac{\omega_{l}}{\left(\epsilon + \beta_{l} \right)^2} \\ \end{array}\end{split}\]

With the WENO3 scheme, one has \(N=2, \; \omega_{1} = 1/3, \; \omega_{2} = 2/3\) and the following closures

\[\begin{split}\begin{array}{l} \beta_{1} = \left(q_{m} - q_{m-1} \right)^2 \\ \beta_{2} = \left(q_{m+1} - q_{m} \right)^2 \\ q_{m + \frac{1}{2}}^{(1)} = -\frac{1}{2} q_{m-1} + \frac{3}{2} q_{m} \\ q_{m + \frac{1}{2}}^{(2)} = \frac{1}{2} q_{m} + \frac{1}{2} q_{m+1} \end{array}\end{split}\]

With the WENO5 scheme, one has \(N=3, \; \omega_{1} = 1/10, \; \omega_{2} = 3/5, \; \omega_{3} = 3/10\) and the following closures

\[\begin{split}\begin{array}{l} \beta_{1} = \frac{13}{12} \left(q_{m-2} - 2 q_{m-1} + q_{m} \right)^2 + \frac{1}{4} \left(q_{m-2} - 4 q_{m-1} + 3 q_{m} \right)^2 \\ \beta_{2} = \frac{13}{12} \left(q_{m-1} - 2 q_{m} + q_{m+1} \right)^2 + \frac{1}{4} \left(q_{m-1} - q_{m+1} \right)^2 \\ \beta_{3} = \frac{13}{12} \left(q_{m} - 2 q_{m+1} + q_{m+2} \right)^2 + \frac{1}{4} \left( 3 q_{m} - 4 q_{m+1} + q_{m+2} \right)^2 \\ q_{m + \frac{1}{2}}^{(1)} = \frac{1}{3} q_{m-2} - \frac{7}{6} q_{m-1} + \frac{11}{6} q_{m} \\ q_{m + \frac{1}{2}}^{(2)} = -\frac{1}{6} q_{m-1} + \frac{5}{6} q_{m} + \frac{1}{3} q_{m+1} \\ q_{m + \frac{1}{2}}^{(3)} = \frac{1}{3} q_{m} + \frac{5}{6} q_{m+1} - \frac{1}{6} q_{m+2} \end{array}\end{split}\]

By default, the WENO3 scheme will be employed for moisture variables if the code is compiled with moisture support. However, users may utilize the WENO scheme for dry scalar variables as well. The scheme for each type is specified by

erf.dryscal_horiz_adv_type   =
erf.dryscal_vert_adv_type    =
erf.moistscal_horiz_adv_type =
erf.moistscal_vert_adv_type  =

Ref: Muñoz-Esparza, D., Sauer, J. A., Jensen, A. A., Xue, L., & Grabowski, W. W. (2022). The FastEddy® resident-GPU accelerated large-eddy simulation framework: Moist dynamics extension, validation and sensitivities of modeling non-precipitating shallow cumulus clouds. Journal of Advances in Modeling Earth Systems, 14, e2021MS002904. doi:10.1029/2021MS002904

Momentum, Thermal, and Scalar Diffusion Contribution to DNS¶

Strain Rate Tensor¶

The schematic below (isomeric view) shows the definition of strain-rate components.

Strain-Rate Components for X-Momentum Equation¶

\[\begin{split}\begin{array}{ll} S_{11,i + \frac{1}{2}} = & \frac{1}{\Delta x}\left( u_{i + 1,j,k} - u_{i,j,k} \right) \\ S_{11,i - \frac{1}{2}} = & \frac{1}{\Delta x}\left( u_{i,j,k} - u_{i - 1,j,k} \right) \\ S_{12,j + \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta y}\left( u_{i,j + 1,k} - u_{i,j,k} \right) + \frac{1}{\Delta x}\left( v_{i,j + 1,k} - v_{i - 1,j + 1,k} \right) \right\rbrack \\ S_{12,j - \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta y}\left( u_{i,j,k} - u_{i,j - 1,k} \right) + \frac{1}{\Delta x}\left( v_{i,j,k} - v_{i - 1,j,k} \right) \right\rbrack \\ S_{13,k + \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta z}\left( u_{i,j,k + 1} - u_{i,j,k} \right) + \frac{1}{\Delta x}\left( w_{i,j,k + 1} - w_{i - 1,j,k + 1} \right) \right\rbrack \\ S_{13,k - \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta z}\left( u_{i,j,k} - u_{i,j,k - 1} \right) + \frac{1}{\Delta x}\left( w_{i,j,k} - w_{i - 1,j,k} \right) \right\rbrack \\ \end{array}\end{split}\]

Strain-Rate Components for Y-Momentum Equation¶

\[\begin{split}\begin{array}{ll} S_{21,i + \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta y}\left( u_{i + 1,j,k} - u_{i + 1,j - 1,k} \right) + \frac{1}{\Delta x}\left( v_{i + 1,j,k} - v_{i,j,k} \right) \right\rbrack \\ S_{21,i - \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta y}\left( u_{i,j,k} - u_{i,j - 1,k} \right) + \frac{1}{\Delta x}\left( v_{i,j,k} - v_{i - 1,j,k} \right) \right\rbrack \\ S_{22,j + \frac{1}{2}} = & \frac{1}{\Delta y}\left( v_{i,j + 1,k} - v_{i,j,k} \right) \\ S_{22,j - \frac{1}{2}} = & \frac{1}{\Delta y}\left( v_{i,j,k} - v_{i,j - 1,k} \right) \\ S_{23,k + \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta z}\left( v_{i,j,k + 1} - v_{i,j,k} \right) + \frac{1}{\Delta y}\left( w_{i,j,k + 1} - w_{i,j - 1,k + 1} \right) \right\rbrack \\ S_{23,k - \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta z}\left( v_{i,j,k} - v_{i,j,k - 1} \right) + \frac{1}{\Delta y}\left( w_{i,j,k} - w_{i,j - 1,k} \right) \right\rbrack \\ \end{array}\end{split}\]

Strain-Rate Components for Z-Momentum Equation¶

\[\begin{split}\begin{array}{ll} S_{31,i + \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta z}\left( u_{i + 1,j,k} - u_{i + 1,j,k - 1} \right) + \frac{1}{\Delta x}\left( w_{i + 1,j,k} - w_{i,j,k} \right) \right\rbrack \\ S_{31,i - \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta z}\left( u_{i,j,k} - u_{i,j,k - 1} \right) + \frac{1}{\Delta x}\left( w_{i,j,k} - w_{i - 1,j,k} \right) \right\rbrack \\ S_{32,j + \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta z}\left( v_{i,j + 1,k} - v_{i,j + 1,k - 1} \right) + \frac{1}{\Delta y}\left( w_{i,j + 1,k} - w_{i,j,k} \right) \right\rbrack \\ S_{32,j - \frac{1}{2}} = & \frac{1}{2}\left\lbrack \frac{1}{\Delta z}\left( v_{i,j,k} - v_{i,j,k - 1} \right) + \frac{1}{\Delta y}\left( w_{i,j,k} - w_{i,j - 1,k} \right) \right\rbrack \\ S_{33,k + \frac{1}{2}} = & \frac{1}{\Delta z}\left( w_{i,j,k + 1} - w_{i,j,k} \right) \\ S_{33,k - \frac{1}{2}} = & \frac{1}{\Delta z}\left( w_{i,j,k} - w_{i,j,k - 1} \right) \\ \end{array}\end{split}\]

Expansion-Rate Tensor¶

Expansion-Rate Components for X-Momentum Equation¶

\[\begin{split}\begin{array}{ll} D_{11,i + \frac{1}{2}} = & \frac{1}{3}\left\lbrack \frac{1}{\Delta x}\left( u_{i + 1,j,k} - u_{i,j,k} \right) + \frac{1}{\Delta y}\left( v_{i,j + 1,k} - v_{i,j,k} \right) + \frac{1}{\Delta z}\left( w_{i,j,k + 1} - w_{i,j,k} \right) \right\rbrack \\ D_{11,i - \frac{1}{2}} = & \frac{1}{3}\left\lbrack \frac{1}{\Delta x}\left( u_{i,j,k} - u_{i-1,j,k} \right) + \frac{1}{\Delta y}\left( v_{i-1,j + 1,k} - v_{i-1,j,k} \right) + \frac{1}{\Delta z}\left( w_{i-1,j,k + 1} - w_{i-1,j,k} \right) \right\rbrack \\ D_{12,j + \frac{1}{2}} = & 0 \\ D_{12,j - \frac{1}{2}} = & 0 \\ D_{13,k + \frac{1}{2}} = & 0 \\ D_{13,k - \frac{1}{2}} = & 0 \\ \end{array}\end{split}\]

Expansion-Rate Components for Y-Momentum Equation¶

\[\begin{split}\begin{array}{ll} D_{21,i + \frac{1}{2}} = & 0 \\ D_{21,i - \frac{1}{2}} = & 0 \\ D_{22,j + \frac{1}{2}} = & \frac{1}{3}\left\lbrack \frac{1}{\Delta x}\left( u_{i + 1,j,k} - u_{i,j,k} \right) + \frac{1}{\Delta y}\left( v_{i,j + 1,k} - v_{i,j,k} \right) + \frac{1}{\Delta z}\left( w_{i,j,k + 1} - w_{i,j,k} \right) \right\rbrack \\ D_{22,j - \frac{1}{2}} = & \frac{1}{3}\left\lbrack \frac{1}{\Delta x}\left( u_{i + 1,j-1,k} - u_{i,j-1,k} \right) + \frac{1}{\Delta y}\left( v_{i,j,k} - v_{i,j-1,k} \right) + \frac{1}{\Delta z}\left( w_{i,j-1,k + 1} - w_{i,j-1,k} \right) \right\rbrack \\ D_{23,k + \frac{1}{2}} = & 0 \\ D_{23,k - \frac{1}{2}} = & 0 \\ \end{array}\end{split}\]

Expansion-Rate Components for Z-Momentum Equation¶

\[\begin{split}\begin{array}{ll} D_{21,i + \frac{1}{2}} = & 0 \\ D_{21,i - \frac{1}{2}} = & 0 \\ D_{22,j + \frac{1}{2}} = & 0 \\ D_{22,j - \frac{1}{2}} = & 0 \\ D_{23,k + \frac{1}{2}} = & \frac{1}{3}\left\lbrack \frac{1}{\Delta x}\left( u_{i + 1,j,k} - u_{i,j,k} \right) + \frac{1}{\Delta y}\left( v_{i,j + 1,k} - v_{i,j,k} \right) + \frac{1}{\Delta z}\left( w_{i,j,k + 1} - w_{i,j,k} \right) \right\rbrack \\ D_{23,k - \frac{1}{2}} = & \frac{1}{3}\left\lbrack \frac{1}{\Delta x}\left( u_{i + 1,j,k-1} - u_{i,j,k-1} \right) + \frac{1}{\Delta y}\left( v_{i,j + 1,k-1} - v_{i,j,k-1} \right) + \frac{1}{\Delta z}\left( w_{i,j,k} - w_{i,j,k-1} \right) \right\rbrack \\ \end{array}\end{split}\]

U Momentum viscous stress divergence¶

\[\begin{split}\begin{align} \left( \rho u \right)_{i,j,k}^{n + 1} = \left( \rho u \right)_{i,j,k}^{n} - & \\ \Delta t & \left\{ \frac{1}{\Delta x}\ \left\lbrack \tau_{11,i + \frac{1}{2}} - \tau_{11,i - \frac{1}{2}} \right\rbrack \right.\ \\ & + \frac{1}{\Delta y}\ \left\lbrack \tau_{12,j + \frac{1}{2}} - \tau_{12,j - \frac{1}{2}} \right\rbrack \\ & + \left. \frac{1}{\Delta z}\ \left\lbrack \tau_{13,k + \frac{1}{2}} - \tau_{13,k - \frac{1}{2}} \right\rbrack \right\} \\ \end{align}\end{split}\]

Note that LES equation has a similar format for computing \(\tau_{11,i + \frac{1}{2}}\), \(\tau_{11,i - \frac{1}{2}}\), \(\tau_{12,j + \frac{1}{2}}\), \(\tau_{12,j - \frac{1}{2}}\), \(\tau_{13,k + \frac{1}{2}}\), and \(\tau_{13,k - \frac{1}{2}}\).

\[\begin{split}\begin{array}{ll} \tau_{ij,m + \frac{1}{2}} = & -\mu_{eff} \ \left\lbrack S_{ij,m + \frac{1}{2}} - D_{ij,m + \frac{1}{2}} \right\rbrack \\ \tau_{ij,m - \frac{1}{2}} = & -\mu_{eff} \ \left\lbrack S_{ij,m - \frac{1}{2}} - D_{ij,m - \frac{1}{2}} \right\rbrack \end{array}\end{split}\]

where \(m = i, j, k\).

\(\mu_{eff} = 2\mu\) if a constant molecular diffusion type is chosen (erf.molec_diff_type = "Constant"). Otherwise (erf.molec_diff_type = "None"), \(\mu_{eff} = 0\).

The nomenclature is similar for other two momentum equations. Note that \(\mu\) is constant in the current implementation and its variation with temperature for low-Mach atmospheric flows has been ignored.

V Momentum viscous stress divergence¶

\[\begin{split}\begin{align} \left( \rho v \right)_{i,j,k}^{n + 1} = \left( \rho v \right)_{i,j,k}^{n} - & \\ \Delta t & \left\{ \frac{1}{\Delta x} \left\lbrack \tau_{21,i + \frac{1}{2}} - \tau_{21,i - \frac{1}{2}} \right\rbrack \right. \\ & + \frac{1}{\Delta y} \left\lbrack \tau_{22,j + \frac{1}{2}} - \tau_{22,j - \frac{1}{2}} \right\rbrack \\ & + \left. \frac{1}{\Delta z} \left\lbrack \tau_{23,k + \frac{1}{2}} - \tau_{23,k - \frac{1}{2}} \right\rbrack \right\} \end{align}\end{split}\]

Note that LES equation has a similar format for computing \(\tau_{21,i + \frac{1}{2}}\), \(\tau_{21,i - \frac{1}{2}}\), \(\tau_{22,j + \frac{1}{2}}\), \(\tau_{22,j - \frac{1}{2}}\), \(\tau_{23,k + \frac{1}{2}}\), and \(\tau_{23,k - \frac{1}{2}}\).

W Momentum viscous stress divergence¶

\[\begin{split}\begin{align} \left( \rho w \right)_{i,j,k}^{n + 1} = \left( \rho w \right)_{i,j,k}^{n} - & \\ \Delta t & \left\{ \frac{1}{\Delta x} \left\lbrack \tau_{31,i + \frac{1}{2}} - \tau_{31,i - \frac{1}{2}} \right\rbrack \right. \\ & + \frac{1}{\Delta y} \left\lbrack \tau_{32,j + \frac{1}{2}} - \tau_{32,j - \frac{1}{2}} \right\rbrack \\ & + \left. \frac{1}{\Delta z} \left\lbrack \tau_{33,k + \frac{1}{2}} - \tau_{33,k - \frac{1}{2}} \right\rbrack \right\} \end{align}\end{split}\]

Note that LES equation has a similar format for computing \(\tau_{31,i + \frac{1}{2}}\), \(\tau_{31,i - \frac{1}{2}}\), \(\tau_{32,j + \frac{1}{2}}\), \(\tau_{32,j - \frac{1}{2}}\), \(\tau_{33,k + \frac{1}{2}}\), and \(\tau_{33,k - \frac{1}{2}}\).

Potential Temperature Diffusion¶

\[\begin{split}\begin{matrix} \left( \rho \theta \right)_{i,j,k}^{n + 1} & = & \left( \rho \theta \right)_{i,j,k}^{n} & + & \Delta t\rho_{i,j,k}\alpha_{T} & \left\{ \frac{1}{{\Delta x}^{2}}\ \left\lbrack \theta_{i + 1,j,k}^{n} - \ {2\theta}_{i,j,k}^{n} + \ \theta_{i - 1,j,k}^{n} \right\rbrack \right.\ \\ & & & & & + \frac{1}{{\Delta y}^{2}}\left\lbrack \theta_{i,j + 1,k}^{n} - \ 2\theta_{i,j,k}^{n} + \ \theta_{i,j - 1,k}^{n} \right\rbrack \\ & & & & & + \left. \frac{1}{{\Delta z}^{2}}\left\lbrack \theta_{i,j,k + 1}^{n} - \ {2\theta}_{i,j,k}^{n} + \ \theta_{i,j,k - 1}^{n} \right\rbrack \right\} \end{matrix}\end{split}\]

Scalar Diffusion¶

\[\begin{split}\begin{matrix} \left( \rho C \right)_{i,j,k}^{n + 1} & = & \left( \rho C \right)_{i,j,k}^{n} & + & \Delta t\rho_{i,j,k}\alpha_{C} & \left\{ \frac{1}{{\Delta x}^{2}}\ \left\lbrack C_{i + 1,j,k}^{n} - \ 2C_{i,j,k}^{n} + \ C_{i - 1,j,k}^{n} \right\rbrack \right.\ \\ & & & & & + \frac{1}{{\Delta y}^{2}}\left\lbrack C_{i,j + 1,k}^{n} - \ 2C_{i,j,k}^{n} + \ C_{i,j - 1,k}^{n} \right\rbrack \\ & & & & & + \left. \frac{1}{{\Delta z}^{2}}\left\lbrack C_{i,j,k + 1}^{n} - \ 2C_{i,j,k}^{n} + \ C_{i,j,k - 1}^{n} \right\rbrack \right\} \end{matrix}\end{split}\]

Note: In WRF, the diffusion coefficients specified in the input file (\(K_h\) and \(K_v\) for horizontal and vertical diffusion) get divided by the Prandtl number for the potential temperature and the scalars. For the momentum, they are used as it is. In ERF, the coefficients specified in the inputs (\(\alpha_T\) and \(\alpha_C\)) are used as it is, and no division by Prandtl number is done.

Momentum, Thermal, and Scalar Diffusion Contribution to LES¶

Strain Rate and Eddy Viscosity¶

The goal is to compute eddy viscosity at the cell centers and interpolated them to the edges. Refer again to the strain rate tensor schematic.

\[\begin{split}\begin{array}{ll} S_{11} = & S_{11i + \frac{1}{2}} \\ S_{22} = & S_{22j + \frac{1}{2}} \\ S_{33} = & S_{33k + \frac{1}{2}} \end{array}\end{split}\]

\[\begin{split}\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} \end{matrix}\end{split}\]

Note that:

\[\begin{split}\begin{array}{ll} S_{12} = & S_{21} \\ S_{13} = & S_{31} \\ S_{23} = & S_{32} \end{array}\end{split}\]

\(K_{i,j,k} = {2\left( C_{S}\ \Delta \right)^{2}\left( {2S}_{mn}S_{mn} \right)}^{\frac{1}{2}} \ \rho_{i,j,k}\), where

\[\begin{split}S_{mn}S_{mn} = S_{11}^{2} + S_{22}^{2} + S_{33}^{2} + S_{12}^{2} + S_{13}^{2} + S_{23}^{2} + S_{21}^{2} + S_{31}^{2} + S_{32}^{2} \\\end{split}\]

Owing to symmetry we need to compute 6 of the 9 tensor components.

\(K_{i,j,k} = 2 \ {\mu_{t}}_{i, j, k}\) is the modeled turbulent viscosity and can be considered analogous to \(2\mu_{i, j, k}\).

\(\mu_{eff} = 2\mu + K = 2\mu + 2\mu_{t}\) if a constant molecular diffusion type is chosen (erf.molec_diff_type = "Constant"). Otherwise (erf.molec_diff_type = "None"), \(\mu_{eff} = K = 2\mu_{t}\).

The interpolated values of eddy-viscosity at the edges are the average of the values at the centers of the 4 cells the edge is part of.

\[\begin{split}\begin{array}{c} K_{i + \frac{1}{2},j - \frac{1}{2},k} = \frac{1}{4}\left\lbrack K_{i,j - 1,k} + K_{i,j,k} + K_{i + 1,j - 1,k} + K_{i + 1,j,k} \right\rbrack \\ K_{i + \frac{1}{2},j + \frac{1}{2},k} = \frac{1}{4}\left\lbrack K_{i,j,k} + K_{i,j + 1,k} + K_{i + 1,j,k} + K_{i + 1,j + 1,k} \right\rbrack \\ K_{i + \frac{1}{2},j,k - \frac{1}{2}} = \frac{1}{4}\left\lbrack K_{i,j,k} + K_{i,j,k - 1} + K_{i + 1,j,k} + K_{i + 1,j,k - 1} \right\rbrack \\ K_{i + \frac{1}{2},j,k + \frac{1}{2}} = \frac{1}{4}\left\lbrack K_{i,j,k + 1} + K_{i,j,k} + K_{i + 1,j,k + 1} + K_{i + 1,j,k} \right\rbrack \\ K_{i,j + \frac{1}{2},k - \frac{1}{2}} = \frac{1}{4}\left\lbrack K_{i,j,k} + K_{i,j,k - 1} + K_{i,j + 1,k} + K_{i,j + 1,k - 1} \right\rbrack \\ K_{i,j + \frac{1}{2},k + \frac{1}{2}} = \frac{1}{4}\left\lbrack K_{i,j,k} + K_{i,j,k + 1} + K_{i,j + 1,k} + K_{i,j + 1,k + 1} \right\rbrack \end{array}\end{split}\]

U Momentum - subfilter stress divergence¶

\[\begin{split}\begin{align} \left( \rho u \right)_{i,j,k}^{n + 1} = \left( \rho u \right)_{i,j,k}^{n} - & \\ \Delta t & \left\{ \frac{1}{\Delta x}\ \left\lbrack \tau_{11,i + \frac{1}{2}} - \tau_{11,i - \frac{1}{2}} \right\rbrack \right.\ \\ & + \frac{1}{\Delta y}\ \left\lbrack \tau_{12,j + \frac{1}{2}} - \tau_{12,j - \frac{1}{2}} \right\rbrack \\ & + \left. \frac{1}{\Delta z}\ \left\lbrack \tau_{13,k + \frac{1}{2}} - \tau_{13,k - \frac{1}{2}} \right\rbrack \right\} \\ \end{align}\end{split}\]

\[\begin{split}\begin{array}{l} \tau_{11,i + \frac{1}{2}} = -K_{i,j,k} \ S_{11,i + \frac{1}{2}} \\ \tau_{11,i - \frac{1}{2}} = -K_{i - 1,j,k} \ S_{11,i - \frac{1}{2}} \\ \tau_{12,j + \frac{1}{2}} = -K_{i - \frac{1}{2},j + \frac{1}{2},k} \ S_{12,j + \frac{1}{2}} \\ \tau_{12,j - \frac{1}{2}} = -K_{i - \frac{1}{2},j - \frac{1}{2},k} \ S_{12,j - \frac{1}{2}} \\ \tau_{13,k + \frac{1}{2}} = -K_{i - \frac{1}{2},j,k + \frac{1}{2}} \ S_{13,k + \frac{1}{2}} \\ \tau_{13,k - \frac{1}{2}} = -K_{i - \frac{1}{2},j,k - \frac{1}{2}} \ S_{13,k - \frac{1}{2}} \end{array}\end{split}\]

V Momentum - subfilter stress divergence¶

\[\begin{split}\begin{align} \left( \rho v \right)_{i,j,k}^{n + 1} = \left( \rho v \right)_{i,j,k}^{n} - & \\ \Delta t & \left\{ \frac{1}{\Delta x} \left\lbrack \tau_{21,i + \frac{1}{2}} - \tau_{21,i - \frac{1}{2}} \right\rbrack \right. \\ & + \frac{1}{\Delta y} \left\lbrack \tau_{22,j + \frac{1}{2}} - \tau_{22,j - \frac{1}{2}} \right\rbrack \\ & + \left. \frac{1}{\Delta z} \left\lbrack \tau_{23,k + \frac{1}{2}} - \tau_{23,k - \frac{1}{2}} \right\rbrack \right\} \end{align}\end{split}\]

\[\begin{split}\begin{array}{ll} \tau_{21,i + \frac{1}{2}} = -K_{i + \frac{1}{2},j - \frac{1}{2},k} \ S_{21,i + \frac{1}{2}} \\ \tau_{21,i - \frac{1}{2}} = -K_{i - \frac{1}{2},j - \frac{1}{2},k} \ S_{21,i - \frac{1}{2}} \\ \tau_{22,j + \frac{1}{2}} = -K_{i,j,k} \ S_{22,j + \frac{1}{2}} \\ \tau_{22,j - \frac{1}{2}} = -K_{i,j - 1,k} \ S_{22,j - \frac{1}{2}} \\ \tau_{23,k + \frac{1}{2}} = -K_{i,j - \frac{1}{2},k + \frac{1}{2}} \ S_{23,k + \frac{1}{2}} \\ \tau_{23,k - \frac{1}{2}} = -K_{i,j - \frac{1}{2}k - \frac{1}{2}} \ S_{23,k - \frac{1}{2}} \end{array}\end{split}\]

W Momentum - subfilter stress divergence¶

\[\begin{split}\begin{align} \left( \rho w \right)_{i,j,k}^{n + 1} = \left( \rho w \right)_{i,j,k}^{n} - & \\ \Delta t & \left\{ \frac{1}{\Delta x} \left\lbrack \tau_{31,i + \frac{1}{2}} - \tau_{31,i - \frac{1}{2}} \right\rbrack \right. \\ & + \frac{1}{\Delta y} \left\lbrack \tau_{32,j + \frac{1}{2}} - \tau_{32,j - \frac{1}{2}} \right\rbrack \\ & + \left. \frac{1}{\Delta z} \left\lbrack \tau_{33,k + \frac{1}{2}} - \tau_{33,k - \frac{1}{2}} \right\rbrack \right\} \end{align}\end{split}\]

\[\begin{split}\begin{array}{ll} \tau_{31,i + \frac{1}{2}} = -K_{i + \frac{1}{2},j,k - \frac{1}{2}} \ S_{31,i + \frac{1}{2}} \\ \tau_{31,i - \frac{1}{2}} = -K_{i - \frac{1}{2},j,k - \frac{1}{2}} \ S_{31,i - \frac{1}{2}} \\ \tau_{32,j + \frac{1}{2}} = -K_{i,j + \frac{1}{2},k - \frac{1}{2}} \ S_{32,j + \frac{1}{2}} \\ \tau_{32,j - \frac{1}{2}} = -K_{i,j - \frac{1}{2},k - \frac{1}{2}} \ S_{32,j - \frac{1}{2}} \\ \tau_{33,k + \frac{1}{2}} = -K_{i,j,k} \ S_{33,k + \frac{1}{2}} \\ \tau_{33,k - \frac{1}{2}} = -K_{i,j, k - 1} \ S_{33,k - \frac{1}{2}} \end{array}\end{split}\]

Energy Conservation- Subgrid heat flux¶

\[\begin{split}\begin{align} \left( \rho \theta \right)_{i,j,k}^{n + 1} = \left( \rho \theta \right)_{i,j,k}^{n} - & \\ & \Delta t \left\{ \frac{1}{\Delta x^2} \left\lbrack K_{i+\frac{1}{2},j,k} (\theta_{i+1,j,k}^{n} - \theta_{i ,j,k}^{n}) - K_{i-\frac{1}{2},j,k} (\theta_{i ,j,k}^{n} - \theta_{i-1,j,k}^{n}) \right\rbrack \right. \\ & + \frac{1}{\Delta y^2} \left\lbrack K_{i,j+\frac{1}{2},k} (\theta_{i,j+1,k}^{n} - \theta_{i,j ,k}^{n}) - K_{i,j-\frac{1}{2},k} (\theta_{i,j ,k}^{n} - \theta_{i,j-1,k}^{n}) \right\rbrack \\ & \left. + \frac{1}{\Delta z^2} \left\lbrack K_{i,j,k-\frac{1}{2}} (\theta_{i,j,k+1}^{n} - \theta_{i,j,k }^{n}) K_{i,j,k-\frac{1}{2}} (\theta_{i,j,k }^{n} - \theta_{i,j,k-1}^{n}) \right\rbrack \right\} \end{align}\end{split}\]

Analogous formulae apply for the subgrid scalar flux.

Prognostic Equation for Subgrid Kinetic Energy¶

\[\begin{split}\begin{align} \left( \rho e \right)_{i,j,k}^{n + 1} = \left( \rho e \right)_{i,j,k}^{n} - & \\ \Delta t & \left\{ \frac{1}{\Delta x}\left\lbrack \left( \rho u \right)_{i+1,j,k}^{n}e_{i+\frac{1}{2},j,k}^{n} - \left( \rho u \right)_{i ,j,k}^{n}e_{i-\frac{1}{2},j,k}^{n} \right\rbrack \right. \\ & + \frac{1}{\Delta y}\left\lbrack \left( \rho v \right)_{i,j+1,k}^{n}e_{i,j+\frac{1}{2},k}^{n} - \left( \rho v \right)_{i,j ,k}^{n}e_{i,j-\frac{1}{2},k}^{n} \right\rbrack \\ & + \frac{1}{\Delta z}\left\lbrack \left( \rho w \right)_{i,j,k+1}^{n}e_{i,j,k+\frac{1}{2}}^{n} - \left( \rho w \right)_{i,j,k }^{n}e_{i,j,k-\frac{1}{2}}^{n} \right\rbrack \\ & + \left. \rho_{i,j,k} K_H \frac{g}{\theta_{i,j,k}} (\frac{\partial\theta}{\partial z})_{i,j,k} - \tau_{mn}\frac{\partial u_{m}}{\partial x_{n}} + (\nabla \cdot (K \nabla e))_{i,j,k} - \rho_{i,j,k} C_{\epsilon} \frac{\left( e_{i,j,k} \right)^{\frac{3}{2}}}{\mathcal{l}} \right\} \end{align}\end{split}\]

where

\[\begin{split}\begin{array}{l} C_{\epsilon} = 0.19 + 0.51\frac{\mathcal{l}}{\Delta s} \\ K_{H} = \left( 1 + 2\frac{\mathcal{l}}{\Delta s} \right)K_{M} \\ K_{M} = 0.1 \mathcal{l} e^{\frac{1}{2}} \end{array}\end{split}\]

For convective or neutral cases, \(\theta^{n}_{i,j,k+1}-\theta^{n}_{i,j,k-1} < 0\), \(\mathcal{l} = \Delta s = \sqrt[3]{\Delta x \Delta y \Delta z}\)

For stable cases, \(\theta^{n}_{i,j,k+1}-\theta^{n}_{i,j,k-1} > 0\),

\[\begin{split}\begin{align} \mathcal{l} & = 0.76\ e^{\frac{1}{2}}\left( \frac{g}{\theta}\frac{\partial\theta}{\partial z} \right) \\ & = 0.76e_{i,j,k}^{\frac{1}{2}}\left\lbrack \frac{g}{\theta_{i,j,k}}\frac{1}{2\Delta z}\left( \theta_{i,j,k + 1}^{n} - \theta_{i,j,k - 1}^{n} \right) \right\rbrack \end{align}\end{split}\]

The second to last term on the right hand side takes the form

\[\begin{split}\begin{align} & \frac{1}{\Delta x^2} \left\lbrack K_{i+\frac{1}{2},j,k} (e_{i+1,j,k}^{n} - e_{i ,j,k}^{n}) - K_{i-\frac{1}{2},j,k} (e_{i ,j,k}^{n} - e_{i-1,j,k}^{n}) \right\rbrack + \\ & \frac{1}{\Delta y^2} \left\lbrack K_{i,j+\frac{1}{2},k} (e_{i,j+1,k}^{n} - e_{i,j ,k}^{n}) - K_{i,j-\frac{1}{2},k} (e_{i,j ,k}^{n} - e_{i,j-1,k}^{n}) \right\rbrack + \\ & \frac{1}{\Delta z^2} \left\lbrack K_{i,j,k-\frac{1}{2}} (e_{i,j,k+1}^{n} - e_{i,j,k }^{n}) - K_{i,j,k-\frac{1}{2}} (e_{i,j,k }^{n} - e_{i,j,k-1}^{n}) \right\rbrack \end{align}\end{split}\]

Finally,

\[\begin{split}\begin{align} \tau_{mn}\frac{\partial u_{m}}{\partial x_{n}} & = K S_{mn}\frac{\partial u_{m}}{\partial x_{n}} \\ & = K S_{mn}S_{mn} \\ & = K (S_{11}^{2} + S_{22}^{2} + S_{33}^{2} + S_{12}^{2} + S_{13}^{2} + S_{23}^{2} + S_{21}^{2} + S_{31}^{2} + S_{32}^{2}) \end{align}\end{split}\]

Contributions from different directions¶

Mesh Refinement¶

ERF allows both static and dynamic mesh refinement, as well as the choice of one-way or two-way coupling. Dynamic refinement is currently only allowed when terrain is not being used.

The refinement ratio is specified by the user at runtime. Refinement is not allowed in the vertical.

Note that any tagged region will be covered by one or more boxes. The user may specify the refinement criteria and/or region to be covered, but not the decomposition of the region into individual grids.

See the Gridding section of the AMReX documentation for details of how individual grids are created.

Coupling Types¶

ERF supports one-way and two-way coupling between levels; this is a run-time input

erf.coupling_type = "OneWay" or "TwoWay"

By one-way coupling, we mean that between each pair of refinement levels, the coarse level communicates data to the fine level to serve as boundary conditions for the time advance of the fine solution. For cell-centered quantities, and face-baced normal momenta on the coarse-fine interface, the coarse data is conservatively interpolated to the fine level.

The interpolated data is utilized to specify ghost cell data (outside of the valid fine region) as well as specified data inside the lateral boundaries of the fine region. See Domain Boundary Conditions for the details of how the relaxation works; when used in the context of mesh refinement we fill the specified values by interpolation from the coarser level rather than reading from the external file. For coarse/fine boundaries, a user may specify the total width of the interior specified (Dirichlet) and relaxation region with erf.cf_width = <Int> (yellow + blue) and analogously the width of the interior specified (Dirichlet) region may be specified with erf.cf_set_width = <Int> (yellow).

Setting erf.cf_set_width = 0 designates that we interpolate the momenta at faces only on the coarse-fine boundary itself; no interior cell-centered data, or momenta inside the fine region, are filled from the coarser level.

By two-way coupling, we mean that in additional to interpolating data from the coarser level to supply boundary conditions for the fine regions, the fine level also communicates data back to the coarse level in two ways:

The fine cell-centered data are conservatively averaged onto the coarse mesh covered by fine mesh.
The fine momenta are conservatively averaged onto the coarse faces covered by fine mesh.
A “reflux” operation is performed for all cell-centered data; this updates values on the coarser level outside of regions covered by the finer level.

We note that when one-way coupling is used, quantities which are advanced in conservation form potentially violate global conservation. Two-way coupling ensures conservation of mass, and of the advective contribution to all scalar updates, but does not account for loss of conservation due to diffusive or source terms.

Domain Boundary Conditions¶

Ideal Domain BCs¶

There are two primary types of physical/domain boundary conditions: those which rely only on the data in the valid regions, and those which rely on externally specified values.

ERF allows users to specify types of boundary condition with keywords in the inputs file. The information for each face is preceded by xlo, xhi, ylo, yhi, zlo, or zhi.

Currently available type of boundary conditions are inflow, outflow, slipwall, noslipwall, symmetry or MOST. (Spelling of the type matters; capitalization does not.)

For example, setting

xlo.type = "Inflow"
xhi.type = "Outflow"
zlo.type = "SlipWall"
zhi.type = "SlipWall"

geometry.is_periodic = 0 1 0

would define a problem with inflow in the low-\(x\) direction, outflow in the high-\(x\) direction, periodic in the \(y\)-direction, and slip wall on the low and high \(y\)-faces, and Note that no keyword is needed for a periodic boundary, here only the specification in geometry.is_periodic is needed.

Each of these types of physical boundary condition has a mapping to a mathematical boundary condition for each type; this is summarized in the table below.

ERF provides the ability to specify a variety of boundary conditions (BCs) in the inputs file. We use the following options preceded by xlo, xhi, ylo, yhi, zlo, and zhi:

Type	Normal vel	Tangential vel	Density	Theta	Scalar
inflow	ext_dir	ext_dir	ext_dir	ext_dir	ext_dir
outflow	foextrap	foextrap	foextrap	foextrap	foextrap
slipwall	ext_dir	foextrap	foextrap	ext_dir/foextrap/neumann	foextrap
noslipwall	ext_dir	ext_dir	foextrap	ext_dir/foextrap/neumann	foextrap
symmetry	reflect_odd	reflect_even	reflect_even	reflect_even	reflect_even
MOST

Here ext_dir, foextrap, and reflect_even refer to AMReX keywords. The ext_dir type refers to an “external Dirichlet” boundary, which means the values must be specified by the user. The foextrap type refers to “first order extrapolation” which sets all the ghost values to the same value in the last valid cell/face (AMReX also has a hoextrap, or “higher order extrapolation” option, which does a linear extrapolation from the two nearest valid values). By contrast, neumann is an ERF-specific boundary type that allows a user to specify a variable gradient. Currently, the neumann BC is only supported for theta to allow for weak capping inversion (\(\partial \theta / \partial z\)) at the top domain.

As an example,

xlo.type                =   "Inflow"
xlo.velocity            =   1. 0.9  0.
xlo.density             =   1.
xlo.theta               =   300.
xlo.scalar              =   2.

sets the boundary condtion 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, xlo.theta = 300. sets the inflow potential temperature, xlo.scalar = 2. sets the inflow value of the advected scalar

The “slipwall” and “noslipwall” types have options for adiabatic vs Dirichlet boundary conditions. If a value for theta is given for a face with type “slipwall” or “noslipwall” then the boundary condition for theta is assumed to be “ext_dir”, i.e. theta is specified on the boundary. If no value is specified then the wall is assumed to be adiabiatc, i.e. there is no temperature flux at the boundary. This is enforced with the “foextrap” designation.

For example

zlo.type  = "NoSlipWall"
zhi.type  = "NoSlipWall"

zlo.theta = 301.0

would designate theta = 301 at the bottom (zlo) boundary, while the top boundary condition would default to a zero gradient (adiabatic) since no value is specified for zhi.theta or zhi.theta_grad. By contrast, thermal inversion may be imposed at the top boundary by providing a specified gradient for theta

zlo.type  = "NoSlipWall"
zhi.type  = "NoSlipWall"

zlo.theta = 301.0
zhi.theta_grad = 1.0

We note that “noslipwall” allows for non-zero tangential velocities to be specified, as in the Couette regression test example, in which we specify

geometry.is_periodic = 1 1 0

zlo.type = "NoSlipWall"
zhi.type = "NoSlipWall"

zlo.velocity    = 0.0 0.0 0.0
zhi.velocity    = 2.0 0.0 0.0

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 MOST Boundaries 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 state data.

Real Domain BCs¶

When using real lateral boundary conditions, time-dependent observation data is read from a file. The observation data is utilized to directly set Dirichlet values on the lateral domain BCs as well as nudge the solution state towards the observation data. The user may specify (in the inputs file) the total width of the interior Dirichlet and relaxation region with erf.real_width = <Int> (yellow + blue) and analogously the width of the interior Dirichlet region may be specified with erf.real_set_width = <Int> (yellow). The real BCs are only imposed for \(\psi = \left\{ \rho; \; \rho \theta; \; \rho q_v; \; u; \; v \right\}\). Due to the staggering of scalars (cell center) and velocities (face center) with an Arakawa C grid, we reduce the relaxation width of the scalars \(\left\{ \rho; \; \rho \theta; \; \rho q_v \right\}\) by 1 to ensure the momentum updates at the last relaxation cell involve a pressure gradient that is computed with relaxed and non-relaxed data.

Lateral boundaries¶

Image taken from Skamarock et al. (2021)

Within the interior Dirichlet cells, the RHS is exactly \(\psi^{n} - \ps^{BDY} / \Delta t\) and, as such, we directly impose this value in the yellow region. Within the relaxation region (blue), the RHS (\(F\)) is given by the following:

\[\begin{split}\begin{align} F &= G + R, \\ R &= \left[ H_{1} \left( \psi^{BDY} - \psi^{\*} \right) - H_{2} \Delta^2 \left( \psi^{BDY} - \psi^{\*} \right) \right] \exp \left(-C_{01} \left(n - {\rm SpecWidth}\right) \right), \\ H_{1} &= \frac{1}{10 \Delta t} \frac{{\rm SpecWidth} + {\rm RelaxWidth} - n}{{\rm RelaxWidth} - 1}, \\ H_{2} &= \frac{1}{50 \Delta t} \frac{{\rm SpecWidth} + {\rm RelaxWidth} - n}{{\rm RelaxWidth} - 1}, \end{align}\end{split}\]

where \(G\) is the RHS of the Navier-Stokes equations, \(\psi^{*}\) is the state variable at the current RK stage, \(\psi^{BDY}\) is temporal interpolation of the observational data, \(C_{01} = -\ln(0.01) / ({\rm RealWidth - SpecWidth})\) is a constant that ensure the exponential blending function obtains a value of 0.01 at the last relaxation cell, and \(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.

\[\frac{dQ}{dt} = \mathrm{RHS} - A\xi^n(Q-Q_\mathrm{target})\]

where RHS are the other right-hand side terms. The parameters to be set by the user are – A is the sponge amplitude, n is the sponge strength and the \(Q_\mathrm{target}\) – the target solution in the sponge. xi is a linear coordinate that is 0 at the beginning of the sponge and 1 at the end. An example of the sponge inputs can be found in Exec/RegTests/Terrain2d_Cylinder.

erf.sponge_strength = 10000.0
erf.use_xlo_sponge_damping = true
erf.xlo_sponge_end = 4.0
erf.use_xhi_sponge_damping = true
erf.xhi_sponge_start = 26.0
erf.use_zhi_sponge_damping = true
erf.zhi_sponge_start = 8.0

erf.sponge_density = 1.2
erf.sponge_x_velocity = 10.0
erf.sponge_y_velocity = 0.0
erf.sponge_z_velocity = 0.0

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 (\(\overline{u^{'}w^{'}}\), \(\overline{v^{'}w^{'}}\), and \(\overline{\theta^{'}w^{'}}\)) are constant with height. \(\Phi_m\) and \(\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:

\[ \begin{align}\begin{aligned}\overline{u^{'}} \overline{w^{'}} = const = -u^{2}_{\star},\\\overline{w^{'}} \overline{\theta^{'}} = const = -u_{\star}\theta_{\star},\\\Phi_{m}(\zeta) = \frac{\kappa z}{u_{\star}} \frac{\partial \overline{u}(z)}{\partial z},\\\Phi_{h}(\zeta) = \frac{\kappa z}{u_{\star}} \frac{\partial \overline{\theta}(z)}{\partial z}\end{aligned}\end{align} \]

where the nondimensional gradients are expressed in terms of the MOST stability parameter, \(\zeta = \frac{z}{L} = -\frac{\kappa z}{u_{\star}^{3}} \frac{g}{\overline{\theta}} \overline{w^{'}\theta^{'}}\), which serves as a surface layer scaling parameter. Here, \(L\) is the Monin-Obukhov length, \(u_{\star}\) is the friction velocity (defined for \(u\) aligned with the wind direction), \(\theta_{\star}\) is the surface layer temperature scale, \(\overline{\theta}\) is the reference virtual potential temperature for the ASL, and \(\kappa\) is the von Karman constant (taken to be \(0.41\)).

Integration of the MOST assumption equations give the classical MOST profiles of mean velocity and potential temperature

\[ \begin{align}\begin{aligned}\overline{u}(z) &= \frac{u_{\star}}{\kappa} \left[ \mathrm{ln} \left(\frac{z}{z_0}\right) - \Psi_m(\zeta)\right],\\\overline{\theta}(z) - \theta_0 &= \frac{\theta_{\star}}{\kappa} \left[ \mathrm{ln}\left(\frac{z}{z_0}\right) - \Psi_{h}(\zeta) \right]\end{aligned}\end{align} \]

where \(\theta_0\) is the surface potential temperature and \(z_0\) is a characteristic roughness height. The integrated similarity functions,

\[ \begin{align}\begin{aligned}\Psi_{m}(\zeta) &= \int_{0} ^{\frac{z}{L}} [1-\Phi_{m}(\zeta)]\frac{d\zeta}{\zeta},\\\Psi_{h}(\zeta) &= \int_{0} ^{\frac{z}{L}} [1-\Phi_{h}(\zeta)]\frac{d\zeta}{\zeta}\end{aligned}\end{align} \]

are calculated analytically from empirical gradient functions \(\Phi_m\) and \(\Phi_h\), which are defined piecewise for stable and unstable values of the stability parameter.

Unstable: \((\zeta < 0)\)

\[ \begin{align}\begin{aligned}\Phi_{m} &= (1-\gamma_{1}\zeta)^{-\frac{1}{4}}, \quad \Psi_{m} = \mathrm{ln}\left[\frac{1}{8}(1+\Phi_{m}^{-2})(1+\Phi_{m}^{-1})^{2}\right]-2\arctan(\Phi_{m}^{-1})+\frac{\pi}{2},\\\Phi_{h} &= \sigma_{\theta}(1-\gamma_{2}\zeta)^{-\frac{1}{2}}, \quad \Psi_{h} = (1+\sigma_{\theta}) \mathrm{ln} \left[\frac{1}{2}(1+\Phi_{h}^{-1}) \right]+(1-\sigma_{\theta}) {\mathrm{ln}} \left[\frac{1}{2}(-1+\Phi_{h}^{-1})\right]\end{aligned}\end{align} \]

Stable: \((\zeta > 0)\)

\[ \begin{align}\begin{aligned}\Phi_{m} &= 1+\beta \zeta, \quad \Psi_{m}=-\beta \zeta,\\\Phi_{h} &= \sigma_{\theta}+\beta \zeta, \quad \Psi_{h}=(1-\sigma_{\theta})\mathrm{ln}(\zeta)-\beta \zeta,\end{aligned}\end{align} \]

where the constants take the values proposed in Dyer, Boundary Layer Meteorology, 1974:

\[\sigma_{\theta}=1, \quad \beta = 5, \quad \gamma_{1}=16, \quad \gamma_{2}=16\]

Inverting the equations above, the MOST stability parameter,

\[\zeta=\frac{z}{L} = -\kappa z \frac{g}{\bar{\theta}} \frac{\theta_{\star}}{u^{2}_{\star}}\]

is determined by the friction velocity

\[u_{\star} = \kappa \overline{u}/[\mathrm{ln}(z/z_0)-\Psi_{m}({z}/{L})]\]

and the characteristic surface layer temperature

\[\theta_{\star} = \kappa (\overline{\theta}-\theta_0)/[\mathrm{ln}(z / z_0)-\Psi_{h}(z/L)]\]

MOST Implementation¶

In ERF, when the MOST boundary condition is applied, velocity and temperature in the ghost cells are set to give stresses that are consistent with the MOST equations laid out above. The code is structured to allow either the surface temperature (\(\theta_0\)) or surface temperature flux (\(\overline{w^{'}\theta^{'}}\)) to be enforced. To apply the MOST boundary, the following algorithm is applied:

Horizontal (planar) averages \(\bar{u}\), \(\bar{v}\) and \(\overline{\theta}\) are computed at a reference height \(z_{ref}\) assumed to be within the surface layer.
Initially, neutral conditions (\(L=\infty, \zeta=0\)) are assumed and used to compute a provisional \(u_{\star}\) using the equation given above. If \(\theta_0\) is specified, the above equation for \(\theta_{\star}\) is applied and then the surface flux is computed \(\overline{w^{'}\theta^{'}} = -u_{\star} \theta_{\star}\). If \(\overline{w^{'}\theta^{'}}\) is specified, \(\theta_{\star}\) is computed as \(-\overline{w^{'}\theta^{'}}/u_{\star}\) and the previous equation is inverted to compute \(\theta_0\).
The stability parameter \(\zeta\) is recomputed using the equation given above based on the provisional values of \(u_{\star}\) and \(\theta_{\star}\).
The previous two steps are repeated iteratively, sequentially updating the values of \(u_{\star}\) and \(\zeta\), until the change in the value of \(u_{\star}\) on each iteration falls below a specified tolerance.
Once the MOST iterations have converged, and the planar average surface flux values are known, the approach from Moeng, Journal of the Atmospheric Sciences, 1984 is applied to consistently compute local surface-normal stress/flux values (e.g., \(\tau_{xz} = - \rho \overline{u^{'}w^{'}}\)):

\[ \begin{align}\begin{aligned}\left. \frac{\tau_{xz}}{\rho} \right|_0 &= u_{\star}^{2} \frac{(u - \bar{u})|\mathbf{\bar{u}}| + \bar{u}\sqrt{u^2 + v^2} }{|\mathbf{\bar{u}}|^2},\\\left. \frac{\tau_{yz}}{\rho} \right|_0 &= u_{\star}^{2} \frac{(v - \bar{v})|\mathbf{\bar{u}}| + \bar{v}\sqrt{u^2 + v^2} }{|\mathbf{\bar{u}}|^2},\\\left. \frac{\tau_{\theta z}}{\rho} \right|_0 &= \theta_\star u_{\star} \frac{|\mathbf{\bar{u}}| ({\theta} - \overline{\theta}) + \sqrt{u^2+v^2} (\overline{\theta} - \theta_0) }{ |\mathbf{\bar{u}}| (\overline{\theta} -\theta_0) } = u_{\star} \kappa \frac{|\mathbf{\bar{u}}| ({\theta} - \overline{\theta}) + \sqrt{u^2+v^2} (\overline{\theta} - \theta_0) }{ |\mathbf{\bar{u}}| [ \mathrm{ln}(z_{ref} / z_0)-\Psi_{h}(z_{ref}/L)] }\end{aligned}\end{align} \]

where \(\bar{u}\), \(\bar{v}\) and \(\overline{\theta}\) are the plane averaged values (at \(z_{ref}\)) of the two horizontal velocity components and the potential temperature, respectively, and \(|\mathbf{\bar{u}}|\) is the plane averaged magnitude of horizontal velocity (plane averaged wind speed). We note a slight variation in the denominator 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) \(K_{m,v}\) and \(K_{\theta,v}\) as follows:

\[ \begin{align}\begin{aligned}\tau_{xz} = K_{m,v} \frac{\partial u}{\partial z}\\\tau_{yz} = K_{m,v} \frac{\partial v}{\partial z}\\\tau_{\theta z} = K_{\theta,v} \frac{\partial \theta}{\partial z}.\end{aligned}\end{align} \]

This implies that, for example, the value set for the conserved \(\rho\theta\) variable in the \(-n\mathrm{th}\) ghost cell is

\[(\rho \theta)_{i,j,-n} = \rho_{i,j,-n} \left[ \frac{(\rho\theta)_{i,j,0}}{\rho_{i,j,0}} - \left. \frac{\tau_{\theta z}}{\rho} \right|_{i,j,0} \frac{\rho_{i,j,0}}{K_{\theta,v,(i,j,0)}} n \Delta z \right]\]

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 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.

MOST Inputs¶

To evaluate the fluxes with MOST, the surface rougness parameter \(z_{0}\) must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity \(u_{\star}\). ERF supports three methods for parameterizing the surface roughness: constant, charnock, and modified_charnock. The latter two methods parameterize \(z_{0} = f(u_{\star})\) and are described in Jimenez & Dudhia, American Meteorological Society, 2018. The rougness calculation method may be specified with

erf.most.roughness_type    = STRING    #Z_0 type (constant, charnock, modified_charnock)

If the charnock method is employed, the \(a\) constant may be specified with erf.most.charnock_constant (defaults to 0.0185). If the modified_charnock method is employed, the depth \(d\) may be specified with erf.most.modified_charnock_depth (defaults to 30 m).

When computing an average \(\overline{\phi}\) for the MOST boundary, where \(\phi\) denotes a generic variable, ERF supports a variety of approaches. Specifically, planar averages and local region averages may be computed with or without time averaging. With each averaging methodology, the query point \(z\) may be determined from the following procedures: specified vertical distance \(z_{ref}\) from the bottom surface, specified \(k_{index}\), or (when employing terrain-fit coordinates) specified normal vector length \(z_{ref}\). The available inputs to the MOST boundary and their associated data types are

erf.most.average_policy    = INT    #POLICY FOR AVERAGING
erf.most.use_normal_vector = BOOL   #USE NORMAL VECTOR W/ TERRAIN?
erf.most.use_interpolation = BOOL   #INTERPOLATE QUERY POINT W/ TERRAIN?
erf.most.time_average      = BOOL   #USE TIME AVERAGING?
erf.most.z0                = FLOAT  #SURFACE ROUGHNESS
erf.most.zref              = FLOAT  #QUERY DISTANCE (HEIGHT OR NORM LENGTH)
erf.most.surf_temp         = FLOAT  #SPECIFIED SURFACE TEMP
erf.most.surf_temp_flux    = FLOAT  #SPECIFIED SURFACE FLUX
erf.most.k_arr_in          = INT    #SPECIFIED K INDEX ARRAY (MAXLEV)
erf.most.radius            = INT    #SPECIFIED REGION RADIUS
erf.most.time_window       = FLOAT  #WINDOW FOR TIME AVG

We now consider two concrete examples. To employ an instantaneous planar average at a specified vertical height above the bottom surface, one would specify:

erf.most.average_policy    = 0
erf.most.use_normal_vector = false
erf.most.time_average      = false
erf.most.z0                = 0.1
erf.most.zref              = 1.0

By contrast, local region averaging would be employed in conjunction with time averaging for the following inputs:

erf.most.average_policy    = 1
erf.most.use_normal_vector = true
erf.most.use_interpolation = true
erf.most.time_average      = true
erf.most.z0                = 0.1
erf.most.zref              = 1.0
erf.most.surf_temp_flux    = 0.0
erf.most.radius            = 1
erf.most.time_window       = 10.0

In the above case, use_normal_vector utilizes the a local surface-normal vector with length \(z_{ref}\) to construct the positions of the query points. Each query point, and surrounding points that are within erf.most.radius from the query point, are interpolated to and averaged; for a radius of 1, 27 points are averaged. The time average is completed by way of an exponential filter function whose peak coincides with the current time step and tail extends backwards in time

\[\frac{1}{\tau} \int_{-\infty}^{0} \exp{\left(t/\tau\right)} \, f(t) \; \rm{d}t.\]

Due to the form of the above integral, it is advantageous to consider \(\tau\) as a multiple of the simulation time step \(\Delta t\), which is specified by erf.most.time_window. As erf.most.time_window is reduced to 0, the exponential filter function tends to a Dirac delta function (prior averages are irrelevant). Increasing erf.most.time_window extends the tail of the exponential and more heavily weights prior averages.

Derived Variables¶

ERF has the ability to created new temporary variables derived from the state variables.

Access to the derived variable is through one of two amrex:AmrLevel functions (which are inherited by ERF)

/**
* \brief Returns a MultiFab containing the derived data for this level.
* The user is responsible for deleting this pointer when done
* with it.  If ngrow>0 the MultiFab is built on the appropriately
* grown BoxArray.
*/
virtual std::unique_ptr<MultiFab> derive (const std::string& name,
                      Real               time,
                      int                ngrow);
/**
* \brief This version of derive() fills the dcomp'th component of mf
* with the derived quantity.
*/
virtual void derive (const std::string& name,
                     Real               time,
                     MultiFab&          mf,
                     int                dcomp);

As an example, pert_prs is a derived variable provided with IAMR, which returns the perturbational pressure field. A multifab filled with the perturbational pressure can be obtained via

std::unique_ptr<MultiFab> pert_pres;
pert_pres = derive(pert_pres, time, ngrow);

Checkpoint / Restart¶

ERF has a standard sort of checkpointing and restarting capability and uses the native AMReX format for reading and writing checkpoints. In the inputs file, the following options control the generation of checkpoint files (which are really directories):

Writing the Checkpoint “Files”¶

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.check_file	prefix for restart files	String	“chk”
erf.check_int	how often (by level-0 time steps) to write restart files	Integer \(> 0\)	-1
erf.check_per	how often in simulation time to write restart files	Real \(> 0\)	-1.0

Restarting¶

Parameter	Definition	Acceptable Values	Default
erf.restart	name of the file (directory) from which to restart files	String	not used if not set

Examples of Usage¶

erf.check_file = chk_run
erf.check_int = 10

means that restart files (really directories) starting with the prefix “chk_run” will be generated every 10 level-0 time steps. The directory names will be chk_run00000, chk_run00010, chk_run00020, etc.
erf.check_per = 5.0

means that restart files (really directories) starting with the prefix “chk_run” will be generated whenever the simulation time passes a multiple of 5.0. The directory names will reflect the integer number of steps which have elapsed.

To restart from chk_run00061,for example, then set

amr.restart = chk_run00061

Plotfiles¶

Controlling PlotFile Generation¶

“Plotfiles” can be written very efficiently in parallel in a native AMReX format or in HDF5. They can also be written in NetCDF.

The following options in the inputs file control the generation of plotfiles. Note that plotfiles can be written at two different frequencies; the names, frequency and content of the two streams are controlled separately.

List of Parameters¶

Parameter	Definition	Acceptable Values	Default
erf.plotfile_type	AMReX, NETCDF or HDF5	“amrex” or “netcdf / “NetCDF” or “hdf5” / “HDF5”	“amrex”
erf.plot_file_1	prefix for plotfiles at first freq.	String	“plt_1_”
erf.plot_file_2	prefix for plotfiles at seoncd freq.	String	“plt_2_”
erf.plot_int_1	how often (by level-0 time steps) to write plot files at first freq.	Integer \(> 0\)	-1
erf.plot_int_2	how often (by level-0 time steps) to write plot files at second freq.	Integer \(> 0\)	-1
erf.plot_per_1	how often in simulation time to write plot files at first freq.	Real \(> 0\)	-1.0
erf.plot_per_2	how often in simulation time to write plot files at second freq.	Real \(> 0\)	-1.0
erf.plot_vars_1	name of variables to include in plotfiles at first freq.	list of names	None
erf.plot_vars_2	name of variables to include in plotfiles at seoncd freq.	list of names	None

Notes¶

The NeTCDF option is only available if ERF has been built with USE_NETCDF enabled.

Examples of Usage¶

erf.plotfile_type = amrex
erf.plot_file_1 = plt_run
erf.plot_int_1 = 10

means that native plot files (actually directories) starting with the prefix “plt_run” will be generated every 10 level-0 time steps. If using amrex format, that directory names will be plt_run00000, plt_run00010, plt_run00020, etc. If using HDF5 format, the names will have “.h5” appended; if using NetCDF format, the names will have “.nc” appended.

In addition, while the amrex plotfiles will contain data at all of the refinement levels, NetCDF files are separated by level.

PlotFile Outputs¶

Plotfiles can include the quantities of several simulation parameters as output. They are summarized in the list below.

Output Options¶

Parameter	Definition
x_velocity	Velocity in x direction
y_velocity	Velocity in y direction
z_velocity	Velocity in z direction
density	Total density
dens_hse	Hydrostatic density
pert_dens	Perturbational density
pressure	Total pressure
pres_hse	Hydrostatic pressure
pert_pres	Perturbational pressure
pres_hse_x	Derivative of hydrostatic pressure in x
pres_hse_y	Derivative of hydrostatic pressure in y
dpdx	Pressure gradient in x direction
dpdy	Pressure gradient in y direction
temp	Temperature
theta	Potential temperature
rhotheta	Density * theta
KE	Kinetic energy
QKE	Turbulent kinetic energy * 2
rhoKE	Density * KE
rhoQKE	Density * QKE
scalar	Scalar magnitude
*vorticity_x	x-component of vorticity
*vorticity_y	y-component of vorticity
*vorticity_z	z-component of vorticity
rhoadv_0	Conserved scalar
soundspeed	Sound speed
z_phys	Terrain height
detJ	Jacobian determinant
mapfac	Map scale factor
Kmv	Vertical Eddy Diffusivity of Momentum
Kmh	Horizontal Eddy Diffusivity of Momentum
Khv	Vertical Eddy Diffusivity of Heat
Khh	Horizontal Eddy Diffusivity of Heat
qt	Nonprecipitating water (qv + qc + qi)
qp	Precipitating water (rain + snow + graupel)
qc	Cloud water mixing ratio
qi	Cloud ice mixing ratio
qv	Water vapor mixing ratio
rhoQt	Density * qt
rhoQp	Density * qp

Examples of Usage¶

In an input file, the user can select parameters to plot by supplying a space-delimited list to erf.plot_vars_1 or erf.plot_vars_2.

erf.plot_vars_1 = option1 option2 option3

Visualization¶

ERF currently generates plotfile in the native AMReX format.

There are several visualization tools that can be used for AMReX plotfiles, specifically ParaView, VisIt and yt.

ParaView¶

The open source visualization package ParaView v5.10 and later can be used to view ERF plotfiles with and without terrain. You can download the paraview executable at https://www.paraview.org/.

To open a plotfile

Run ParaView v5.10, then select “File” \(\rightarrow\) “Open”.
Navigate to your run directory, and select either a single plotfile or a set of plotfiles. Open multiple plotfile at once by selecting plt.. Paraview will load the plotfiles as a time series. ParaView will ask you about the file type – choose “AMReX/BoxLib Grid Reader”.
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 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. Change to the required number of AMR level before clicking “Apply”.
For “Representation” select “Surface”.
For “Coloring” select the variable you chose above.
To add planes, near the top left you will see a cube icon with a green plane slicing through it. If you hover your mouse over it, it will say “Slice”. Click that button.
You can play with the Plane Parameters to define a plane of data to view, as shown in 4.

_images/ParaView.png — : Plotfile image generated with ParaView¶

VisIt¶

AMReX data can also be visualized by VisIt, an open source visualization and analysis software. To follow along with this example, first build and run the first heat equation tutorial code.

Next, download VisIt from https://visit-dav.github.io/visit-website/ and install. To open a single plotfile, run VisIt, then select “File” \(\rightarrow\) “Open file …”, then select the Header file associated the the plotfile of interest (e.g., plt00000/Header). Assuming you ran the simulation in 2D, here are instructions for making a simple plot:

To view the data, select “Add” \(\rightarrow\) “Pseudocolor” \(\rightarrow\) “phi”, and then select “Draw”.
To view the grid structure (not particularly interesting yet, but when we add AMR it will be), select “Add” \(\rightarrow\) “Subset” \(\rightarrow\) “levels”. Then double-click the text “Subset - levels”, enable the “Wireframe” option, select “Apply”, select “Dismiss”, and then select “Draw”.
To save the image, select “File” \(\rightarrow\) “Set save options”, then customize the image format to your liking, then click “Save”.

Your image should look similar to the left side of 2.

: 2D (left) and 3D (right) images generated using VisIt.¶

In 3D, you must apply the “Operators” \(\rightarrow\) “Slicing” \(\rightarrow\) “ThreeSlice”, with the “ThreeSlice operator attribute” set to x=0.25, y=0.25, and z=0.25. You can left-click and drag over the image to rotate the image to generate something similar to right side of 2.

To make a movie, you must first create a text file named movie.visit with a list of the Header files for the individual frames. This can most easily be done using the command:

~/amrex/Tutorials/Basic/HeatEquation_EX1_C> ls -1 plt*/Header | tee movie.visit
plt00000/Header
plt01000/Header
plt02000/Header
plt03000/Header
plt04000/Header
plt05000/Header
plt06000/Header
plt07000/Header
plt08000/Header
plt09000/Header
plt10000/Header

The next step is to run VisIt, select “File” \(\rightarrow\) “Open file…”, then select movie.visit. Create an image to your liking and press the “play” button on the VCR-like control panel to preview all the frames. To save the movie, choose “File” \(\rightarrow\) “Save movie …”, and follow the on-screen instructions.

Caveat:

The Visit reader determines “Cycle” from the name of the plotfile (directory), specifically from the integer that follows the string “plt” in the plotfile name.

So … if you call it plt00100 or myplt00100 or this_is_my_plt00100 then it will correctly recognize and print Cycle: 100.

If you call it plt00100_old it will also correctly recognize and print Cycle: 100

But, if you do not have “plt” followed immediately by the number, e.g. you name it pltx00100, then VisIt will not be able to correctly recognize and print the value for “Cycle”. (It will still read and display the data itself.)

yt¶

yt, an open source Python package available at https://yt-project.org/, can be used for analyzing and visualizing mesh and particle data generated by AMReX codes. Some of the AMReX developers are also yt project members. Below we describe how to use on both a local workstation, as well as at the NERSC HPC facility for high-throughput visualization of large data sets.

Note - AMReX datasets require yt version 3.4 or greater.

Using on a local workstation¶

Running yt on a local system generally provides good interactivity, but limited performance. Consequently, this configuration is best when doing exploratory visualization (e.g., experimenting with camera angles, lighting, and color schemes) of small data sets.

To use yt on an AMReX plot file, first start a Jupyter notebook or an IPython kernel, and import the yt module:

In [1]: import yt

In [2]: print(yt.__version__)
3.4-dev

Next, load a plot file; in this example we use a plot file from the Nyx cosmology application:

In [3]: ds = yt.load("plt00401")
yt : [INFO     ] 2017-05-23 10:03:56,182 Parameters: current_time              = 0.00605694344696544
yt : [INFO     ] 2017-05-23 10:03:56,182 Parameters: domain_dimensions         = [128 128 128]
yt : [INFO     ] 2017-05-23 10:03:56,182 Parameters: domain_left_edge          = [ 0.  0.  0.]
yt : [INFO     ] 2017-05-23 10:03:56,183 Parameters: domain_right_edge         = [ 14.24501  14.24501  14.24501]

In [4]: ds.field_list
Out[4]:
[('DM', 'particle_mass'),
 ('DM', 'particle_position_x'),
 ('DM', 'particle_position_y'),
 ('DM', 'particle_position_z'),
 ('DM', 'particle_velocity_x'),
 ('DM', 'particle_velocity_y'),
 ('DM', 'particle_velocity_z'),
 ('all', 'particle_mass'),
 ('all', 'particle_position_x'),
 ('all', 'particle_position_y'),
 ('all', 'particle_position_z'),
 ('all', 'particle_velocity_x'),
 ('all', 'particle_velocity_y'),
 ('all', 'particle_velocity_z'),
 ('boxlib', 'density'),
 ('boxlib', 'particle_mass_density')]

From here one can make slice plots, 3-D volume renderings, etc. An example of the slice plot feature is shown below:

In [9]: slc = yt.SlicePlot(ds, "z", "density")
yt : [INFO     ] 2017-05-23 10:08:25,358 xlim = 0.000000 14.245010
yt : [INFO     ] 2017-05-23 10:08:25,358 ylim = 0.000000 14.245010
yt : [INFO     ] 2017-05-23 10:08:25,359 xlim = 0.000000 14.245010
yt : [INFO     ] 2017-05-23 10:08:25,359 ylim = 0.000000 14.245010

In [10]: slc.show()

In [11]: slc.save()
yt : [INFO     ] 2017-05-23 10:08:34,021 Saving plot plt00401_Slice_z_density.png
Out[11]: ['plt00401_Slice_z_density.png']

The resulting image is 5. One can also make volume renderings with ; an example is show below:

_images/yt_Nyx_density_slice.png — : Slice plot of \(128^3\) Nyx simulation using yt.¶

In [12]: sc = yt.create_scene(ds, field="density", lens_type="perspective")

In [13]: source = sc[0]

In [14]: source.tfh.set_bounds((1e8, 1e15))

In [15]: source.tfh.set_log(True)

In [16]: source.tfh.grey_opacity = True

In [17]: sc.show()
<Scene Object>:
Sources:
    source_00: <Volume Source>:YTRegion (plt00401): , center=[  1.09888770e+25   1.09888770e+25   1.09888770e+25] cm, left_edge=[ 0.  0.  0.] cm, right_edge=[  2.19777540e+25   2.19777540e+25   2.19777540e+25] cm transfer_function:None
Camera:
    <Camera Object>:
    position:[ 14.24501  14.24501  14.24501] code_length
    focus:[ 7.122505  7.122505  7.122505] code_length
    north_vector:[ 0.81649658 -0.40824829 -0.40824829]
    width:[ 21.367515  21.367515  21.367515] code_length
    light:None
    resolution:(512, 512)
Lens: <Lens Object>:
    lens_type:perspective
    viewpoint:[ 0.95423473  0.95423473  0.95423473] code_length

In [19]: sc.save()
yt : [INFO     ] 2017-05-23 10:15:07,825 Rendering scene (Can take a while).
yt : [INFO     ] 2017-05-23 10:15:07,825 Creating volume
yt : [INFO     ] 2017-05-23 10:15:07,996 Creating transfer function
yt : [INFO     ] 2017-05-23 10:15:07,997 Calculating data bounds. This may take a while.
Set the TransferFunctionHelper.bounds to avoid this.
yt : [INFO     ] 2017-05-23 10:15:16,471 Saving render plt00401_Render_density.png

The output of this is 6.

_images/yt_Nyx_density_vol_rend.png — Volume rendering of \(128^3\) Nyx simulation using yt. This corresponds to the same plot file used to generate the slice plot in 5.¶

Coupling To AMR-Wind¶

The simplest form of coupling is one-way, file-based coupling. With this approach, an ERF simulation is run and output files with relevant simulation data are stored at regular time intervals. Once the ERF simulation is complete, an AMR-Wind simulation begins, using the ERF output files as time-dependent boundary conditions or forcing.

1D File-based coupling¶

The 1D file-based coupling scheme is based on the Time-Height Profile Assimilation Technique proposed by Allaerts et al. In this approach, the time evolution of the vertical profile of the horizontal velocity components and potential temperature at a particular position from the mesoscale simulation is saved, and a forcing term is added to the microscale simulation to indirectly assimilate the mesoscale data at the microscale. The output files generated by ERF are NetCDF files of a format compatible with the AMR-Wind implementation of this form of coupling.

We note that the solution variables in ERF are

\[(\rho, \rho u, \rho v, \rho w, \rho \theta)\]

while those in AMR-Wind are

\[(u, v, w, T, p)\]

We convert the variables from conservative to primitive form prior to saving the output files.

The following parameters control the generation of 1D column output files in ERF. To use this capability, the user must install the NetCDF library and link to it at compile time, which can be done by setting USE_NETCDF = TRUE in the GNUMakefile for ERF.

Parameter	Definition	Acceptable Values	Default
erf.output_1d_column	whether to output for coupling to AMR-Wind	0 or 1	0
erf.column_file_name	prefix for output files	String	“column_data.nc”
erf.column_interval	how often (by level-0 time steps) to output data files	Integer \(> 0\)	-1
erf.column_per	how often (by simulation time) to write output data files	Real \(> 0\)	-1.0
erf.column_loc_x	x-coordinate where vertical profile will be extracted	prob_lo(0) <= x <= prob_hi(0)	0.0
erf.column_loc_y	y-coordinate where vertical profile will be extracted	prob_lo(1) <= y <= prob_hi(1)	0.0

You should specify either erf.output_int or erf.output_per, but not both.

2D File-based coupling¶

For 2D file-based coupling to AMR-Wind, vertical planes are saved from an ERF simulation to define inflow boundary conditions for AMR-Wind. These files are in native-AMReX format and contain the AMReX data structure known as a BndryRegister. These files are specifically meant to be read by AMR-Wind, which is also AMReX-based, or ERF itself.

To generate the boundary data files for each of the four planes, the following (or similar) should be added to the input file:

erf.output_bndry_planes = 1
erf.bndry_output_planes_interval = 2
erf.bndry_output_start_time = 0.0
erf.bndry_output_planes_file = "BndryFiles"
erf.bndry_output_var_names = temperature velocity density
erf.bndry_output_box_lo = 256. 256.
erf.bndry_output_box_hi = 768. 768.

The above inputs will output boundary planes of data at x=xlo, x=xhi, y=ylo and y=yhi, where xlo and ylo are defined by bndry_output_box_lo and xhi and yhi are defined by bndry_output_box_hi. Note that the selected lo and hi values must be interior to the first and last cell centers on the coarse grid (i.e., at least one half cell width from the boundary). Except if the domain is periodic in that direction, any location up to and including the domain boundary may be chosen. In this case the variables that are written are temperature, velocity and density, and they are written every 2 coarse time steps starting at bndry_output_start_time which is 0 in this case.

We also have the functionality in ERF to read in these types of files; for this one would add the following (or similar) line to the inputs file:

erf.input_bndry_planes = 1
erf.bndry_file = "BndryFiles"
erf.bndry_input_var_names = density temperature velocity

When run with these inputs, ERF will read in the time sequence of files contained in the folder BndryFiles, and perform time interpolation as necessary. The only assumption about the times associated with the files is that the start and end times of the current simulation lie in the time period covered by the files in BndryFiles. Within BndryFiles there is an ascii file time.dat which contains the (originating) timesteps and physical times associated with each of the files.

It is assumed at this point that the physical domain of the simulation reading the files is exactly the physical domain specified by bndry_output_box_lo and bndry_output_box_hi when the files were written. If not, ERF will abort with an error message.

We note that the boundary plane data will only be used on faces identified in the inputs file as inflow faces, i.e. if we specific inflow/outflow in the x-direction, and periodic in the y-direction, as below, then only the “xlo” boundary data from BndryFiles will actually be used.

geometry.is_periodic = 0 1 0

xlo.type = "Inflow"
xhi.type = "Outflow"

ERF vs WRF¶

The following comparison is based on the WRF Version 4 Technical Report, titled “A Description of the Advancd Research WRF Model Version 4”

Similarities¶

Equations: both ERF and WRF solve the fully-compressible, Eulerian nonhydrostatic equations, and conserve dry air mass and scalar mass. ERF does not have a hydrostatic option.

Prognostic Variables: velocity components (u,v,w); perturbation moist potential temperature. Optionally, turbulent kinetic energy and any number of scalars such as water vapor mixing ratio, rain/snow mixing ratio, cloud water / ice mixing ratio.

Horizontal grid: both ERF and WRF use Arakawa C-grid staggering.

Time Integration: Time-split integration using 3rd-order Runge-Kutta scheme with smaller time step for acoustic and gravity wave modes. Variable time step capability.

Spatial Discretization: 2nd- to 6th-order advection options in horizontal and vertical. In addition, several different WENO schemes are available for scalar variables other than density and potential temperature.

Turbulent Mixing: Sub-grid scale turbulence formulation. Vertically implicit acoustic step off-centering.

Diffusion: In WRF, the diffusion coefficients specified in the input file (\(K_h\) and \(K_v\) for horizontal and vertical diffusion) get divided by the Prandtl number for the potential temperature and the scalars. For the momentum, they are used as it is. In ERF, the coefficients specified in the inputs (\(\alpha_T\) and \(\alpha_C\)) are used as it is, and no division by Prandtl number is done.

Initial conditions: both ERF and WRF have the ability to initialize problems from 3D “real” data (output of real.exe), “ideal” data (output of ideal.exe) and from 1D input soundings.

Lateral boundary conditions: Periodic, open, symmetric and specified (in wrfbdy* files).

Bottom boundary conditions: Frictional or free-slip

Earth’s Rotation: Coriolis terms in ERF controlled by run-time input flag

Mapping to Sphere: ERF supports the use of map scale factors for isotropic projections (read in from wrfinput files).

Nesting: One-way or two-way. Multiple levels and integer ratios.

Key Differences¶

Vertical coordinates: Unlike WRF, ERF uses a terrain-following height-based vertical coordinate, with vertical grid stretching permitted.

Time Integration: ERF supports using a 3rd-order Runge-Kutta scheme with no substepping as alternative to RK3 with acoustic substepping.

Initial conditions: ERF has an additional mode of “custom” initialization in which the user writes the initialization routine.

ERF does not have the capability for global simulation

Verification¶

The following tests are used to verify the correct behavior of different algorithmic components of ERF.

Scalar Advection¶

Here we present spatial and temporal convergence studies for simple scalar advection with a uniform velocity field. The initial data has constant density and pressure, constant velocity \(u=10\) in the x-direction, and a scalar initialized with profile \(cos(\pi x)\) in a domain that is 2 units wide and periodic in the lateral directions with slip walls on top and bottom. The simulation is run for 10 periods, i.e. until time \(t=2.0\)

The first study, shown below on the left, tests the horizontal centered/upwind advection stencils for second through sixth order. The study on the right tests the WENO 3rd and 5th order stencils, with and without the smoothing contributions in the stencil. In all of these cases, the time step was held fixed at \(\Delta t = 0.0000078125\) to ensure that the spatial error dominates the temporal error.

Convergence studies of spatial error¶

Spatial convergence study (centered/upwind)	Spatial convergence study (WENO)

The second study tests the temporal accuracy by first setting \(\Delta t = 0.0005\) and \(\Delta x = 1/8\), then reducing both \(\Delta t\) and \(\Delta x\) by a factor of two, keeping the ratio of \(\Delta t\) to \(\Delta x\) constant. These tests were run with the 6th order accurate spatial stencil so that the temporal error dominated the spatial error. Here we recover the expected 3rd order accuracy of the RK3 scheme.

Convergence study of temporal error¶

Temporal convergence study

Nonlinear Density Current¶

The density current problem tests the effects of gravity and the behavior at a slip wall.

A detailed description of the problem and a comparison of solutions using a number of different codes can be found in the Straka 1993 paper

Potential temperature perturbation at 600s and 900s¶

Perturbational potential temperature at t = 600s	Perturbational potential temperature at t = 900s

Ekman Spiral¶

The Ekman spiral problem tests the computation of the stress term internally and at no-slip walls, as well as Coriolis and geostrophic forcing.

A description of the problem, including the exact solution, can be found at Ekman Spiral Description

The steady solution is shown below, as well as a log-log plot showing the error scaling as \(O(\Delta z^2)\).

Flow profile and Error¶

Flow profiles	Convergence study

Potential flow over a semi-cylinder¶

The potential flow over a semi-cylinder problem tests the terrain feature in two dimensions and the effectiveness of sponge zones in preventing spurious reflections. This is a classic text book problem which has an exact result. The density is constant and the streamwise velocity is 10 m/s, and sponge zones are used on both the streamwise boundaries as well as the top boundary. The bottom wall is an inviscid, slip wall. The schematic of the computational domain and the comparison of the steady state velocity profiles with the exact solution at two different horizontal and vertical locations are shown below. This case runs without any diffusion.

Potential flow over a hemisphere¶

The potential flow over a hemisphere problem tests the terrain feature in three dimensions and the effectiveness of sponge zones in preventing spurious reflections. This is a classic text book problem which has an exact solution. The density is constant and the streamwise velocity is 10 m/s, and sponge zones are used on both the streamwise boundaries as well as the top boundary. The bottom wall is an inviscid, slip wall. The schematic of the computational domain and the comparison of the steady state velocity profiles with the exact solution at two different horizontal locations are shown below. A small amount of physical diffusion was added to keep the simulation stable in this case. Hence the wall velocities in the plot below do not match exactly to the inviscid solution. But away from the walls, excellent quantitative agreement is observed.

Stokes second problem - viscous flow over an oscillating flat plate¶

The Stokes second problem - viscous flow over an oscillating flat plate in quiescent initial conditions, tests the viscous terms in the governing equations with terrain and the feature for specifying a custom terrain velocity. This is a classic text book problem which has an analytical solution. The bottom wall is specified as a flat terrain and all other boundaries are outflow. The contours of horizontal velocity and the comparison of the numerical and exact solutions at a given time are shown in the figure below.

Dry bubble and moist bubble rise simulations¶

Benchmark simulations of dry and moist bubble rises in Bryan and Fritsch are reproduced. The test case consists of a warm bubble rising in dry and moist conditions. The potential temperature perturbation and the vertical velocity are compared as shown in the figures below. Excellent quantitative comparisons are obtained. The dry and moist bubble cases are in Exec/RegTests/Bubble with the corresponding input files inputs_BF02_dry_bubble and inputs_BF02_moist_bubble.

GPU weak scaling¶

The plot shows weak scaling of the ABL application with the Smagorinsky LES model using A100 GPUs on the Perlmutter system at NERSC. The domain size is amr.n_cell = 128 128 512 for a single GPU; this is progressively scaled up to 2048 1024 512 for 128 GPUs. This test uses all 4 GPUs per node with GPU-aware MPI communication and runs 100 time steps.

The times shown includes both initialization and time evolution. All I/O and diagnostic calculations are turned off.

Regression Tests¶

There are currently 26 tests which are run as part of every PR. The CI tests use cmake and are based on the version of AMReX in the ERF submodule.

In addition there are nightly tests that use GNUMake and use the current development branch of AMReX.

Results from the nightly CPU tests can be found here: CPU tests

Results from the nightly GPU tests can be found here: GPU tests

The following problems are currently tested in the CI:

Test	nx ny nz	xbc	ybc	zbc	Ext	Other
Bubble_Density_Current	256 4 64	Symmetry Outflow	Periodic	SlipWall SlipWall	None	moist bubble
CouetteFlow	32 4 16	Periodic	Periodic	SlipWall SlipWall	None	inhomogeneous bc at zhi
DensityCurrent	256 4 64	Symmetry Outflow	Periodic	SlipWall SlipWall	None	+gravity
DensityCurrent_detJ2	256 4 64	Symmetry Outflow Outflow	Periodic	SlipWall SlipWall SlipWall	None	use_terrain = true uses zlevels detJ = 2 everywhere
DensityCurrent_detJ2_nosub	256 4 64	Symmetry Outflow	Periodic	SlipWall SlipWall	None	use_terrain = true uses zlevels detJ = 2 everywhere no substepping
DensityCurrent_detJ2_MT	256 4 64	Symmetry Outflow Outflow	Periodic	SlipWall SlipWall SlipWall	None	use_terrain = true uses zlevels detJ = 2 everywhere terrain_type = Moving
EkmanSpiral	4 4 400	Periodic	Periodic	NoSlipWall SlipWall	Geo	+Coriolis +gravity
IsentropicVortexAdvecting	48 48 4	Periodic	Periodic	SlipWall SlipWall	None
IsentropicVortexStationary	48 48 4	Periodic	Periodic	SlipWall SlipWall	None
MSF_NoSub_IsentropicVortexAdv	48 48 4	Periodic	Periodic	SlipWall SlipWall	None	tests map factors without substepping
MSF_Sub_IsentropicVortexAdv	48 48 4	Periodic	Periodic	SlipWall SlipWall	None	tests map factors with substepping
PoiseuilleFlow	32 4 16	Periodic	Periodic	NoSlipWall NoSlipWall	GradP in x
RayleighDamping	64 4 64	Periodic	Periodic	SlipWall SlipWall	None	Rayleigh damping
ScalarAdvectionUniformU	64 64 4	Periodic	Periodic	SlipWall SlipWall	None
ScalarAdvectionShearedU	64 4 64	Periodic	Periodic	SlipWall SlipWall	None
ScalarAdvDiff_order2	32 32 32	Periodic	Periodic	SlipWall SlipWall	None	advection + diffusion “Centered_2nd”
ScalarAdvDiff_order3	32 32 32	Periodic	Periodic	SlipWall SlipWall	None	advection + diffusion “Upwind_3rd”
ScalarAdvDiff_order4	32 32 32	Periodic	Periodic	SlipWall SlipWall	None	advection + diffusion “Centered_4th”
ScalarAdvDiff_order5	32 32 32	Periodic	Periodic	SlipWall SlipWall	None	advection + diffusion “Upwind_5th”
ScalarAdvDiff_order6	32 32 32	Periodic	Periodic	SlipWall SlipWall	None	advection + diffusion “Centered_6th”
ScalarDiffusionGaussian	16 16 16	Periodic	Periodic	SlipWall SlipWall	None
ScalarDiffusionSine	16 16 4	Periodic	Periodic	SlipWall SlipWall	None
TaylorGreenAdvecting	16 16 16	Periodic	Periodic	SlipWall SlipWall	None
TaylorGreenAdvectingDiffusing	16 16 16	Periodic	Periodic	SlipWall SlipWall	None

while the following tests are run nightly:

Test	nx ny nz	xbc	ybc	zbc	Ext	Other
ABL-Deardorff	64 64 64	Periodic	Periodic	NoSlipWall SlipWall	None	LES
ABL-Smag	64 64 64	Periodic	Periodic	NoSlipWall SlipWall	None	LES
CouetteFlow-x	32 4 16	Periodic	Periodic	NoSlipWall NoSlipWall	None	inhomogeneous bc at zhi
CouetteFlow-y	4 32 16	Periodic	Periodic	NoSlipWall NoSlipWall	None	inhomogeneous bc at zhi
DensityCurrent	256 4 64	Symmetry Outflow	Periodic	SlipWall SlipWall	None	+gravity
EkmanSpiral	4 4 400	Periodic	Periodic	NoSlipWall SlipWall	Geo	+Coriolis +gravity
EkmanSpiral_restart	4 4 400	Periodic	Periodic	NoSlipWall SlipWall	Geo	restart test
IsentropicVortexAdvecting	48 48 4	Periodic	Periodic	SlipWall SlipWall	None
IsentropicVortexStationary	48 48 4	Periodic	Periodic	SlipWall SlipWall	None
PoiseuilleFlow-x	32 4 16	Periodic	Periodic	NoSlipWall NoSlipWall	GradP in x
PoiseuilleFlow-y	4 32 16	Periodic	Periodic	NoSlipWall NoSlipWall	GradP in y
ScalarAdvecDiffDoubleDen	32 32 32	Periodic	Periodic	SlipWall SlipWall	None	Density = 2
ScalarAdvDiffInflowOutflow	32 32 32	Inflow Outflow	Periodic	SlipWall SlipWall	None
ScalarAdvecDiffUniformU	32 32 32	Periodic	Periodic	SlipWall SlipWall	None
ScalarAdvecUniformU	64 64 4	Periodic	Periodic	SlipWall SlipWall	None
ScalarAdvecShearedU	64 4 64	Periodic	Periodic	SlipWall SlipWall	None
ScalarAdvecUniformU	64 64 4	Periodic	Periodic	SlipWall SlipWall	None
ScalarDiffusionGaussian	64 64 64	Periodic	Periodic	SlipWall SlipWall	None
ScalarDiffusionSine	64 64 4	Periodic	Periodic	SlipWall SlipWall	None
TaylorGreenAdvecting	64 64 64	Periodic	Periodic	SlipWall SlipWall	None
TaylorGreenAdvDiffDoubleDen	64 64 64	Periodic	Periodic	SlipWall SlipWall	None	Density = 2

More details about the CI tests are given below.

Scalar Advection by Uniform Flow in XY Plane¶

This tests scalar advection with periodic boundaries in the lateral directions and slip walls at low and high z.

Test Location: Tests/test_files/ScalarAdvectionUniformU

Problem Location: Exec/RegTests/ScalarAdvDiff

X-Y slice of a 2-d cylindrical blob in a uniform velocity field (10,5,0)¶

Scalar concentration at t=0.	Scalar concentration at 20 steps.

Scalar Advection by Sheared Flow¶

This tests scalar advection with periodic boundaries in the lateral directions and slip walls at low and high z.

Test Location: Tests/test_files/ScalarAdvectionShearedU

Problem Location: Exec/RegTests/ScalarAdvDiff

X-Z slice of a 2-d cylindrical blob in a uniform shearing velocity field (8 log( (z+z0)/z0 ) / log ( (zref+z0)/z0 ) with z0 = 0.1 and zref = 80 in a triply periodic domain 8x8x8¶

Scalar concentration at t=0.	Scalar concentration at 80 steps

Scalar Diffusion: Sphere of Scalar¶

This tests scalar diffusion with periodic boundaries in the lateral directions and slip walls at low and high z.

Test Location: Tests/test_files/ScalarDiffusionGaussian

Problem Location: Exec/RegTests/ScalarAdvDiff

Diffusion of a spherical blob of scalar¶

Scalar concentration at t=0.	Scalar concentration at 20 steps (t = 0.01).

Scalar Diffusion: Sinusoidal Variation of Scalar¶

This tests scalar diffusion with periodic boundaries in the lateral directions and slip walls at low and high z.

Test Location: Tests/test_files/ScalarDiffusionSine

Problem Location: Exec/RegTests/ScalarAdvDiff

Diffusion of a scalar initialized as sin(x)¶

Scalar concentration at t=0.	Scalar concentration at 20 steps (t = 0.2).

Scalar Advection/Diffusion by Uniform Flow With Different Spatial Orders¶

This tests scalar advection and diffusion with periodic boundaries in the lateral directions and slip walls at low and high z.

Test Location (for 2nd order): Tests/test_files/ScalarAdvDiff_order2

Problem Location: Exec/RegTests/ScalarAdvDiff

Advection and diffusion of a spherical blob in a uniform velocity field (100,0,0)¶

Scalar concentration at t=0.	Scalar concentration at 20 steps (t = 0.01).

Rayleigh Damping¶

This tests Rayleigh damping. The problem is initialized as in the shear flow case, then Rayleigh damping is applied with a target mean profile of (2,1,0).

Test Location: Tests/test_files/RayleighDamping

Problem Location: Exec/RegTests/ScalarAdvDiff

Isentropic Vortex: Stationary¶

This tests advection of an isentropic vortex with triply periodic boundaries.

Test Location: Tests/test_files/IsentropicVortexStationary

Problem Location: Exec/RegTests/IsentropicVortex

Isentropic Vortex: Advecting¶

This tests advection of an isentropic vortex with triply periodic boundaries.

Test Location: Tests/test_files/IsentropicVortexAdvecting

Problem Location: Exec/RegTests/IsentropicVortex

Taylor Green Vortex: Advection¶

This tests advection and diffusion with triply periodic boundaries.

Test Location: Tests/test_files/TaylorGreenAdvecting

Problem Location: Exec/RegTests/TaylorGreenVortex

Taylor Green Vortex: Advection and Diffusion¶

This tests advection and diffusion with triply periodic boundaries.

Test Location: Tests/test_files/TaylorGreenAdvectingDiffusing

Problem Location: Exec/RegTests/TaylorGreenVortex

Scalar concentration¶

Flow field at t=0.	Flow field at 10 steps (t = 1.6).

Couette Flow¶

This tests Couette flow in a channel. The domain is periodic in the x- and y-directions, and has NoSlipWall bc’s on the low-z and high-z faces. At the high-z boundary the velocity is specified to be \(U = (2,0,0)\). The steady solution for this problem is \(U = (z/8,0,0)\) in the domain which is 16 units high in z.

Test Location: Tests/test_files/CouetteFlow

Problem Location: `Exec/RegTests/CouetteFlow`_

Poiseuille Flow¶

This tests Poiseuille flow in a channel. The domain is periodic in the x- and y-directions, and has NoSlipWall bc’s on the low-z and high-z faces. We initialize the solution with the steady parabolic profile \(U = (1-z^2,0,0)\) in the domain which runs from -1. to 1. in z. The viscosity is specified to be 0.1 and the imposed pressure gradient is \(Gp = (-0.2,0,0)\).

Test Location: Tests/test_files/PoiseuilleFlow

Problem Location: `Exec/RegTests/PoiseuilleFlow`_

Nonlinear Density Current¶

The density current problem tests the effects of gravity and the behavior at a slip wall.

See Verification for more information.

Test Location: Tests/test_files/DensityCurrent

Problem Location: `Exec/RegTests/DensityCurrent`_

Ekman Spiral¶

The Ekman spiral problem tests the computation of the stress term internally and at no-slip walls, as well as Coriolis and geostrophic forcing.

See Verification for more information.

Test Location: Tests/test_files/EkmanSpiral

Problem Location: Exec/RegTests/EkmanSpiral

Applications and Requirements¶

Last update: 2021-01-14

Overview¶

The Energy Research and Forecasting (ERF) model addresses a key component of the wind energy simulation supply chain, the downscaling of the energy contained within mesoscale atmospheric flows into the microscale wind plant environment, where turbines generate electricity. The efficient, robust and accurate downscaling of mesoscale flow information, which defines not only the energy available for conversion to electricity, but also many characteristics influencing the extraction of that energy, is crucial for the design and optimization of wind plant layout and operation in realistic, complicated operating conditions, including the more challenging complex terrain and offshore settings defining the majority of future deployment.

Much of the functionality ERF will provide is obtained today from weather prediction models such as the Weather Research and Forecasting (WRF) model. However, neither WRF nor similar models were designed to prioritize accurate prediction of low-level winds, or to be easily configurable for wind energy applications, focusing instead on larger-scale weather parameters such as storm dynamics, precipitation and temperature. While among such models WRF is most amenable to, and has been the most widely adopted into wind energy work flows, these and other limitations of its software architecture and design, including poor load balancing and inefficient data transfer in downscaling configurations, no ability to decompose in the vertical direction which limits parallelizability, and complications arising from its pressure-based vertical coordinate, render WRF inadequate as a long-term solution for wind energy applications.

Most importantly, WRF’s software architecture cannot effectively utilize Graphics Processing Units (GPUs), a requirement both for the increased computational expense of emerging applications, as well as for use of the next generation of high-performance computing (HPC) assets slated for rollout first at the DOE national laboratories, and eventually defining the majority of available scientific computing hardware. Moreover, no current code development efforts pursuing GPU-compatible architectures are addressing the downscaling functionality required of wind energy applications, with the other efforts pursuing either numerical weather prediction (NWP) on global-scale grids, or focusing on microscale simulation without including the critical mesoscale component defining the scales of energy input.

Beyond bridging the scale gap between mesoscale and microscale energy flow, ERF will provide functionality not contained within other mesoscale or microscale simulation codes, specifically the offshore and complex terrain environments that define key challenges for future wind energy expansion. To provide the required functionality, ERF will include advanced physical process models for mesoscale-to-microscale coupling, while seamlessly interfacing with external solvers, such as wave and sea-state models, and the ExaWind microscale wind plant simulation code. Coupling ERF with ExaWind is especially critical, as ExaWind relies on a larger-scale model such as ERF for energy inflow and boundary condition data, in order to provide accurate wind plant optimization guidance in complicated settings and conditions.

In short, ERF will provide a modern, flexible, and efficient GPU-capable software framework to supply critical atmospheric and environmental drivers of energy availability to the microscale wind plant environment, thereby enabling the significant advances in siting, design and operation required to support continued industry expansion into increasingly challenging environments and high penetration scenarios.

ERF User Base¶

The ERF model is being designed for use by moderately skilled to advanced practitioners of computational fluid dynamics (CFD) codes such as OpenFOAM and WRF, working within the wind energy industry, national laboratories and university research groups, in applications involving numerical simulations of atmospheric flows, and the interactions of those flows with wind turbines, wind plants, and multiple interacting plants in regions of dense development. Beyond its applications, an additional goal of ERF is to serve as a bridge for users and developers of existing CFD and NWP codes to transition away from older legacy codes into a modern software architecture and programming paradigm that efficiently utilizes the next generation of GPU-accelerated HPC hardware, while providing opportunities for expanded applicability to modern wind energy research challenges that require the implementation and integration of new computational capabilities. ERF also targets user-developers who seek to contribute new code back to the ERF code base to improve and expand it, as has occurred within other open-source codes such as WRF and OpenFOAM.

ERF Applications¶

ERF is being designed for wind energy applications requiring simulation of the atmospheric flows that provide the energy available for conversion to power within wind plants. ERF will focus on capturing at high fidelity the atmospheric flow quantities of primary interest for wind turbine and plant design, operation, and optimization, at multiple time and space scales, and in variable operating conditions. The primary applications and phenomena that ERF will address, and how ERF will improve upon the present state of the art, are described below.

1. Resource Characterization¶

A key application targeted by ERF is wind resource characterization, the assessment of the potential of a given site to produce power over time. Today’s commonly applied resource characterization approaches used in both mesoscale and microscale settings are necessarily of restricted fidelity due to the computational expense of higher-fidelity techniques exceeding industry resources in typical workflows. The computational burden of high fidelity will be significantly mitigated using ERF’s more efficient code architecture, algorithms, and ability to use GPU-accelerated HPC hardware, enabling improved resource characterization at all scales.

At the largest mesoscale applications envisioned for ERF (domains of hundreds to thousands of km on a side; grid spacings of one to ten km), increased throughput will enable several advantages. One is the ability to perform larger ensembles for more robust quantification of uncertainty. A second is the ability to use finer grid spacings over the same geographical footprint in areas requiring increased resolution of landscape or surface features. A third will be the ability to use higher-fidelity physical process models, including the three-dimensional planetary boundary layer (PBL) scheme and machine-learning-based surface layer flux models, each developed within the portfolio of projects supported by the Wind Energy Technologies Office (WETO). ERF will also support the implementation of more advanced offshore wave and sea-state models, including new capabilities developed under the WETO-supported Wind Forecast Improvement Project 3, which is focused on improving computational approaches for offshore environments. All of these new capabilities will enable much more accurate simulation and characterization of wind-energy-relevant flow phenomena within the atmospheric boundary layer (ABL) than is currently achievable within any existing code base.

At the smallest microscale applications envisioned for ERF (involving domains of a few to tens of km on a side, with grid spacings of a few to tens of m), higher computational throughput will enable improved characterization of the turbulence environment, as well as the impacts of smaller-scale terrain or surface heterogeneities on flow quantities of interest. More accurate treatment of complex terrain effects via immersed boundary methods, as well as the ability to integrate resolved wave and sea-state forcing for offshore applications will significantly enhance microscale resource assessment in these settings.

In addition to providing improved resource characterization in both mesoscale and microscale applications, ERF will enable much more efficient mesoscale-to-microscale coupling via efficient dynamic downscaling, interfacing the microscale turbulence field with the mesoscale forcing that drives it. Such multiscale coupling is especially critical in settings involving complex meteorology and landscape characteristics that supply forcing at scales larger than can be encompassed within even a very large single-domain microscale setup.

2. Forensics¶

Another application ERF will support is the ability to simulate unique meteorological events of importance, such as those leading to damage or some other outcome for which improved understanding is desired. ERF will enable enhanced forensics abilities via higher-resolution, and higher-fidelity treatment of relevant physical processes impacting the flow, coupled with the flexibility to either ingest larger-scale forcing datasets from forecast models or analysis products, or to set up idealized process-level simulations with controlled forcing.

3. Wind Plant Inflow¶

A primary use case for ERF is the downscaling of mesoscale atmospheric flows to microscale grid spacing, for which all of the relevant scales of motion, including turbulence, are sufficiently resolved to specify turbine-airflow interactions. While this information can be used to estimate resulting power, fatigue loading, and other data, ERF is also being designed to couple directly with the ExaWind microscale wind plant simulation code, within which turbine performance, loading, and controls models of various fidelity can operate directly within the ERF-generated inflow. These coupled ERF-ExaWind simulations will provide unprecedented levels of full-spectrum fidelity, information required to understand and optimize wind plant performance in general, complicated operating conditions and environments.

An additional aspect of wind plant inflow is the impact of entire wind plants on both their own inflow, via blockage effects, as well as on downstream plants via wind plant wakes, gravity waves, or other atmospheric disturbances that large wind plants generate. These issues are of particular importance in areas of dense development, and require a larger simulation footprint than is practical within even a very large single-resolution microscale domain. In addition, difficulties simulating gravity waves using the incompressible solvers that form the basis of many microscale wind plant simulation codes are ameliorated using a fully compressible solver such as ERF will employ. The multiple-resolution capability of ERF, coupled with the incorporation of wind plant wake models applicable at both mesoscale and microscale grid spacings, will provide a flexible framework to better understand wind plant interactions, regional wind power generation, and the regional integration of wind-generated power.

4. Offshore Development¶

A defining challenge of future wind energy development is the offshore environment, which presents unique operating conditions that require the creation of suitable simulation tools to understand and design for. Among the unique offshore conditions impacting wind energy are swell, wave and sea-state variability that impact the low-level atmospheric flow, hence the available energy. Sea surface temperature and roughness variability can also influence submesoscale motions that are important in offshore environments. The large thermal inertia of water can also support persistent static stability regimes with strong and long-lived impacts on flow and wake propagation, for example, via synoptic-scale advection of air masses with different thermal properties over the water, and due to variability in sea-surface temperature due to the existence of currents or bathymetric influences.

Improved parameterizations to represent these unique features of offshore environments at various scales will be implemented into ERF, along with abilities to explicitly specify wave characteristics in large-eddy simulation (LES) domains containing sufficient mesh resolution to capture the wave shape and impacts of moving wave surfaces on the flow. IBMs may provide a pathway to efficiently implement resolved wave impacts into LES domains as well.

5. Impacts of Complex Terrain¶

ERF will be designed to improve the representation of complex terrain and its impacts on the flow, including gravity flows, gravity waves, mountain-valley circulations, and coastal jets, in mesoscale simulations, relative to other mesoscale models. Improved complex terrain capabilities will be incorporated via the use of higher-order numerical stencils for the evaluation of horizontal derivatives over moderately steep terrain, and immersed boundary methods (IBMs) for very steep terrain. These approaches will reduce numerical errors while extending ERF to much steeper slopes than the standard WRF model and similar codes can simulate. IBMs will also permit use of higher resolution, less smoothed terrain than with WRF, which will improve simulation fidelity in complex terrain, and thus improve the local wind accuracy around turbines and plants. IBMs can also stabilize numerical solutions over steep terrain, even if not strictly required, allowing for larger model time steps and therefore accelerating execution. ERF’s use of a vertical coordinate with a fixed height will lead to much more efficient use of IBMs in ERF than in WRF, where the changing heights require new interpolations and projections at every time step.

6. Impacts of Low-Level Jets¶

An important meteorological feature defining the energy resources in many geographic locations, including the US central great plains and offshore regions, is the low-level jet (LLJ), a narrow ribbon of fast moving air that occurs within the lowest several hundred meters above the surface. While LLJs provide a rich energy resource, LLJs characteristically contain strong sheer, veer and intermittent atmospheric wave and turbulence activity, all of which can increase fatigue loading. Moreover, details of their height and strength, as well as the timing of their onset and dissolution, which impact the integration of power produced, present numerous challenges to development within such regions. More efficient downscaling and higher-fidelity mesoscale and microscale turbulence models will provide enhanced understanding of LLJ impacts on wind power applications.

7. Impact of Clouds and Precipitation¶

Clouds are important modulators of the atmospheric flow, impacting turbulence intensity via shading of the surface, which also influences boundary layer growth and the vertical transfer of momentum, leading to changes in mean wind speed, shear and veer across the turbine rotor swept area. In the offshore environment, radiative cooling from liquid water in the stratocumulus field that often surmounts the marine ABL can drive increased turbulence within the flow below the cloud deck. Clouds also produce various species of precipitation which impact turbine and plant performance, including glaze and rime ice that reduce aerodynamic performance, raindrops that can accelerate leading edge erosion, and graupel and hail which can be particularly damaging. The greater computational expense of cloud microphysics models that can accurately depict these processes, as well as running these schemes at the fine mesh scales required, have hindered understanding and predictive capabilities for the impacts of clouds and precipitation on wind power generation thus far. ERF’s more efficient solution framework will facilitate addressing these gaps. Higher-fidelity representation of cloud impacts on surface shading will also improve ERF’s applicability to solar and hybrid plant operation.

8. Impacts of Thunderstorms and Tropical Disturbances¶

Thunderstorms and tropical disturbances are multiscale phenomena with potentially high impacts on wind plant operation and reliability due to strong winds, high turbulence levels and thus gustiness, large raindrops and hail, and, in offshore environments, associated surface waves that can cause significant damage to turbine platforms. Current atmospheric simulation codes such as WRF provide some capability to investigate storm impacts on wind plants, but simulations are too expensive for routine application. ERF will provide a superior framework for investigating the impacts of strong storms on both onshore and especially offshore environments, where integration with high-fidelity wave and sea-state models will significantly extend predictive capabilities for ABL flow and potentially relevant sea-state characteristics.

9. Regional Integration and Hybrid Plants¶

The integration of wind-generated electricity at high penetration scenarios is critical to the robust and efficient operation of the grid, and the ability to safely reduce baseload generation from fossil fuel-based plants. Integration relies on accurate predictions of both power generation and demand, each of which is impacted by atmospheric parameters computed by ERF. With a regional mesoscale footprint, ERF will be able to predict the regional distributions of atmospheric quantities impacting both production and demand. ERF can also improve solar forecasting, for many of the reasons it enhances wind predictions (greater throughput, more efficient downscaling, and higher-fidelity physics parameterizations), thereby facilitating integration of wind with another rapidly developing weather-dependent generation source. Moreover, the integration of solar and wind with data analytics from the grid will facilitate exploration of integration with energy storage, transmission, and grid services. As such, ERF will facilitate the development and operation of hybrid plants which balance multiple sources of production, storage and grid services, helping to transform meteorologically-dependent, intermittent sources of energy into robust, dispatchable power generation and delivery.

10. Development of Machine Learning and Artificial Intelligence Methods¶

Machine learning (ML) and artificial intelligence (AI) approaches represent significant opportunities to improve understanding of physical phenomena, enhance the fidelity of numerical simulations, improve the computational efficiency of such simulation codes, and to develop faster running reduced-order models for applications requiring increased throughput in industry workflows. ERF will address development of ML and AI methods by providing the ability to generate simulation-based training datasets via its integration of high-fidelity physical process parameterizations with efficient multiscale simulation. Coupling of ERF with the ExaWind code will provide unique datasets enabling the examination of relationships between mesoscale flow features and turbine-level response that can bypass the most expensive portion of full-physics numerical simulations, the downscaling of flow into and simulation within the microscale domain.

ERF’s Role within the Spectrum of Geophysical and Wind Energy Applications¶

These above listed activities and phenomena define the key applications envisioned for the ERF model. However, there are two other flow simulation regimes of relevance to wind-energy that ERF is not intended to address. The first of these is weather forecasting. While ERF could, in principle, be extended to capture larger-scales of meteorological forcing, efforts are already underway elsewhere to create next-generation numerical weather NWP systems, operating at global scales, to capture the largest scales impacting weather system evolution, while also being designed to utilize GPUs for enhanced speed and efficiency. ERF will leverage these concurrent developments by interfacing with a data preprocessor to ingest forecast and analysis fields produced by these new larger-scale models. ERF will focus on the efficient downscaling of those solutions to footprints of relevance to wind energy applications, capturing the associated finer-scale mesoscale and turbulence features, as well as wind plant interactions, along the way.

The second application ERF will not address is very fine-scale microscale simulation and interaction with resolved turbine components. As is the case with NWP, other modeling tools, including the ExaWind code, are being developed to support those activities, with ERF functioning primarily as the provider of inflow and boundary information, downscaled from mesoscale or NWP scales, to those fine-scale application domains, through robust model coupling interfaces. ERF will also be designed to upscale information from finer-scale offline-coupled simulations back into its domain(s) for improved fidelity.

ERF Features and Requirements¶

Below is a list of the features that the ERF code must provide, and requirements of the code to support those features, in order to satisfy the above described user base and applications.

1. Excellent Performance on Both CPU- and GPU-Based HPC Platforms¶

ERF must be able to run efficiently on both CPU- and GPU-based HPC platforms. This flexibility is required to support enhanced utilization of existing HPC architectures for which significant industry investments have been made, to serve as a vehicle to transition those users to next-generation platforms and programming paradigms, and to support current applications using GPU-accelerated hardware being rolled out today at leadership computing facilities (LCFs). To meet these use cases, ERF must compile and run on a variety of platforms, but also must be configurable for optimal performance on LCFs, including coupling with ExaWind tools on those platforms. Key metrics to assess adequate performance include superior scaling up to tens of thousands of cores on LCF systems, with several levels of mesh refinement, and solution accuracy that meets or exceeds that of legacy codes such as WRF and OpenFoam in similar applications. In addition to LCF machines at DOE labs, integration with new disk storage approaches (e.g., burst buffer at NERSC) should also be explored. Other platforms that would be desirable for ERF to utilize include emerging GPU-based small sized clusters, commodity desktops and laptops with GPU cards, and cloud resources, which are increasingly coming to replace industry-owned HPC resources at many wind energy companies.

2. C++ Base Code with Ability to Incorporate FORTRAN¶

For optimal utilization of GPU-based hardware, the ERF source code must be written in C++. However, ERF should also be able to incorporate legacy Fortran source code from other models. While Fortran code incorporated into the C++ code base will not result in optimal performance, it will allow for the rapid expansion of ERF’s capabilities, while providing a pathway to facilitate adoption of ERF by users and developers familiar with Fortran programming and legacy codes. Future development of ERF, including potential community development, can target the rewriting of desired Fortran modules into C++ for enhanced performance.

3. Configurability for Mesoscale, Microscale and Multiscale Simulations¶

ERF must provide flexible configurability that supports mesoscale or microscale simulations, each using a single mesh level, and for multiscale simulations, with as many levels of mesh refinement as required to span a given scale range, from horizontal grid spacings as large as several kilometers to as fine as several meters. It would be preferable to be able to use integer mesh refinement ratios from two to approximately ten or so in order to support the ability to avoid certain grid spacings depending on the application. The code must also support vertical refinement, however due to the nature of the vertical coordinate and required applications, vertical refinement should enable arbitrary height levels to be specified on all meshes.

Atmospheric downscaling applications will focus primarily on increased resolution over particular geographical areas, rather than tracking flow features that move in time, and therefore static refinement over rectangular volumes is sufficient for the near term (although eventually adaptive refinement will be useful to track features such as storms or turbine wakes). One-way coupling at the refinement interfaces, for which the fine mesh only receives information from, but does not transmit information back to the parent mesh, is acceptable at the beginning. Two-way nesting, for which the fine mesh provides information back to the parent mesh, is eventually needed to improve simulation also on those coarser domains.

The downscaling information exchange at mesh interfaces can follow procedures used in other codes, such as WRF, in which the parent mesh executes one time step, after which variables at the beginning and end of each parent time step are interpolated linearly in time to the refined-mesh time step to support integration of the finer mesh solution. For two-way coupling, the fine-mesh solution can be averaged to the coarse mesh after advancing to the coarse-mesh time step, where it can either replace the solution on that mesh, or provide a target for relaxation of the coarser-mesh solution.

4. Initial and Boundary Condition Preprocessing for Real-Data Simulations¶

To facilitate adoption of ERF for its primary intended workflows, those involving mesoscale energy simulation and its downscaling, ERF must support straightforward methods to prepare initial conditions and forcing for its lateral and surface boundaries. This process requires interpolation in time and space from the locations of state variables within the native datasets, to the locations of ERF’s mesh (or meshes in downscaling setups). The input data must also be transferred onto the ERF model’s map projection. Adopting some of the functionality of the WRF Preprocessing System (WPS) would provide an excellent pathway, as WPS has the ability to read in data from multiple larger-scale forecast and analysis products, and to project those data onto a number of standard map projections. As the WRF model supports global simulation, as well as configurability over arbitrary locations, including equatorial and polar latitudes, it contains several map projections. With ERF’s focus on mid-latitude locations, a Lambert Conformal projection would be ideal at the initial development phase. ERF can also initially focus on global analysis and forecast models such as NCEP’s Global Forecast System (GFS), as well as for North American applications, NOAA’s High Resolution Rapid Refresh (HRRR) model.

5. Initial and Boundary Condition Support for Idealized Simulations¶

To serve a variety of process-level applications and interaction studies, ERF must provide support for idealized setups, easily configurable into arbitrary (rectangular) domains with simplified, user-defined specification of input (vertical distributions of state variables), forcing (horizontal pressure gradient, geostrophic wind, and wind profile assimilation methods), and support for idealized boundary conditions, including open and periodic lateral boundaries, and appropriate top and surface boundary conditions such as free slip, no flux and no normal flow. Radiative fluxes at the domain top should also be easy to specify for use with low model tops, which is commonly done to reduce cost in high-resolution simulations under applicable atmospheric conditions.

6. Compressible Nonhydrostatic Dynamic Core¶

To simulate atmospheric forcing arising from mesoscale processes, a fully compressible, non-hydrostatic equation set appropriate for dynamics involving vertical density variability and well-resolved vertical motions must be used. While the WRF model uses a pressure-based vertical coordinate system for which the height above the surface of model grid points change in time, ERF should adopt a system for which the heights remain constant over time, simplifying applicability to wind energy use cases, including the coupling with other codes such as ExaWind. The equation set and vertical coordinate used by the COSMO (Consortium for Small Scale Modeling) model would serve as an excellent template for ERF.

7. Second-Order Finite Difference Spatial Discretization¶

The ERF equation set will require a discretization strategy and numerical solution procedures amenable to optimization for different mesh spacings and applications. For ease of implementation and familiarity with users of other code bases, as well as ability to incorporate modules from other codes such as WRF, a finite-difference spatial discretization strategy with second-order accuracy should be employed for ERF. Options for higher-order spatial differences can be included as well, however those methods may not scale as well on GPU-based hardware.

8. Fully Explicit and Mixed Implicit-Explicit Time Discretizations¶

For time advancement, a fully explicit method is required for general applications, however the flexibility to implement options for implicit treatment in the vertical, or all three directions, should be included. Different meshes must be able to use different time steps (selectable by the user based upon the spatial refinement ratio). Substepping in time for acoustic mode propagation should be included for the advancement of the pressure or density field in mesoscale domains where a time step sufficient to resolve those modes would be untenable.

9. Application-Relevant Physical Process Parameterizations¶

To support mesoscale, microscale and downscaling simulations, ERF must contain physical process parameterizations appropriate to all relevant scales and processes, including:

Monin-Obukhov Similarity Theory (MOST) surface stress boundary condition
Wave damping treatments for the upper domain
Subgrid turbulence closure for large-eddy simulation
Subgrid turbulence closure for mesoscale simulation
Surface energy budget parameterization for surface energy and momentum fluxes
Vegetation canopy model for improved flow over tall vegetation
Shortwave and longwave radiative transfer to capture solar/diurnal forcing and cloud-induced radiation impacts
Cloud microphysics to capture cloud-radiative forcing and precipitation
Offshore wave, sea-state and ocean current models appropriate for various spatial scales
Wind plant wake models for microscale and mesoscale applications

Following the WRF model, the physical process parameterizations should be callable on user-adjustable time steps (for faster simulation execution), and if a Runge-Kutta time advancement scheme is chosen, computed on the first predictor step of the sequence.

Hindcasting and forensics applications will require methods to incorporate forcing data from either analysis datasets (“analysis” nudging) or from observations (“observation” nudging), using Newtonian relaxation, or spectral methods.

Many of ERF’s required parameterizations and other capabilities can be taken from existing models such as WRF and implemented into the ERF source code. Another pathway is to interface ERF with the community physics package (CPP), a repository of common physics modules developed for integration with various community models such as WRF, MPAS and FV3. As described previously, providing an ability to compile modules written in Fortran into the ERF executable will facilitate the incorporation of existing codes from other models, and from the CPP, as well as facilitating the transition of users of older legacy codes and modules to ERF.

10. Robust and Efficient Interfacing of ERF with Other Code Bases¶

ERF must robustly and efficiently interface with other codes, such as ExaWind, modules within the CPP, wave, sea-state, and wind turbine aerodynamics and load models to support ERF’s applications. The interfacing must allow for the efficient exchange of required variable values or forcing terms at the domain boundaries or overlap regions of the coupled codes.

11. Efficient and Flexible Selection of Output and Internal Diagnostics¶

ERF should incorporate efficient parallel input/output (I/O) strategies to support runs over large processor counts, and frequent output, including options for asynchronous I/O to overlap I/O and computation. For variables on staggered grids, projection to cell-centered locations is desirable. For variables decomposed into perturbation and base states, reconstruction of the entire physical field before outputting is desirable to reduce file sizes.

To further reduce output file sizes and time spent writing output fields, ERF must provide abilities to select output for specific variables over domain subvolumes and slices of arbitrary size and orientation, and to compute diagnostic quantities such as spectra and Reynolds stresses, either over time or arbitrary spatial directions and volumes, during code execution.

ERF should adopt NetCDF climate and forecast (CF)-compliant metadata or a sufficiently close proximity thereof for output files to enable the use of common plotting and analysis software like xarray in Python, Ferret, NCL, etc.

12. Flexibility and Extensibility¶

For the ERF code to remain viable on evolving computational hardware while maintaining abilities to address evolving wind energy applications, the code architecture must provide high levels of flexibility and extensibility. These attributes must encompass data and memory flow and management, as required to utilize future HPC hardware, ability to incorporate more efficient and accurate numerical solution methods and physical process models, and ability to interface with other codes, including not only those related to geophysical process and wind energy design, but other related applications such as grid, transmission, storage and others.

13. Detailed Documentation and Test Cases¶

To attract new users and encourage development by the user community, ERF must contain detailed documentation and test cases describing configurability and pathways for extensibility. The documentation should take the form of a users’ guide for a high-level understanding of how to use the code for various applications, a detailed technical volume describing the equations, discretization strategy, solution methods, and data management, and extensive comments within the source code that describe the function of specific segments. The WRF model provides an excellent example of a documentation and code structure that facilitates adoption by new users, while encouraging community contributed extensions to its capabilities.

14. Balance between Performance, Robustness and Usability¶

ERF development must strike an appropriate balance between achievement of optimal performance on priority HPC platforms and other considerations, including robustness of the code to faithfully compile and run on different platforms, tractability for use by a diverse set of practitioners in a diverse set of applications, and maintainability to ensure a manageable workload for continued utility.

ERF Code Design¶

Development of the ERF source code requires a detailed strategy for the implementation of specific components and capabilities, a vision for overall code topology to provide the required flexibility and extensibility, protocols for code development consistency and documentation, and robust testing to establish code performance and ensure continuity of performance and functionality under continuing development.

1. Development of Appropriate Code Development Protocols¶

To ensure a robust record of code development history, ERF will be developed and managed within a GIT commit structure, with protocols for description and verification of implementations following the ExaWind project.

Documentation will be required on three levels:

Unit-level information within the source code to aid future users and developers
Higher-level design and implementation information in accompanying documentation on a readthedocs.org portal
A users’ guide detailing code functionality and describing accompanying test cases to demonstrate that functionality to new users; to learn by going through examples.
Shared analysis scripts to ensure consistency of results across the team

Development of the code will follow modern software project paradigms, including

Configurability of software packages

Different combinations of components constituting different applications

E.g., single- versus multi-resolution, real data versus idealized, microscale only versus mesoscale or multiscale (requiring physics packages), stand-alone versus coupled with other codes, …

Assessment of compile-time versus run-time configurability to guide optimization

Encapsulation
- Related functionality and data that can be grouped as a single class are organized into encapsulated code units that can have multiple alternative implementations.
- Code libraries to provide services such as discretization, data management, and orchestration of data movement for parallelization as well as I/O.
- Physics solvers are largely ignorant of the details managed by the framework (dynamic core and I/O), and operate in “plug-and-play” mode within the framework to enable streamlining for specific applications, or to perform ensembles for which configuration, dynamics and physics options can be modified at compile or run time
Separation of concerns so that, e.g., infrastructure and physics solvers development do not intrude into each other’s space
Hierarchy of parallelism within and across code units

2. Development of Robust Testing Strategies¶

Robust testing at multiple levels is required for verification, validation, and characterization of code performance and simulation accuracy.

At the smallest unit level, test cases for unit-level commits will be provided by the originators of those code units and included in a master suite of unit tests which must exhibit acceptable performance during future code development.
At intermediate aggregation levels for which functional units group together to provide a moderately complex capability, those sub-model aggregations must exhibit acceptable performance against suitable test cases, under future code development.
At the highest aggregation levels for which functional groups combine together to provide a complex capability, those whole-model aggregations must exhibit acceptable performance against suitable aggregate-level test cases under future development.

Acceptable performance should consist of bit-for-bit agreement for the addition of unrelated code components, or components that operate on data management and movements, but are not expected to alter values.

Higher-level tests should consist of performance against standards, such as analytical or manufactured solutions, as well as convergence under spatial or temporal refinement, depending on the sub-model component. The highest-level or “whole-model” tests should consist of evaluation against data from field campaigns or previous appropriate whole-model results, such as obtained from WRF or OpenFoam.

These standards should be articulated at a sufficient level of detail to guide community contributors of requirements for code additions or improvements.

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.

1. Governing Equation Set¶

Details of ERF’s governing equation set are currently still being formulated, however it will closely follow the formulation used by the COSMO model, with the following attributes:

Fully compressible non-hydrostatic formulation
Terrain-following computational mesh
Fixed-height vertical coordinate
Prognostic variables (three dimensional winds, a conserved temperature variable (potential or moist potential temperature, and pressure) cast in perturbation form, relative to hydrostatic base state.

This formulation does not conserve mass, but errors are negligible over intended timescales of simulations (hours to days).

2. Spatial Discretization¶

Spatial discretization of the ERF governing equation set will follow a second-order finite difference strategy, on an Arakawa C-grid. Horizontal grid spacing will be equal in each direction and constant within each mesh refinement level with respect to adjustments necessary for mapping factors. The vertical coordinate will support user specification of heights on each refinement level. Higher-order difference formulations, including upwinding schemes, will also be provided if needed and resources allow. A Lambert conformal map projection will be used.

Immersed boundary methods will also be implemented for a variety of applications, including reduced numerical errors over complex terrain, extension to very steep terrain up to and including vertical cliff walls, incorporation of embedded turbine and structure geometries, and potentially as a method to incorporate resolved wave impacts on atmospheric flows. Various methods developed within the WRF model can be straightforwardly extracted from that code and implemented into ERF, within which the ability to maintain constant heights of the model vertical coordinate will greatly enhance code performance when in use.

3. Temporal Discretization¶

The base temporal discretization strategy for ERF will be an explicit third-order Runge-Kutta predictor-corrector method, similar to that used in the WRF code and others. The scheme will permit users to specify time steps for each mesh refinement level. For improved performance, smaller time steps for acoustic modes within each Runge-Kutta time step will be enabled for applications with grid spacings large enough to improve performance. Implicit methods will also be investigated for further performance enhancements. Implicit methods for the vertical direction only will be explored for cases in which the vertical grid spacing is small relative to the horizontal, requiring a prohibitively small model time step, as well as to stabilize the code with a larger time step. Advantages of also using implicit methods in the horizontal direction relative to the increased memory requirements of such formulations will also be explored.

4. Numerical Solution Methods¶

ERF should be formulated to support integration of a variety of numerical solution methods for different applications, and to incorporate future techniques that can provide superior performance or other advantages. The team will prioritize use of solvers available within the PETSc libraries that have already been incorporated into the mesh refinement framework that is also being used for ERF.

6. Boundary Conditions¶

List of top, bottom and lateral boundary condition implementations provided here.

The highest priority needed for initial code testing and to simplify additional development includes capability to support ABL flow, including free-slip boundaries at the top and bottom, with no flow normal to the domain top and surface, and no fluxes of energy or constituents, periodic lateral boundaries in each direction, and wave damping at the model top.

Once the basic code is functional, the ability to read time dependent boundary conditions from input files using modified WPS code will be implemented to enable real-world simulations.

7. Physical Process Parameterizations¶

List of physical process models described here. Priority development includes:

MOST surface stress boundary condition for surface momentum fluxes
Subgrid turbulence flux model for large-eddy simulation
Subgrid turbulence flux model (PBL scheme) for mesoscale simulation
Hooks to the CPP for radiation, cloud, surface, and boundary layer parameterizations

8. Model Coupling¶

Explore pathways to interface both ERF and ExaWind within the same AMReX platform.

9. Input/Output¶

Short-term goal: Utilize native AMReX output which can be read by ParaView and VisIt. Add capability to write to NetCDF in a mostly-CF-compliant format, enabling integration with Xarray in Python, and potentially other plotting utilities.

Long-term goal: Implement capability to configure I/O at run time based, for example, on a YAML file that is separate from the configuration file used to run the code. The separate I/O file would support:

Adding and removing specific variables from the output
Changing output frequency for different variables
Use of multiple output files with different sets of variables and output frequencies