ERF_ComputeDiffusivityMYNN25.cpp File Reference

#include "ERF_SurfaceLayer.H"
#include "ERF_DirectionSelector.H"
#include "ERF_Diffusion.H"
#include "ERF_Constants.H"
#include "ERF_TurbStruct.H"
#include "ERF_PBLModels.H"

Include dependency graph for ERF_ComputeDiffusivityMYNN25.cpp:

Macros
#define	EXTRA_MYNN25_CHECKS   0

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

Macro Definition Documentation

EXTRA_MYNN25_CHECKS

#define EXTRA_MYNN25_CHECKS   0

Function Documentation

ComputeDiffusivityMYNN25()

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

Compute eddy diffusivities using the MYNN PBL model.

Parameters

[in]	xvel	X-velocity MultiFab.
[in]	yvel	Y-velocity MultiFab.
[in]	cons_in	Conserved variables MultiFab.
[out]	eddyViscosity	MultiFab for storing computed eddy diffusivities.
[in]	geom	Grid geometry.
[in]	turbChoice	Turbulence model options and configuration.
[in]	SurfLayer	Surface layer data.
[in]	use_terrain_fitted_coords	Flag to use terrain-fitted coordinates.
[in]	use_moisture	Flag to include moisture.
[in]	level	Level index.
[in]	bc_ptr	Boundary condition record pointers.
[in]	z_phys_nd	Nodal physical height.
[in]	z_phys_cc	Cell-centered physical height.
[in]	moisture_indices	Moisture component indices.

{
    const StratType pbl_strat = turbChoice.pbl_strat_type;
    auto mynn     = turbChoice.pbl_mynn;
    auto level2   = turbChoice.pbl_mynn_level2;
 
    Real Lt_alpha = mynn.Lt_alpha;
 
    // Optional boundary-layer-depth taper of l_T; needs a diagnosed PBL height.
    const bool use_Lt_taper = (mynn.Lt_taper_exp > zero);
    const Real Lt_taper_exp = mynn.Lt_taper_exp;
    if (use_Lt_taper && !SurfLayer->computes_pblh()) {
        Abort("erf.pbl_mynn_Lt_taper_exp > 0 needs a diagnosed PBL height; "
              "set erf.most.pblh_calc = MYNN25");
    }
 
    // Dirichlet flags to switch derivative stencil
    bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
    bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].hi(2) == ERFBCType::ext_dir) );
    bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
    bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].hi(2) == ERFBCType::ext_dir) );
    bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
    bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir) );
 
    // Epsilon
    Real eps = std::numeric_limits<Real>::epsilon();
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for ( MFIter mfi(eddyViscosity,TileNoZ()); mfi.isValid(); ++mfi) {
 
        const Box& bx = mfi.tilebox();
        const Array4<Real const>& cell_data = cons_in.array(mfi);
        const Array4<Real      >& K_turb    = eddyViscosity.array(mfi);
        const Array4<Real const>& uvel      = xvel.array(mfi);
        const Array4<Real const>& vvel      = yvel.array(mfi);
 
        // Compute some quantities that are constant in each column
        // Sbox is shrunk to only include the interior of the domain in the vertical direction to compute integrals
        // Box includes one ghost cell in each direction
        const Box& dbx = geom.Domain();
        AMREX_ALWAYS_ASSERT(bx.smallEnd(2) == dbx.smallEnd(2) && bx.bigEnd(2) == dbx.bigEnd(2));
 
        const GeometryData gdata = geom.data();
 
        const Box xybx = makeSlab(bx,2,0);
        FArrayBox qturb(bx,1,The_Async_Arena());
        FArrayBox qintegral(xybx,2,The_Async_Arena());
        qintegral.setVal<RunOn::Device>(zero);
        const Array4<Real> qint = qintegral.array();
        const Array4<Real> qvel = qturb.array();
 
        // vertical integrals to compute lengthscale
        if (use_terrain_fitted_coords) {
            const Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
            const auto invCellSize = geom.InvCellSizeArray();
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                // q^2 / 2 is the TKE
                qvel(i,j,k) = std::sqrt(two * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp));
                AMREX_ALWAYS_ASSERT_WITH_MESSAGE(qvel(i,j,k) > zero, "KE must have a positive value");
 
                const Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr);
                const Real dz   = Compute_h_zeta_AtCellCenter(i,j,k,invCellSize,z_nd_arr);
                Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz);
                Gpu::Atomic::Add(&qint(i,j,0,1),      qvel(i,j,k)*dz);
            });
        } else {
            ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                // q^2 / 2 is the TKE
                qvel(i,j,k) = std::sqrt(two * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp));
                AMREX_ALWAYS_ASSERT_WITH_MESSAGE(qvel(i,j,k) > zero, "KE must have a positive value");
 
                // Not multiplying by dz: it's constant and would fall out when we divide qint0/qint1 anyway
 
                const Real Zval = (k + myhalf)*gdata.CellSize(2);
                Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k));
                Gpu::Atomic::Add(&qint(i,j,0,1),      qvel(i,j,k));
            });
        }
 
        int izmin = geom.Domain().smallEnd(2);
        int izmax = geom.Domain().bigEnd(2);
 
        // Spatially varying MOST
        Real d_kappa   = KAPPA;
        Real d_gravity = CONST_GRAV;
 
        const auto& t_mean_mf = SurfLayer->get_mac_avg(level,5); // theta_v
        const auto& q_mean_mf = SurfLayer->get_mac_avg(level,4); // q_v
        const auto& u_star_mf = SurfLayer->get_u_star(level);
        const auto& t_star_mf = SurfLayer->get_t_star(level);
        const auto& q_star_mf = SurfLayer->get_q_star(level);
 
        const auto& tm_arr     = t_mean_mf->const_array(mfi);
        const auto& qm_arr     = q_mean_mf->const_array(mfi);
        const auto& u_star_arr = u_star_mf->const_array(mfi);
        const auto& t_star_arr = t_star_mf->const_array(mfi);
        const auto& q_star_arr = (use_moisture) ? q_star_mf->const_array(mfi) : Array4<Real>{};
 
        const auto& pblh_mf  = SurfLayer->get_pblh(level);
        const auto  pblh_arr = (use_Lt_taper) ? pblh_mf->const_array(mfi)
                                              : Array4<const Real>{};
 
        const Array4<Real const> z_nd_arr = z_phys_nd->const_array(mfi);
        const PBLDerivativeDzInv_T pbl_derivative_dz_inv{z_phys_cc->const_array(mfi)};
 
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            // Compute some partial derivatives that we will need (second order)
            // U and V derivatives are interpolated to account for staggered grid
            Real dthetavdz, dudz, dvdz;
            ComputeVerticalDerivativesPBL(i, j, k,
                                          uvel, vvel, cell_data, izmin, izmax, pbl_derivative_dz_inv(i,j,k),
                                          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,
                                          dthetavdz, dudz, dvdz,
                                          moisture_indices, pbl_strat);
 
            // Spatially varying MOST
            Real theta0 = tm_arr(i,j,0);
            Real qv0    = qm_arr(i,j,0);
            Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0);
            Real surface_latent_heat{0};
            if (use_moisture) {
                // Compute buoyancy flux (Stull Eqn. 4.4.5d)
                surface_latent_heat = -u_star_arr(i,j,0) * q_star_arr(i,j,0);
                surface_heat_flux *= (one + epsv*qv0);
                surface_heat_flux += epsv * theta0 * surface_latent_heat;
            }
 
            Real l_obukhov;
            if (std::abs(surface_heat_flux) > eps && u_star_arr(i,j,0) > eps) {
                l_obukhov = -( theta0 * u_star_arr(i,j,0)*u_star_arr(i,j,0)*u_star_arr(i,j,0) )
                           / ( d_kappa * d_gravity * surface_heat_flux );
            } else {
                l_obukhov = std::numeric_limits<Real>::max();
            }
 
            // Surface-layer length scale (NN09, Eqn. 53)
            AMREX_ASSERT(l_obukhov != 0);
            const Real zval = use_terrain_fitted_coords ? Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr) :
                                                          (k + myhalf)*gdata.CellSize(2);
            const Real zeta = zval/l_obukhov;
            Real l_S;
            if (zeta >= one) {
                l_S = KAPPA*zval/Real(3.7);
            } else if (zeta >= 0) {
                l_S = KAPPA*zval/(one + Real(2.7) * zeta);
            } else {
                l_S = KAPPA*zval*std::pow(one - Real(100.0) * zeta, Real(0.2));
            }
 
            // At the first cell center the vertical gradients computed above are
            // built from the foextrap (cons) and hoextrap (velocity) ghost values
            // that zlo.type = surface_layer installs, and those carry no surface
            // information: 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). Replace them with the MOST
            // profile gradients that the surface layer is imposing on this cell,
            //   |dU/dz|    = u_*      phi_m(zeta) / (kappa z)
            //   dthetav/dz = thetav_* phi_h(zeta) / (kappa z)
            // which is the same profile the stress and heat flux were derived from.
            if (k == izmin) {
                PBLSurfaceLayerGradient sl;
                sl.u_star  = u_star_arr(i,j,0);
                sl.tstar_v = ComputeVirtualTStarPBL(t_star_arr(i,j,0),
                                                    (use_moisture) ? q_star_arr(i,j,0) : zero,
                                                    theta0, qv0, use_moisture);
                sl.zval    = zval;
                sl.zeta    = zeta;
                ApplySurfaceLayerGradientsPBL(sl, dthetavdz, dudz, dvdz);
            }
 
            // ABL-depth length scale (NN09, Eqn. 54)
            Real l_T;
            if (qint(i,j,0,1) > zero) {
                l_T = Lt_alpha*qint(i,j,0,0)/qint(i,j,0,1);
            } else {
                l_T = std::numeric_limits<Real>::max();
            }
 
            // l_T is also the ABL-depth scale that the convective branch of
            // l_B below is built from; keep that one untapered, since there it
            // stands for the depth of the layer rather than a mixing length.
            const Real l_T_abl = l_T;
 
            // Optional taper toward the top of the boundary layer.  The factor
            // is floored rather than allowed to reach zero so that the 1/l_T
            // term in the master length scale stays finite.
            //
            // NOTE: above z = h the factor clamps at that floor, so the master
            // length scale is driven to ~0 there.  That is harmless when the
            // turbulence is confined to the boundary layer, but it would
            // suppress genuine elevated shear-driven mixing; hence the taper
            // is opt-in rather than a default.
            if (use_Lt_taper) {
                const Real pblh = pblh_arr(i,j,0);
                if (pblh > zero) {
                    Real f = one - zval/pblh;
                    f = (f < Real(1.0e-4)) ? Real(1.0e-4) : f;
                    l_T *= std::pow(f, Lt_taper_exp);
                }
            }
 
            // Buoyancy length scale (NN09, Eqn. 55)
            Real l_B;
            if (dthetavdz > zero) {
                Real N_brunt_vaisala = std::sqrt(CONST_GRAV/theta0 * dthetavdz);
                if (zeta < zero) {
                    Real qc = CONST_GRAV/theta0 * surface_heat_flux * l_T_abl; // velocity scale
                    qc = std::pow(qc,one/three);
                    l_B = (one + Real(5.0)*std::sqrt(qc/(N_brunt_vaisala * l_T_abl))) * qvel(i,j,k)/N_brunt_vaisala;
                } else {
                    l_B = qvel(i,j,k) / N_brunt_vaisala;
                }
            } else {
                l_B = std::numeric_limits<Real>::max();
            }
 
            // Master length scale
            Real Lm;
            if (mynn.config == MYNNConfigType::Chen2021) {
                Lm = std::pow(one/(l_S*l_S) + one/(l_T*l_T) + one/(l_B*l_B), -myhalf);
            } else {
                // NN09, Eqn 52
                Lm = one / (one/l_S + one/l_T + one/l_B);
            }
 
            // Calculate nondimensional production terms
            Real shearProd  = dudz*dudz + dvdz*dvdz;
            Real buoyProd   = -(CONST_GRAV/theta0) * dthetavdz;
            Real L2_over_q2 = Lm*Lm/(qvel(i,j,k)*qvel(i,j,k));
            Real GM         = L2_over_q2 * shearProd;
            Real GH         = L2_over_q2 * buoyProd;
 
            // Equilibrium (Level-2) q calculation follows NN09, Appendix A
            Real Rf  = level2.calc_Rf(GM, GH);
            Real SM2 = level2.calc_SM(Rf);
            Real qe2 = mynn.B1 * Lm*Lm * SM2 * (one-Rf) * shearProd;
            Real qe  = (qe2 < zero) ? zero : amrex::max(std::sqrt(qe2),eps);
 
            // Level 2 limiting introduced by Helfand and Labraga 1988 (NN09, Eqn. 42)
            Real alphac  = (qvel(i,j,k) >= qe) ? one : qvel(i,j,k) / qe;
//#if EXTRA_MYNN25_CHECKS
#if 0
            // VERY verbose diagnostic
            Real lGM = std::copysign(std::max(std::fabs(GM),level2.eps),GM);
            Real Ri  = -GH/lGM;
            if (alphac < one) {
                AllPrint() << "Level 2 limiter at " << IntVect(i,j,k) << " :"
                    << " ustar= " << u_star_arr(i,j,0)
                    << " alphac= " << alphac
                    << " Ri,SM2,SH2= " << Ri << " " << SM2 << " " << level2.calc_SH(Rf)
                    << std::endl;
            }
#endif
 
            // Level Real(2.5) stability functions
            Real SM, SH, SQ;
            mynn.calc_stability_funcs(SM,SH,SQ,GM,GH,alphac);
 
            // Clip SM, SH following WRF
            SM = amrex::min(amrex::max(SM, mynn.SMmin), mynn.SMmax);
            SH = amrex::min(amrex::max(SH, mynn.SHmin), mynn.SHmax);
            SQ = amrex::min(amrex::max(SQ, mynn.SQmin), mynn.SQmax);
#if EXTRA_MYNN25_CHECKS
            if (SM == mynn.SMmin) {
                Warning("SM clipped at min val");
            } else if (SM == mynn.SMmax) {
                Warning("SM clipped at max val");
            }
            if (SH == mynn.SHmin) {
                Warning("SH clipped at min val");
            } else if (SH == mynn.SHmax) {
                Warning("SH clipped at max val");
            }
#endif
 
            // Finally, compute the eddy viscosity/diffusivities
            const Real rho = cell_data(i,j,k,Rho_comp);
            K_turb(i,j,k,EddyDiff::Mom_v)   = rho * Lm * qvel(i,j,k) * SM;
            K_turb(i,j,k,EddyDiff::Theta_v) = rho * Lm * qvel(i,j,k) * SH;
            K_turb(i,j,k,EddyDiff::KE_v)    = rho * Lm * qvel(i,j,k) * SQ;
 
            // TODO: implement partial-condensation scheme?
            // Currently, implementation matches NN09 without rain (i.e.,
            // the liquid water potential temperature is equal to the
            // potential temperature.
 
            // NN09 gives the total water content flux; this assumes that
            // all the species have the same eddy diffusivity.  eddyDiffs is
            // zeroed once at allocation and nothing else writes Q_v, so this
            // must be unconditional: skipping it leaves the moisture variables
            // with a surface flux and no vertical turbulent transport.
            K_turb(i,j,k,EddyDiff::Q_v) = rho * Lm * qvel(i,j,k) * SH;
 
            K_turb(i,j,k,EddyDiff::Turb_lengthscale) = Lm;
        });
    }
}

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:

Here is the caller graph for this function: