ERF_PBLModels.H File Reference

#include <functional>
#include <limits>
#include <memory>
#include <AMReX_MultiFab.H>
#include <AMReX_Periodicity.H>
#include "ERF_Constants.H"
#include "ERF_DataStruct.H"
#include "ERF_MOSTUtils.H"

Include dependency graph for ERF_PBLModels.H:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	PBLDerivativeDzInv
struct	PBLSurfaceLayerGradient
struct	PBLDerivativeDzInv_N
	Functor for inverse vertical spacings with constant grid spacing. More...
struct	PBLDerivativeDzInv_S
	Functor for inverse vertical spacings for stretched grids using a spacing array. More...
struct	PBLDerivativeDzInv_T
	Functor for inverse vertical spacings for terrain-following grids using cell-center heights. More...

Functions
void	ComputeDiffusivityMYJ (double dt, const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const std::unique_ptr< amrex::MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices)
void	ComputeDiffusivityMYNN25 (const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const std::unique_ptr< amrex::MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices)
void	ComputeDiffusivityMYNNEDMF (const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const std::unique_ptr< amrex::MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices)
void	ComputeDiffusivityYSU (const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const std::unique_ptr< amrex::MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices)
void	ComputeDiffusivityMRF (const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const std::unique_ptr< amrex::MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices, const amrex::MultiFab *terrain_blank=nullptr)
void	ComputeDiffusivityYSUNew (const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const std::unique_ptr< amrex::MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices, const amrex::MultiFab *qheating_rates=nullptr, const amrex::MultiFab *terrain_blank=nullptr)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	GetPBLStratVar (int i, int j, int k, const StratType strat, const amrex::Array4< const amrex::Real > &cell_data, const MoistureComponentIndices &moisture_indices)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void	ComputeVerticalDerivativesPBL (int i, int j, int k, const amrex::Array4< const amrex::Real > &uvel, const amrex::Array4< const amrex::Real > &vvel, const amrex::Array4< const amrex::Real > &cell_data, const int izmin, const int izmax, const PBLDerivativeDzInv &dz_inv, const bool c_ext_dir_on_zlo, const bool c_ext_dir_on_zhi, const bool u_ext_dir_on_zlo, const bool u_ext_dir_on_zhi, const bool v_ext_dir_on_zlo, const bool v_ext_dir_on_zhi, amrex::Real &dthetadz, amrex::Real &dudz, amrex::Real &dvdz, const MoistureComponentIndices &moisture_indices, const StratType strat=StratType::thetav)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	ComputeVirtualTStarPBL (amrex::Real t_star, amrex::Real q_star, amrex::Real theta, amrex::Real qv, bool use_moisture)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void	ApplySurfaceLayerGradientsPBL (const PBLSurfaceLayerGradient &sl, amrex::Real &dthetadz, amrex::Real &dudz, amrex::Real &dvdz)
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real	ComputeQKESourceTerms (int i, int j, int k, const amrex::Array4< const amrex::Real > &uvel, const amrex::Array4< const amrex::Real > &vvel, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< const amrex::Real > &K_turb, const amrex::Box &domain, const PBLDerivativeDzInv dz_inv, amrex::Real pbl_mynn_B1_l, const amrex::Real theta_mean, const MoistureComponentIndices &moisture_indices, bool c_ext_dir_on_zlo, bool c_ext_dir_on_zhi, bool u_ext_dir_on_zlo, bool u_ext_dir_on_zhi, bool v_ext_dir_on_zlo, bool v_ext_dir_on_zhi, const amrex::Array4< const amrex::Real > &tau13, const amrex::Array4< const amrex::Real > &tau23, const amrex::Array4< const amrex::Real > &hfx_z, const amrex::Array4< const amrex::Real > &qfx1_z, bool surface_layer_on_zlo)
void	ApplyPBLHSmoothing (amrex::FArrayBox &pblh_fab, const amrex::Box &xybx_valid, const amrex::Real weight, const int passes, const amrex::Box &domain, const amrex::Periodicity &periodicity)
	Apply spatial smoothing to PBLH field using 5-point stencil. More...
std::unique_ptr< amrex::MultiFab >	WidenForPBLH (const amrex::MultiFab &src, const int ng, const std::function< void(amrex::MultiFab &)> &fill)

Function Documentation

ApplyPBLHSmoothing()

void ApplyPBLHSmoothing ( amrex::FArrayBox &  pblh_fab, const amrex::Box &  xybx_valid, const amrex::Real  weight, const int  passes, const amrex::Box &  domain, const amrex::Periodicity &  periodicity  )

inline

Apply spatial smoothing to PBLH field using 5-point stencil.

Applies a simple 5-point spatial averaging to reduce grid-to-grid noise in the diagnosed PBLH field. This is standard post-processing practice for discrete Rib-crossing detection as described in Seibert et al. (2000): "Review and intercomparison of operational methods for the determination of the mixing height", Atmospheric Environment, 34, 1001-1027.

Stencil: PBLH_smooth(i,j) = w * PBLH(i,j) + (1-w)/4 * [PBLH(i+-1,j) + PBLH(i,j+-1)]

The stencil reads one column outside the region it writes, so the result on a cell depends on columns that a single tile – or a single box – does not own. The caller therefore names the region it wants smoothed (xybx_valid, normally the planar tile box) inside a work array carrying passes columns of halo around it, filled from the same inputs. Pass p writes the region grown by passes-1-p and reads the one grown by passes-p, so each pass consumes one column of halo and the last lands exactly on xybx_valid. Every read stays inside the array, and the result does not depend on how the domain is divided into boxes or tiles.

At a non-periodic domain edge there is no column to read, so the stencil reuses the cell's own value (reflective/Neumann). Across a periodic edge the halo holds the wrapped columns and is read directly, which is why the periodicity has to be passed in: clamping there would fold an edge into a domain that has none.

Parameters

[in,out]	pblh_fab	FArrayBox holding the 2D PBLH field (component 0), sized to contain xybx_valid grown by passes in x and y
[in]	xybx_valid	Region to smooth (x,y indices only)
[in]	weight	Center cell weight in stencil (must be in [0,1])
[in]	passes	Number of smoothing iterations
[in]	domain	Domain box, for the reflective edge condition
[in]	periodicity	Domain periodicity, so periodic edges read the halo

{
    if (passes <= 0) { return; }
 
    // The caller must supply the halo the passes will consume; see the note above.
    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
        pblh_fab.box().contains(amrex::grow(xybx_valid, amrex::IntVect(passes,passes,0))),
        "ApplyPBLHSmoothing: the PBLH work array needs `passes` columns of halo "
        "around the region being smoothed");
 
    // Use a temporary to hold intermediate results
    amrex::FArrayBox pblh_temp(pblh_fab.box(), 1, amrex::The_Async_Arena());
    auto pblh = pblh_fab.array();
    auto pblh_tmp = pblh_temp.array();
 
    const auto& dom_lo = amrex::lbound(domain);
    const auto& dom_hi = amrex::ubound(domain);
 
    const bool per_x = periodicity.isPeriodic(0);
    const bool per_y = periodicity.isPeriodic(1);
 
    const amrex::Real wt_side = (one - weight) / four; // Weight for each of 4 neighbors
 
    // Apply smoothing passes
    for (int pass = 0; pass < passes; ++pass) {
        const int halo = passes - 1 - pass;
        const amrex::Box bx = amrex::grow(xybx_valid, amrex::IntVect(halo,halo,0));
 
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
        {
            // Reflective (Neumann) at a non-periodic domain edge: reuse the cell's
            // own value. Everywhere else -- including across a periodic edge, where
            // the halo holds the wrapped column -- read the neighbour.
            const int i_xlo = (per_x || i > dom_lo.x) ? (i - 1) : i;
            const int i_xhi = (per_x || i < dom_hi.x) ? (i + 1) : i;
            const int j_ylo = (per_y || j > dom_lo.y) ? (j - 1) : j;
            const int j_yhi = (per_y || j < dom_hi.y) ? (j + 1) : j;
 
            // 5-point stencil
            pblh_tmp(i, j, 0) = weight * pblh(i, j, 0)
                              + wt_side * (pblh(i_xlo, j, 0) + pblh(i_xhi, j, 0) +
                                          pblh(i, j_ylo, 0) + pblh(i, j_yhi, 0));
        });
 
        // Copy back only what this pass computed: outside bx the array still holds
        // the previous pass, which no later pass reads.
        pblh_fab.copy<amrex::RunOn::Device>(pblh_temp, bx, 0, bx, 0, 1);
    }
}

Referenced by ComputeDiffusivityMRF(), and ComputeDiffusivityYSUNew().

Here is the call graph for this function:

Here is the caller graph for this function:

ApplySurfaceLayerGradientsPBL()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void ApplySurfaceLayerGradientsPBL	(	const PBLSurfaceLayerGradient &	sl,
		amrex::Real &	dthetadz,
		amrex::Real &	dudz,
		amrex::Real &	dvdz
	)

Overwrite the vertical gradients in the first cell above the surface with the MOST profile gradients, |dU/dz| = u_* phi_m(zeta) / (kappa z) dthetav/dz = thetav_* phi_h(zeta) / (kappa z)

With zlo.type = surface_layer the state carries foextrap and the velocities hoextrap at the lower boundary, so the centered stencil in ComputeVerticalDerivativesPBL reads ghost values that hold no surface information there: the potential-temperature gradient comes out half the one-sided value and the velocity gradient has no dependence on the friction velocity at all (ERF #4037). The MOST profile used here is the one the surface stress and heat flux were themselves derived from.

Every PBL scheme that consumes these gradients uses only the shear magnitude dudz^2 + dvdz^2, so that magnitude is carried in dudz and dvdz is left at zero rather than being split by wind direction.

Parameters

[in]	sl	Monin-Obukhov inputs for this column
[out]	dthetadz	vertical derivative of virtual potential temperature
[out]	dudz	vertical derivative of x-velocity (carries the shear magnitude)
[out]	dvdz	vertical derivative of y-velocity (set to zero)

{
    const similarity_funs sfuns{};
    const amrex::Real inv_kappa_z = one / (KAPPA * sl.zval);
    dudz     = sl.u_star  * sfuns.calc_phi_m(sl.zeta) * inv_kappa_z;
    dvdz     = zero;
    dthetadz = sl.tstar_v * sfuns.calc_phi_h(sl.zeta) * inv_kappa_z;
}

Referenced by ComputeDiffusivityMYJ(), ComputeDiffusivityMYNN25(), ComputeDiffusivityMYNNEDMF(), ComputeDiffusivityYSU(), and ComputeDiffusivityYSUNew().

Here is the caller graph for this function:

ComputeDiffusivityMRF()

void ComputeDiffusivityMRF	(	const amrex::MultiFab &	xvel,
		const amrex::MultiFab &	yvel,
		const amrex::MultiFab &	cons_in,
		amrex::MultiFab &	eddyViscosity,
		const amrex::Geometry &	geom,
		const TurbChoice &	turbChoice,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		bool	use_terrain_fitted_coords,
		bool	use_moisture,
		int	level,
		const amrex::BCRec *	bc_ptr,
		bool	,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_cc,
		const MoistureComponentIndices &	moisture_indices,
		const amrex::MultiFab *	terrain_blank = nullptr
	)

Compute eddy diffusivities of momentum (eddy viscosity) and heat using the Medium Range Forecast (MRF) boundary layer scheme based on Hong and Pan (1996). Reference: Hong, S. Y., and H.-L. Pan, 1996: Nonlocal Boundary Layer Vertical Diffusion in a Medium-Range Forecast Model. Monthly Weather Review, 124, 2322-2339. https://doi.org/10.1175/1520-0493(1996)124<2322:NBLVDI>2.0.CO;2

Parameters

[in]	xvel	Velocity in x-dir
[in]	yvel	Velocity in y-dir
[in]	cons_in	Cell center conserved quantities
[out]	eddyViscosity	Holds turbulent viscosity
[in]	geom	Problem geometry
[in]	turbChoice	Container with turbulence parameters
[in]	SurfLayer	Pointer to Monin-Obukhov class if instantiated
[in]	use_moisture	If we have microphysics enabled
[in]	level	Current level
[in]	bc_ptr	Pointer to array with boundary condition info
[in]	vert_only	Only compute vertical eddy diffusivities
[in]	z_phys_nd	Physical location of grid nodes, if terrain (or grid stretching) is enabled
[in]	z_phys_cc	Physical location of cell centers

ComputeDiffusivityMYJ()

void ComputeDiffusivityMYJ	(	double	dt,
		const amrex::MultiFab &	xvel,
		const amrex::MultiFab &	yvel,
		amrex::MultiFab &	cons_in,
		amrex::MultiFab &	eddyViscosity,
		const amrex::Geometry &	geom,
		const TurbChoice &	turbChoice,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		bool	use_terrain_fitted_coords,
		bool	use_moisture,
		int	level,
		const amrex::BCRec *	bc_ptr,
		bool	,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_cc,
		const MoistureComponentIndices &	moisture_indices
	)

Compute eddy diffusivities of momentum (eddy viscosity) and heat using the Mellor-Yamada-Janjić (MYJ) boundary layer scheme. Reference: Janjić, Z. I., 1994: The Step-Mountain Eta Coordinate Model: Further Developments of the Convection, Viscous Sublayer, and Turbulence Closure Schemes. Monthly Weather Review, 122, 927–945. https://doi.org/10.1175/1520-0493(1994)122<0927:TSMECM>2.0.CO;2

Parameters

[in]	xvel	Velocity in x-dir
[in]	yvel	Velocity in y-dir
[in]	cons_in	Cell center conserved quantities
[out]	eddyViscosity	Holds turbulent viscosity
[in]	geom	Problem geometry
[in]	turbChoice	Container with turbulence parameters
[in]	SurfLayer	Pointer to Monin-Obukhov class if instantiated
[in]	use_moisture	If we have microphysics enabled
[in]	level	Current level
[in]	bc_ptr	Pointer to array with boundary condition info
[in]	vert_only	Only compute vertical eddy diffusivities
[in]	z_phys_nd	Physical location of grid nodes, if terrain (or grid stretching) is enabled
[in]	z_phys_cc	Physical location of cell centers

ComputeDiffusivityMYNN25()

void ComputeDiffusivityMYNN25	(	const amrex::MultiFab &	xvel,
		const amrex::MultiFab &	yvel,
		const amrex::MultiFab &	cons_in,
		amrex::MultiFab &	eddyViscosity,
		const amrex::Geometry &	geom,
		const TurbChoice &	turbChoice,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		bool	use_terrain_fitted_coords,
		bool	use_moisture,
		int	level,
		const amrex::BCRec *	bc_ptr,
		bool	,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_cc,
		const MoistureComponentIndices &	moisture_indices
	)

Compute eddy diffusivities of momentum (eddy viscosity) and heat using the Mellor-Yamada-Nakanishi-Niino Level amrex::Real(2.5) PBL scheme

Parameters

[in]	xvel	Velocity in x-dir
[in]	yvel	Velocity in y-dir
[in]	cons_in	Cell center conserved quantities
[out]	eddyViscosity	Holds turbulent viscosity
[in]	geom	Problem geometry
[in]	turbChoice	Container with turbulence parameters
[in]	SurfLayer	Pointer to Monin-Obukhov class if instantiated
[in]	use_moisture	If we have microphysics enabled
[in]	level	Current level
[in]	bc_ptr	Pointer to array with boundary condition info
[in]	vert_only	Only compute vertical eddy diffusivities
[in]	z_phys_nd	Physical location of grid nodes, if terrain (or grid stretching) is enabled
[in]	z_phys_cc	Physical location of cell centers

ComputeDiffusivityMYNNEDMF()

void ComputeDiffusivityMYNNEDMF	(	const amrex::MultiFab &	xvel,
		const amrex::MultiFab &	yvel,
		const amrex::MultiFab &	cons_in,
		amrex::MultiFab &	eddyViscosity,
		const amrex::Geometry &	geom,
		const TurbChoice &	turbChoice,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		bool	use_terrain_fitted_coords,
		bool	use_moisture,
		int	level,
		const amrex::BCRec *	bc_ptr,
		bool	,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_cc,
		const MoistureComponentIndices &	moisture_indices
	)

Compute eddy diffusivities of momentum (eddy viscosity) and heat using the Mellor-Yamada-Nakanishi-Niino Level amrex::Real(2.5) PBL scheme

Parameters

[in]	xvel	Velocity in x-dir
[in]	yvel	Velocity in y-dir
[in]	cons_in	Cell center conserved quantities
[out]	eddyViscosity	Holds turbulent viscosity
[in]	geom	Problem geometry
[in]	turbChoice	Container with turbulence parameters
[in]	SurfLayer	Pointer to Monin-Obukhov class if instantiated
[in]	use_moisture	If we have microphysics enabled
[in]	level	Current level
[in]	bc_ptr	Pointer to array with boundary condition info
[in]	vert_only	Only compute vertical eddy diffusivities
[in]	z_phys_nd	Physical location of grid nodes, if terrain (or grid stretching) is enabled
[in]	z_phys_cc	Physical location of cell centers

ComputeDiffusivityYSU()

void ComputeDiffusivityYSU	(	const amrex::MultiFab &	xvel,
		const amrex::MultiFab &	yvel,
		const amrex::MultiFab &	cons_in,
		amrex::MultiFab &	eddyViscosity,
		const amrex::Geometry &	geom,
		const TurbChoice &	turbChoice,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		bool	use_terrain_fitted_coords,
		bool	use_moisture,
		int	level,
		const amrex::BCRec *	bc_ptr,
		bool	,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_cc,
		const MoistureComponentIndices &	moisture_indices
	)

Compute eddy diffusivities of momentum (eddy viscosity) and heat using the Yonsei University PBL scheme

Parameters

[in]	xvel	Velocity in x-dir
[in]	yvel	Velocity in y-dir
[in]	cons_in	Cell center conserved quantities
[out]	eddyViscosity	Holds turbulent viscosity
[in]	geom	Problem geometry
[in]	turbChoice	Container with turbulence parameters
[in]	SurfLayer	Pointer to Monin-Obukhov class if instantiated
[in]	use_moisture	If we have microphysics enabled
[in]	level	Current level
[in]	bc_ptr	Pointer to array with boundary condition info
[in]	vert_only	Only compute vertical eddy diffusivities
[in]	z_phys_nd	Physical location of grid nodes, if terrain (or grid stretching) is enabled
[in]	z_phys_cc	Physical location of cell centers

ComputeDiffusivityYSUNew()

void ComputeDiffusivityYSUNew	(	const amrex::MultiFab &	xvel,
		const amrex::MultiFab &	yvel,
		const amrex::MultiFab &	cons_in,
		amrex::MultiFab &	eddyViscosity,
		const amrex::Geometry &	geom,
		const TurbChoice &	turbChoice,
		std::unique_ptr< SurfaceLayer > &	SurfLayer,
		bool	use_terrain_fitted_coords,
		bool	use_moisture,
		int	level,
		const amrex::BCRec *	bc_ptr,
		bool	,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_nd,
		const std::unique_ptr< amrex::MultiFab > &	z_phys_cc,
		const MoistureComponentIndices &	moisture_indices,
		const amrex::MultiFab *	qheating_rates = nullptr,
		const amrex::MultiFab *	terrain_blank = nullptr
	)

Compute eddy diffusivities of momentum (eddy viscosity) and heat using the WRF Yonsei University (YSU) PBL scheme.

Based on: Hong, Noh & Dudhia (2006), MWR, https://doi.org/10.1175/MWR3250.1 Hong (2010), QJRMS, https://doi.org/10.1002/qj.665

Features: three-pass bulk-Richardson PBLH diagnosis, nonlocal countergradient fluxes, explicit entrainment at PBL top, and grid-adaptive Richardson-number mixing in free atmosphere.

Parameters

[in]	xvel	x-direction velocity
[in]	yvel	y-direction velocity
[in]	cons_in	cell-center conserved quantities
[out]	eddyViscosity	turbulent eddy viscosity
[in]	geom	problem geometry
[in]	turbChoice	turbulence parameter container
[in]	SurfLayer	Monin-Obukhov surface layer object
[in]	use_terrain_fitted_coords	flag for terrain coordinates
[in]	use_moisture	flag for moisture active
[in]	level	AMR level
[in]	bc_ptr	boundary condition array
[in]	vert_only	compute vertical diffusivity only
[in]	z_phys_nd	physical node heights (terrain)
[in]	z_phys_cc	physical cell-center heights
[in]	moisture_indices	component index mapping for moisture variables
[in]	qheating_rates	radiation heating rates (SW, LW components); optional

ComputeQKESourceTerms()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeQKESourceTerms	(	int	i,
		int	j,
		int	k,
		const amrex::Array4< const amrex::Real > &	uvel,
		const amrex::Array4< const amrex::Real > &	vvel,
		const amrex::Array4< const amrex::Real > &	cell_data,
		const amrex::Array4< const amrex::Real > &	cell_prim,
		const amrex::Array4< const amrex::Real > &	K_turb,
		const amrex::Box &	domain,
		const PBLDerivativeDzInv	dz_inv,
		amrex::Real	pbl_mynn_B1_l,
		const amrex::Real	theta_mean,
		const MoistureComponentIndices &	moisture_indices,
		bool	c_ext_dir_on_zlo,
		bool	c_ext_dir_on_zhi,
		bool	u_ext_dir_on_zlo,
		bool	u_ext_dir_on_zhi,
		bool	v_ext_dir_on_zlo,
		bool	v_ext_dir_on_zhi,
		const amrex::Array4< const amrex::Real > &	tau13,
		const amrex::Array4< const amrex::Real > &	tau23,
		const amrex::Array4< const amrex::Real > &	hfx_z,
		const amrex::Array4< const amrex::Real > &	qfx1_z,
		bool	surface_layer_on_zlo
	)

Function for computing the QKE source terms (NN09, Eqn. 5).

Parameters

[in]	u	velocity in x-dir
[in]	v	velocity in y-dir
[in]	cell_data	conserved cell center vars
[in]	cell_prim	primitive cell center vars
[in]	K_turb	turbulent viscosity
[in]	domain	box of the whole domain
[in]	dz_inv	inverse spacings for interior and Dirichlet-boundary differences
[in]	pbl_mynn_B1_l	a parameter
[in]	theta_mean	average theta
[in]	moisture_indices	moisture variable component indices
[in]	tau13	rho<u'w'> on z-faces; the surface-layer stress at the zlo face
[in]	tau23	rho<v'w'> on z-faces; the surface-layer stress at the zlo face
[in]	hfx_z	rho<w'theta'> on z-faces; the surface-layer heat flux at the zlo face
[in]	qfx1_z	rho<w'qv'> on z-faces; the surface-layer moisture flux at the zlo face
[in]	surface_layer_on_zlo	a SurfaceLayer supplies the fluxes through the zlo face

{
    amrex::Real source_term = zero;
 
    int izmin = domain.smallEnd(2);
    int izmax = domain.bigEnd(2);
 
    const amrex::Real rho = cell_data(i,j,k,Rho_comp);
 
    // Notes:
    // - We transport TKE = myhalf*QKE rather than QKE, so the RHS terms do not
    //   have a factor of two.
    // - Transport terms due to turbulence and pressure are included when
    //   DiffusionSrcForState_* is called from ERF_slow_rhs_post.
    // - Eddy diffusivities are updated at the beginning of each time step only.
 
    const bool use_surface_fluxes = ( surface_layer_on_zlo && (k == izmin) &&
                                      tau13 && tau23 && hfx_z );
 
    if (use_surface_fluxes)
    {
        // The first cell center lies inside the unresolved surface layer. With
        // zlo.type = surface_layer the state carries foextrap and the velocities
        // hoextrap at the lower boundary, so the stencil used below reads ghost
        // values that hold no surface information: the resolved gradients there
        // are not a discretization of anything (ERF #4037). Build the production
        // terms from the fluxes the model actually carries instead:
        //
        //   tau13, tau23  = rho<u'w'>, rho<v'w'>      (z-faces, staggered in x,y)
        //   hfx_z, qfx1_z = rho<w'theta'>, rho<w'qv'> (z-faces)
        //
        // Two different heights are wanted here, so the two are kept apart:
        //   - u_*, thetav_* and L are surface-layer parameters and have no
        //     cell-centered counterpart, so they come from the zlo face alone;
        //   - the production terms belong at the cell center, so their flux
        //     factors average the two z-faces bounding this cell.
        // Only the profile shape, phi(zeta_1)/(kappa z_1), ties the two together.
 
        const amrex::Real theta = cell_data(i,j,k,RhoTheta_comp) / rho;
        const amrex::Real qv    = (moisture_indices.qv >= 0) ?
                                  cell_data(i,j,k,moisture_indices.qv) / rho : zero;
 
        // ---- Surface-layer parameters, on the zlo face ----
 
        // |tau|/rho = u_*^2
        const amrex::Real t13_s = myhalf * ( tau13(i,j,k) + tau13(i+1,j,k) );
        const amrex::Real t23_s = myhalf * ( tau23(i,j,k) + tau23(i,j+1,k) );
        const amrex::Real tau_s = std::sqrt(t13_s*t13_s + t23_s*t23_s);
        const amrex::Real ustar2 = tau_s / rho;
        const amrex::Real ustar  = std::sqrt(ustar2);
 
        // Surface virtual heat flux rho<w'theta_v'> (Stull Eqn. 4.4.5d)
        amrex::Real hfxv_s = hfx_z(i,j,k) * (one + epsv*qv);
        if (qfx1_z) { hfxv_s += epsv * theta * qfx1_z(i,j,k); }
 
        // Height of this cell center above the surface. dz_inv.dirichlet_lo is by
        // construction the inverse of the boundary-to-first-center distance, so it
        // carries the right height for uniform, stretched and terrain-fitted grids.
        const amrex::Real z1 = one / dz_inv.dirichlet_lo;
 
        // Obukhov length, from the surface fluxes. A vanishing stress (a building
        // face, say) leaves the shear production zero whatever zeta is, so the
        // neutral value is the safe one to carry into phi_m there.
        const amrex::Real thetav = theta * (one + epsv*qv);
        const amrex::Real eps    = std::numeric_limits<amrex::Real>::epsilon();
        amrex::Real zeta = zero;
        if (std::abs(hfxv_s) > eps && ustar > eps) {
            const amrex::Real l_obukhov = -( thetav * rho * ustar2 * ustar )
                                        / ( KAPPA * CONST_GRAV * hfxv_s );
            zeta = z1 / l_obukhov;
        }
 
        // phi_m, the companion of the psi_m that SurfaceLayer uses
        const similarity_funs sfuns{};
        const amrex::Real phi_m = sfuns.calc_phi_m(zeta);
 
        // ---- Production terms, at the cell center ----
 
        // Cell-centered stress: average the two z-faces bounding this cell
        const amrex::Real t13_c = fourth * ( tau13(i,j,k  ) + tau13(i+1,j,k  )
                                           + tau13(i,j,k+1) + tau13(i+1,j,k+1) );
        const amrex::Real t23_c = fourth * ( tau23(i,j  ,k) + tau23(i,j  ,k+1)
                                           + tau23(i,j+1,k) + tau23(i,j+1,k+1) );
 
        // Cell-centered virtual heat flux. Both faces are from this stage: the
        // QKE sources are added after RhoQ1 is diffused (see AddTurbKESources).
        amrex::Real hfxv_c = myhalf * ( hfx_z(i,j,k) + hfx_z(i,j,k+1) ) * (one + epsv*qv);
        if (qfx1_z) {
            hfxv_c += epsv * theta * myhalf * ( qfx1_z(i,j,k) + qfx1_z(i,j,k+1) );
        }
 
        // Shear production, -<u'w'> dU/dz - <v'w'> dV/dz. MOST sets the direction
        // of the mean shear anti-parallel to the surface stress and its magnitude
        // to u_* phi_m/(kappa z_1), so the dot product is the cell-centered stress
        // projected onto the surface stress. This reduces to |tau| u_* phi_m /
        // (kappa z_1) when the stress does not turn across the cell.
        if (tau_s > eps) {
            const amrex::Real tau_proj = (t13_c*t13_s + t23_c*t23_s) / tau_s;
            source_term += std::max(tau_proj,amrex::Real(0.)) * ustar * phi_m / (KAPPA * z1);
        }
 
        // Buoyancy production, which needs no gradient at all
        source_term += (CONST_GRAV/theta_mean) * hfxv_c;
    }
    else
    {
        // Compute some relevant derivatives
        amrex::Real dthetadz, dudz, dvdz;
        ComputeVerticalDerivativesPBL(i, j, k,
                                      uvel, vvel, cell_data, izmin, izmax, dz_inv,
                                      c_ext_dir_on_zlo, c_ext_dir_on_zhi,
                                      u_ext_dir_on_zlo, u_ext_dir_on_zhi,
                                      v_ext_dir_on_zlo, v_ext_dir_on_zhi,
                                      dthetadz, dudz, dvdz,
                                      moisture_indices);
 
        // Second-order turbulent fluxes, e.g.:
        // -<uw> = L q SM dU/dz  (NN09, Eqn. 18)
        //       = Kmv/rho dU/dz
 
        // Shear Production
        source_term += K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz);
 
        // Buoyancy Production
        source_term -= (CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz;
    }
 
    // Dissipation (NN09, Eqn. 12)
    amrex::Real qke = two * cell_prim(i,j,k,PrimKE_comp);
    if (std::abs(qke) > zero) {
        source_term -= rho * std::pow(qke,amrex::Real(1.5)) /
                       (pbl_mynn_B1_l * K_turb(i,j,k,EddyDiff::Turb_lengthscale));
    }
 
    return source_term;
}

Referenced by AddTurbKESources().

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeVerticalDerivativesPBL()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void ComputeVerticalDerivativesPBL	(	int	i,
		int	j,
		int	k,
		const amrex::Array4< const amrex::Real > &	uvel,
		const amrex::Array4< const amrex::Real > &	vvel,
		const amrex::Array4< const amrex::Real > &	cell_data,
		const int	izmin,
		const int	izmax,
		const PBLDerivativeDzInv &	dz_inv,
		const bool	c_ext_dir_on_zlo,
		const bool	c_ext_dir_on_zhi,
		const bool	u_ext_dir_on_zlo,
		const bool	u_ext_dir_on_zhi,
		const bool	v_ext_dir_on_zlo,
		const bool	v_ext_dir_on_zhi,
		amrex::Real &	dthetadz,
		amrex::Real &	dudz,
		amrex::Real &	dvdz,
		const MoistureComponentIndices &	moisture_indices,
		const StratType	strat = StratType::thetav
	)

{
    // ext_dir ghost values hold the value at the physical boundary, which is
    // half a ghost-center spacing from the adjacent cell center. Use a
    // two-point, first-order boundary-to-center difference with inverse
    // distance 2/dz.
    if ( k==izmax && c_ext_dir_on_zhi ) {
        dthetadz = ( GetPBLStratVar(i,j,k+1,strat,cell_data,moisture_indices)
                    -GetPBLStratVar(i,j,k  ,strat,cell_data,moisture_indices) )*dz_inv.dirichlet_hi;
    } else if ( k==izmin && c_ext_dir_on_zlo ) {
        dthetadz = ( GetPBLStratVar(i,j,k  ,strat,cell_data,moisture_indices)
                    -GetPBLStratVar(i,j,k-1,strat,cell_data,moisture_indices) )*dz_inv.dirichlet_lo;
    } else {
        dthetadz = ( GetPBLStratVar(i,j,k+1,strat,cell_data,moisture_indices)
                    -GetPBLStratVar(i,j,k-1,strat,cell_data,moisture_indices) )*dz_inv.interior;
    }
 
    if ( k==izmax && u_ext_dir_on_zhi ) {
        dudz = myhalf*( uvel(i  ,j,k+1) - uvel(i  ,j,k)
                       +uvel(i+1,j,k+1) - uvel(i+1,j,k) )*dz_inv.dirichlet_hi;
    } else if ( k==izmin && u_ext_dir_on_zlo ) {
        dudz = myhalf*( uvel(i  ,j,k) - uvel(i  ,j,k-1)
                       +uvel(i+1,j,k) - uvel(i+1,j,k-1) )*dz_inv.dirichlet_lo;
    } else {
        dudz = myhalf*( uvel(i  ,j,k+1) - uvel(i  ,j,k-1)
                       +uvel(i+1,j,k+1) - uvel(i+1,j,k-1) )*dz_inv.interior;
    }
 
    if ( k==izmax && v_ext_dir_on_zhi ) {
        dvdz = myhalf*( vvel(i,j  ,k+1) - vvel(i,j  ,k)
                       +vvel(i,j+1,k+1) - vvel(i,j+1,k) )*dz_inv.dirichlet_hi;
    } else if ( k==izmin && v_ext_dir_on_zlo ) {
        dvdz = myhalf*( vvel(i,j  ,k) - vvel(i,j  ,k-1)
                       +vvel(i,j+1,k) - vvel(i,j+1,k-1) )*dz_inv.dirichlet_lo;
    } else {
        dvdz = myhalf*( vvel(i,j  ,k+1) - vvel(i,j  ,k-1)
                       +vvel(i,j+1,k+1) - vvel(i,j+1,k-1) )*dz_inv.interior;
    }
}

Referenced by ComputeDiffusivityMRF(), ComputeDiffusivityMYJ(), ComputeDiffusivityMYNN25(), ComputeDiffusivityMYNNEDMF(), ComputeDiffusivityYSU(), ComputeDiffusivityYSUNew(), and ComputeQKESourceTerms().

Here is the call graph for this function:

Here is the caller graph for this function:

ComputeVirtualTStarPBL()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeVirtualTStarPBL	(	amrex::Real	t_star,
		amrex::Real	q_star,
		amrex::Real	theta,
		amrex::Real	qv,
		bool	use_moisture
	)

Virtual temperature scale, thetav_* = theta_*(1 + epsv qv) + epsv theta q_*. The surface virtual heat flux is <w'thetav'> = -u_* thetav_*.

Parameters

[in]	t_star	temperature scale
[in]	q_star	moisture scale
[in]	theta	potential temperature
[in]	qv	water vapor mixing ratio
[in]	use_moisture	whether the moisture term is active

{
    amrex::Real tstar_v = t_star * (one + epsv*qv);
    if (use_moisture) { tstar_v += epsv * theta * q_star; }
    return tstar_v;
}

Referenced by ComputeDiffusivityMYJ(), ComputeDiffusivityMYNN25(), ComputeDiffusivityMYNNEDMF(), ComputeDiffusivityYSU(), and ComputeDiffusivityYSUNew().

Here is the caller graph for this function:

GetPBLStratVar()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetPBLStratVar	(	int	i,
		int	j,
		int	k,
		const StratType	strat,
		const amrex::Array4< const amrex::Real > &	cell_data,
		const MoistureComponentIndices &	moisture_indices
	)

Compute vertical derivatives of potential temperature and velocity for PBL models.

Parameters

	i	x-index
	j	y-index
	k	z-index
[in]	uvel	x-velocity
[in]	vvel	y-velocity
[in]	cell_data	conserved cell-center variables
[in]	izmin	lower z-index limit
[in]	izmax	upper z-index limit
[in]	dz_inv	inverse physical distance spacings
[in]	c_ext_dir_on_zlo	Dirichlet condition on z-low face for cell data
[in]	c_ext_dir_on_zhi	Dirichlet condition on z-high face for cell data
[in]	u_ext_dir_on_zlo	Dirichlet condition on z-low face for u-velocity
[in]	u_ext_dir_on_zhi	Dirichlet condition on z-high face for u-velocity
[in]	v_ext_dir_on_zlo	Dirichlet condition on z-low face for v-velocity
[in]	v_ext_dir_on_zhi	Dirichlet condition on z-high face for v-velocity
[out]	dthetadz	computed vertical derivative of potential temperature
[out]	dudz	computed vertical derivative of x-velocity
[out]	dvdz	computed vertical derivative of y-velocity
[in]	moisture_indices	moisture variable component indices Thermodynamic variable whose vertical gradient the PBL schemes use for buoyancy, selected by erf.thermal_stratification.
[in]	strat	stratification variable to return

{
    if (strat == StratType::theta) {
        return cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp);
    } else if (strat == StratType::thetal) {
        return GetThetal(i,j,k,cell_data,moisture_indices);
    }
    return GetThetav(i,j,k,cell_data,moisture_indices);
}

Referenced by ComputeVerticalDerivativesPBL().

Here is the call graph for this function:

Here is the caller graph for this function:

WidenForPBLH()

std::unique_ptr<amrex::MultiFab> WidenForPBLH ( const amrex::MultiFab &  src, const int  ng, const std::function< void(amrex::MultiFab &)> &  fill  )

inline

Copy of a field read by the PBLH passes, widened to ng horizontal ghost cells.

With erf.pblh_smoothing_passes = n the PBLH passes read their inputs n columns outside the tile (see ApplyPBLHSmoothing), but some of those inputs carry fewer ghost cells: the surface-layer fields have one on the bottom surface. The PBL schemes read such fields from a copy made here instead.

The valid cells and the field's own ghost cells are copied unchanged, so a run that needs no more halo than the field has gives the same answer as reading the field itself. fill then fills the new ghost cells that a box of the level covers, inside the domain or across a periodic edge, the same way the field's owner fills its own (for a surface-layer field, SurfaceLayer::fill_planar_boundary). Ghost cells no box covers – outside a non-periodic domain edge, or next to a coarser level – take the value of the nearest cell of the field's original halo. The smoothing never reads a column outside a non-periodic edge; next to a coarser level this is a zeroth-order extrapolation of the surface scales, while the state there is the interpolated one.

Parameters

[in]	src	Field to widen
[in]	ng	Ghost cells wanted in x and y; z keeps those of src
[in]	fill	Callable taking the widened MultiFab& and filling its ghost cells from the valid cells of the other boxes of the level

fill is a std::function, not a template parameter: nvcc does not accept the device lambda below inside a function template instantiated with a lambda local to the caller.

{
    const amrex::IntVect ng_src = src.nGrowVect();
    const amrex::IntVect ng_dst(ng, ng, ng_src[2]);
    const int ncomp = src.nComp();
 
    // A value no field holds, marking the cells nothing below fills
    constexpr amrex::Real unset = std::numeric_limits<amrex::Real>::max();
 
    auto dst = std::make_unique<amrex::MultiFab>(src.boxArray(), src.DistributionMap(),
                                                 ncomp, ng_dst);
    dst->setVal(unset);
    amrex::MultiFab::Copy(*dst, src, 0, 0, ncomp, ng_src);
    fill(*dst);
 
    for (amrex::MFIter mfi(*dst); mfi.isValid(); ++mfi) {
        // The cells holding src's own data; never written below, so they can be read
        const amrex::Box have = amrex::grow(mfi.validbox(), ng_src);
        const auto lo = amrex::lbound(have);
        const auto hi = amrex::ubound(have);
        auto a = dst->array(mfi);
        amrex::ParallelFor(mfi.fabbox(), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            if (have.contains(i, j, k)) { return; }
            const int ic = amrex::min(amrex::max(i, lo.x), hi.x);
            const int jc = amrex::min(amrex::max(j, lo.y), hi.y);
            for (int n = 0; n < ncomp; ++n) {
                if (a(i, j, k, n) == unset) { a(i, j, k, n) = a(ic, jc, k, n); }
            }
        });
    }
    return dst;
}

Referenced by ComputeDiffusivityMRF(), and ComputeDiffusivityYSUNew().

Here is the call graph for this function:

Here is the caller graph for this function: