27 const bool use_terrain = (z_phys_nd !=
nullptr);
28 const bool use_most = (most !=
nullptr);
44 Real eps = std::numeric_limits<Real>::epsilon();
47 #pragma omp parallel if (Gpu::notInLaunchRegion())
49 for ( MFIter mfi(eddyViscosity,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
51 const Box &bx = mfi.growntilebox(1);
52 const Array4<Real const>& cell_data = cons_in.array(mfi);
53 const Array4<Real >& K_turb = eddyViscosity.array(mfi);
54 const Array4<Real const>& uvel =
xvel.array(mfi);
55 const Array4<Real const>& vvel =
yvel.array(mfi);
60 const Box &dbx = geom.Domain();
61 Box sbx(bx.smallEnd(), bx.bigEnd());
63 AMREX_ALWAYS_ASSERT(sbx.smallEnd(2) == dbx.smallEnd(2) && sbx.bigEnd(2) == dbx.bigEnd(2));
65 const GeometryData gdata = geom.data();
67 const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
68 FArrayBox qintegral(xybx,2);
69 qintegral.setVal<RunOn::Device>(0.0);
70 FArrayBox qturb(bx,1); FArrayBox qturb_old(bx,1);
71 const Array4<Real> qint = qintegral.array();
72 const Array4<Real> qvel = qturb.array();
76 const Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
77 const auto invCellSize = geom.InvCellSizeArray();
78 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
80 qvel(i,j,k) = std::sqrt(2.0 * cell_data(i,j,k,
RhoKE_comp) / cell_data(i,j,k,
Rho_comp));
81 AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0,
"KE must have a positive value");
83 Real fac = (sbx.contains(i,j,k)) ? 1.0 : 0.0;
86 Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz*fac);
87 Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k)*dz*fac);
90 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
92 qvel(i,j,k) = std::sqrt(2.0 * cell_data(i,j,k,
RhoKE_comp) / cell_data(i,j,k,
Rho_comp));
93 AMREX_ASSERT_WITH_MESSAGE(qvel(i,j,k) > 0.0,
"KE must have a positive value");
97 Real fac = (sbx.contains(i,j,k)) ? 1.0 : 0.0;
98 const Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
99 Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*fac);
100 Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k)*fac);
104 Real dz_inv = geom.InvCellSize(2);
105 const auto& dxInv = geom.InvCellSizeArray();
106 int izmin = geom.Domain().smallEnd(2);
107 int izmax = geom.Domain().bigEnd(2);
110 Real d_kappa =
KAPPA;
113 const auto& t_mean_mf = most->get_mac_avg(level,4);
114 const auto& q_mean_mf = most->get_mac_avg(level,3);
115 const auto& u_star_mf = most->get_u_star(level);
116 const auto& t_star_mf = most->get_t_star(level);
117 const auto& q_star_mf = most->get_q_star(level);
119 const auto& tm_arr = t_mean_mf->const_array(mfi);
120 const auto& qm_arr = q_mean_mf->const_array(mfi);
121 const auto& u_star_arr = u_star_mf->const_array(mfi);
122 const auto& t_star_arr = t_star_mf->const_array(mfi);
123 const auto& q_star_arr = (use_moisture) ? q_star_mf->const_array(mfi) : Array4<Real>{};
125 const Array4<Real const> z_nd_arr = use_terrain ? z_phys_nd->const_array(mfi) : Array4<Real>{};
127 ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k) noexcept
137 Real dthetadz, dudz, dvdz;
139 uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta,
140 c_ext_dir_on_zlo, c_ext_dir_on_zhi,
141 u_ext_dir_on_zlo, u_ext_dir_on_zhi,
142 v_ext_dir_on_zlo, v_ext_dir_on_zhi,
143 dthetadz, dudz, dvdz,
144 RhoQv_comp, RhoQc_comp, RhoQr_comp, use_most);
147 Real theta0 = tm_arr(i,j,0);
148 Real qv0 = qm_arr(i,j,0);
149 Real surface_heat_flux = -u_star_arr(i,j,0) * t_star_arr(i,j,0);
150 Real surface_latent_heat{0};
153 surface_latent_heat = -u_star_arr(i,j,0) * q_star_arr(i,j,0);
154 surface_heat_flux *= (1.0 + 0.61*qv0);
155 surface_heat_flux += 0.61 * theta0 * surface_latent_heat;
159 if (std::abs(surface_heat_flux) > eps) {
160 l_obukhov = -( theta0 * u_star_arr(i,j,0)*u_star_arr(i,j,0)*u_star_arr(i,j,0) )
161 / ( d_kappa * d_gravity * surface_heat_flux );
163 l_obukhov = std::numeric_limits<Real>::max();
167 AMREX_ASSERT(l_obukhov != 0);
168 int lk = amrex::max(k,0);
170 : gdata.ProbLo(2) + (lk + 0.5)*gdata.CellSize(2);
171 const Real zeta = zval/l_obukhov;
174 l_S =
KAPPA*zval/3.7;
175 }
else if (zeta >= 0) {
176 l_S =
KAPPA*zval/(1+2.7*zeta);
178 l_S =
KAPPA*zval*std::pow(1.0 - 100.0 * zeta, 0.2);
183 if (qint(i,j,0,1) > 0.0) {
184 l_T = Lt_alpha*qint(i,j,0,0)/qint(i,j,0,1);
186 l_T = std::numeric_limits<Real>::max();
192 Real N_brunt_vaisala = std::sqrt(
CONST_GRAV/theta0 * dthetadz);
195 qc = std::pow(
qc,1.0/3.0);
196 l_B = (1.0 + 5.0*std::sqrt(
qc/(N_brunt_vaisala * l_T))) * qvel(i,j,k)/N_brunt_vaisala;
198 l_B = qvel(i,j,k) / N_brunt_vaisala;
201 l_B = std::numeric_limits<Real>::max();
207 Lm = std::pow(1.0/(l_S*l_S) + 1.0/(l_T*l_T) + 1.0/(l_B*l_B), -0.5);
210 Lm = 1.0 / (1.0/l_S + 1.0/l_T + 1.0/l_B);
214 Real shearProd = dudz*dudz + dvdz*dvdz;
215 Real buoyProd = -(
CONST_GRAV/theta0) * dthetadz;
216 Real L2_over_q2 = Lm*Lm/(qvel(i,j,k)*qvel(i,j,k));
217 Real GM = L2_over_q2 * shearProd;
218 Real GH = L2_over_q2 * buoyProd;
221 Real Rf = level2.calc_Rf(GM, GH);
222 Real SM2 = level2.calc_SM(Rf);
223 Real qe2 = mynn.B1*Lm*Lm*SM2*(1.0-Rf)*shearProd;
224 Real qe = (qe2 < 0.0) ? 0.0 : std::sqrt(qe2);
227 Real alphac = (qvel(i,j,k) > qe) ? 1.0 : qvel(i,j,k) / (qe + eps);
231 mynn.calc_stability_funcs(SM,SH,SQ,GM,GH,alphac);
234 SM = amrex::min(amrex::max(SM,mynn.SMmin), mynn.SMmax);
235 SH = amrex::min(amrex::max(SH,mynn.SHmin), mynn.SHmax);
250 if (mynn.diffuse_moistvars) {
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
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 int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp, const bool use_most=true)
Definition: ERF_PBLModels.H:82
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:39
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:352
@ yvel_bc
Definition: ERF_IndexDefines.H:89
@ cons_bc
Definition: ERF_IndexDefines.H:76
@ xvel_bc
Definition: ERF_IndexDefines.H:88
@ ext_dir
Definition: ERF_IndexDefines.H:180
@ Theta_v
Definition: ERF_IndexDefines.H:157
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:161
@ Q_v
Definition: ERF_IndexDefines.H:160
@ Mom_v
Definition: ERF_IndexDefines.H:156
@ KE_v
Definition: ERF_IndexDefines.H:158
@ rho
Definition: ERF_Kessler.H:30
@ qc
Definition: ERF_SatAdj.H:36
@ xvel
Definition: ERF_IndexDefines.H:130
@ yvel
Definition: ERF_IndexDefines.H:131
MYNNLevel25 pbl_mynn
Definition: ERF_TurbStruct.H:203
MYNNLevel2 pbl_mynn_level2
Definition: ERF_TurbStruct.H:204