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

#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 int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp) const
 

Static Public Attributes

static constexpr amrex::Real thetamin_height = 200.0
 
static constexpr amrex::Real theta_incr_land = 1.25
 
static constexpr amrex::Real theta_incr_water = 1.0
 
static constexpr amrex::Real sbl_lim = 200.0
 
static constexpr amrex::Real sbl_damp = 400.0
 

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 int  RhoQv_comp,
const int  RhoQc_comp,
const int  RhoQr_comp 
) const
inline
30  {
31 #if 0
32  // NOTE: Cannot use ReduceToPlane because it clips the box to the
33  // validbox only, i.e., lateral ghost cells aren't updated
34  int dir = 2; // z
35  auto const& cons_arrs = cons.const_arrays();
36  auto thetav_min = amrex::ReduceToPlane<amrex::ReduceOpMin,amrex::Real>(dir, bxlow, cons,
37  [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) -> amrex::Real
38  {
39  return Thetav(i,j,k,cons_arrs[box_no],RhoQv_comp,RhoQc_comp,RhoQr_comp);
40  });
41 #endif
42 
43  // Create 2D multifabs like pblh
44  auto const& ba = pblh->boxArray();
45  auto const& dm = pblh->DistributionMap();
46  auto const& ng = pblh->nGrowVect();
47 
48  amrex::MultiFab min_thetav(ba,dm,1,ng);
49  min_thetav.setVal(1.E34);
50 
51  amrex::MultiFab pblh_tke(ba,dm,1,ng);
52  pblh_tke.setVal(0);
53 
54  pblh->setVal(0);
55 
56  // Now, loop over columns...
57  for (amrex::MFIter mfi(cons,TileNoZ()); mfi.isValid(); ++mfi)
58  {
59  const amrex::Box& domain = geom.Domain();
60  amrex::Box gtbx = mfi.growntilebox();
61  gtbx.setSmall(2,domain.smallEnd(2)); // don't loop over ghost cells
62  gtbx.setBig(2,domain.bigEnd(2)); // in z
63 
64  auto min_thv_arr = min_thetav.array(mfi);
65  auto pblh_arr = pblh->array(mfi);
66  auto pblh_tke_arr = pblh_tke.array(mfi);
67 
68  const auto cons_arr = cons.const_array(mfi);
69  const auto lmask_arr = (lmask) ? lmask->const_array(mfi) : amrex::Array4<int> {};
70 
71  // -----------------------------------------------------
72  // WITH terrain/grid stretching
73  // -----------------------------------------------------
74  if (z_phys_cc)
75  {
76  const auto zphys_arr = z_phys_cc->const_array(mfi);
77 
78  // Need to sort out ghost cell differences (z_phys_cc has ng=1)
79  int imin = lbound(zphys_arr).x;
80  int jmin = lbound(zphys_arr).y;
81  int imax = ubound(zphys_arr).x;
82  int jmax = ubound(zphys_arr).y;
83 
84  // Find minimum thetav in the surface layer (this updates
85  // ghost cells, too)
86  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
87  {
88  int ii = amrex::max(amrex::min(i,imax),imin);
89  int jj = amrex::max(amrex::min(j,jmax),jmin);
90 
91  if (zphys_arr(ii,jj,k) < thetamin_height) {
92  amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
93  if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv;
94  }
95  });
96 
97  // This depends on TileNoZ and k increasing monotonically
98  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
99  {
100  int ii = amrex::max(amrex::min(i,imax),imin);
101  int jj = amrex::max(amrex::min(j,jmax),jmin);
102 
103  if (pblh_arr(i,j,0) == 0)
104  {
105  //
106  // Find PBL height based on thetav increase (best for CBLs)
107  //
108  amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
109  amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
110 
111  int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
112  if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land)
113  && (thv < min_thv_arr(i,j,0) + theta_incr_land))
114  {
115  // Interpolate to get lowest height where theta = min_theta + theta_incr
116  pblh_arr(i,j,0) = zphys_arr(ii,jj,k)
117  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
118  * (min_thv_arr(i,j,0) + theta_incr_land - thv);
119  }
120  else if (!is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_water)
121  && (thv < min_thv_arr(i,j,0) + theta_incr_water))
122  {
123  // Interpolate to get lowest height where theta = min_theta + theta_incr
124  pblh_arr(i,j,0) = zphys_arr(ii,jj,k)
125  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(thv1-thv)
126  * (min_thv_arr(i,j,0) + theta_incr_water - thv);
127  }
128  }
129  if (pblh_tke_arr(i,j,0) == 0)
130  {
131  //
132  // Find PBL height based on TKE (for SBLs only)
133  //
134  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
135  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
136  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
137  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
138  amrex::Real TKEeps = 0.05 * maxtke;
139  TKEeps = amrex::max(TKEeps, 0.02); // min val from WRF
140 
141  if ((tke1 <= TKEeps) && (tke > TKEeps))
142  {
143  // Interpolate to get lowest height where TKE -> 0
144  pblh_tke_arr(i,j,0) = zphys_arr(ii,jj,k)
145  + (zphys_arr(ii,jj,k+1)-zphys_arr(ii,jj,k))/(tke1-tke)
146  * (TKEeps - tke);
147  }
148  }
149  });
150  }
151  else
152  // -----------------------------------------------------
153  // NO terrain
154  // -----------------------------------------------------
155  {
156  const amrex::Real dz_no_terrain = geom.CellSize(2);
157 
158  // Find minimum thetav in the surface layer (this updates
159  // ghost cells, too)
160  // - box size is known a priori
161  int kmax = static_cast<int>(thetamin_height / dz_no_terrain);
162  AMREX_ASSERT(kmax > 0);
163  amrex::Box gtbxlow = gtbx;
164  gtbxlow.setBig(2,kmax);
165 
166  ParallelFor(gtbxlow, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
167  {
168  amrex::Real thv = Thetav(i,j,k,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
169  if (min_thv_arr(i,j,0) > thv) min_thv_arr(i,j,0) = thv;
170  });
171 
172  // This depends on TileNoZ and k increasing monotonically
173  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
174  {
175  if (pblh_arr(i,j,0) == 0)
176  {
177  //
178  // Find PBL height based on thetav increase (best for CBLs)
179  //
180  amrex::Real thv = Thetav(i,j,k ,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
181  amrex::Real thv1 = Thetav(i,j,k+1,cons_arr,RhoQv_comp,RhoQc_comp,RhoQr_comp);
182 
183  int is_land = (lmask_arr) ? lmask_arr(i,j,0) : 1;
184  if (is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_land)
185  && (thv < min_thv_arr(i,j,0) + theta_incr_land))
186  {
187  // Interpolate to get lowest height where theta = min_theta + theta_incr
188  pblh_arr(i,j,0) = (k+0.5)*dz_no_terrain
189  + dz_no_terrain/(thv1-thv)
190  * (min_thv_arr(i,j,0) + theta_incr_land - thv);
191  }
192  else if (!is_land && (thv1 >= min_thv_arr(i,j,0) + theta_incr_water)
193  && (thv < min_thv_arr(i,j,0) + theta_incr_water))
194  {
195  // Interpolate to get lowest height where theta = min_theta + theta_incr
196  pblh_arr(i,j,0) = (k+0.5)*dz_no_terrain
197  + dz_no_terrain/(thv1-thv)
198  * (min_thv_arr(i,j,0) + theta_incr_water - thv);
199  }
200  }
201  if (pblh_tke_arr(i,j,0) == 0)
202  {
203  //
204  // Find PBL height based on TKE (for SBLs only)
205  //
206  amrex::Real tke = cons_arr(i,j,k ,RhoKE_comp) / cons_arr(i,j,k ,Rho_comp);
207  amrex::Real tke1 = cons_arr(i,j,k+1,RhoKE_comp) / cons_arr(i,j,k+1,Rho_comp);
208  amrex::Real maxtke = cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
209  // - threshold is 5% of max TKE (Kosovic & Curry 2000, JAS)
210  amrex::Real TKEeps = 0.05 * maxtke;
211  TKEeps = amrex::max(TKEeps, 0.02); // min val from WRF
212 
213  if ((tke1 <= TKEeps) && (tke > TKEeps))
214  {
215  // Interpolate to get lowest height where TKE -> 0
216  pblh_tke_arr(i,j,0) = (k+0.5)*dz_no_terrain
217  + dz_no_terrain/(tke1-tke) * (TKEeps - tke);
218  }
219  }
220  });
221  }
222  }// MFIter
223 
224  //
225  // Calculate hybrid PBL height
226  //
227  for (amrex::MFIter mfi(*pblh); mfi.isValid(); ++mfi)
228  {
229  const auto cons_arr = cons.const_array(mfi);
230  auto pblh_tke_arr = pblh_tke.array(mfi);
231  auto pblh_arr = pblh->array(mfi);
232 
233  amrex::Box gtbx = mfi.growntilebox();
234  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept
235  {
236  //
237  // Clip PBLH_TKE to more realistic values
238  //
239  // Note from WRF MYNN-EDMF: TKE-based PBLH can be very large in cells
240  // with convective precipiation (> 8km!), so an artificial limit is
241  // imposed to not let PBLH_TKE exceed the theta_v-based PBL height
242  // +/- 350 m. This has no impact on 98-99% of the domain, but is the
243  // simplest patch that adequately addresses these extremely large
244  // PBLHs.
245  amrex::Real zi = pblh_arr(i,j,0);
246  pblh_tke_arr(i,j,0) = amrex::max(
247  amrex::min(pblh_tke_arr(i,j,0), zi+350.),
248  amrex::max(zi-350., 10.));
249 
250  //
251  // Finally, blend between the two PBLH estimates
252  //
253  amrex::Real maxqke = 2.0 * cons_arr(i,j,0 ,RhoKE_comp) / cons_arr(i,j,0 ,Rho_comp);
254  if (maxqke > 0.05) {
255  amrex::Real wt = 0.5*std::tanh((zi - sbl_lim)/sbl_damp) + 0.5;
256  pblh_arr(i,j,0) = (1.-wt)*pblh_tke_arr(i,j,0) + wt*pblh_arr(i,j,0);
257  }
258  });
259  }//MFIter
260  }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
Thetav
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Thetav(int i, int j, int k, const amrex::Array4< const amrex::Real > &cell_data, const int RhoQv_comp, const int RhoQc_comp, const int RhoQr_comp)
Definition: ERF_Thetav.H:10
TileNoZ
AMREX_FORCE_INLINE amrex::IntVect TileNoZ()
Definition: ERF_TileNoZ.H:11
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:129
MYNNPBLH::sbl_lim
static constexpr amrex::Real sbl_lim
Definition: ERF_PBLHeight.H:265
MYNNPBLH::theta_incr_land
static constexpr amrex::Real theta_incr_land
Definition: ERF_PBLHeight.H:263
MYNNPBLH::thetamin_height
static constexpr amrex::Real thetamin_height
Definition: ERF_PBLHeight.H:262
MYNNPBLH::theta_incr_water
static constexpr amrex::Real theta_incr_water
Definition: ERF_PBLHeight.H:264
MYNNPBLH::sbl_damp
static constexpr amrex::Real sbl_damp
Definition: ERF_PBLHeight.H:266
Here is the call graph for this function:

Member Data Documentation

◆ sbl_damp

constexpr amrex::Real MYNNPBLH::sbl_damp = 400.0
staticconstexpr

Referenced by compute_pblh().

◆ sbl_lim

constexpr amrex::Real MYNNPBLH::sbl_lim = 200.0
staticconstexpr

Referenced by compute_pblh().

◆ theta_incr_land

constexpr amrex::Real MYNNPBLH::theta_incr_land = 1.25
staticconstexpr

Referenced by compute_pblh().

◆ theta_incr_water

constexpr amrex::Real MYNNPBLH::theta_incr_water = 1.0
staticconstexpr

Referenced by compute_pblh().

◆ thetamin_height

constexpr amrex::Real MYNNPBLH::thetamin_height = 200.0
staticconstexpr

Referenced by compute_pblh().

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