ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Public Member Functions | Static Public Attributes | List of all members
MYNNPBLH Struct Reference

Diagnostic utility for the planetary boundary layer height. More...

#include <ERF_PBLHeight.H>

Collaboration diagram for MYNNPBLH:
Collaboration graph
[legend]

Public Member Functions

AMREX_GPU_HOST AMREX_FORCE_INLINE void compute_pblh (const amrex::Geometry &geom, const amrex::MultiFab *z_phys_cc, amrex::MultiFab *pblh, const amrex::MultiFab &cons, const amrex::iMultiFab *lmask, const MoistureComponentIndices &moisture_indices) const
 

Static Public Attributes

static constexpr amrex::Real thetamin_height = amrex::Real(200.0)
 Height below which minimum theta-v is determined [m]. More...
 
static constexpr amrex::Real theta_incr_land = amrex::Real(1.25)
 Theta increase determining the capping inversion height over land [K]. More...
 
static constexpr amrex::Real theta_incr_water = one
 Theta increase determining the capping inversion height over water [K]. More...
 
static constexpr amrex::Real sbl_lim = amrex::Real(200.0)
 Upper limit of the stable boundary layer height [m]. More...
 
static constexpr amrex::Real sbl_damp = amrex::Real(400.0)
 Transition length for blending [m]. More...
 

Detailed Description

Diagnostic utility for the planetary boundary layer height.

Member Function Documentation

◆ compute_pblh()

AMREX_GPU_HOST AMREX_FORCE_INLINE void MYNNPBLH::compute_pblh ( const amrex::Geometry &  geom,
const amrex::MultiFab *  z_phys_cc,
amrex::MultiFab *  pblh,
const amrex::MultiFab &  cons,
const amrex::iMultiFab *  lmask,
const MoistureComponentIndices &  moisture_indices 
) const
inline
31  {
32 #if 0
33  // NOTE: Cannot use ReduceToPlane because it clips the box to the
34  // validbox only, i.e., lateral ghost cells aren't updated
35  int dir = 2; // z
36  auto const& cons_arrs = cons.const_arrays();
37  auto thetav_min = amrex::ReduceToPlane<amrex::ReduceOpMin,amrex::Real>(dir, bxlow, cons,
38  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) -> amrex::Real
39  {
40  return GetThetav(i,j,k,cons_arrs[box_no],moisture_indices);
41  });
42 #endif
43 
44  // Create 2D multifabs like pblh
45  auto const& ba = pblh->boxArray();
46  auto const& dm = pblh->DistributionMap();
47  auto const& ng = pblh->nGrowVect();
48 
49  amrex::MultiFab min_thetav(ba,dm,1,ng);
50  min_thetav.setVal(bogus_large_value);
51 
52  amrex::MultiFab pblh_tke(ba,dm,1,ng);
53  pblh_tke.setVal(zero);
54 
55  pblh->setVal(zero);
56 
57  // Now, loop over columns...
58  for (amrex::MFIter mfi(cons,TileNoZ()); mfi.isValid(); ++mfi)
59  {
60  const amrex::Box& domain = geom.Domain();
61 
62  // NOTE: this MFIter is over cons, which carries more ghost cells than the
63  // 2D fields we write here (pblh, pblh_tke, min_thetav all have the
64  // ghost cells of pblh). Grow by *their* ghost vector, not by the
65  // default, or we write past the end of their allocations.
66  amrex::Box gtbx = mfi.growntilebox(ng);
67  gtbx.setSmall(2,domain.smallEnd(2)); // don't loop over ghost cells
68  gtbx.setBig(2,domain.bigEnd(2)); // in z
69 
70  // Columns are parallelized over; the k loops below are serialized within
71  // each thread. Everything in this routine is a per-column scan -- the
72  // minimum theta_v, and the lowest k at which each PBLH criterion is met --
73  // so parallelizing over k would both race on the 2D output and lose the
74  // "lowest k wins" semantics.
75  amrex::Box gtbx2d = gtbx; gtbx2d.setRange(2,0);
76  const int klo = gtbx.smallEnd(2);
77  const int khi = gtbx.bigEnd(2);
78 
79  auto min_thv_arr = min_thetav.array(mfi);
80  auto pblh_arr = pblh->array(mfi);
81  auto pblh_tke_arr = pblh_tke.array(mfi);
82 
83  const auto cons_arr = cons.const_array(mfi);
84  const auto lmask_arr = (lmask) ? lmask->const_array(mfi) : amrex::Array4<int> {};
85 
86  // -----------------------------------------------------
87  // WITH terrain/grid stretching
88  // -----------------------------------------------------
89  if (z_phys_cc)
90  {
91  const auto zphys_arr = z_phys_cc->const_array(mfi);
92 
93  // Need to sort out ghost cell differences (z_phys_cc has ng=1)
94  int imin = lbound(zphys_arr).x;
95  int jmin = lbound(zphys_arr).y;
96  int imax = ubound(zphys_arr).x;
97  int jmax = ubound(zphys_arr).y;
98 
99  // Find minimum thetav in the surface layer (this updates
100  // ghost cells, too)
101  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
102  {
103  int ii = amrex::max(amrex::min(i,imax),imin);
104  int jj = amrex::max(amrex::min(j,jmax),jmin);
105 
106  amrex::Real min_thv = min_thv_arr(i,j,0);
107  for (int k(klo); k <= khi; ++k) {
108  if (zphys_arr(ii,jj,k) < thetamin_height) {
109  amrex::Real thv = GetThetav(i, j, k, cons_arr, moisture_indices);
110  min_thv = amrex::min(min_thv, thv);
111  }
112  }
113  min_thv_arr(i,j,0) = min_thv;
114  });
115 
116  // This depends on TileNoZ and k increasing monotonically
117  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
118  {
119  int ii = amrex::max(amrex::min(i,imax),imin);
120  int jj = amrex::max(amrex::min(j,jmax),jmin);
121 
122  const amrex::Real min_thv = min_thv_arr(i,j,0);
123  const int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
124 
125  amrex::Real zi = pblh_arr(i,j,0);
126  amrex::Real zi_tke = pblh_tke_arr(i,j,0);
127 
128  for (int k(klo); k <= khi; ++k)
129  {
130  if (zi == 0)
131  {
132  //
133  // Find PBL height based on thetav increase (best for CBLs)
134  //
135  amrex::Real thv = GetThetav(i, j, k , cons_arr, moisture_indices);
136  amrex::Real thv1 = GetThetav(i, j, k+1, cons_arr, moisture_indices);
137 
138  if (is_land && (thv1 >= min_thv + theta_incr_land)
139  && (thv < min_thv + theta_incr_land))
140  {
141  // Interpolate to get lowest height where theta = min_theta + theta_incr
142  zi = zphys_arr(ii,jj,k)
143  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
144  * (min_thv + theta_incr_land - thv);
145  }
146  else if (!is_land && (thv1 >= min_thv + theta_incr_water)
147  && (thv < min_thv + theta_incr_water))
148  {
149  // Interpolate to get lowest height where theta = min_theta + theta_incr
150  zi = zphys_arr(ii,jj,k)
151  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
152  * (min_thv + theta_incr_water - thv);
153  }
154  }
155  if (zi_tke == 0)
156  {
157  //
158  // Find PBL height based on TKE (for SBLs only)
159  //
160  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
161  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
162  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
163  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
164  amrex::Real TKEeps = amrex::Real(0.05) * maxtke;
165  TKEeps = amrex::max(TKEeps, amrex::Real(0.02)); // min val from WRF
166 
167  if ((tke1 <= TKEeps) && (tke > TKEeps))
168  {
169  // Interpolate to get lowest height where TKE -> 0
170  zi_tke = zphys_arr(ii,jj,k)
171  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(tke1-tke)
172  * (TKEeps - tke);
173  }
174  }
175  // Both criteria are satisfied; the remaining levels are no-ops
176  if ((zi != 0) && (zi_tke != 0)) { break; }
177  }
178 
179  pblh_arr(i,j,0) = zi;
180  pblh_tke_arr(i,j,0) = zi_tke;
181  });
182  }
183  else
184  // -----------------------------------------------------
185  // NO terrain
186  // -----------------------------------------------------
187  {
188  const amrex::Real dz_no_terrain = geom.CellSize(2);
189 
190  // Find minimum thetav in the surface layer (this updates
191  // ghost cells, too)
192  // - box size is known a priori
193  int kmax = static_cast<int>(thetamin_height / dz_no_terrain);
194  AMREX_ASSERT(kmax > 0);
195  const int khi_low = amrex::min(kmax, khi);
196 
197  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
198  {
199  amrex::Real min_thv = min_thv_arr(i,j,0);
200  for (int k(klo); k <= khi_low; ++k) {
201  amrex::Real thv = GetThetav(i, j, k, cons_arr, moisture_indices);
202  min_thv = amrex::min(min_thv, thv);
203  }
204  min_thv_arr(i,j,0) = min_thv;
205  });
206 
207  // This depends on TileNoZ and k increasing monotonically
208  ParallelFor(gtbx2d, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
209  {
210  const amrex::Real min_thv = min_thv_arr(i,j,0);
211  const int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
212 
213  amrex::Real zi = pblh_arr(i,j,0);
214  amrex::Real zi_tke = pblh_tke_arr(i,j,0);
215 
216  for (int k(klo); k <= khi; ++k)
217  {
218  if (zi == 0)
219  {
220  //
221  // Find PBL height based on thetav increase (best for CBLs)
222  //
223  amrex::Real thv = GetThetav(i, j, k , cons_arr, moisture_indices);
224  amrex::Real thv1 = GetThetav(i, j, k+1, cons_arr, moisture_indices);
225 
226  if (is_land && (thv1 >= min_thv + theta_incr_land)
227  && (thv < min_thv + theta_incr_land))
228  {
229  // Interpolate to get lowest height where theta = min_theta + theta_incr
230  zi = (k+myhalf)*dz_no_terrain
231  + dz_no_terrain/(thv1-thv)
232  * (min_thv + theta_incr_land - thv);
233  }
234  else if (!is_land && (thv1 >= min_thv + theta_incr_water)
235  && (thv < min_thv + theta_incr_water))
236  {
237  // Interpolate to get lowest height where theta = min_theta + theta_incr
238  zi = (k+myhalf)*dz_no_terrain
239  + dz_no_terrain/(thv1-thv)
240  * (min_thv + theta_incr_water - thv);
241  }
242  }
243  if (zi_tke == 0)
244  {
245  //
246  // Find PBL height based on TKE (for SBLs only)
247  //
248  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
249  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
250  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
251  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
252  amrex::Real TKEeps = amrex::Real(0.05) * maxtke;
253  TKEeps = amrex::max(TKEeps, amrex::Real(0.02)); // min val from WRF
254 
255  if ((tke1 <= TKEeps) && (tke > TKEeps))
256  {
257  // Interpolate to get lowest height where TKE -> 0
258  zi_tke = (k+myhalf)*dz_no_terrain
259  + dz_no_terrain/(tke1-tke) * (TKEeps - tke);
260  }
261  }
262  // Both criteria are satisfied; the remaining levels are no-ops
263  if ((zi != 0) && (zi_tke != 0)) { break; }
264  }
265 
266  pblh_arr(i,j,0) = zi;
267  pblh_tke_arr(i,j,0) = zi_tke;
268  });
269  }
270  }// MFIter
271 
272  //
273  // Calculate hybrid PBL height
274  //
275  for (amrex::MFIter mfi(*pblh); mfi.isValid(); ++mfi)
276  {
277  const auto cons_arr = cons.const_array(mfi);
278  auto pblh_tke_arr = pblh_tke.array(mfi);
279  auto pblh_arr = pblh->array(mfi);
280 
281  amrex::Box gtbx = mfi.growntilebox();
282  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
283  {
284  //
285  // Clip PBLH_TKE to more realistic values
286  //
287  // Note from WRF MYNN-EDMF: TKE-based PBLH can be very large in cells
288  // with convective precipitation (> 8km!), so an artificial limit is
289  // imposed to not let PBLH_TKE exceed the theta_v-based PBL height
290  // +/- 350 m. This has no impact on 98-99% of the domain, but is the
291  // simplest patch that adequately addresses these extremely large
292  // PBLHs.
293  amrex::Real zi = pblh_arr(i,j,0);
294  pblh_tke_arr(i,j,0) = amrex::max(
295  amrex::min(pblh_tke_arr(i,j,0), zi+amrex::Real(350.)),
296  amrex::max(zi-amrex::Real(350.), amrex::Real(10.)));
297 
298  //
299  // Finally, blend between the two PBLH estimates
300  //
301  amrex::Real maxqke = two * cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
302  if (maxqke > amrex::Real(0.05)) {
303  amrex::Real wt = myhalf*std::tanh((zi - sbl_lim)/sbl_damp) + myhalf;
304  pblh_arr(i,j,0) = (one-wt)*pblh_tke_arr(i,j,0) + wt*pblh_arr(i,j,0);
305  }
306  });
307  }//MFIter
308  }
two
constexpr amrex::Real two
Definition: ERF_Constants.H:10
bogus_large_value
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:26
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:41
khi
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
GetThetav
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetav(const int &i, const int &j, const int &k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:72
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
MicVar_Morr::ng
@ ng
Definition: ERF_Morrison.H:49
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:214
WDM6SedNodeScratch::zi
@ zi
Definition: ERF_AdvanceWDM6.cpp:274
MYNNPBLH::sbl_lim
static constexpr amrex::Real sbl_lim
Upper limit of the stable boundary layer height [m].
Definition: ERF_PBLHeight.H:325
MYNNPBLH::sbl_damp
static constexpr amrex::Real sbl_damp
Transition length for blending [m].
Definition: ERF_PBLHeight.H:326
MYNNPBLH::theta_incr_water
static constexpr amrex::Real theta_incr_water
Theta increase determining the capping inversion height over water [K].
Definition: ERF_PBLHeight.H:324
MYNNPBLH::thetamin_height
static constexpr amrex::Real thetamin_height
Height below which minimum theta-v is determined [m].
Definition: ERF_PBLHeight.H:322
MYNNPBLH::theta_incr_land
static constexpr amrex::Real theta_incr_land
Theta increase determining the capping inversion height over land [K].
Definition: ERF_PBLHeight.H:323
Here is the call graph for this function:

Member Data Documentation

◆ sbl_damp

MYNNPBLH::sbl_damp = amrex::Real(400.0)
staticconstexpr

Transition length for blending [m].

Referenced by compute_pblh().

◆ sbl_lim

MYNNPBLH::sbl_lim = amrex::Real(200.0)
staticconstexpr

Upper limit of the stable boundary layer height [m].

Referenced by compute_pblh().

◆ theta_incr_land

MYNNPBLH::theta_incr_land = amrex::Real(1.25)
staticconstexpr

Theta increase determining the capping inversion height over land [K].

Referenced by compute_pblh().

◆ theta_incr_water

MYNNPBLH::theta_incr_water = one
staticconstexpr

Theta increase determining the capping inversion height over water [K].

Referenced by compute_pblh().

◆ thetamin_height

MYNNPBLH::thetamin_height = amrex::Real(200.0)
staticconstexpr

Height below which minimum theta-v is determined [m].

Referenced by compute_pblh().

The documentation for this struct was generated from the following file: