ERF_ComputeDiffusivityYSU.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 "ERF_TileNoZ.H"
#include <cmath>

Include dependency graph for ERF_ComputeDiffusivityYSU.cpp:

Functions
void	ComputeDiffusivityYSU (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)

Function Documentation

ComputeDiffusivityYSU()

void ComputeDiffusivityYSU	(	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 diffusivity using the YSU PBL scheme.

Parameters

[in]	xvel	x-velocity field
[in]	yvel	y-velocity field
[in]	cons_in	Conservative state
[out]	eddyViscosity	Computed eddy viscosity coefficients
[in]	geom	Geometry used for spacing and domain info
[in]	turbChoice	Turbulence model options
[in]	SurfLayer	Surface layer model data
[in]	use_terrain_fitted_coords	Flag to use terrain-fitted coordinates
[in]	level	AMR level
[in]	bc_ptr	Boundary condition records
[in]	z_phys_nd	Nodal physical heights
[in]	z_phys_cc	Cell-centered physical heights
[in]	moisture_indices	Indices for moisture components

{
    const StratType pbl_strat = turbChoice.pbl_strat_type;
    /*
      YSU PBL initially introduced by S.-Y. Hong, Y. Noh, and J. Dudhia, MWR, 2006 [HND06]
 
      Further Modifications from S.-Y. Hong, Q. J. R. Meteorol. Soc., 2010 [H10]
 
      Implementation follows WRF as of early 2024 with some simplifications
    */
 
    const Real most_zref = SurfLayer->get_zref_min(level);
 
    // Require that MOST zref is 10 m so we get the wind speed at 10 m from most
    bool invalid_zref = false;
    if (use_terrain_fitted_coords) {
        invalid_zref = most_zref != Real(10.0);
    } else {
        // zref gets reset to nearest cell center, so assert that zref is in the same cell as the 10m point
        Real dz = geom.CellSize(2);
        invalid_zref = int((most_zref - myhalf*dz)/dz) != int((Real(10.0) - myhalf*dz)/dz);
    }
    if (invalid_zref) {
        Print() << "most_zref = " << most_zref << std::endl;
        Abort("MOST Zref must be 10m for YSU PBL scheme");
    }
 
    // The PBL height for the SurfaceLayer (the w* correction and diagnostics), as MRF
    // and YSUNew hand theirs over
    MultiFab pblh_mf(eddyViscosity.boxArray(), eddyViscosity.DistributionMap(), 1, 0);
 
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    // TileNoZ, not TilingIfNotGPU: every iterate must span the full column (asserted
    // below), and the CPU default tile size would split any box of 16+ cells in z.
    for ( MFIter mfi(eddyViscosity,TileNoZ()); mfi.isValid(); ++mfi) {
 
        // Pull out the box we're working on, make sure it covers full domain in z-direction
        const Box &bx = mfi.growntilebox(1);
        const Box &dbx = geom.Domain();
        Box sbx(bx.smallEnd(), bx.bigEnd());
        sbx.grow(2,-1);
        AMREX_ALWAYS_ASSERT(sbx.smallEnd(2) == dbx.smallEnd(2) && sbx.bigEnd(2) == dbx.bigEnd(2));
 
        // Get some data in arrays
        const auto& cell_data = cons_in.const_array(mfi);
        const auto& uvel = xvel.const_array(mfi);
        const auto& vvel = yvel.const_array(mfi);
 
        const auto& z0_arr        = SurfLayer->get_z0(level)->const_array(mfi);
        const auto& ws10av_arr    = SurfLayer->get_mac_avg(level,6)->const_array(mfi);
        const auto& t10av_arr     = SurfLayer->get_mac_avg(level,3)->const_array(mfi);
        const auto& t_surf_arr    = SurfLayer->get_t_surf(level)->const_array(mfi);
        const auto& over_land_arr = (SurfLayer->get_lmask(level)) ? SurfLayer->get_lmask(level)->const_array(mfi) :
                                                                  Array4<int> {};
        const Array4<Real const> z_nd_arr = z_phys_nd->array(mfi);
        const PBLDerivativeDzInv_T pbl_derivative_dz_inv{z_phys_cc->const_array(mfi)};
 
        // create flattened boxes to store PBL height
        const GeometryData gdata = geom.data();
        const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
        FArrayBox pbl_height(xybx,1,The_Async_Arena());
        IArrayBox pbl_index(xybx,1,The_Async_Arena());
        const auto& pblh_arr = pbl_height.array();
        const auto& pbli_arr = pbl_index.array();
 
        // -- Diagnose PBL height - starting out assuming non-moist --
        // loop is only over i,j in order to find height at each x,y
        const Real f0 = turbChoice.pbl_ysu_coriolis_freq;
        const bool force_over_water = turbChoice.pbl_ysu_force_over_water;
        const Real land_Ribcr = turbChoice.pbl_ysu_land_Ribcr;
        const Real unst_Ribcr = turbChoice.pbl_ysu_unst_Ribcr;
        ParallelFor(xybx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
        {
            // Reconstruct a surface bulk Richardson number from the surface layer model
            // In WRF, this value is supplied to YSU by the surface layer model
            const Real t_surf = t_surf_arr(i,j,0);
            const Real t_layer = t10av_arr(i,j,0);
            const Real ws_layer = ws10av_arr(i,j,0);
            // The squared wind speed in the bulk Richardson numbers is floored at 1 m^2/s^2,
            // as in WRF's YSU (and YSUNew), so a calm column gives a finite number
            const Real ws2_floor = Real(1.0);
            const Real Rib_layer = CONST_GRAV * most_zref / amrex::max(ws_layer*ws_layer, ws2_floor)
                                 * (t_layer - t_surf)/(t_layer);
 
            // For now, we only support stable boundary layers
            if (Rib_layer < unst_Ribcr) {
                Abort("For now, YSU PBL only supports stable conditions");
            }
 
            // TODO: unstable BLs
 
            // PBL Height: Stable Conditions
            Real Rib_cr;
            bool over_land = (over_land_arr) ? over_land_arr(i,j,0) : 1;
            if (over_land && !force_over_water) {
                Rib_cr = land_Ribcr;
            } else { // over water
                // Velocity at z=10 m comes from MOST -> currently the average using whatever averaging MOST uses.
                // TODO: Revisit this calculation with local ws10?
                const Real z0 = z0_arr(i,j,0);
                const Real Rossby = ws_layer/(f0*z0);
                Rib_cr = min(Real(0.16)*std::pow(Real(1.0e-7)*Rossby,-Real(0.18)),Real(0.3)); // Note: upper bound in WRF code, but not H10 paper
            }
 
            bool above_critical = false;
            int kpbl = 0;
            // Rib_dn is overwritten on the first pass of the loop; it starts equal to
            // Rib_up only so that a column with no k = 1 does not read it unset.
            Real Rib_up = Rib_layer, Rib_dn = Rib_layer;
            const Real base_theta = cell_data(i,j,0,RhoTheta_comp) / cell_data(i,j,0,Rho_comp);
            while (!above_critical and bx.contains(i,j,kpbl+1)) {
                kpbl += 1;
                const Real zval = use_terrain_fitted_coords
                                ? Compute_Zrel_AtCellCenter(i,j,kpbl,z_nd_arr)
                                : (kpbl + myhalf)*gdata.CellSize(2);
                const Real ws2_level = fourth*( (uvel(i,j,kpbl)+uvel(i+1,j  ,kpbl))*(uvel(i,j,kpbl)+uvel(i+1,j  ,kpbl))
                                              + (vvel(i,j,kpbl)+vvel(i  ,j+1,kpbl))*(vvel(i,j,kpbl)+vvel(i  ,j+1,kpbl)) );
                const Real theta = cell_data(i,j,kpbl,RhoTheta_comp) / cell_data(i,j,kpbl,Rho_comp);
                Rib_dn = Rib_up;
                Rib_up = (theta-base_theta)/base_theta * CONST_GRAV * zval / amrex::max(ws2_level, ws2_floor);
                above_critical = Rib_up >= Rib_cr;
            }
 
            Real interp_fact;
            if (Rib_dn >= Rib_cr) {
                interp_fact = zero;
            } else if (Rib_up <= Rib_cr)
                interp_fact = one;
            else {
                interp_fact = (Rib_cr - Rib_dn) / (Rib_up - Rib_dn);
            }
 
            const Real zval_up = use_terrain_fitted_coords
                               ? Compute_Zrel_AtCellCenter(i,j,kpbl,z_nd_arr)
                               : (kpbl + myhalf)*gdata.CellSize(2);
            const Real zval_dn = use_terrain_fitted_coords
                               ? Compute_Zrel_AtCellCenter(i,j,kpbl-1,z_nd_arr)
                               : (kpbl-1 + myhalf)*gdata.CellSize(2);
            pblh_arr(i,j,0) = zval_dn + interp_fact*(zval_up-zval_dn);
 
            const Real zval_0 = use_terrain_fitted_coords
                              ? Compute_Zrel_AtCellCenter(i,j,0,z_nd_arr)
                              : (myhalf)*gdata.CellSize(2);
            const Real zval_1 = use_terrain_fitted_coords
                              ? Compute_Zrel_AtCellCenter(i,j,1,z_nd_arr)
                              : (Real(1.5))*gdata.CellSize(2);
            if (pblh_arr(i,j,0) < myhalf*(zval_0+zval_1) ) {
                kpbl = 0;
            }
            pbli_arr(i,j,0) = kpbl;
        });
 
        // -- Compute nonlocal/countergradient mixing parameters --
        // Not included for stable so nothing to do until unstable treatment is added
 
        // -- Compute entrainment parameters --
        // 0 for stable so nothing to do?
 
        // -- Compute diffusion coefficients --
 
        const auto& u_star_arr = SurfLayer->get_u_star(level)->const_array(mfi);
        const auto& t_star_arr = SurfLayer->get_t_star(level)->const_array(mfi);
        const auto  q_star_arr = (use_moisture) ?
                                 SurfLayer->get_q_star(level)->const_array(mfi) : Array4<const Real>{};
        const auto& l_obuk_arr = SurfLayer->get_olen(level)->const_array(mfi);
        const Array4<Real      > &K_turb = eddyViscosity.array(mfi);
 
        // 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) );
 
        const auto& dxInv = geom.InvCellSizeArray();
        const Real dz_inv = geom.InvCellSize(2);
        const int izmin = geom.Domain().smallEnd(2);
        const int izmax = geom.Domain().bigEnd(2);
 
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            const Real zval = use_terrain_fitted_coords
                            ? Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr)
                            : (k + myhalf)*gdata.CellSize(2);
            const Real rho = cell_data(i,j,k,Rho_comp);
            const Real met_h_zeta = use_terrain_fitted_coords ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr) : one;
            const Real dz_terrain = met_h_zeta/dz_inv;
            const Real small  = Real(1.0e-10); // floor on |L| (same as MRF and YSUNew) and on phi_term
            const Real l_obuk = std::copysign(std::max(std::fabs(l_obuk_arr(i,j,0)),small),l_obuk_arr(i,j,0));
            if (k < pbli_arr(i,j,0)) {
                // -- Compute diffusion coefficients within PBL
                constexpr Real zfacmin = Real(1e-8); // value from WRF
                constexpr Real pfac    = Real(2.);  // profile exponent
                const Real zfac = std::min(std::max(amrex::Real(1) - zval / pblh_arr(i,j,0), zfacmin), amrex::Real(1));
                // stable only
                const Real phi_term = std::max(amrex::Real(1) + amrex::Real(5) * zval / l_obuk, small); // phi_term appears in WRF but not papers
                const Real wscalek = std::max(u_star_arr(i,j,0) / phi_term, Real(0.001)); // Real(0.001) limit appears in WRF but not papers
                K_turb(i,j,k,EddyDiff::Mom_v)   = rho * wscalek * KAPPA * zval * std::pow(zfac, pfac);
                K_turb(i,j,k,EddyDiff::Theta_v) = K_turb(i,j,k,EddyDiff::Mom_v);
            } else {
                // -- Compute coefficients in free stream above PBL
                constexpr Real lam0 = Real(30.0);
                constexpr Real min_richardson = -Real(100.0);
                constexpr Real prandtl_max = Real(4.0);
                Real dthetadz, 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,
                                              dthetadz, dudz, dvdz, moisture_indices, pbl_strat);
 
                // This branch is the free atmosphere above the PBL, so it only
                // reaches the first cell when the PBL index collapses to it. The
                // resolved gradients are unusable there (ERF #4037), so fall back
                // on the MOST profile; see ApplySurfaceLayerGradientsPBL.
                if (k == izmin) {
                    const Real theta_k = cell_data(i,j,k,RhoTheta_comp) / rho;
                    const Real qv_k    = (moisture_indices.qv >= 0) ?
                                         cell_data(i,j,k,moisture_indices.qv) / rho : zero;
                    PBLSurfaceLayerGradient sl;
                    sl.u_star  = u_star_arr(i,j,0);
                    sl.tstar_v = ComputeVirtualTStarPBL(t_star_arr(i,j,0),
                                                        (q_star_arr) ? q_star_arr(i,j,0) : zero,
                                                        theta_k, qv_k, use_moisture);
                    sl.zval    = zval;
                    sl.zeta    = zval / l_obuk;
                    ApplySurfaceLayerGradientsPBL(sl, dthetadz, dudz, dvdz);
                }
 
                const Real shear_squared = dudz*dudz + dvdz*dvdz + Real(1.0e-9); // Real(1.0e-9) from WRF to avoid divide by zero
                const Real theta = cell_data(i,j,k,RhoTheta_comp) / cell_data(i,j,k,Rho_comp);
                Real richardson = CONST_GRAV / theta * dthetadz / shear_squared;
                const Real lambdadz = std::min(std::max(Real(0.1)*dz_terrain , lam0), Real(300.0)); // in WRF, H10 paper just says use lam0
                const Real lengthscale = lambdadz * KAPPA * zval / (lambdadz + KAPPA * zval);
                const Real turbfact = lengthscale * lengthscale * std::sqrt(shear_squared);
 
                if (richardson < 0) {
                    richardson = max(richardson, min_richardson);
                    Real sqrt_richardson = std::sqrt(-richardson);
                    K_turb(i,j,k,EddyDiff::Mom_v) = rho * turbfact * (one - Real(8.0) * richardson / (one + Real(1.746) * sqrt_richardson));
                    K_turb(i,j,k,EddyDiff::Theta_v) = rho * turbfact * (one - Real(8.0) * richardson / (one + Real(1.286) * sqrt_richardson));
                } else {
                    const Real oneplus5ri = one + Real(5.0) * richardson;
                    K_turb(i,j,k,EddyDiff::Theta_v) = rho * turbfact / (oneplus5ri * oneplus5ri);
                    const Real prandtl = std::min(one+Real(2.1)*richardson, prandtl_max); // limit from WRF
                    K_turb(i,j,k,EddyDiff::Mom_v) = K_turb(i,j,k,EddyDiff::Theta_v) * prandtl;
                }
            }
 
            // limit both diffusion coefficients - from WRF, not documented in papers
            constexpr Real ckz = Real(0.001);
            constexpr Real Kmax = Real(1000.0);
            const Real rhoKmin = ckz * dz_terrain * rho;
            const Real rhoKmax = rho * Kmax;
            K_turb(i,j,k,EddyDiff::Mom_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Mom_v) ,rhoKmax), rhoKmin);
            K_turb(i,j,k,EddyDiff::Theta_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Theta_v) ,rhoKmax), rhoKmin);
            K_turb(i,j,k,EddyDiff::Turb_lengthscale) = pblh_arr(i,j,0);
        });
 
        // HACK set bottom ghost cell to 1st cell
        ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
        {
            if (k==-1) {
                K_turb(i,j,k,EddyDiff::Mom_v) = K_turb(i,j,0,EddyDiff::Mom_v);
                K_turb(i,j,k,EddyDiff::Theta_v) = K_turb(i,j,0,EddyDiff::Theta_v);
            }
        });
 
        // Under TileNoZ the tile spans the full column, and pblh_mf has no ghost cells,
        // so the tiles together fill every cell
        {
            auto pblh_out = pblh_mf.array(mfi);
            ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
            {
                pblh_out(i,j,k) = pblh_arr(i,j,0);
            });
        }
    } // mfi
 
    SurfLayer->set_pblh(level, pblh_mf);
}

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:

Here is the caller graph for this function: