36 #pragma omp parallel if (Gpu::notInLaunchRegion())
38 for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
41 const Box &bx = mfi.growntilebox(1);
42 const Box &dbx = geom.Domain();
43 Box sbx(bx.smallEnd(), bx.bigEnd());
45 AMREX_ALWAYS_ASSERT(sbx.smallEnd(2) == dbx.smallEnd(2) && sbx.bigEnd(2) == dbx.bigEnd(2));
48 const auto& cell_data = cons_in.const_array(mfi);
49 const auto& uvel =
xvel.const_array(mfi);
50 const auto& vvel =
yvel.const_array(mfi);
52 const auto& z0_arr = SurfLayer->get_z0(level)->const_array();
53 const auto& ws10av_arr = SurfLayer->get_mac_avg(level,5)->const_array(mfi);
54 const auto& t10av_arr = SurfLayer->get_mac_avg(level,2)->const_array(mfi);
55 const auto& t_surf_arr = SurfLayer->get_t_surf(level)->const_array(mfi);
56 const auto& over_land_arr = (SurfLayer->get_lmask(level)) ? SurfLayer->get_lmask(level)->const_array(mfi) :
58 const Array4<Real const> z_nd_arr = z_phys_nd->array(mfi);
59 const Real most_zref = SurfLayer->get_zref();
62 bool invalid_zref =
false;
63 if (use_terrain_fitted_coords) {
64 invalid_zref = most_zref != 10.0;
67 Real dz = geom.CellSize(2);
68 invalid_zref = int((most_zref - 0.5*dz)/dz) != int((10.0 - 0.5*dz)/dz);
71 Print() <<
"most_zref = " << most_zref << std::endl;
72 Abort(
"MOST Zref must be 10m for YSU PBL scheme");
76 const GeometryData gdata = geom.data();
77 const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
78 FArrayBox pbl_height(xybx,1);
79 IArrayBox pbl_index(xybx,1);
80 const auto& pblh_arr = pbl_height.array();
81 const auto& pbli_arr = pbl_index.array();
89 ParallelFor(xybx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int) noexcept
93 const Real t_surf = t_surf_arr(i,j,0);
94 const Real t_layer = t10av_arr(i,j,0);
95 const Real ws_layer = ws10av_arr(i,j,0);
96 const Real Rib_layer =
CONST_GRAV * most_zref / (ws_layer*ws_layer) * (t_layer - t_surf)/(t_layer);
99 if (Rib_layer < unst_Ribcr) {
100 Abort(
"For now, YSU PBL only supports stable conditions");
107 bool over_land = (over_land_arr) ? over_land_arr(i,j,0) : 1;
108 if (over_land && !force_over_water) {
113 const Real z0 = z0_arr(i,j,0);
114 const Real Rossby = ws_layer/(f0*z0);
115 Rib_cr = min(0.16*std::pow(1.0e-7*Rossby,-0.18),0.3);
118 bool above_critical =
false;
120 Real Rib_up = Rib_layer, Rib_dn;
122 while (!above_critical and bx.contains(i,j,kpbl+1)) {
124 const Real zval = use_terrain_fitted_coords ?
126 const Real ws2_level = 0.25*( (uvel(i,j,kpbl)+uvel(i+1,j ,kpbl))*(uvel(i,j,kpbl)+uvel(i+1,j ,kpbl))
127 + (vvel(i,j,kpbl)+vvel(i ,j+1,kpbl))*(vvel(i,j,kpbl)+vvel(i ,j+1,kpbl)) );
130 Rib_up = (
theta-base_theta)/base_theta *
CONST_GRAV * zval / ws2_level;
131 above_critical = Rib_up >= Rib_cr;
135 if (Rib_dn >= Rib_cr) {
137 }
else if (Rib_up <= Rib_cr)
140 interp_fact = (Rib_cr - Rib_dn) / (Rib_up - Rib_dn);
143 const Real zval_up = use_terrain_fitted_coords ?
145 const Real zval_dn = use_terrain_fitted_coords ?
147 pblh_arr(i,j,0) = zval_dn + interp_fact*(zval_up-zval_dn);
149 const Real zval_0 = use_terrain_fitted_coords ?
151 const Real zval_1 = use_terrain_fitted_coords ?
153 if (pblh_arr(i,j,0) < 0.5*(zval_0+zval_1) ) {
156 pbli_arr(i,j,0) = kpbl;
167 const auto& u_star_arr = SurfLayer->get_u_star(level)->const_array(mfi);
168 const auto& l_obuk_arr = SurfLayer->get_olen(level)->const_array(mfi);
169 const Array4<Real > &K_turb = eddyViscosity.array(mfi);
179 const auto& dxInv = geom.InvCellSizeArray();
180 const Real dz_inv = geom.InvCellSize(2);
181 const int izmin = geom.Domain().smallEnd(2);
182 const int izmax = geom.Domain().bigEnd(2);
184 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
186 const Real zval = use_terrain_fitted_coords ?
190 const Real dz_terrain = met_h_zeta/dz_inv;
191 if (k < pbli_arr(i,j,0)) {
193 constexpr Real zfacmin = 1e-8;
194 constexpr Real phifac = 8.0;
195 constexpr Real wstar3 = 0.0;
196 constexpr Real pfac = 2.0;
197 const Real zfac = std::min(std::max(1 - zval / pblh_arr(i,j,0), zfacmin ), 1.0);
199 const Real ust3 = u_star_arr(i,j,0) * u_star_arr(i,j,0) * u_star_arr(i,j,0);
200 Real wscalek = ust3 + phifac *
KAPPA * wstar3 * (1.0 - zfac);
201 wscalek = std::pow(wscalek, 1.0/3.0);
203 const Real phi_term = 1 + 5 * zval / l_obuk_arr(i,j,0);
204 wscalek = std::max(u_star_arr(i,j,0) / phi_term, 0.001);
209 constexpr Real lam0 = 30.0;
210 constexpr Real min_richardson = -100.0;
211 constexpr Real prandtl_max = 4.0;
212 Real dthetadz, dudz, dvdz;
214 uvel, vvel, cell_data, izmin, izmax, 1.0/dz_terrain,
215 c_ext_dir_on_zlo, c_ext_dir_on_zhi,
216 u_ext_dir_on_zlo, u_ext_dir_on_zhi,
217 v_ext_dir_on_zlo, v_ext_dir_on_zhi,
218 dthetadz, dudz, dvdz, moisture_indices);
219 const Real shear_squared = dudz*dudz + dvdz*dvdz + 1.0e-9;
222 const Real lambdadz = std::min(std::max(0.1*dz_terrain , lam0), 300.0);
223 const Real lengthscale = lambdadz *
KAPPA * zval / (lambdadz +
KAPPA * zval);
224 const Real turbfact = lengthscale * lengthscale * std::sqrt(shear_squared);
226 if (richardson < 0) {
227 richardson = max(richardson, min_richardson);
228 Real sqrt_richardson = std::sqrt(-richardson);
229 K_turb(i,j,k,
EddyDiff::Mom_v) =
rho * turbfact * (1.0 - 8.0 * richardson / (1.0 + 1.746 * sqrt_richardson));
230 K_turb(i,j,k,
EddyDiff::Theta_v) =
rho * turbfact * (1.0 - 8.0 * richardson / (1.0 + 1.286 * sqrt_richardson));
232 const Real oneplus5ri = 1.0 + 5.0 * richardson;
234 const Real
prandtl = std::min(1.0+2.1*richardson, prandtl_max);
240 constexpr Real ckz = 0.001;
241 constexpr Real Kmax = 1000.0;
242 const Real rhoKmin = ckz * dz_terrain *
rho;
243 const Real rhoKmax =
rho * Kmax;
250 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
TurbChoice turbChoice
Definition: ERF_DiffSetup.H:2
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
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 amrex::Real &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)
Definition: ERF_PBLModels.H:147
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real Compute_h_zeta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:47
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_Zrel_AtCellCenter(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:390
@ yvel_bc
Definition: ERF_IndexDefines.H:88
@ cons_bc
Definition: ERF_IndexDefines.H:76
@ xvel_bc
Definition: ERF_IndexDefines.H:87
@ ext_dir
Definition: ERF_IndexDefines.H:209
@ Theta_v
Definition: ERF_IndexDefines.H:176
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:180
@ Mom_v
Definition: ERF_IndexDefines.H:175
@ theta
Definition: ERF_MM5.H:20
@ rho
Definition: ERF_Kessler.H:22
@ xvel
Definition: ERF_IndexDefines.H:141
@ yvel
Definition: ERF_IndexDefines.H:142
real(c_double), parameter prandtl
Definition: ERF_module_model_constants.F90:88
amrex::Real pbl_ysu_land_Ribcr
Definition: ERF_TurbStruct.H:395
amrex::Real pbl_ysu_coriolis_freq
Definition: ERF_TurbStruct.H:386
bool pbl_ysu_force_over_water
Definition: ERF_TurbStruct.H:392
amrex::Real pbl_ysu_unst_Ribcr
Definition: ERF_TurbStruct.H:397