MSGWaM

MSGWaM

3D transient implementation of MSGWaM.

See also

• PinCFlow.MSGWaM.BoundaryRays

• PinCFlow.MSGWaM.RayUpdate

• PinCFlow.MSGWaM.MeanFlowEffect

External links

1. Muraschko et al. (2014)

2. Boeloeni et al. (2016)

3. Wilhelm et al. (2018)

4. Wei et al. (2019)

5. Jochum et al. (2025)

Interpolation

Module for interpolating mean-flow quantities to ray-volume positions.

Provides functions that find the grid points closest to a given ray-volume position and perform trilinear interpolation of mean-flow quantities.

See also

• PinCFlow.Types

• PinCFlow.Update

DChiDX

Singleton for dispatch to interpolation of $\partial \chi_\mathrm{b} / \partial x$.

DChiDY

Singleton for dispatch to interpolation of $\partial \chi_\mathrm{b} / \partial y$.

DChiDZ

Singleton for dispatch to interpolation of $\partial \chi_\mathrm{b} / \partial z$.

DN2DZ

Singleton for dispatch to interpolation of $\partial N^2 / \partial z$.

DUDX

Singleton for dispatch to interpolation of $\partial u_\mathrm{b} / \partial x$.

DUDY

Singleton for dispatch to interpolation of $\partial u_\mathrm{b} / \partial y$.

DUDZ

Singleton for dispatch to interpolation of $\partial u_\mathrm{b} / \partial z$.

DVDX

Singleton for dispatch to interpolation of $\partial v_\mathrm{b} / \partial x$.

DVDY

Singleton for dispatch to interpolation of $\partial v_\mathrm{b} / \partial y$.

DVDZ

Singleton for dispatch to interpolation of $\partial v_\mathrm{b} / \partial z$.

N2

Singleton for dispatch to interpolation of $N^2$.

compute_derivatives(
    state::State,
    i::Integer,
    j::Integer,
    kd::Integer,
    ku::Integer,
    phitype::DUDX,
)::NTuple{2, <:AbstractFloat}

Compute and return the zonal derivative of the zonal wind ($\partial u_\mathrm{b} / \partial x$) at $\left(i, j, k_\mathrm{D}\right)$ and $\left(i, j, k_\mathrm{U}\right)$.

The derivative is given by

\[\left(\frac{\partial u_\mathrm{b}}{\partial x}\right) = \frac{u_{\mathrm{b}, i + 1 / 2} - u_{\mathrm{b}, i - 1 / 2}}{\Delta \widehat{x}} + G^{13} \frac{u_{\mathrm{b}, k + 1} - u_{\mathrm{b}, k - 1}}{2 \Delta \widehat{z}}.\]

compute_derivatives(
    state::State,
    i::Integer,
    j::Integer,
    kd::Integer,
    ku::Integer,
    phitype::DUDY,
)::NTuple{2, <:AbstractFloat}

Compute and return the meridional derivative of the zonal wind ($\partial u_\mathrm{b} / \partial y$) at $\left(i + 1 / 2, j + 1 / 2, k_\mathrm{D}\right)$ and $\left(i + 1 / 2, j + 1 / 2, k_\mathrm{U}\right)$.

The derivative is given by

\[\begin{align*} \left(\frac{\partial u_\mathrm{b}}{\partial y}\right)_{i + 1 / 2, j + 1 / 2} & = \frac{u_{\mathrm{b}, i + 1 / 2, j + 1} - u_{\mathrm{b}, i + 1 / 2}}{\Delta \widehat{y}} + G_{i + 1 / 2, j + 1 / 2}^{23} \frac{u_{\mathrm{b}, i + 1 / 2, j + 1 / 2, k + 1} - u_{\mathrm{b}, i + 1 / 2, j + 1 / 2, k - 1}}{2 \Delta \widehat{z}}. \end{align*}\]

compute_derivatives(
    state::State,
    i::Integer,
    j::Integer,
    kd::Integer,
    ku::Integer,
    phitype::DUDZ,
)::NTuple{2, <:AbstractFloat}

Compute and return the vertical derivative of the zonal wind ($\partial u_\mathrm{b} / \partial z$) at $\left(i + 1 / 2, j, k_\mathrm{D} + 1 / 2\right)$ and $\left(i + 1 / 2, j, k_\mathrm{U} + 1 / 2\right)$.

The derivative is given by

\[\left(\frac{\partial u_\mathrm{b}}{\partial z}\right)_{i + 1 / 2, k + 1 / 2} = \frac{u_{\mathrm{b}, i + 1 / 2, k + 1} - u_{\mathrm{b}, i + 1 / 2}}{J_{i + 1 / 2, k + 1 / 2} \Delta \widehat{z}}.\]

At grid points beyond the vertical boundaries, it is set to zero.

compute_derivatives(
    state::State,
    i::Integer,
    j::Integer,
    kd::Integer,
    ku::Integer,
    phitype::DVDX,
)::NTuple{2, <:AbstractFloat}

Compute and return the zonal derivative of the meridional wind ($\partial v_\mathrm{b} / \partial x$) at $\left(i + 1 / 2, j + 1 / 2, k_\mathrm{D}\right)$ and $\left(i + 1 / 2, j + 1 / 2, k_\mathrm{U}\right)$.

The derivative is given by

\[\left(\frac{\partial v_\mathrm{b}}{\partial x}\right)_{i + 1 / 2, j + 1 / 2} = \frac{v_{\mathrm{b}, i + 1, j + 1 / 2} - v_{\mathrm{b}, j + 1 / 2}}{\Delta \widehat{x}} + G_{i + 1 / 2, j + 1 / 2}^{13} \frac{v_{\mathrm{b}, i + 1 / 2, j + 1 / 2, k + 1} - v_{\mathrm{b}, i + 1 / 2, j + 1 / 2, k - 1}}{2 \Delta \widehat{z}}.\]

compute_derivatives(
    state::State,
    i::Integer,
    j::Integer,
    kd::Integer,
    ku::Integer,
    phitype::DVDY,
)::NTuple{2, <:AbstractFloat}

Compute and return the meridional derivative of the meridional wind ($\partial v_\mathrm{b} / \partial y$) at $\left(i, j, k_\mathrm{D}\right)$ and $\left(i, j, k_\mathrm{U}\right)$.

The derivative is given by

\[\left(\frac{\partial v_\mathrm{b}}{\partial y}\right) = \frac{v_{\mathrm{b}, j + 1 / 2} - v_{\mathrm{b}, j - 1 / 2}}{\Delta \widehat{y}} + G^{23} \frac{v_{\mathrm{b}, k + 1} - v_{\mathrm{b}, k - 1} }{2 \Delta \widehat{z}}.\]

compute_derivatives(
    state::State,
    i::Integer,
    j::Integer,
    kd::Integer,
    ku::Integer,
    phitype::DVDZ,
)::NTuple{2, <:AbstractFloat}

Compute and return the vertical derivative of the meridional wind ($\partial v_\mathrm{b} / \partial z$) at $\left(i, j + 1 / 2, k_\mathrm{D} + 1 / 2\right)$ and $\left(i, j + 1 / 2, k_\mathrm{U} + 1 / 2\right)$.

The derivative is given by

\[\left(\frac{\partial v_\mathrm{b}}{\partial z}\right)_{j + 1 / 2, k + 1 / 2} = \frac{v_{\mathrm{b}, j + 1 / 2, k + 1} - v_{\mathrm{b}, j + 1 / 2}}{J_{j + 1 / 2, k + 1 / 2} \Delta \widehat{z}}.\]

At grid points beyond the vertical boundaries, it is set to zero.

compute_derivatives(
    state::State,
    i::Integer,
    j::Integer,
    kd::Integer,
    ku::Integer,
    phitype::DChiDX,
)::NTuple{2, <:AbstractFloat}

Compute and return the zonal derivative of the tracer field ($\partial \chi_\mathrm{b} / \partial x$) at $\left(i + 1 / 2, j, k_\mathrm{D}\right)$ and $\left(i + 1 / 2, j, k_\mathrm{U}\right)$.

The derivative is given by

\[\left(\frac{\partial \chi_\mathrm{b}}{\partial x}\right)_{i + 1 / 2} = \frac{\chi_{\mathrm{b}, i + 1} - \chi_\mathrm{b}}{\Delta \widehat{x}} + G_{i + 1 / 2}^{13} \frac{\chi_{\mathrm{b}, i + 1 / 2, k + 1} - \chi_{\mathrm{b}, i + 1 / 2, k - 1}}{2 \Delta \widehat{z}}.\]

compute_derivatives(
    state::State,
    i::Integer,
    j::Integer,
    kd::Integer,
    ku::Integer,
    phitype::DChiDY,
)::NTuple{2, <:AbstractFloat}

Compute and return the meridional derivative of the tracer field ($\partial \chi_\mathrm{b} / \partial y$) at $\left(i, j + 1 / 2, k_\mathrm{D}\right)$ and $\left(i, j + 1 / 2, k_\mathrm{U}\right)$.

The derivative is given by

\[\left(\frac{\partial \chi_\mathrm{b}}{\partial y}\right)_{j + 1 / 2} = \frac{\chi_{\mathrm{b}, j + 1} - \chi_\mathrm{b}}{\Delta \widehat{y}} + G_{j + 1 / 2}^{23} \frac{\chi_{\mathrm{b}, j + 1 / 2, k + 1} - \chi_{\mathrm{b}, j + 1 / 2, k - 1}}{2 \Delta \widehat{z}}.\]

compute_derivatives(
    state::State,
    i::Integer,
    j::Integer,
    kd::Integer,
    ku::Integer,
    phitype::DChiDZ,
)::NTuple{2, <:AbstractFloat}

Compute and return the vertical derivative of the tracer field ($\partial \chi_\mathrm{b} / \partial z$) at $\left(i, j, k_\mathrm{D} + 1 / 2\right)$ and $\left(i, j, k_\mathrm{U} + 1 / 2\right)$.

The derivative is given by

\[\left(\frac{\partial \chi_\mathrm{b}}{\partial z}\right)_{k + 1 / 2} = \frac{\chi_{\mathrm{b}, k + 1} - \chi_\mathrm{b}}{J_{k + 1 / 2} \Delta \widehat{z}}.\]

At grid points beyond the vertical boundaries, it is set to zero.

Arguments

• state: Model state.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• kd: Lower vertical grid-cell index.

• ku: Upper vertical grid-cell index.

• phitype: Type of derivative to compute.

get_next_half_level(
    i::Integer,
    j::Integer,
    z::AbstractFloat,
    state::State;
    dkd::Integer = 0,
    dku::Integer = 0,
)::Integer

Determine and return the index of the next half-level above z at the horizontal position $\left(i, j\right)$.

This method is heavily used for interpolation to ray-volume positions. To ensure that the vertical boundary conditions are met and no out-of-bounds errors occur, the following constraints are set.

• In MPI processes at the lower boundary of the domain, the computed index has the lower bound state.domain.k0. In other processes, an error is thrown if it is below 1 + dkd.

• In MPI processes at the upper boundary of the domain, the computed index has the upper bound state.domain.k1. In other processes, an error is thrown if it is above state.domain.nzz - dku.

In case an error is thrown, the parameter wkb_cfl_number of the discretization namelist should be set to a smaller value.

Arguments

• i: Zonal index.

• j: Meridional index.

• z: Vertical position.

• state: Model state.

Keywords

• dkd: Number of levels needed below.

• dku: Number of levels needed above.

get_next_level(
    i::Integer,
    j::Integer,
    z::AbstractFloat,
    state::State;
    dkd::Integer = 0,
    dku::Integer = 0,
)::Integer

Determine and return the index of the next level above z at the horizontal position $\left(i, j\right)$.

This method is heavily used for interpolation to ray-volume positions. To ensure that the vertical boundary conditions are met and no out-of-bounds errors occur, the following constraints are set.

• In MPI processes at the lower boundary of the domain, the computed index has the lower bound state.domain.k0. In other processes, an error is thrown if it is below 1 + dkd.

• In MPI processes at the upper boundary of the domain, the computed index has the upper bound state.domain.k1 + 1. In other processes, an error is thrown if it is above state.domain.nzz - dku.

In case an error is thrown, the parameter wkb_cfl_number of the discretization namelist should be set to a smaller value.

Arguments

• i: Zonal index.

• j: Meridional index.

• z: Vertical position.

• state: Model state.

Keywords

• dkd: Number of levels needed below.

• dku: Number of levels needed above.

interpolate(
    state::State;
    philbd::AbstractFloat = NaN,
    philbu::AbstractFloat = NaN,
    philfd::AbstractFloat = NaN,
    philfu::AbstractFloat = NaN,
    phirbd::AbstractFloat = NaN,
    phirbu::AbstractFloat = NaN,
    phirfd::AbstractFloat = NaN,
    phirfu::AbstractFloat = NaN,
    zlbd::AbstractFloat = NaN,
    zlbu::AbstractFloat = NaN,
    zlfd::AbstractFloat = NaN,
    zlfu::AbstractFloat = NaN,
    zrbd::AbstractFloat = NaN,
    zrbu::AbstractFloat = NaN,
    zrfd::AbstractFloat = NaN,
    zrfu::AbstractFloat = NaN,
    zlc::AbstractFloat = NaN,
    yb::AbstractFloat = NaN,
    yf::AbstractFloat = NaN,
    ylc::AbstractFloat = NaN,
    xl::AbstractFloat = NaN,
    xr::AbstractFloat = NaN,
    xlc::AbstractFloat = NaN,
)::AbstractFloat

Perform trilinear interpolation to (xlc, ylc, zlc), with values from eight surrounding grid points (two zonal positions, two meridional positions and eight vertical positions), and return the result.

Out of the eight grid points, four each are assumed to be to the left, to the right, behind, in front of, below and above the location of interest. Due to the grid being terrain-following, this includes eight different vertical positions, but only two zonal and two meridional positions. This is handled by performing successive linear interpolations, where the vertical position is interpolated along with the field of interest.

The exact algorithm is as follows.

1. Interpolation in $x$:

   \[\begin{align*} \psi_\mathrm{BD} & = f_x \psi_\mathrm{LBD} + (1 - f_x) \psi_\mathrm{RBD},\\ \psi_\mathrm{BU} & = f_x \psi_\mathrm{LBU} + (1 - f_x) \psi_\mathrm{RBU},\\ \psi_\mathrm{FD} & = f_x \psi_\mathrm{LFD} + (1 - f_x) \psi_\mathrm{RFD},\\ \psi_\mathrm{FU} & = f_x \psi_\mathrm{LFU} + (1 - f_x) \psi_\mathrm{RFU} \end{align*}\]

2. Interpolation in $y$:

   \[\begin{align*} \psi_\mathrm{D} & = f_y \psi_\mathrm{BD} + (1 - f_y) \psi_\mathrm{FD},\\ \psi_\mathrm{U} & = f_y \psi_\mathrm{BU} + (1 - f_y) \psi_\mathrm{FU} \end{align*}\]

3. Interpolation in $z$:

   \[\phi_\mathrm{C} = f_z \phi_\mathrm{D} + (1 - f_z) \phi_\mathrm{U}\]

Therein, $\mathrm{L}$, $\mathrm{R}$, $\mathrm{B}$, $\mathrm{F}$, $\mathrm{D}$ and $\mathrm{U}$ denote grid points to the left, to the right, forward, backward, downward and upward of the location of interest (denoted by $\mathrm{C}$), respectively, $\psi = \left(\phi, z\right)$ and

\[f_\alpha = \begin{cases} 0 & \mathrm{if} \quad \alpha_\beta = \alpha_\gamma,\\ 1 & \mathrm{if} \quad \alpha_\mathrm{C} < \alpha_\beta,\\ \frac{\alpha_\gamma - \alpha_\mathrm{C}}{\alpha_\gamma - \alpha_\beta} & \mathrm{if} \quad \alpha_\beta \leq \alpha_\mathrm{C} \leq \alpha_\gamma,\\ 0 & \mathrm{if} \quad \alpha_\gamma < \alpha_\mathrm{C}, \end{cases}\]

where $\left(\alpha, \beta, \gamma\right) \in \left\{\left(x, \mathrm{L}, \mathrm{R}\right), \left(y, \mathrm{B}, \mathrm{F}\right), \left(z, \mathrm{D}, \mathrm{U}\right)\right\}$.

Due to their large number, the positions and values are given as keyword arguments with the default value NaN, so that their order does not matter and calls with missing arguments are easy to detect.

Arguments

• state: Model state.

Keywords

• philbd: Value at the point to the left, behind and below.

• philbu: Value at the point to the left, behind and above.

• philfd: Value at the point to the left, in front and below.

• philfu: Value at the point to the left, in front and above.

• phirbd: Value at the point to the right, behind and below.

• phirbu: Value at the point to the right, behind and above.

• phirfd: Value at the point to the right, in front and below.

• phirfu: Value at the point to the right, in front and above.

• zlbd: Vertical coordinate of the point to the left, behind and below.

• zlbu: Vertical coordinate of the point to the left, behind and above.

• zlfd: Vertical coordinate of the point to the left, in front and below.

• zlfu: Vertical coordinate of the point to the left, in front and above.

• zrbd: Vertical coordinate of the point to the right, behind and below.

• zrbu: Vertical coordinate of the point to the right, behind and above.

• zrfd: Vertical coordinate of the point to the right, in front and below.

• zrfu: Vertical coordinate of the point to the right, in front and above.

• zlc: Vertical position of interest.

• yb: Meridional coordinate of the points behind.

• yf: Meridional coordinate of the points in front.

• ylc: Meridional position of interest.

• xl: Zonal coordinate of the points to the left.

• xr: Zonal coordinate of the points to the right.

• xlc: Zonal position of interest.

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::U,
)::AbstractFloat

Interpolate the zonal wind ($u_\mathrm{b}$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x} + \Delta \widehat{x} / 2$ and $\widehat{y}$ that are closest to xlc and ylc, respectively. For each of these four horizontal positions, it then determines the two points in $z$ that are closest to zlc. The resulting eight grid points are used to interpolate $u_\mathrm{b}$ to the location of interest, using interpolate.

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::V,
)::AbstractFloat

Interpolate the meridional wind ($v_\mathrm{b}$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x}$ and $\widehat{y} + \Delta \widehat{y} / 2$ that are closest to xlc and ylc, respectively. The steps that follow are analogous to those in the method for the zonal wind ($u_\mathrm{b}$).

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::W,
)::AbstractFloat

Interpolate the vertical wind ($w_\mathrm{b}$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x}$ and $\widehat{y}$ that are closest to xlc and ylc, respectively. For each of these four horizontal positions, it then determines the two points in $z + J \Delta \widehat{z} / 2$ that are closest to zlc. The resulting eight grid points are used to interpolate $w_\mathrm{b}$ to the location of interest, using compute_vertical_wind and interpolate. At grid points beyond the vertical boundaries, the values used in the interpolation are replaced with zeros.

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::DUDX,
)::AbstractFloat

Interpolate the zonal derivative of the zonal wind ($\partial u_\mathrm{b} / \partial x$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x}$ and $\widehat{y}$ that are closest to xlc and ylc, respectively. For each of these four horizontal positions, it then determines the two points in $z$ that are closest to zlc. The resulting eight grid points are used to interpolate $\partial u_\mathrm{b} / \partial x$ to the location of interest, using compute_derivatives and interpolate.

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::DUDY,
)::AbstractFloat

Interpolate the meridional derivative of the zonal wind ($\partial u_\mathrm{b} / \partial y$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x} + \Delta \widehat{x} / 2$ and $\widehat{y} + \Delta \widehat{y} / 2$ that are closest to xlc and ylc, respectively. For each of these four horizontal positions, it then determines the two points in $z$ that are closest to zlc. The resulting eight grid points are used to interpolate $\partial u_\mathrm{b} / \partial y$ to the location of interest, using compute_derivatives and interpolate.

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::DUDZ,
)::AbstractFloat

Interpolate the vertical derivative of the zonal wind ($\partial u_\mathrm{b} / \partial z$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x} + \Delta \widehat{x} / 2$ and $\widehat{y}$ that are closest to xlc and ylc, respectively. For each of these four horizontal positions, it then determines the two points in $z + J \Delta \widehat{z} / 2$ that are closest to zlc. The resulting eight grid points are used to interpolate $\partial u_\mathrm{b} / \partial z$ to the location of interest, using compute_derivatives and interpolate.

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::DVDX,
)::AbstractFloat

Interpolate the zonal derivative of the meridional wind ($\partial v_\mathrm{b} / \partial x$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x} + \Delta \widehat{x} / 2$ and $\widehat{y} + \Delta \widehat{y} / 2$ that are closest to xlc and ylc, respectively. The steps that follow are analogous to those in the method for the meridional derivative of the zonal wind ($\partial u_\mathrm{b} / \partial y$).

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::DVDY,
)::AbstractFloat

Interpolate the meridional derivative of the meridional wind ($\partial v_\mathrm{b} / \partial y$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x}$ and $\widehat{y}$ that are closest to xlc and ylc, respectively. The steps that follow are analogous to those in the method for the zonal derivative of the zonal wind ($\partial u_\mathrm{b} / \partial x$).

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::DVDZ,
)::AbstractFloat

Interpolate the vertical derivative of the meridional wind ($\partial v_\mathrm{b} / \partial z$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x}$ and $\widehat{y} + \Delta \widehat{y} / 2$ that are closest to xlc and ylc, respectively. The steps that follow are analogous to those in the method for the vertical derivative of the zonal wind ($\partial u_\mathrm{b} / \partial z$).

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::DChiDX,
)::AbstractFloat

Interpolate the zonal derivative of the tracer mixing ratio ($\partial \chi_\mathrm{b} / \partial x$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x} + \Delta \widehat{x} / 2$ and $\widehat{y}$ that are closest to xlc and ylc, respectively. For each of these four horizontal positions, it then determines the two points in $z$ that are closest to zlc. The resulting eight grid points are used to interpolate $\partial \chi_\mathrm{b} / \partial x$ to the location of interest, using compute_derivatives and interpolate.

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::DChiDY,
)::AbstractFloat

Interpolate the meridional derivative of the tracer mixing ratio ($\partial \chi_\mathrm{b} / \partial y$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x}$ and $\widehat{y} + \Delta \widehat{y} / 2$ that are closest to xlc and ylc, respectively. For each of these four horizontal positions, it then determines the two points in $z$ that are closest to zlc. The resulting eight grid points are used to interpolate $\partial \chi_\mathrm{b} / \partial y$ to the location of interest, using compute_derivatives and interpolate.

interpolate_mean_flow(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
    phitype::DChiDZ,
)::AbstractFloat

Interpolate the vertical derivative of the tracer mixing ratio ($\partial \chi_\mathrm{b} / \partial z$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x}$ and $\widehat{y}$ that are closest to xlc and ylc, respectively. For each of these four horizontal positions, it then determines the two points in $z + J \Delta \widehat{z} / 2$ that are closest to zlc. The resulting eight grid points are used to interpolate $\partial \chi_\mathrm{b} / \partial z$ to the location of interest, using compute_derivatives and interpolate.

Arguments

• xlc: Zonal position of interest.

• ylc: Meridional position of interest.

• zlc: Vertical position of interest.

• state: Model state.

• phitype: Mean-flow quantity to interpolate.

See also

• PinCFlow.MSGWaM.Interpolation.get_next_level

• PinCFlow.MSGWaM.Interpolation.interpolate

• PinCFlow.MSGWaM.Interpolation.get_next_half_level

• PinCFlow.Update.compute_vertical_wind

• PinCFlow.MSGWaM.Interpolation.compute_derivatives

interpolate_sponge(
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    state::State,
)::AbstractFloat

Interpolate the Rayleigh-damping coefficient of the LHS sponge ($\alpha_\mathrm{R}$) to (xlc, ylc, zlc), using a trilinear-interpolation algorithm, and return the result.

This method first determines the two points in $\widehat{x}$ and $\widehat{y}$ that are closest to xlc and ylc, respectively. For each of these four horizontal positions, it then determines the two points in $z$ that are closest to zlc. The resulting eight grid points are used to interpolate $\alpha_\mathrm{R}$ to the location of interest, using interpolate.

Arguments

• xlc: Zonal position of interest.

• ylc: Meridional position of interest.

• zlc: Vertical position of interest.

• state: Model state.

See also

• PinCFlow.MSGWaM.Interpolation.get_next_level

• PinCFlow.MSGWaM.Interpolation.interpolate

interpolate_stratification(
    zlc::AbstractFloat,
    state::State,
    strtype::N2,
)::AbstractFloat

Interpolate the squared buoyancy frequency ($N^2$) to zlc and return the result.

This method first determines the two points in $z$ that are closest to zlc. As horizontal position, it uses (i0, j0), which is arbitrary, since $N^2$ has no horizontal dependence. Subsequently, simple linear interpolation is performed to find $N^2$ at zlc.

interpolate_stratification(
    zlc::AbstractFloat,
    state::State,
    strtype::DN2DZ,
)::AbstractFloat

Interpolate the vertical derivative of the squared buoyancy frequency ($\partial N^2 / \partial z$) to zlc and return the result.

This method first determines the two points in $z + J \Delta \widehat{z} / 2$ that are closest to zlc. As for $N^2$, (i0, j0) is used as the horizontal position, and simple linear interpolation is performed to find $\partial N^2 / \partial z$ at zlc.

Arguments

• zlc: Vertical position of interest.

• state: Model state.

• strtype: Stratification quantity to interpolate.

See also

• PinCFlow.MSGWaM.Interpolation.get_next_level

• PinCFlow.MSGWaM.Interpolation.get_next_half_level

RayOperations

Module for various ray-volume operations needed throughout MSGWaM.

See also

• PinCFlow.Types

• PinCFlow.MSGWaM.Interpolation

check_rays(state::State)

Check if all ray volumes are assigned to the correct grid cells.

Arguments

• state: Model state.

compute_intrinsic_frequency(
    state::State,
    r::Integer,
    i::Integer,
    j::Integer,
    k::Integer,
)::AbstractFloat

Return the intrinsic frequency of the ray volume specified by $\left(r, i, j, k\right)$.

The intrinsic frequency is calculated from the dispersion relation

\[\widehat{\omega}_r = \sigma \sqrt{\frac{N_r^2 \left(k_r^2 + l_r^2\right) + f^2 m_r^2}{\left|\boldsymbol{k}_r\right|^2}},\]

where $N_r^2$ is the squared buoyancy frequency interpolated to the ray volume's vertical position and $\sigma$ is the frequency branch.

Arguments

• state: Model state

• r: Ray-volume index.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• k: Vertical grid-cell index.

See also

• PinCFlow.MSGWaM.Interpolation.interpolate_stratification

compute_merge_index(
    wnr::AbstractFloat,
    wnr_min_p::AbstractFloat,
    wnr_max_p::AbstractFloat,
    wnr_min_n::AbstractFloat,
    wnr_max_n::AbstractFloat,
    dwnr_mrg_p::AbstractFloat,
    dwnr_mrg_n::AbstractFloat,
    nray::Integer,
)::Integer

Return the index of the wavenumber wnr on a 1D spectral grid specified by the other arguments.

This method is used by PinCFlow.MSGWaM.RayUpdate.merge_rays! to sort ray volumes into spectral bins.

Arguments

• wnr: Wavenumber value.

• wnr_min_p: Minimum positive wavenumber.

• wnr_max_p: Maximum positive wavenumber.

• wnr_min_n: Minimum negative wavenumber.

• wnr_max_n: Maximum negative wavenumber.

• dwnr_mrg_p: Logarithmic spacing of discrete positive wavenumbers.

• dwnr_mrg_n: Logarithmic spacing of discrete negative wavenumbers.

• nray: Number of spectral grid points.

compute_saturation_integrals(
    state::State,
    i::Integer,
    j::Integer,
    k::Integer,
)::NTuple{2, <:AbstractFloat}

Compute and return the two spectral integrals $S_1$ and $S_2$, as needed by the saturation scheme (in the grid cell $\left(i, j, k\right)$).

Computes the sums

\[\begin{align*} S_1 & \approx \sum\limits_r \left(m_r \left|b_{\mathrm{w}, r}\right|\right)^2 f_r,\\ S_2 & \approx \sum\limits_r \left(m_r \left|b_{\mathrm{w}, r}\right| \left|\boldsymbol{k}_r\right|\right)^2 f_r, \end{align*}\]

where

\[f_r = \max \left(1, \frac{\Delta x_r}{\Delta \widehat{x}}\right) \max \left(1, \frac{\Delta y_r}{\Delta \widehat{y}}\right) \max \left(1, \frac{\Delta z_r}{J \Delta \widehat{z}}\right)\]

is the maximum grid-cell fraction that can be covered by each ray volume (with $\left(\Delta x_r, \Delta y_r, \Delta z_r\right)$ being the ray-volume extents in physical space) and

\[\left|b_{\mathrm{w}, r}\right|^2 = \frac{2}{\overline{\rho}} \frac{N_r^4 \left(k_r^2 + l_r^2\right)}{\widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2} \mathcal{N}_r \Delta k_r \Delta l_r \Delta m_r\]

is the squared gravity-wave amplitude of the buoyancy. Therein, $N_r^2$ is the squared buoyancy frequency interpolated to the ray-volume position (using interpolate_stratification) and $\left(\Delta k_r, \Delta l_r, \Delta m_r\right)$ are the ray-volume extents in spectral space.

Arguments

• state: Model state.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• k: Vertical grid-cell index.

See also

• PinCFlow.MSGWaM.Interpolation.get_next_half_level

• PinCFlow.MSGWaM.RayOperations.compute_intrinsic_frequency

• PinCFlow.MSGWaM.Interpolation.interpolate_stratification

compute_spectral_bounds(
    wavenumbers::AbstractVector{<:AbstractFloat},
)::NTuple{4, <:AbstractFloat}

Compute the minima and maxima of positive and negative entries in wavenumbers and return them (in the order positive minimum, positive maximum, negative minimum, negative maximum).

This method is used by PinCFlow.MSGWaM.RayUpdate.merge_rays! to create spectral bins.

Arguments

• wavenumbers: Vector of wavenumbers in the considered spectral dimension.

compute_wave_action_integral(
    merge_mode::ConstantWaveAction,
    nr::AbstractFloat,
    omegar::AbstractFloat,
    fxk::AbstractFloat,
    fyl::AbstractFloat,
    fzm::AbstractFloat,
)

Return the wave action obtained by multiplying the given phase-space wave-action density with the given phase-space volume.

This method is used to implement conservation of wave action in ray-volume merging.

compute_wave_action_integral(
    merge_mode::ConstantWaveEnergy,
    nr::AbstractFloat,
    omegar::AbstractFloat,
    fxk::AbstractFloat,
    fyl::AbstractFloat,
    fzm::AbstractFloat,
)::AbstractFloat

Return the wave energy obtained by multiplying the given phase-space wave-action density with the given intrinsic frequency and phase-space volume.

This method is used to implement conservation of wave energy in ray-volume merging.

Arguments

• merge_mode: Merging strategy.

• nr: Phase-space wave-action density.

• omegar: Intrinsic frequency.

• fxk: Phase space factor in the $x$-$k$ subspace.

• fyl: Phase space factor in the $y$-$l$ subspace.

• fzm: Phase space factor in the $z$-$m$ subspace.

copy_rays!(
    rays::Rays,
    r::Pair{<:Integer, <:Integer},
    i::Pair{<:Integer, <:Integer},
    j::Pair{<:Integer, <:Integer},
    k::Pair{<:Integer, <:Integer},
)

Copy all properties of the ray volume specified by the first components of the index pairs to that specified by the second components.

Arguments

• rays: Collection of ray-volume-property arrays.

• r: Ray-volume indices.

• i: Zonal grid-cell indices.

• j: Meridional grid-cell indices.

• k: Vertical grid-cell indices.

get_physical_extent(
    rays::Rays,
    r::Integer,
    i::Integer,
    j::Integer,
    k::Integer,
)::NTuple{3, <:AbstractFloat}

Return the physical extents of the ray volume specified by $\left(r, i, j, k\right)$ as the tuple (dxr, dyr, dzr).

Arguments

• rays: Collection of ray-volume-property arrays.

• r: Ray-volume index.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• k: Vertical grid-cell index.

get_physical_position(
    rays::Rays,
    r::Integer,
    i::Integer,
    j::Integer,
    k::Integer,
)::NTuple{3, <:AbstractFloat}

Return the physical position of the ray volume specified by $\left(r, i, j, k\right)$ as the tuple (xr, yr, zr).

Arguments

• rays: Collection of ray-volume-property arrays.

• r: Ray-volume index.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• k: Vertical grid-cell index.

get_spectral_extent(
    rays::Rays,
    r::Integer,
    i::Integer,
    j::Integer,
    k::Integer,
)::NTuple{3, <:AbstractFloat}

Return the spectral extents of the ray volume specified by $\left(r, i, j, k\right)$ as the tuple (dkr, dlr, dmr).

Arguments

• rays: Collection of ray-volume-property arrays.

• r: Ray-volume index.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• k: Vertical grid-cell index.

get_spectral_position(
    rays::Rays,
    r::Integer,
    i::Integer,
    j::Integer,
    k::Integer,
)::NTuple{3, <:AbstractFloat}

Return the spectral position of the ray volume specified by $\left(r, i, j, k\right)$ as the tuple (kr, lr, mr).

Arguments

• rays: Collection of ray-volume-property arrays.

• r: Ray-volume index.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• k: Vertical grid-cell index.

get_surfaces(
    rays::Rays,
    r::Integer,
    i::Integer,
    j::Integer,
    k::Integer,
)::NTuple{3, <:AbstractFloat}

Compute phase-space surface areas of the ray volume specified by $\left(r, i, j, k\right)$ and return them as the tuple (axk, ayl, azm).

Arguments

• rays: Collection of ray-volume-property arrays.

• r: Ray-volume index.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• k: Vertical grid-cell index.

remove_rays!(state::State)

Remove gaps (i.e. zero-wave-action ray volumes between nonzero-wave-action ray volumes) in the ray-volume arrays.

In each grid cell, this method moves all ray volumes as far to the front of the arrays possible and updates nray accordingly, so that every ray volume in the range 1:nray[i, j, k] has nonzero wave action.

Arguments

• state: Model state.

See also

• PinCFlow.MSGWaM.RayOperations.copy_rays!

update_merged_rays!(
    merge_mode::AbstractMergeMode,
    merged_rays::MergedRays,
    bin::Integer,
    xr::AbstractFloat,
    dxr::AbstractFloat,
    yr::AbstractFloat,
    dyr::AbstractFloat,
    zr::AbstractFloat,
    dzr::AbstractFloat,
    kr::AbstractFloat,
    dkr::AbstractFloat,
    lr::AbstractFloat,
    dlr::AbstractFloat,
    mr::AbstractFloat,
    dmr::AbstractFloat,
    fxk::AbstractFloat,
    fyl::AbstractFloat,
    fzm::AbstractFloat,
    nr::AbstractFloat,
    omegar::AbstractFloat,
)

Update the fields of merged_rays at bin such that they contain the outermost bounds and total wave action/energy of all contributing ray volumes.

This method is used to compute the properties of merged ray volumes. It is called for every old ray volume that contributes to the new, merged volume and updates the outermost bounds in physical and spectral space, as well as the total wave action/energy, accordingly.

Arguments

• merge_mode: Merging strategy.

• merged_rays: Properties of merged ray volumes.

• bin: Index of the bin to update.

• xr: Position of the old ray volume in $x$.

• dxr: Extent of the old ray volume in $x$.

• yr: Position of the old ray volume in $y$.

• dyr: Extent of the old ray volume in $y$.

• zr: Position of the old ray volume in $z$.

• dzr: Extent of the old ray volume in $z$.

• kr: Position of the old ray volume in $k$.

• dkr: Extent of the old ray volume in $k$.

• lr: Position of the old ray volume in $l$.

• dlr: Extent of the old ray volume in $k$.

• mr: Position of the old ray volume in $m$.

• dmr: Extent of the old ray volume in $m$.

• fxk: Phase-space factor of the old ray volume in $x$-$k$ subspace.

• fyl: Phase-space factor of the old ray volume in $y$-$l$ subspace.

• fzm: Phase-space factor of the old ray volume in $z$-$m$ subspace.

• nr: Phase-space wave-action density of the old ray volume.

• omegar: Intrinsic frequency of the old ray volume.

See also

• PinCFlow.MSGWaM.RayOperations.compute_wave_action_integral

RaySources

Module for ray-volume sources.

See also

• PinCFlow.Types

• PinCFlow.MSGWaM.RayOperations

activate_orographic_source!(
    state::State,
    omi_ini::AbstractArray{<:AbstractFloat, 4},
    wnk_ini::AbstractArray{<:AbstractFloat, 4},
    wnl_ini::AbstractArray{<:AbstractFloat, 4},
    wnm_ini::AbstractArray{<:AbstractFloat, 4},
    wad_ini::AbstractArray{<:AbstractFloat, 4},
)

Compute ray-volume properties in the launch layer (i.e. at k = k0 - 1) for the initialization of MSGWaM.

Sets the launch-layer values of arrays for initial ray-volume properties (intrinsic frequencies, wavenumbers and wave-action densities). For this purpose, the horizontal components of the resolved wind, the background density and the squared buoyancy frequency are vertically averaged between the surface and an approximation for the summits of the unresolved orography. The vertical averages are then used to compute a non-dimensionalized mountain wave amplitude, from which an approximate reduction of the generated wave amplitude due to blocking is inferred (see below). Afterwards, the ray-volume properties are obtained by calling compute_orographic_mode with the correspondingly scaled mode of the orographic spectrum and the vertical averages as arguments.

activate_orographic_source!(state::State)

Launch ray volumes that represent unresolved orographic gravity waves.

In each column of MPI processes at the lower boundary, this method first computes vertical averages of the horizontal components of the resolved wind, the background density and the squared buoyancy frequency between $h_\mathrm{b}$ (the surface) and $h_\mathrm{b} + \Delta h$ (an approximation for the summits of the unresolved orography, with $\Delta h = \sum_\alpha \left|h_{\mathrm{w}, \alpha}\right|$). The vertical averages are then used to compute a non-dimensionalized mountain wave amplitude, from which an approximate reduction of the generated wave amplitude due to blocking, as well as the depth of the blocked layer, is inferred. A loop over the spectral modes of the unresolved orography follows, in which the properties of each mode are computed, using compute_orographic_mode with the scaled mode of the orographic spectrum and vertical averages as arguments, and corresponding ray volumes are launched at k = k0 - 1.

The parameterization of blocking is built around the non-dimensionalized mountain wave amplitude, or Long number,

\[\mathrm{Lo} = \frac{N_h \Delta h}{\left|\boldsymbol{u}_h\right|},\]

where $N_h$ is the square root of the vertically averaged squared buoyancy frequency and $\boldsymbol{u}_h$ is the vertically averaged resolved horizontal wind. This number is used to estimate the depth of the blocked layer as

\[\Delta z_\mathrm{B} = 2 \Delta h \max \left(0, \frac{\mathrm{Lo} - C}{\mathrm{Lo}}\right),\]

where $C$ is a critical value represented by the model parameter state.namelists.wkb.long_threshold. The corresponding scaling of the orographic spectrum is given by

\[r \left(\mathrm{Lo}\right) = \frac{2 \Delta h - \Delta z_\mathrm{B}}{2 \Delta h} = \min \left(1, \frac{C}{\mathrm{Lo}}\right),\]

so that $\Delta z_\mathrm{B} = 2 \Delta h \left(1 - r\right)$. In addition to the reduction of the mountain-wave amplitude, the present blocked-layer scheme adds a blocked-flow drag to the mean-flow impact. This is implemented in PinCFlow.MSGWaM.MeanFlowEffect.apply_blocked_layer_scheme!.

The launch algorithm distinguishes between the following situations (regarding previously launched ray volumes).

1. There is no ray volume with nonzero wave-action density. A new ray volume is launched.

2. There is a ray volume with nonzero wave-action density that has at least partially passed through the lower boundary. The ray volume is either clipped or extended, such that its lower edge coincides with the surface, and the part below the surface is discarded. Then, it is assigned to the first model layer k0, i.e. its indices are changed from (r, i, j, k0 - 1) to (rray, i, j, k0), where rray is the new last ray-volume index at (i, j, k0). Finally, a new ray volume is launched.

3. There is a ray volume with nonzero wave-action density, which has not yet crossed the lower boundary. It is replaced with a new one.

Arguments

• state: Model state.

• omi_ini: Array for intrinsic frequencies.

• wnk_ini: Array for zonal wavenumbers.

• wnl_ini: Array for meridional wavenumbers.

• wnm_ini: Array for vertical wavenumbers.

• wad_ini: Array for wave-action densities.

See also

• PinCFlow.MSGWaM.RaySources.compute_orographic_mode

• PinCFlow.MSGWaM.RayOperations.copy_rays!

compute_orographic_mode(
    displm::AbstractFloat,
    wnk::AbstractFloat,
    wnl::AbstractFloat,
    uavg::AbstractFloat,
    vavg::AbstractFloat,
    rhoavg::AbstractFloat,
    bvsavg::AbstractFloat,
    fc::AbstractFloat,
    branch::Integer,
)::NTuple{5, <:AbstractFloat}

Compute and return the properties of an orographic-gravity-wave mode.

Calculates the intrinsic frequency, wavenumbers and wave-action density of a gravity-wave mode generated by flow over topography, using linear mountain-wave theory. The implemented formulas are as follows (where this method's arguments are represented by $h_{\mathrm{w}, \alpha}$, $k_{h, \alpha}$, $l_{h, \alpha}$, $u_h$, $v_h$, $\overline{\rho}_h$, $N_h^2$, $f$ and $\sigma$).

• Intrinsic frequency:

  \[\widehat{\omega}_\alpha = \sigma \left|- k_{h, \alpha} u_h - l_{h, \alpha} v_h\right|\]

• Horizontal wavenumbers:

  \[\begin{align*} k_\alpha & = \sigma \mathrm{sgn} \left(- k_{h, \alpha} u_h - l_{h, \alpha} v_h\right) k_{h, \alpha},\\ l_\alpha & = \sigma \mathrm{sgn} \left(- k_{h, \alpha} u_h - l_{h, \alpha} v_h\right) l_{h, \alpha} \end{align*}\]

• Vertical wavenumber:

  \[m_\alpha = - \sigma \sqrt{\frac{\left(k_\alpha^2 + l_\alpha^2\right) \left(N_h^2 - \widehat{\omega}_\alpha^2\right)}{\widehat{\omega}_\alpha^2 - f^2}}\]

• Wave-action density:

  \[\mathcal{A}_\alpha = \frac{\overline{\rho}_h}{2} \frac{\widehat{\omega}_\alpha \left|\boldsymbol{k}_\alpha\right|^2}{k_\alpha^2 + l_\alpha^2} \left|h_{\mathrm{w}, \alpha}\right|^2\]

If the squared intrinsic frequency is smaller than the squared Coriolis parameter or larger than the squared buoyancy frequency (and thus outside of the gravity-wave spectrum), the vertical wavenumber and wave-action density are set to zero.

Arguments

• displm: Topographic displacement.

• wnk: Zonal wavenumber from the topographic spectrum.

• wnl: Meridional wavenumber from the topographic spectrum.

• uavg: Vertically averaged resolved zonal wind.

• vavg: Vertically averaged resolved meridional wind.

• rhoavg: Vertically averaged background density.

• bvsavg: Vertically averaged squared buoyancy frequency.

• fc: Coriolis parameter.

• branch: Frequency branch.

BoundaryRays

Module for enforcing boundary conditions for ray volumes.

Provides functions for configurations that are serial or parallel in any dimension of physical space. Assumes periodicity in the horizontal and solid-wall boundaries in the vertical.

See also

• PinCFlow.Types

• PinCFlow.MPIOperations

• PinCFlow.Boundaries

• PinCFlow.MSGWaM.RayOperations

set_boundary_rays!(state::State)

Enforce boundary conditions for ray volumes by dispatching to a WKB-mode-specific method.

set_boundary_rays!(state::State, wkb_mode::NoWKB)

Return for non-WKB configurations.

set_boundary_rays!(state::State, wkb_mode::SteadyState)

Enforce horizontal boundary conditions for "ray volumes" in steady-state mode.

Zonal (meridional) boundary conditions are only enforced if state.namelists.domain.x_size > 1 (state.namelists.domain.y_size > 1).

set_boundary_rays!(state::State, wkb_mode::Union{SingleColumn, MultiColumn})

Enforce horizontal and vertical boundary conditions for ray volumes in single-column or multi-column mode.

Zonal (meridional) boundary conditions are only enforced if state.namelists.domain.x_size > 1 (state.namelists.domain.y_size > 1).

Arguments

• state: Model state.

• wkb_mode: Approximations used by MSGWaM.

See also

• PinCFlow.MSGWaM.BoundaryRays.set_zonal_boundary_rays!

• PinCFlow.MSGWaM.BoundaryRays.set_meridional_boundary_rays!

• PinCFlow.MSGWaM.BoundaryRays.set_vertical_boundary_rays!

set_meridional_boundary_rays!(state::State)

Set the ray volumes at the meridional boundaries.

This method first enforces meridional boundary conditions for state.wkb.nray (by applying set_meridional_boundaries_of_field! to it) and then sets the corresponding boundary ray volumes, assuming periodicity. If the domain is parallelized in $\widehat{y}$, ray volumes are communicated between MPI processes, using set_meridional_halo_rays!. At the meridional boundaries of the domain, the $y$-coordinates of ray volumes are adjusted such that shifting works properly.

Arguments

• state: Model state.

See also

• PinCFlow.Boundaries.set_meridional_boundaries_of_field!

• PinCFlow.MSGWaM.BoundaryRays.set_meridional_halo_rays!

• PinCFlow.MSGWaM.RayOperations.copy_rays!

set_meridional_halo_rays!(state::State)

Exchange ray volumes in meridional halo cells.

Performs bidirectional MPI communication between backward and forward neighbor processes. The number of communicated ray volumes is determined from the maximum counts at the backward and forward boundaries of the MPI subdomains.

Arguments

• state: Model state.

set_vertical_boundary_rays!(state::State)

Enforce vertical boundary conditions for ray volumes.

If the domain is parallelized in $\widehat{z}$, ray-volume counts and the ray volumes themselves are first communicated between MPI processes, using set_vertical_halos_of_field! and set_vertical_halo_rays!, respectively. The vertical boundary conditions are then enforced by cutting (removing) ray volumes that have partially (fully) crossed the upper boundary and reflecting ray volumes (by adjusting the vertical position and wavenumber) that have at least partially crossed the lower boundary from above.

Arguments

• state: Model state.

See also

• PinCFlow.MPIOperations.set_vertical_halos_of_field!

• PinCFlow.MSGWaM.BoundaryRays.set_vertical_halo_rays!

• PinCFlow.MSGWaM.RayOperations.copy_rays!

set_vertical_halo_rays!(state::State)

Exchange ray volumes in vertical halo cells.

Performs MPI communication between downward and upward neighbor processes. The number of communicated ray volumes is determined from the maximum counts at the downward and upward boundaries of the MPI subdomains. Solid walls are assumed at the vertical boundaries of the full domain. The corresponding ghost-cell ray volumes are not changed.

Arguments

• state: Model state.

set_zonal_boundary_rays!(state::State)

Set the ray volumes at the zonal boundaries.

This method first enforces zonal boundary conditions for state.wkb.nray (by applying set_zonal_boundaries_of_field! to it) and then sets the corresponding boundary ray volumes, assuming periodicity. If the domain is parallelized in $\widehat{x}$, ray volumes are communicated between MPI processes, using set_zonal_halo_rays!. At the zonal boundaries of the domain, the $x$-coordinates of ray volumes are adjusted such that shifting works properly.

Arguments

• state: Model state.

See also

• PinCFlow.Boundaries.set_zonal_boundaries_of_field!

• PinCFlow.MSGWaM.BoundaryRays.set_zonal_halo_rays!

• PinCFlow.MSGWaM.RayOperations.copy_rays!

set_zonal_halo_rays!(state::State)

Exchange ray volumes in zonal halo cells.

Performs bidirectional MPI communication between left and right neighbor processes. The number of communicated ray volumes is determined from the maximum counts at the left and right boundaries of the MPI subdomains.

Arguments

• state: Model state.

RayUpdate

Module for the integration of the ray equations.

In addition to ray-volume initialization and propagation, functions for tracking ray volumes on the model grid and controlling their count, as well as a saturation scheme for capturing wave breaking, are provided.

See also

• PinCFlow.Types

• PinCFlow.MSGWaM.BoundaryRays

• PinCFlow.MSGWaM.Interpolation

• PinCFlow.MSGWaM.RayOperations

• PinCFlow.MSGWaM.RaySources

X

Singleton for dispatch to operations in $x$-direction.

XYZ

Singleton for dispatch to operations in all directions.

XZ

Singleton for dispatch to operations in $x$- and $z$-direction.

Y

Singleton for dispatch to operations in $y$-direction.

YZ

Singleton for dispatch to operations in $y$- and $z$-direction.

Z

Singleton for dispatch to operations in $z$-direction.

apply_saturation_scheme!(state::State, dt::AbstractFloat)

Apply the saturation scheme by dispatching to a WKB-mode-specific method.

apply_saturation_scheme!(
    state::State,
    dt::AbstractFloat,
    wkb_mode::Union{NoWKB, SteadyState},
)

Return for configurations without WKB / with steady-state WKB.

In steady-state mode, saturation is handled by PinCFlow.MSGWaM.RayUpdate.propagate_rays!.

apply_saturation_scheme!(
    state::State,
    dt::AbstractFloat,
    wkb_mode::Union{SingleColumn, MultiColumn},
)

Apply the saturation scheme.

Saturation is assumed to occur when the static-instability criterion

\[\sum\limits_r \left(m_r \left|b_{\mathrm{w}, r}\right|\right)^2 f_r \geq \alpha_\mathrm{s}^2 N^4\]

is locally fulfilled (i.e. within a grid cell). Therein, $\left|b_{\mathrm{w}, r}\right|^2$ is the squared gravity-wave amplitude of the buoyancy, $f_r$ is the maximum grid-cell fraction each ray volume can cover and $\alpha_\mathrm{s}$ is a saturation coefficient that represents the uncertainties of the criterion. The phase-space wave-action density is then reduced in accordance with

\[\frac{\Delta \left|b_{\mathrm{w}, r}\right|^2}{\Delta t} = - 2 K \left|\boldsymbol{k}_r\right|^2 \left|b_{\mathrm{w}, r}\right|^2,\]

which is based on the assumption that wave breaking leads to turbulent fluxes that may be parameterized with a flux-gradient ansatz. The turbulent viscosity and diffusivity

\[K = \left[2 \Delta t\sum\limits_r \left(m_r \left|b_{\mathrm{w}, r}\right| \left|\boldsymbol{k}_r\right|\right)^2 f_r\right]^{- 1} \max \left[0, \sum_r \left(m_r \left|b_{\mathrm{w}, r}\right|\right)^2 f_r - \alpha_\mathrm{s}^2 N^4\right]\]

is such that wave action is reduced exactly to the saturation threshold. The two sums involved in this scheme (discretizations of spectral integrals) are computed with compute_spectral_integrals.

Arguments

• state: Model state.

• dt: Time step.

• wkb_mode: Approximations used by MSGWaM.

See also

• PinCFlow.MSGWaM.RayOperations.compute_saturation_integrals

• PinCFlow.MSGWaM.Interpolation.interpolate_stratification

• PinCFlow.MSGWaM.RayOperations.remove_rays!

initialize_rays!(state::State)

Complete the initialization of MSGWaM by dispatching to a WKB-mode-specific method.

initialize_rays!(state::State, wkb_mode::NoWKB)

Return for non-WKB configurations.

initialize_rays!(
    state::State,
    wkb_mode::Union{SteadyState, SingleColumn, MultiColumn},
)

Complete the initialization of MSGWaM.

In each grid cell, wave_modes wave modes are computed, using state.namelists.wkb.initial_wave_field, as well as activate_orographic_source! for mountain waves. For each of these modes, nrx * nry * nrz * nrk * nrl * nrm ray volumes are then defined such that they evenly divide the volume one would get for nrx = nry = nrz = nrk = nrl = nrm = 1 (the parameters are taken from state.namelists.wkb). Finally, the maximum group velocities are determined for the corresponding CFL condition that is used in the computation of the time step.

Arguments

• state: Model state.

• wkb_mode: Approximations used by MSGWaM.

See also

• PinCFlow.MSGWaM.RaySources.activate_orographic_source!

• PinCFlow.MSGWaM.Interpolation.interpolate_stratification

• PinCFlow.MSGWaM.Interpolation.interpolate_mean_flow

merge_rays!(state::State)

Merge ray volumes by dispatching to a WKB-mode-specific method.

merge_rays!(state::State, wkb_mode::Union{NoWKB, SteadyState})

Return for configurations without WKB / with steady-state WKB.

merge_rays!(state::State, wkb_mode::Union{SingleColumn, MultiColumn})

Merge ray volumes in grid cells in which their count exceeds a threshold.

This method checks in each grid cell if the number of ray volumes exceeds a maximum that was determined from namelist parameters (state.wkb.nray_max). If it does, the ray volumes in that cell are merged such that the new count is smaller or equal to the threshold. This is done by binning them on a spectral grid with logarithmic spacing, defined from the minima and maxima of the contributing negative and positive wavenumbers in all spectral dimensions. The merging is performed such that the bounds of the new ray volumes coincide with the outermost bounds of the old ray volumes and wave action (or wave energy, depending on the merging strategy) is conserved.

Arguments

• state: Model state.

• wkb_mode: Approximations used by MSGWaM.

See also

• PinCFlow.MSGWaM.RayOperations.compute_spectral_bounds

• PinCFlow.MSGWaM.RayOperations.get_physical_position

• PinCFlow.MSGWaM.RayOperations.get_spectral_position

• PinCFlow.MSGWaM.RayOperations.get_physical_extent

• PinCFlow.MSGWaM.RayOperations.get_spectral_extent

• PinCFlow.MSGWaM.RayOperations.get_surfaces

• PinCFlow.MSGWaM.RayOperations.compute_intrinsic_frequency

• PinCFlow.MSGWaM.RayOperations.compute_merge_index

• PinCFlow.MSGWaM.RayOperations.update_merged_rays!

• PinCFlow.MSGWaM.RayOperations.compute_wave_action_integral

propagate_rays!(state::State, dt::AbstractFloat, rkstage::Integer)

Integrate the wave-action-density and ray equations by dispatching to a WKB-mode-specific method.

propagate_rays!(
    state::State,
    dt::AbstractFloat,
    rkstage::Integer,
    wkb_mode::NoWKB,
)

Return for non-WKB configurations.

propagate_rays!(
    state::State,
    dt::AbstractFloat,
    rkstage::Integer,
    wkb_mode::Union{SingleColumn, MultiColumn},
)

Integrate the wave-action-density and ray equations derived from 1D or 3D transient WKB theory.

The updates of the RK tendencies for the phase-space position of each ray volume are given by

\[\begin{align*} q_r^x & \rightarrow \Delta t \left(u_{\mathrm{b}, r} + k_r \frac{N_r^2 - \widehat{\omega}_r^2}{\widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2}\right) + \alpha_\mathrm{RK} q_r^x,\\ q_r^y & \rightarrow \Delta t \left(v_{\mathrm{b}, r} + l_r \frac{N_r^2 - \widehat{\omega}_r^2}{\widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2}\right) + \alpha_\mathrm{RK} q_r^y,\\ q_r^z & \rightarrow - \Delta t \frac{m_r \left(\widehat{\omega}_r^2 - f^2\right)}{\widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2} + \alpha_\mathrm{RK} q_r^z,\\ q_r^k & \rightarrow - \Delta t \left[k_r \left(\frac{\partial u_\mathrm{b}}{\partial x}\right)_r - l_r \left(\frac{\partial v_\mathrm{b}}{\partial x}\right)_r\right] + \alpha_\mathrm{RK} q_r^k,\\ q_r^l & \rightarrow - \Delta t \left[k_r \left(\frac{\partial u_\mathrm{b}}{\partial y}\right)_r - l_r \left(\frac{\partial v_\mathrm{b}}{\partial y}\right)_r\right] + \alpha_\mathrm{RK} q_r^l,\\ q_r^m & \rightarrow - \Delta t \left[k_r \left(\frac{\partial u_\mathrm{b}}{\partial z}\right)_r - l_r \left(\frac{\partial v_\mathrm{b}}{\partial z}\right)_r - \frac{k_r^2 + l_r^2}{2 \widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2} \left(\frac{\partial N^2}{\partial z}\right)_r\right] + \alpha_\mathrm{RK} q_r^m \end{align*}\]

and the position update is

\[\begin{align*} x_r & = x_r + \beta_\mathrm{RK} q_r^x, & y_r & \rightarrow y_r + \beta_\mathrm{RK} q_r^y, & z_r & \rightarrow z_r + \beta_\mathrm{RK} q_r^z,\\ k_r & \rightarrow k_r + \beta_\mathrm{RK} q_r^k, & l_r & \rightarrow l_r + \beta_\mathrm{RK} q_r^l, & m_r & \rightarrow m_r + \beta_\mathrm{RK} q_r^m, \end{align*}\]

where the subscript $r$ indicates either a ray-volume property or a mean-flow property interpolated to the ray-volume position, via interpolate_mean_flow and interpolate_stratification. In addition, MSGWaM updates the ray-volume extents, following

\[\begin{align*} q_r^{\Delta x} & \rightarrow \Delta t \left(u_{\mathrm{b}, r, +} - u_{\mathrm{b}, r, -}\right) + \alpha_\mathrm{RK} q_r^{\Delta x}, & \Delta x_r & \rightarrow \Delta x_r + \beta_\mathrm{RK} q_r^{\Delta x},\\ q_r^{\Delta y} & \rightarrow \Delta t \left(v_{\mathrm{b}, r, +} - v_{\mathrm{b}, r, -}\right) + \alpha_\mathrm{RK} q_r^{\Delta y}, & \Delta y_r & \rightarrow \Delta y_r + \beta_\mathrm{RK} q_r^{\Delta y},\\ q_r^{\Delta z} & \rightarrow \Delta t \left(c_{\mathrm{g} z, r, +} - c_{\mathrm{g} z, r, -}\right) + \alpha_\mathrm{RK} q_r^{\Delta z}, & \Delta z_r & \rightarrow \Delta z_r + \beta_\mathrm{RK} q_r^{\Delta z}, \end{align*}\]

where $u_{\mathrm{b}, r, \pm}$ is the interpolation of $u_\mathrm{b}$ to $x_{r, \pm} = x_r \pm \Delta x_r / 2$ (from before the position update) and $v_{\mathrm{b}, r, \pm}$ is the equivalent for $v_\mathrm{b}$ in $y$-direction. In the computation of $c_{\mathrm{g} z, r, \pm}$, the intrinsic frequency and squared buoyancy frequency are interpolated to $z_{r, \pm} = z_r \pm \Delta z_r / 2$ (also from before the position update). The update of the spectral ray-volume extents uses the fact that the surfaces in the $x$-$k$, $y$-$l$ and $z$-$m$ subspaces are conserved. Finally, the update of the phase-space wave-action density reads

\[\mathcal{N}_r \rightarrow \left(1 + 2 \alpha_{\mathrm{R}, r} f_\mathrm{RK} \Delta t\right)^{- 1} \mathcal{N}_r,\]

where $\alpha_{\mathrm{R}, r}$ is the interpolation of the Rayleigh-damping coefficient to the updated ray-volume position, obtained from interpolate_sponge.

The group velocities that are calculated for the propagation in physical space are also used to determine the maxima needed for the WKB-CFL condition used in the time-step computation.

propagate_rays!(
    state::State,
    dt::AbstractFloat,
    rkstage::Integer,
    wkb_mode::SteadyState,
)

Update the vertical wavenumber and wave-action density, using steady-state WKB theory.

In steady-state mode, the ray volumes are stationary in physical space. In mountain-wave simulations, this method first updates the ray volumes in the launch layer by calling activate_orographic_source!. Subsequently, it performs a vertical sweep to update all other ray volumes. Therein, the vertical wavenumber is set to

\[m_r = - \sigma \sqrt{\frac{\left(k_r^2 + l_r^2\right) \left(N_r^2 - \widehat{\omega}_r^2\right)}{\widehat{\omega}_r^2 - f^2}},\]

where $N_r^2$ is the squared buoyancy frequency interpolated to the ray-volume position (with interpolate_stratification) and $\widehat{\omega}_r = - k_r u_\mathrm{b} - l_r v_\mathrm{b}$ (in the case of mountain waves, for which $\omega_r = 0$). The new wave-action-density field is obtained by integrating

\[\frac{\partial}{\partial z} \left(c_{\mathrm{g} z, r} \mathcal{A}_r\right) = - 2 \alpha_{\mathrm{R}, r} \mathcal{A}_r - 2 K \left|\boldsymbol{k}_r\right|^2 \mathcal{A}_r,\]

where $\alpha_{\mathrm{R}, r}$ is the Rayleigh-damping coefficient interpolated to the ray-volume position (using interpolate_sponge) and

\[K = \left[2 \sum\limits_r \frac{J \Delta \widehat{z}}{c_{\mathrm{g}, z, r}} \left(m_r \left|b_{\mathrm{w}, r}\right| \left|\boldsymbol{k}_r\right|\right)^2 f_r\right]^{- 1} \max \left[0, \sum_r \left(m_r \left|b_{\mathrm{w}, r}\right|\right)^2 f_r - \alpha_\mathrm{s}^2 N^4\right]\]

is the turbulent viscosity and diffusivity due to wave breaking (see PinCFlow.MSGWaM.RayUpdate.apply_saturation_scheme! for more details). After the first right-hand-side term has been integrated with the implicit step

\[\mathcal{A}_r = \left[1 + \frac{2 \alpha_{\mathrm{R}, r}}{c_{\mathrm{g} z, r}} \left(z_r - z_{r, k - 1}\right)\right]^{- 1} \frac{c_{\mathrm{g} z, r, k - 1}}{c_{\mathrm{g} z, r}} \mathcal{A}_{r, k - 1},\]

the second term is integrated with the pseudo-time step $J \Delta \widehat{z} / c_{\mathrm{g} z, r}$, which corresponds to the substitution $\mathcal{A}_r \rightarrow \left(1 - 2 J \Delta \widehat{z} / c_{\mathrm{g} z, r} K \left|\boldsymbol{k}_r\right|^2\right) \mathcal{A}_r$.

If the domain is parallelized in the vertical, the integration in vertical subdomains is performed sequentially, with one-way communication providing boundary conditions.

Arguments

• state: Model state.

• dt: Time step.

• rkstage: Runge-Kutta-stage index.

• wkb_mode: Approximations used by MSGWaM.

See also

• PinCFlow.MSGWaM.RayOperations.get_physical_position

• PinCFlow.MSGWaM.RayOperations.get_spectral_position

• PinCFlow.MSGWaM.RayOperations.get_physical_extent

• PinCFlow.MSGWaM.RayOperations.get_spectral_extent

• PinCFlow.MSGWaM.Interpolation.interpolate_stratification

• PinCFlow.MSGWaM.Interpolation.interpolate_mean_flow

• PinCFlow.MSGWaM.Interpolation.interpolate_sponge

• PinCFlow.MSGWaM.RaySources.activate_orographic_source!

• PinCFlow.MSGWaM.RayOperations.copy_rays!

shift_rays!(state::State)

Shift the array positions of ray volumes such that they are attributed to the correct grid cells by dispatching to a WKB-mode-specific method.

shift_rays!(state::State, wkb_mode::Union{NoWKB, SteadyState})

Return for configurations without WKB / with steady-state WKB.

shift_rays!(state::State, wkb_mode::SingleColumn)

Shift the vertical array positions of ray volumes such that they are attributed to the correct grid cells.

This method enforces the vertical boundary conditions (via set_vertical_boundary_rays!), checks if ray volumes need to be shifted and, if they do, copies them to the correct grid cells and marks them for removal (by dispatching to the appropriate method). A second call of set_vertical_boundary_rays! ensures that ray volumes that have moved across MPI processes are included in the appropriate halo cells. Finally, the gaps that were created by marking ray volumes for removal are filled (via remove_rays!).

shift_rays!(state::State, wkb_mode::MultiColumn)

Shift the array positions of ray volumes such that they are attributed to the correct grid cells.

For each dimension in physical space (with more than one grid point), this method performs the corresponding equivalent of the algorithm that is implemented in the method for single-column mode.

shift_rays!(state::State, direction::X)

For each ray volume, check if it is attributed to the correct position in $\widehat{x}$ and, if it is not, create a copy that is and mark the original for removal.

Ray volumes that should be attributed to a halo cell are marked for removal but not copied, since the copies are created from the corresponding halo cell in the adjacent MPI process.

shift_rays!(state::State, direction::Y)

For each ray volume, check if it is attributed to the correct position in $\widehat{y}$ and, if it is not, create a copy that is and mark the original for removal.

Ray volumes in halo cells are treated in the same way as in the method for shifting in $\widehat{x}$.

shift_rays!(state::State, direction::Z)

For each ray volume, check if it is attributed to the correct position in $\widehat{z}$ and, if it is not, create a copy that is and mark the original for removal.

Ray volumes in halo cells are treated in the same way as in the methods for shifting in $\widehat{x}$ and $\widehat{z}$.

Arguments

• state: Model state.

• wkb_mode: Approximations used by MSGWaM.

• direction: Shift direction.

See also

• PinCFlow.MSGWaM.BoundaryRays.set_zonal_boundary_rays!

• PinCFlow.MSGWaM.BoundaryRays.set_meridional_boundary_rays!

• PinCFlow.MSGWaM.BoundaryRays.set_vertical_boundary_rays!

• PinCFlow.MSGWaM.RayOperations.remove_rays!

• PinCFlow.MSGWaM.RayOperations.check_rays

• PinCFlow.MSGWaM.RayOperations.copy_rays!

• PinCFlow.MSGWaM.Interpolation.get_next_half_level

split_rays!(state::State)

Split ray volumes that have become larger than the local grid cell by dispatching to a WKB-mode-specific method.

split_rays!(state::State, wkb_mode::Union{NoWKB, SteadyState})

Return for configurations without WKB / with steady-state WKB.

split_rays!(state::State, wkb_mode::SingleColumn)

Split ray volumes which have a vertical extent larger than the local vertical grid spacing.

split_rays!(state::State, wkb_mode::MultiColumn)

In each dimension of physical space, split ray volumes which have an extent larger than the local grid spacing.

The splitting is performed sequentially, such that a ray volume with extents that are all between once and twice as large as allowed is split into exactly eight smaller ray volumes (all of which have the same size).

split_rays!(i::Integer, j::Integer, k::Integer, state::State, axis::X)

In the grid cell specified by $\left(i, j, k\right)$, split ray volumes with $\Delta x_r > \Delta \widehat{x}$.

The number of splits is the result of ceiling division of $\Delta x_r$ by $\Delta \widehat{x}$. Each split is carried out by adjusting the position and extent of the ray volume, copying it and changing the position of the copy appropriately.

split_rays!(i::Integer, j::Integer, k::Integer, state::State, axis::Y)

In the grid cell specified by $\left(i, j, k\right)$, split ray volumes with $\Delta y_r > \Delta \widehat{y}$.

The splitting is analogous to that in $\widehat{x}$.

split_rays!(i::Integer, j::Integer, k::Integer, state::State, axis::Z)

In the grid cell specified by $\left(i, j, k\right)$, split ray volumes with $\Delta z_r > J_{\min} \Delta \widehat{z}$, with $J_{\min}$ being the minimum value of the Jacobian in all grid cells that are at least partially covered by the ray volume (at its true horizontal position on the grid).

The splitting is analogous to that in $\widehat{x}$ and $\widehat{y}$.

Arguments

• state: Model state.

• wkb_mode: Approximations used by MSGWaM.

• i: Grid-cell index in $\widehat{x}$-direction

• j: Grid-cell index in $\widehat{y}$-direction

• k: Grid-cell index in $\widehat{z}$-direction

• axis: Axis perpendicular to the split.

See also

• PinCFlow.MSGWaM.RayOperations.copy_rays!

• PinCFlow.MSGWaM.Interpolation.get_next_half_level

MeanFlowEffect

Module for computing the mean-flow effect of gravity waves.

Provides functions that compute mean-flow tendencies by integrating ray-volume properties in spectral space and mapping the result to physical grid cells. Also provides two filters for smoothing the tendencies, as well as a simple blocked-layer scheme that includes a blocked-flow drag in mountain-wave simulations.

See also

• PinCFlow.Types

• PinCFlow.Boundaries

• PinCFlow.MSGWaM.Interpolation

• PinCFlow.MSGWaM.RayUpdate

UChi

Singleton for dispatch to calculation of zonal gravity-wave-tracer fluxes.

VChi

Singleton for dispatch to calculation of meridional gravity-wave-tracer fluxes.

WChi

Singleton for dispatch to calculation of vertical gravity-wave-tracer fluxes.

apply_blocked_layer_scheme!(state::State)

Compute the blocked-flow drag and adjust the mean-flow impact accordingly.

The blocked-flow drag is given by

\[\left(\frac{\partial \boldsymbol{u}_\mathrm{b}}{\partial t}\right)_\mathrm{B} = - D \frac{\left|\boldsymbol{k}_h\right| \left|\boldsymbol{u}_\mathrm{p}\right|}{2 \pi} \boldsymbol{u}_\mathrm{p},\]

where $\boldsymbol{u}_\mathrm{b} = \left(u_\mathrm{b}, v_\mathrm{b}, 0\right)^\mathrm{T}$ is the resolved horizontal wind, $D$ is a dimensionless drag coefficient (represented by state.namelists.wkb.drag_coefficient),

\[\boldsymbol{k}_h = \frac{\sum_\alpha \left|h_{\mathrm{w}, \alpha}\right| \boldsymbol{k}_{h, \alpha}}{\sum_\alpha \left|h_{\mathrm{w}, \alpha}\right|}\]

is a weighted average of the orography's wavevectors $\boldsymbol{k}_{h, \alpha}$ (with $h_{\mathrm{w}, \alpha}$ being the corresponding spectral modes) and

\[\boldsymbol{u}_\mathrm{p} = \frac{\left(\boldsymbol{u}_\mathrm{b} \cdot \boldsymbol{k}_h\right) \boldsymbol{k}_h}{\left|\boldsymbol{k}_h\right|^2}\]

is the projection of $\boldsymbol{u}_\mathrm{b}$ onto $\boldsymbol{k}_h$. This drag replaces the drag due to gravity waves below $z_\mathrm{B}$, the upper edge of the blocked layer that has been determined by PinCFlow.MSGWaM.RaySources.activate_orographic_source!. In grid cells that contain this upper edge, blocking and gravity waves both contribute to the total drag, weighted by the corresponding grid-cell fractions. The gravity-wave heating is treated similarly, with the "blocking contribution" being zero.

Arguments

• state: Model state.

apply_shapiro_filter!(
    output::AbstractVector{<:AbstractFloat},
    input::AbstractVector{<:AbstractFloat},
    scope::UnitRange{<:Integer},
    order::Val{1},
)

Apply the first-order Shapiro filter to input.

The elements of output are given by

\[\widetilde{\phi}_i = \frac{1}{4} \left(\phi_{i - 1} + 2 \phi_i + \phi_{i + 1}\right),\]

where $\phi_i$ are the elements of input.

apply_shapiro_filter!(
    output::AbstractVector{<:AbstractFloat},
    input::AbstractVector{<:AbstractFloat},
    scope::UnitRange{<:Integer},
    order::Val{2},
)

Apply the second-order Shapiro filter to input.

The elements of output are given by

\[\widetilde{\phi}_i = \frac{1}{16} \left(- \phi_{i - 2} + 4 \phi_{i - 1} + 10 \phi_i + 4 \phi_{i + 1} - \phi_{i + 2}\right),\]

where $\phi_i$ are the elements of input.

apply_shapiro_filter!(
    output::AbstractVector{<:AbstractFloat},
    input::AbstractVector{<:AbstractFloat},
    scope::UnitRange{<:Integer},
    order::Val{3},
)

Apply the third-order Shapiro filter to input.

The elements of output are given by

\[\widetilde{\phi}_i = \frac{1}{64} \left(\phi_{i - 3} - 6 \phi_{i - 2} + 15 \phi_{i - 1} + 44 \phi_i + 15 \phi_{i + 1} - 6 \phi_{i + 2} + \phi_{i + 3}\right),\]

where $\phi_i$ are the elements of input.

apply_shapiro_filter!(
    output::AbstractVector{<:AbstractFloat},
    input::AbstractVector{<:AbstractFloat},
    scope::UnitRange{<:Integer},
    order::Val{4},
)

Apply the fourth-order Shapiro filter to input.

The elements of output are given by

\[\widetilde{\phi}_i = \frac{1}{256} \left(- \phi_{i - 4} + 8 \phi_{i - 3} - 28 \phi_{i - 2} + 56 \phi_{i - 1} + 186 \phi_i + 56 \phi_{i + 1} - 28 \phi_{i + 2} + 8 \phi_{i + 3} - \phi_{i + 4}\right),\]

where $\phi_i$ are the elements of input.

Arguments

• output: Filtered output vector.

• input: Input vector.

• scope: Index range.

• order: Order of the Shapiro filter.

compute_gw_integrals!(state::State)

Compute the gravity-wave integrals needed for the computation of the mean-flow impact by dispatching to a WKB-mode-specific method.

compute_gw_integrals!(state::State, wkb_mode::MultiColumn)

Compute the gravity-wave integrals needed for the computation of the mean-flow impact in multi-column mode.

This method computes the sums

\[\begin{align*} \overline{\rho} \left\langle u' u' \right\rangle & = \overline{\rho} \sum_{r, \lambda, \mu, \nu} \left(F u_\mathrm{w} u_\mathrm{w}^*\right)_{r, i + \lambda, j + \mu, k + \nu},\\ \overline{\rho} \left\langle u' v' \right\rangle & = \overline{\rho} \sum_{r, \lambda, \mu, \nu} \left(F u_\mathrm{w} v_\mathrm{w}^*\right)_{r, i + \lambda, j + \mu, k + \nu},\\ \overline{\rho} \left\langle u' w' \right\rangle & = \overline{\rho} \sum_{r, \lambda, \mu, \nu} \left(F u_\mathrm{w} w_\mathrm{w}^*\right)_{r, i + \lambda, j + \mu, k + \nu},\\ \overline{\rho} \left\langle v' v' \right\rangle & = \overline{\rho} \sum_{r, \lambda, \mu, \nu} \left(F v_\mathrm{w} v_\mathrm{w}^*\right)_{r, i + \lambda, j + \mu, k + \nu},\\ \overline{\rho} \left\langle v' w' \right\rangle & = \overline{\rho} \sum_{r, \lambda, \mu, \nu} \left(F v_\mathrm{w} w_\mathrm{w}^*\right)_{r, i + \lambda, j + \mu, k + \nu},\\ \left\langle \theta' u' \right\rangle & = \sum_{r, \lambda, \mu, \nu} \left(F \theta_\mathrm{w} u_\mathrm{w}^*\right)_{r, i + \lambda, j + \mu, k + \nu},\\ \left\langle \theta' v' \right\rangle & = \sum_{r, \lambda, \mu, \nu} \left(F \theta_\mathrm{w} v_\mathrm{w}^*\right)_{r, i + \lambda, j + \mu, k + \nu},\\ \mathcal{E} & = \sum_{r, \lambda, \mu, \nu} \left(F \mathcal{A} \widehat{\omega}\right)_{r, i + \lambda, j + \mu, k + \nu}. \end{align*}\]

Therein, $\left(\lambda, \mu, \nu\right)$ are index shifts to ray volumes that are at least partially within the grid cell at $\left(i, j, k\right)$, $F_{r, i + \lambda, j + \mu, k + \nu}$ are the corresponding ray volume fractions and $\left(u_\mathrm{w}, v_\mathrm{w}, w_\mathrm{w}, \theta_\mathrm{w}\right)_{r, i + \lambda, j + \mu, k + \nu}$ are the wave amplitudes of the wind and the potential temperature. The computation is based on the relations

\[\begin{align*} \overline{\rho} u_{\mathrm{w}, r} u_{\mathrm{w}, r}^* & = \left(k_r \widehat{c}_{\mathrm{g} x, r} - \mathrm{sgn} \left(\left|f\right|\right) \frac{k_r \widehat{c}_{\mathrm{g} x, r} + l_r \widehat{c}_{\mathrm{g} y, r}}{1 - \left(\widehat{\omega}_r / f\right)^2}\right) \mathcal{A}_r,\\ \overline{\rho} u_{\mathrm{w}, r} v_{\mathrm{w}, r}^* & = l_r \widehat{c}_{\mathrm{g} x, r} \mathcal{A}_r,\\ \overline{\rho} u_{\mathrm{w}, r} w_{\mathrm{w}, r}^* & = \frac{k_r \widehat{c}_{\mathrm{g} z, r}}{1 - \left(f / \widehat{\omega}_r\right)^2} \mathcal{A}_r,\\ \overline{\rho} v_{\mathrm{w}, r} v_{\mathrm{w}, r}^* & = \left(l_r \widehat{c}_{\mathrm{g} y, r} - \mathrm{sgn} \left(\left|f\right|\right) \frac{k_r \widehat{c}_{\mathrm{g} x, r} + l_r \widehat{c}_{\mathrm{g} y, r}}{1 - \left(\widehat{\omega}_r / f\right)^2}\right) \mathcal{A}_r,\\ \overline{\rho} v_{\mathrm{w}, r} w_{\mathrm{w}, r}^* & = \frac{l_r \widehat{c}_{\mathrm{g} z, r}}{1 - \left(f / \widehat{\omega}_r\right)^2} \mathcal{A}_r,\\ \theta_{\mathrm{w}, r} u_{\mathrm{w}, r}^* & = \frac{f \overline{\theta}}{g \overline{\rho}} \frac{l_r m_r N_r^2}{\widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2} \mathcal{A}_r,\\ \theta_{\mathrm{w}, r} v_{\mathrm{w}, r}^* & = - \frac{f \overline{\theta}}{g \overline{\rho}} \frac{k_r m_r N_r^2}{\widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2} \mathcal{A}_r, \end{align*}\]

where $N_r^2$ is the squared buoyancy frequency interpolated to the ray-volume position. The components of the intrinsic group velocity are given by

\[\begin{align*} \widehat{c}_{\mathrm{g} x, r} & = \frac{k_r \left(N_r^2 - \widehat{\omega}_r^2\right)}{\widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2},\\ \widehat{c}_{\mathrm{g} y, r} & = \frac{l_r \left(N_r^2 - \widehat{\omega}_r^2\right)}{\widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2},\\ \widehat{c}_{\mathrm{g} z, r} & = - \frac{m_r \left(\widehat{\omega}_r^2 - f^2\right)}{\widehat{\omega}_r \left|\boldsymbol{k}_r\right|^2}. \end{align*}\]

compute_gw_integrals!(state::State, wkb_mode::SingleColumn)

Compute the gravity-wave integrals needed for the computation of the mean-flow impact in single-column mode.

This method computes $\overline{\rho} \left\langle u' w' \right\rangle$, $\overline{\rho} \left\langle v' w' \right\rangle$, $\left\langle \theta' u' \right\rangle$, $\left\langle \theta' v' \right\rangle$ and $\mathcal{E}$ (see above for details).

compute_gw_integrals!(state::State, wkb_mode::SteadyState)

Compute the gravity-wave integrals needed for the computation of the mean-flow impact in steady-state mode.

This method computes the sums $\overline{\rho} \left\langle u' w' \right\rangle$ and $\overline{\rho} \left\langle v' w' \right\rangle$ (see above for details). In contrast to the multi-column and single-column modes, the steady-state mode uses the pseudo-momentum approximation

\[\begin{align*} \overline{\rho} u_{\mathrm{w}, r} w_{\mathrm{w}, r}^* & = k_r \widehat{c}_{\mathrm{g} z, r} \mathcal{A}_r,\\ \overline{\rho} v_{\mathrm{w}, r} w_{\mathrm{w}, r}^* & = l_r \widehat{c}_{\mathrm{g} z, r} \mathcal{A}_r. \end{align*}\]

Arguments

• state::State: Model state.

• wkb_mode: Approximations used by MSGWaM.

See also

• PinCFlow.MSGWaM.Interpolation.interpolate_stratification

• PinCFlow.MSGWaM.MeanFlowEffect.compute_horizontal_cell_indices

• PinCFlow.MSGWaM.Interpolation.get_next_half_level

compute_gw_tendencies!(state::State)

Compute the gravity-wave impact on the momentum and mass-weighted potential temperature.

Calculates the tendencies that are to be added to the equations for momentum and mass-weighted potential temperature. These are given by

\[\begin{align*} \left(\frac{\partial \rho_\mathrm{b} u_\mathrm{b}}{\partial t}\right)_\mathrm{w} & = - \frac{\rho_\mathrm{b}}{\overline{\rho}} \left[\frac{\left(\overline{\rho} \left\langle u' u' \right\rangle\right)_{i + 1} - \left(\overline{\rho} \left\langle u' u' \right\rangle\right)_{i - 1}}{2 \Delta \widehat{x}} + G^{13} \frac{\left(\overline{\rho} \left\langle u' u' \right\rangle\right)_{k + 1} - \left(\overline{\rho} \left\langle u' u' \right\rangle\right)_{k - 1}}{2 \Delta \widehat{z}}\right.\\ & \qquad \qquad + \frac{\left(\overline{\rho} \left\langle u' v' \right\rangle\right)_{j + 1} - \left(\overline{\rho} \left\langle u' v' \right\rangle\right)_{j - 1}}{2 \Delta \widehat{y}} + G^{23} \frac{\left(\overline{\rho} \left\langle u' v' \right\rangle\right)_{k + 1} - \left(\overline{\rho} \left\langle u' v' \right\rangle\right)_{k - 1}}{2 \Delta \widehat{z}}\\ & \qquad \qquad + \left.\frac{\left(\overline{\rho} \left\langle u' w' \right\rangle\right)_{k + 1} - \left(\overline{\rho} \left\langle u' w' \right\rangle\right)_{k - 1}}{2 J \Delta \widehat{z}}\right] - \rho_\mathrm{b} \frac{f}{\overline{\theta}} \left\langle \theta' v' \right\rangle,\\ \left(\frac{\partial \rho_\mathrm{b} v_\mathrm{b}}{\partial t}\right)_\mathrm{w} & = - \frac{\rho_\mathrm{b}}{\overline{\rho}} \left[\frac{\left(\overline{\rho} \left\langle v' u' \right\rangle\right)_{i + 1} - \left(\overline{\rho} \left\langle v' u' \right\rangle\right)_{i - 1}}{2 \Delta \widehat{x}} + G^{13} \frac{\left(\overline{\rho} \left\langle v' u' \right\rangle\right)_{k + 1} - \left(\overline{\rho} \left\langle v' u' \right\rangle\right)_{k - 1}}{2 \Delta \widehat{z}}\right.\\ & \qquad \qquad + \frac{\left(\overline{\rho} \left\langle v' v' \right\rangle\right)_{j + 1} - \left(\overline{\rho} \left\langle v' v' \right\rangle\right)_{j - 1}}{2 \Delta \widehat{y}} + G^{23} \frac{\left(\overline{\rho} \left\langle v' v' \right\rangle\right)_{k + 1} - \left(\overline{\rho} \left\langle v' v' \right\rangle\right)_{k - 1}}{2 \Delta \widehat{z}}\\ & \qquad \qquad + \left.\frac{\left(\overline{\rho} \left\langle v' w' \right\rangle\right)_{k + 1} - \left(\overline{\rho} \left\langle v' w' \right\rangle\right)_{k - 1}}{2 J \Delta \widehat{z}}\right] + \rho_\mathrm{b} \frac{f}{\overline{\theta}} \left\langle \theta' u' \right\rangle,\\ \left(\frac{\partial \rho_\mathrm{b} \widehat{w}_\mathrm{b}}{\partial t}\right)_\mathrm{w} & = G^{13} \left(\frac{\partial \rho_\mathrm{b} u_\mathrm{b}}{\partial t}\right)_\mathrm{w} + G^{23} \left(\frac{\partial \rho_\mathrm{b} v_\mathrm{b}}{\partial t}\right)_\mathrm{w},\\ \left(\frac{\partial P_\mathrm{b}}{\partial t}\right)_\mathrm{w} & = - \rho_\mathrm{b} \left[\frac{\left(\overline{\rho} \left\langle \theta' u' \right\rangle\right)_{i + 1} - \left(\overline{\rho} \left\langle \theta' u' \right\rangle\right)_{i - 1}}{2 \Delta \widehat{x}} + G^{13} \frac{\left(\overline{\rho} \left\langle \theta' u' \right\rangle\right)_{k + 1} - \left(\overline{\rho} \left\langle \theta' u' \right\rangle\right)_{k - 1}}{2 \Delta \widehat{z}}\right.\\ & \qquad \qquad + \left.\frac{\left(\overline{\rho} \left\langle \theta' v' \right\rangle\right)_{j + 1} - \left(\overline{\rho} \left\langle \theta' v' \right\rangle\right)_{j - 1}}{2 \Delta \widehat{y}} + G^{23} \frac{\left(\overline{\rho} \left\langle \theta' v' \right\rangle\right)_{k + 1} - \left(\overline{\rho} \left\langle \theta' v' \right\rangle\right)_{k - 1}}{2 \Delta \widehat{z}}\right], \end{align*}\]

where $\left(u_\mathrm{b}, v_\mathrm{b}, \widehat{w}_\mathrm{b}\right)$ are the components of the transformed (i.e. terrain-following) resolved wind, $\rho_\mathrm{b}$ is the resolved density (including the reference part $\overline{\rho}$) and $P_\mathrm{b}$ is the resolved mass-weighted potential temperature. For a documentation of the fluxes, see PinCFlow.MSGWaM.MeanFlowEffect.compute_gw_integrals!. Below state.namelists.wkb.impact_altitude, all tendencies are set to zero.

Arguments

• state::State: Model state.

compute_horizontal_cell_indices(
    state::State,
    xr::AbstractFloat,
    yr::AbstractFloat,
    dxr::AbstractFloat,
    dyr::AbstractFloat,
)::NTuple{4, <:Integer}

From the given horizontal ray-volume position and extent, determine the indices of the grid cells that contain the ray-volume edges and return them (in the order left, right, backward and forward).

Arguments

• state: Model state.

• xr: Ray-volume position in $x$.

• yr: Ray-volume position in $y$.

• dxr: Ray-volume extent in $x$.

• dyr: Ray-volume extent in $y$.

compute_leading_order_tracer_fluxes!(
    state::State,
    tracer_setup::NoTracer,
    fc::AbstractFloat,
    omir::AbstractFloat,
    wnrk::AbstractFloat,
    wnrl::AbstractFloat,
    wnrm::AbstractFloat,
    wadr::AbstractFloat,
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    i::Integer,
    j::Integer,
    k::Integer,
)

Return for configurations without tracer transport.

compute_leading_order_tracer_fluxes!(
    state::State,
    tracer_setup::TracerOn,
    fc::AbstractFloat,
    omir::AbstractFloat,
    wnrk::AbstractFloat,
    wnrl::AbstractFloat,
    wnrm::AbstractFloat,
    wadr::AbstractFloat,
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    i::Integer,
    j::Integer,
    k::Integer,
)

Compute the leading-order gravity-wave tracer fluxes at $(i, j, k)$.

The zonal, meridional, and vertical fluxes are given by

\[\begin{align*} \overline{\rho}\left\langle u' \chi'\right\rangle & = \overline{\rho}\sum_{r, \lambda,\mu,\nu}\left(F u_\mathrm{w}\chi^*_\mathrm{w}\right)_{r, i + \lambda, j + \mu, k + \nu},\\ \overline{\rho}\left\langle v' \chi'\right\rangle & = \overline{\rho}\sum_{r, \lambda, \mu, \nu} \left(F v_\mathrm{w} \chi^*_\mathrm{w}\right)_{r, i + \lambda, j + \mu, k + \nu},\\ \overline{\rho}\left\langle w' \chi'\right\rangle & = \overline{\rho}\sum_{r, \lambda, \mu, \nu} \left(F w_\mathrm{w} \chi^*_\mathrm{w}\right)_{r, i + \lambda, j + \mu, k + \nu}. \end{align*}\]

Arguments:

• state: Model state.

• tracer_setup: General tracer-transport configuration.

• fc: Coriolis parameter.

• omir: Gravity-wave intrinsic frequency.

• wnrk: Zonal wavenumber.

• wnrl: Meridional wavenumber.

• wnrm: Vertical wavenumber.

• wadr: Phase-space wave-action density.

• xlc: Zonal location of the ray-volume.

• ylc: Meridional location of the ray-volume.

• zlc: Vertical location of the ray-volume.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• k: Vertical grid-cell index.

See also:

• PinCFlow.MSGWaM.MeanFlowEffect.leading_order_tracer_fluxes

compute_leading_order_tracer_forcing!(
    state::State,
    i::Integer,
    j::Integer,
    k::Integer,
    tracer_setup::TracerOn,
)

Compute and return the leading-order tracer forcing at $\left(i, j, k\right)$.

compute_leading_order_tracer_forcing!(
    state::State,
    i::Integer,
    j::Integer,
    k::Integer,
    tracer_setup::NoTracer,
)

Return for configurations without tracer transport.

Arguments

• state: Model state.

• i: Zonal grid-cell index.

• j: Meridional grid-cell index.

• k: Vertical grid-cell index.

• tracer_setup: General tracer-transport configuration.

compute_mean_flow_effect!(state::State)

Calculate the mean-flow impact of unresolved gravity waves by dispatching to a WKB-mode-specific method.

compute_mean_flow_effect!(state::State, wkb_mode::NoWKB)

Return for non-WKB configurations.

compute_mean_flow_effect!(
    state::State,
    wkb_mode::Union{SteadyState, SingleColumn, MultiColumn},
)

Calculate the mean-flow impact of unresolved gravity waves.

This method first computes several spectral integrals (using compute_gw_integrals!), most of which represent gravity-wave fluxes. After the boundary conditions for these have been enforced (using set_boundaries!), the corresponding tendencies are calculated (using compute_gw_tendencies!). These also have boundary conditions that need to be enforced (once again using set_boundaries!) before they are smoothed to remove small-scale features that may occur due to a coarse ray-volume distribution (using smooth_gw_tendencies!). Afterwards, if MSGWaM parameterizes mountain waves, the tendencies are adjusted to account for the formation of blocked layers (using apply_blocked_layer_scheme!), before the boundary conditions are enforced again.

Arguments

• state: Model state.

• wkb_mode: Approximations used by MSGWaM.

See also

• PinCFlow.MSGWaM.MeanFlowEffect.compute_gw_integrals!

• PinCFlow.Boundaries.set_boundaries!

• PinCFlow.MSGWaM.MeanFlowEffect.compute_gw_tendencies!

• PinCFlow.MSGWaM.MeanFlowEffect.smooth_gw_tendencies!

• PinCFlow.MSGWaM.MeanFlowEffect.apply_blocked_layer_scheme!

leading_order_tracer_fluxes(
    state::State,
    fc::AbstractFloat,
    omir::AbstractFloat,
    wnrk::AbstractFloat,
    wnrl::AbstractFloat,
    wnrm::AbstractFloat,
    wadr::AbstractFloat,
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    direction::UChi,
)::AbstractFloat

Compute and return the contribution of a ray volume located at (xlc,ylc,zlc) to the zonal leading-order gravity-wave tracer fluxes $\overline{\rho}\langle u'\chi'\rangle$.

The flux-contributions are given by

\[\overline{\rho} u_{\mathrm{w}, r}\chi^*_{\mathrm{w}, r} = \frac{f}{\widehat{\omega}_r} \frac{m_r}{\left|\boldsymbol{k}_r\right|^2} \mathcal{A}_r \left[l_r \left(\frac{\partial \chi_\mathrm{b}}{\partial z}\right)_r - m_r \left(\frac{\partial \chi_\mathrm{b}}{\partial y}\right)_r\right].\]

leading_order_tracer_fluxes(
    state::State,
    fc::AbstractFloat,
    omir::AbstractFloat,
    wnrk::AbstractFloat,
    wnrl::AbstractFloat,
    wnrm::AbstractFloat,
    wadr::AbstractFloat,
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    direction::VChi,
)::AbstractFloat

Compute and return the contribution of a ray volume located at (xlc,ylc,zlc) to the meridional leading-order gravity-wave tracer fluxes $\overline{\rho}\langle v'\chi'\rangle$.

The fluxes are given by

\[\overline{\rho} v_{\mathrm{w}, r} \chi^*_{\mathrm{w}, r} = \frac{f}{\widehat{\omega}_{r}} \frac{m_r}{\left|\boldsymbol{k}_{r}\right|^2} \mathcal{A}_r \left[m_r \left(\frac{\partial \chi_\mathrm{b}}{\partial x}\right)_r - k_r \left(\frac{\partial \chi_\mathrm{b}}{\partial z}\right)_r\right].\]

leading_order_tracer_fluxes(
    state::State,
    fc::AbstractFloat,
    omir::AbstractFloat,
    wnrk::AbstractFloat,
    wnrl::AbstractFloat,
    wnrm::AbstractFloat,
    wadr::AbstractFloat,
    xlc::AbstractFloat,
    ylc::AbstractFloat,
    zlc::AbstractFloat,
    direction::WChi,
)::AbstractFloat

Compute and return the contribution of a ray volume located at (xlc,ylc,zlc) to the vertical leading-order gravity-wave tracer fluxes $\overline{\rho}\langle w'\chi'\rangle$.

The fluxes are given by

\[\overline{\rho} w_{\mathrm{w}, r} \chi^*_{\mathrm{w}, r} = \frac{f}{\widehat{\omega}_r} \frac{m_r}{\left|\boldsymbol{k}_r\right|^2} \mathcal{A}_r \left[k_r \left(\frac{\partial \chi_\mathrm{b}}{\partial y}\right)_r - l_r \left(\frac{\partial \chi_\mathrm{b}}{\partial x}\right)_r\right].\]

set_tracer_field_zero!(state)

Reset the gravity-wave-induced tracer fluxes to zero by dispatching over tracer configurations.

set_tracer_field_zero!(state::State, tracer_setup::NoTracer)

Return for configurations without tracer transport.

set_tracer_field_zero!(state::State, tracer_setup::TracerOn)

Set the gravity-wave-induced tracer fluxes to zero.

Arguments:

• state: Model state.

• tracer_setup: General tracer-transport configuration.

smooth_gw_tendencies!(state::State)

Apply spatial smoothing to gravity-wave tendency fields by dispatching to a method specific for the chosen filter (state.namelists.wkb.filter_type) and dimensionality of the domain.

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Box,
    direction::XYZ,
)

Apply a 3D box filter to smooth in all spatial directions.

Applies the moving average

\[\widetilde{\phi}_{i, j, k} = \left(2 N_\mathrm{s} + 1\right)^{- 3} \sum\limits_{\lambda = i - N_\mathrm{s}}^{i + N_\mathrm{s}} \sum\limits_{\mu = j - N_\mathrm{s}}^{j + N_\mathrm{s}} \sum\limits_{\nu = k - N_\mathrm{s}}^{k + N_\mathrm{s}} \phi_{\lambda, \mu, \nu},\]

where $N_\mathrm{s}$ is the order of the filter (state.namelists.wkb.filter_order).

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Box,
    direction::XZ,
)

Apply a 2D box filter to smooth in $\widehat{x}$ and $\widehat{z}$.

Applies the moving average

\[\widetilde{\phi}_{i, j, k} = \left(2 N_\mathrm{s} + 1\right)^{- 2} \sum\limits_{\lambda = i - N_\mathrm{s}}^{i + N_\mathrm{s}} \sum\limits_{\nu = k - N_\mathrm{s}}^{k + N_\mathrm{s}} \phi_{\lambda, j, \nu},\]

where $N_\mathrm{s}$ is the order of the filter (state.namelists.wkb.filter_order).

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Box,
    direction::YZ,
)

Apply a 2D box filter to smooth in $\widehat{y}$ and $\widehat{z}$.

Applies the moving average

\[\widetilde{\phi}_{i, j, k} = \left(2 N_\mathrm{s} + 1\right)^{- 2} \sum\limits_{\mu = j - N_\mathrm{s}}^{j + N_\mathrm{s}} \sum\limits_{\nu = k - N_\mathrm{s}}^{k + N_\mathrm{s}} \phi_{i, \mu, \nu},\]

where $N_\mathrm{s}$ is the order of the filter (state.namelists.wkb.filter_order).

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Box,
    direction::Z,
)

Apply a 1D box filter to smooth in $\widehat{z}$.

Applies the moving average

\[\widetilde{\phi}_{i, j, k} = \left(2 N_\mathrm{s} + 1\right)^{- 1} \sum\limits_{\nu = k - N_\mathrm{s}}^{k + N_\mathrm{s}} \phi_{i, j, \nu},\]

where $N_\mathrm{s}$ is the order of the filter (state.namelists.wkb.filter_order).

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Shapiro,
    direction::XYZ,
)

Apply a 3D Shapiro filter to smooth in all spatial directions.

A 1D Shapiro filter is applied sequentially in $\widehat{x}$, $\widehat{y}$ and $\widehat{z}$.

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Shapiro,
    direction::XZ,
)

Apply a 2D Shapiro filter to smooth in $\widehat{x}$ and $\widehat{z}$.

A 1D Shapiro filter is applied sequentially in $\widehat{x}$ and $\widehat{z}$.

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Shapiro,
    direction::YZ,
)

Apply a 2D Shapiro filter to smooth in $\widehat{y}$ and $\widehat{z}$.

A 1D Shapiro filter is applied sequentially in $\widehat{y}$ and $\widehat{z}$.

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Shapiro,
    direction::Z,
)

Apply a 1D Shapiro filter to smooth in $\widehat{z}$.

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Shapiro,
    direction::Y,
)

Apply a 1D Shapiro filter to smooth in $\widehat{y}$.

smooth_gw_tendencies!(
    output::AbstractArray{<:AbstractFloat, 3},
    state::State,
    filter_type::Shapiro,
    direction::X,
)

Apply a 1D Shapiro filter to smooth in $\widehat{x}$.

smooth_gw_tendencies!(state::State, tracer_setup::TracerOn)

Apply smoothing to tracer tendencies.

smooth_gw_tendencies!(state::State, tracer_setup::NoTracer)

Return for configurations without tracer transport.

Arguments

• state: Model state.

• output: Field to smooth.

• filter_type: Filter type.

• direction: Directions to smooth in.

• tracer_setup: General tracer-transport configuration.

See also

• PinCFlow.MSGWaM.MeanFlowEffect.apply_shapiro_filter!