ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_MetgridUtils.H
Go to the documentation of this file.
1 /**
2  * \file ERF_MetgridUtils.H
3  */
4 #ifndef ERF_METGRIDUTIL_H_
5 #define ERF_METGRIDUTIL_H_
6 
7 #include <cmath>
8 #include <limits>
9 
10 #include <ERF.H>
11 #include <ERF_EOS.H>
12 #include <ERF_Utils.H>
13 #include <ERF_ProbCommon.H>
14 #include <ERF_HSEUtils.H>
15 #include <ERF_BaseStateParams.H>
16 
17 /**
18  * Read metgrid NetCDF data and populate ERF initialization FABs.
19  */
20 void
21 read_from_metgrid (int lev,
22  int itime,
23  const amrex::Box& domain,
24  const std::string& fname,
25  std::string& NC_dateTime,
26  double& NC_epochTime,
27  int& flag_psfc,
28  int& flag_msf,
29  int& flag_sst,
30  int& flag_tsk,
31  int& flag_lmask,
32  int& NC_nx,
33  int& NC_ny,
34  amrex::Real& NC_dx,
35  amrex::Real& NC_dy,
36  amrex::FArrayBox& NC_xvel_fab,
37  amrex::FArrayBox& NC_yvel_fab,
38  amrex::FArrayBox& NC_temp_fab,
39  amrex::FArrayBox& NC_rhum_fab,
40  amrex::FArrayBox& NC_pres_fab,
41  amrex::FArrayBox& NC_ght_fab,
42  amrex::FArrayBox& NC_hgt_fab,
43  amrex::FArrayBox& NC_psfc_fab,
44  amrex::FArrayBox& NC_msfu_fab,
45  amrex::FArrayBox& NC_msfv_fab,
46  amrex::FArrayBox& NC_msfm_fab,
47  amrex::FArrayBox& NC_sst_fab,
48  amrex::FArrayBox& NC_tsk_fab,
49  amrex::FArrayBox& NC_LAT_fab,
50  amrex::FArrayBox& NC_LON_fab,
51  amrex::IArrayBox& NC_lmask_iab,
52  amrex::Geometry& geom);
53 
54 /**
55  * Initialize nodal terrain heights from metgrid surface-height data.
56  */
57 void
58 init_terrain_from_metgrid (const bool& avg_grid_faces_to_nodes,
59  amrex::Geometry& geom,
60  amrex::MultiFab* z_phys_nd,
61  amrex::FArrayBox& NC_hgt_fab);
62 
63 /**
64  * Initialize ERF state, velocity, and boundary data from metgrid fields.
65  */
66 void
67 init_state_from_metgrid (const int lev,
68  const int itime,
69  const bool use_moisture,
70  const bool interp_theta,
71  const bool metgrid_debug_quiescent,
72  const bool metgrid_debug_isothermal,
73  const bool metgrid_debug_dry,
74  const bool metgrid_basic_linear,
75  const bool metgrid_use_below_sfc,
76  const bool metgrid_use_sfc,
77  const bool metgrid_retain_sfc,
78  const amrex::Real metgrid_proximity,
79  const int metgrid_order,
80  const int metgrid_metgrid_force_sfc_k,
81  const amrex::Real l_rdOcp,
82  amrex::Box& tbxc,
83  amrex::Box& tbxu,
84  amrex::Box& tbxv,
85  amrex::Box& tbxw,
86  amrex::FArrayBox& state_fab,
87  amrex::FArrayBox& x_vel_fab,
88  amrex::FArrayBox& y_vel_fab,
89  amrex::FArrayBox& z_vel_fab,
90  amrex::FArrayBox& z_phys_nd_fab,
91  const amrex::FArrayBox& NC_ght_fab,
92  const amrex::FArrayBox& NC_xvel_fab,
93  const amrex::FArrayBox& NC_yvel_fab,
94  const amrex::FArrayBox& NC_temp_fab,
95  const amrex::FArrayBox& NC_rhum_fab,
96  const amrex::FArrayBox& NC_pres_fab,
97  amrex::FArrayBox& tmp_src_fab,
98  amrex::FArrayBox& tmp_dst_fab,
99  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_xlo,
100  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_xhi,
101  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_ylo,
102  amrex::Vector<amrex::Vector<amrex::FArrayBox>>& fabs_for_bcs_yhi,
103  const amrex::Array4<const int>& mask_c_arr,
104  const amrex::Array4<const int>& mask_u_arr,
105  const amrex::Array4<const int>& mask_v_arr);
106 
107 /**
108  * Initialize map-scale factors from metgrid fields or defaults.
109  */
110 void
111 init_msfs_from_metgrid (const bool metgrid_debug_msf,
112  amrex::FArrayBox& msfu_fab,
113  amrex::FArrayBox& msfv_fab,
114  amrex::FArrayBox& msfm_fab,
115  const int& flag_msf,
116  amrex::FArrayBox& NC_MSFU_fab,
117  amrex::FArrayBox& NC_MSFV_fab,
118  amrex::FArrayBox& NC_MSFM_fab);
119 
120 
121 /**
122  * Initialize hydrostatic base-state quantities from metgrid data.
123  */
124 void
125 init_base_state_from_metgrid (const bool use_moisture,
126  const bool metgrid_debug_psfc,
127  const amrex::Real l_rdOcp,
128  const amrex::Box& domain,
129  const amrex::Box& valid_bx,
130  const int& flag_psfc,
131  amrex::FArrayBox& state_fab,
132  amrex::FArrayBox& r_hse_fab,
133  amrex::FArrayBox& p_hse_fab,
134  amrex::FArrayBox& pi_hse_fab,
135  amrex::FArrayBox& th_hse_fab,
136  amrex::FArrayBox& qv_hse_fab,
137  amrex::FArrayBox& z_phys_nd_fab,
138  amrex::FArrayBox& z_phys_cc_fab,
139  const amrex::FArrayBox& NC_psfc_fab,
140  const BaseStateParams& bsp);
141 
142 /**
143  * Interpolate a value with a Lagrange polynomial stencil.
144  */
145 AMREX_FORCE_INLINE
146 AMREX_GPU_DEVICE
147 void
148 lagrange_interp (const int& order,
149  amrex::Real* x,
150  amrex::Real* y,
151  amrex::Real& new_x,
152  amrex::Real& new_y)
153 {
154  // Interpolation using Lagrange polynomials.
155  // P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x)
156  // where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn)
157  // ---------------------------------------------
158  // (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn)
159  amrex::Real Px = zero;
160  for (int i=0; i <= order; i++) {
161  amrex::Real n = one;
162  amrex::Real d = one;
163  for (int k=0; k <= order; k++) {
164  if (k == i) continue;
165  n *= new_x-x[k];
166  d *= x[i]-x[k];
167  }
168  if (d != zero) {
169  Px += y[i]*n/d;
170  }
171  }
172  new_y = Px;
173 }
174 
175 /**
176  * Prepare interpolation stencils and apply Lagrange interpolation or extrapolation.
177  *
178  * Both orig_npts and new_npts are point counts, not last indices; valid indices
179  * into the orig_* arrays are 0 through orig_npts-1, and into the new_* arrays
180  * are 0 through new_npts-1.
181  */
182 AMREX_FORCE_INLINE
183 AMREX_GPU_DEVICE
184 void
185 lagrange_setup (char var_type,
186  const bool& exp_interp,
187  const int& orig_npts,
188  const int& new_npts,
189  const int& order,
190  const int& i,
191  const int& j,
192  amrex::Real* orig_x_z,
193  amrex::Real* orig_x_p,
194  amrex::Real* orig_y,
195  amrex::Real* new_x_z,
196  amrex::Real* new_x_p,
197  amrex::Real* new_y)
198 {
199 
200  amrex::Real CRC_const1 = amrex::Real(11880.516); // m
201  amrex::Real CRC_const2 = amrex::Real(0.1902632);
202  amrex::Real CRC_const3 = amrex::Real(0.0065); // K km-1
203 
204  amrex::ignore_unused(i,j);
205 #ifndef AMREX_USE_GPU
206  bool debug = false;
207 #endif
208 
209  for (int new_k=0; new_k < new_npts; new_k++) {
210 #ifndef AMREX_USE_GPU
211  if (debug) amrex::Print() << "new_k=" << new_k;
212 #endif
213  // Determine if interpolating or extrapolating.
214  // If interpolating, find bounding x values and store the indices.
215  bool extrapolating = true;
216  int kl, kr;
217  for (int ko=0; ko < orig_npts-1; ko++) {
218  amrex::Real a = new_x_z[new_k]-orig_x_z[ko];
219  amrex::Real b = new_x_z[new_k]-orig_x_z[ko+1];
220  if (a*b <= zero) {
221  kl = ko;
222  kr = ko+1;
223  extrapolating = false;
224  break;
225  }
226  }
227 
228  if (extrapolating) {
229  if (var_type == 'T') {
230  // Assume a standard atmosphere -amrex::Real(6.5) K km-1 lapse rate.
231  // Comparable to the WRF default, t_extrap_type=two
232  amrex::Real depth_of_extrap_in_p = new_x_p[new_k]-orig_x_p[0];
233  amrex::Real avg_of_extrap_p = myhalf*(new_x_p[new_k]+orig_x_p[0]);
234  amrex::Real temp_extrap_starting_point = orig_y[0]*std::pow(orig_x_p[0]/amrex::Real(100000.0), RdoCp);
235  amrex::Real dZdP = CRC_const1*CRC_const2*std::pow(avg_of_extrap_p/amrex::Real(100.0), CRC_const2-one);
236  amrex::Real dZ = dZdP*(depth_of_extrap_in_p/amrex::Real(100.0));
237  amrex::Real dT = dZ*CRC_const3;
238  new_y[new_k] = (temp_extrap_starting_point+dT)*std::pow(amrex::Real(100000.0)/new_x_p[new_k], RdoCp);
239  } else {
240  // Use a constant value below ground.
241  // Comparable to the WRF default, extrap_type=two
242  new_y[new_k] = orig_y[0];
243  }
244  continue;
245  }
246 
247  if (order%2 != 0) {
248  if ((kl-((order+1)/2-1) >= 0) && (kr+((order+1)/2-1) <= orig_npts-1)) {
249  // Odd order interpolation.
250  int ksta = kl-(((order+1)/2)-1);
251  int kend = ksta+order;
252 #ifndef AMREX_USE_GPU
253  int ksize = kend-ksta;
254  if (debug) amrex::Print() << " (1a) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
255 #endif
256  amrex::Real new_x;
257  amrex::GpuArray<amrex::Real,9> orig_x_sub;
258  amrex::Real* orig_x_sub_p = orig_x_sub.data();
259  if (exp_interp) {
260  new_x = new_x_p[new_k];
261  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
262  } else {
263  new_x = new_x_z[new_k];
264  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
265  }
266  amrex::GpuArray<amrex::Real,9> orig_y_sub;
267  amrex::Real* orig_y_sub_p = orig_y_sub.data();
268  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
269 #ifndef AMREX_USE_GPU
270  if (debug) {
271  amrex::Print() << " orig_x_sub_p = [";
272  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
273  amrex::Print() << "]" << std::endl;
274  amrex::Print() << " orig_y_sub_p = [";
275  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
276  amrex::Print() << "]" << std::endl;
277  }
278 #endif
279  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
280  } else {
281  // Interpolation stencil is too big due to proximity to the surface or model top.
282  // We have bounding points so resort to a simple linear interpolation.
283  int ksta = kl;
284  int kend = kr;
285 #ifndef AMREX_USE_GPU
286  int ksize = kend-ksta+1;
287  if (debug) amrex::Print() << " (1b) order=" << order << " new_x=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
288 #endif
289  amrex::Real new_x;
290  amrex::GpuArray<amrex::Real,2> orig_x_sub;
291  amrex::GpuArray<amrex::Real,2> orig_y_sub;
292  amrex::Real* orig_x_sub_p = orig_x_sub.data();
293  amrex::Real* orig_y_sub_p = orig_y_sub.data();
294  if (exp_interp) {
295  new_x = new_x_p[new_k];
296  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
297  } else {
298  new_x = new_x_z[new_k];
299  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
300  }
301  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
302 #ifndef AMREX_USE_GPU
303  if (debug) {
304  amrex::Print() << " orig_x_sub_p = [";
305  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
306  amrex::Print() << "]" << std::endl;
307  amrex::Print() << " orig_y_sub_p = [";
308  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
309  amrex::Print() << "]" << std::endl;
310  }
311 #endif
312  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
313  }
314  } else if (order%2 == 0) {
315  if ((kl-(order/2) >= 0) && (kr+order/2 <= orig_npts-1)) {
316  // Even order interpolation.
317  // Interpolate twice. Offset the range by 1 the 2nd time. Average the results.
318  amrex::Real new_y_l, new_y_r;
319  {
320  int ksta = kl-(order/2-1);
321  int kend = ksta+order;
322 #ifndef AMREX_USE_GPU
323  int ksize = kend-ksta;
324  if (debug) amrex::Print() << " (2a) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
325 #endif
326  amrex::Real new_x;
327  amrex::GpuArray<amrex::Real,10> orig_x_sub;
328  amrex::GpuArray<amrex::Real,10> orig_y_sub;
329  amrex::Real* orig_x_sub_p = orig_x_sub.data();
330  amrex::Real* orig_y_sub_p = orig_y_sub.data();
331  if (exp_interp) {
332  new_x = new_x_p[new_k];
333  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
334  } else {
335  new_x = new_x_z[new_k];
336  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
337  }
338  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
339 #ifndef AMREX_USE_GPU
340  if (debug) {
341  amrex::Print() << " orig_x_sub_p = [";
342  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
343  amrex::Print() << "]" << std::endl;
344  amrex::Print() << " orig_y_sub_p = [";
345  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
346  amrex::Print() << "]" << std::endl;
347  }
348 #endif
349  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_l);
350  }
351  {
352  int ksta = kl-order/2;
353  int kend = ksta+order;
354 #ifndef AMREX_USE_GPU
355  int ksize = kend-ksta;
356  if (debug) amrex::Print() << "new_k=" << new_k << " (2b) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
357 #endif
358  amrex::Real new_x;
359  amrex::GpuArray<amrex::Real,10> orig_x_sub;
360  amrex::GpuArray<amrex::Real,10> orig_y_sub;
361  amrex::Real* orig_x_sub_p = orig_x_sub.data();
362  amrex::Real* orig_y_sub_p = orig_y_sub.data();
363  if (exp_interp) {
364  new_x = new_x_p[new_k];
365  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
366  } else {
367  new_x = new_x_z[new_k];
368  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
369  }
370  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
371 #ifndef AMREX_USE_GPU
372  if (debug) {
373  amrex::Print() << " orig_x_sub_p = [";
374  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
375  amrex::Print() << "]" << std::endl;
376  amrex::Print() << " orig_y_sub_p = [";
377  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
378  amrex::Print() << "]" << std::endl;
379  }
380 #endif
381  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_r);
382  }
383  new_y[new_k] = myhalf*(new_y_l+new_y_r);
384  } else if ((kl-(order/2-1) >= 0) && (kr+order/2 <= orig_npts-1)) {
385  int ksta = kl-(order/2-1);
386  int kend = ksta+order;
387 #ifndef AMREX_USE_GPU
388  int ksize = kend-ksta;
389  if (debug) amrex::Print() << " (3) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
390 #endif
391  amrex::Real new_x;
392  amrex::GpuArray<amrex::Real,10> orig_x_sub;
393  amrex::GpuArray<amrex::Real,10> orig_y_sub;
394  amrex::Real* orig_x_sub_p = orig_x_sub.data();
395  amrex::Real* orig_y_sub_p = orig_y_sub.data();
396  if (exp_interp) {
397  new_x = new_x_p[new_k];
398  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
399  } else {
400  new_x = new_x_z[new_k];
401  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
402  }
403  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
404 #ifndef AMREX_USE_GPU
405  if (debug) {
406  amrex::Print() << " orig_x_sub_p = [";
407  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
408  amrex::Print() << "]" << std::endl;
409  amrex::Print() << " orig_y_sub_p = [";
410  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
411  amrex::Print() << "]" << std::endl;
412  }
413 #endif
414  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
415  } else if ((kl-order/2 >= 0) && (kr+order/2-1 <= orig_npts-1)) {
416  int ksta = kl-order/2;
417  int kend = ksta+order;
418 #ifndef AMREX_USE_GPU
419  int ksize = kend-ksta;
420  if (debug) amrex::Print() << " (4) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
421 #endif
422  amrex::Real new_x;
423  amrex::GpuArray<amrex::Real,10> orig_x_sub;
424  amrex::GpuArray<amrex::Real,10> orig_y_sub;
425  amrex::Real* orig_x_sub_p = orig_x_sub.data();
426  amrex::Real* orig_y_sub_p = orig_y_sub.data();
427  if (exp_interp) {
428  new_x = new_x_p[new_k];
429  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
430  } else {
431  new_x = new_x_z[new_k];
432  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
433  }
434  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
435 #ifndef AMREX_USE_GPU
436  if (debug) {
437  amrex::Print() << " orig_x_sub_p = [";
438  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
439  amrex::Print() << "]" << std::endl;
440  amrex::Print() << " orig_y_sub_p = [";
441  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
442  amrex::Print() << "]" << std::endl;
443  }
444 #endif
445  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
446  } else {
447  // Linear interpolation.
448  int ksta = kl;
449  int kend = kr;
450 #ifndef AMREX_USE_GPU
451  int ksize = kend-ksta+1;
452  if (debug) amrex::Print() << " (5) order=" << order << " new_x=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
453 #endif
454  amrex::Real new_x;
455  amrex::GpuArray<amrex::Real,2> orig_x_sub;
456  amrex::GpuArray<amrex::Real,2> orig_y_sub;
457  amrex::Real* orig_x_sub_p = orig_x_sub.data();
458  amrex::Real* orig_y_sub_p = orig_y_sub.data();
459  if (exp_interp) {
460  new_x = new_x_p[new_k];
461  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
462  } else {
463  new_x = new_x_z[new_k];
464  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
465  }
466  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
467 #ifndef AMREX_USE_GPU
468  if (debug) {
469  amrex::Print() << " orig_x_sub_p = [";
470  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
471  amrex::Print() << "]" << std::endl;
472  amrex::Print() << " orig_y_sub_p = [";
473  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
474  amrex::Print() << "]" << std::endl;
475  }
476 #endif
477  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
478  }
479  }
480 #ifndef AMREX_USE_GPU
481  if (debug) amrex::Print() << " new_y[" << new_k << "]=" << new_y[new_k] << std::endl;
482 #endif
483  }
484 }
485 
486 /**
487  * Compute isothermal pressure at height z.
488  */
489 AMREX_FORCE_INLINE
490 AMREX_GPU_DEVICE
491 void
492 calc_p_isothermal (const amrex::Real& z,
493  amrex::Real& p)
494 {
495  p = p_0*std::exp(-CONST_GRAV*z/(amrex::Real(290.0)*R_d));
496 }
497 
498 /**
499  * Interpolate one metgrid column onto ERF vertical levels with quality control.
500  */
501 AMREX_FORCE_INLINE
502 AMREX_GPU_DEVICE
503 void
504 interpolate_column_metgrid (const bool& metgrid_use_below_sfc,
505  const bool& metgrid_use_sfc,
506  const bool& exp_interp,
507  const bool& metgrid_retain_sfc,
508  const amrex::Real& metgrid_proximity,
509  const int& metgrid_order,
510  const int& metgrid_force_sfc_k,
511  const int& i,
512  const int& j,
513  const int& kmax,
514  const int& src_comp,
515  const int& dst_comp,
516  char var_type,
517  char stag,
518  const amrex::Array4<amrex::Real const>& orig_z_full,
519  const amrex::Array4<amrex::Real const>& orig_data,
520  const amrex::Array4<amrex::Real const>& new_z_full,
521  const amrex::Array4<amrex::Real>& new_data_full)
522 {
523  // Here we closely follow WRF's vert_interp from
524  // dyn_em/module_initialize_real.F, although changes have been
525  // made to accommodate interpolation relative to height instead of
526  // pressure.
527  int imin_orig = amrex::lbound(amrex::Box(orig_data)).x;
528  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
529  int jmin_orig = amrex::lbound(amrex::Box(orig_data)).y;
530  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
531  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
532  int kmax_new = kmax + 1;
533 
534  AMREX_ASSERT(kmax_orig < 256);
535  AMREX_ASSERT(kmax_new < 256);
536 
537  amrex::GpuArray<amrex::Real,256> new_z;
538  amrex::GpuArray<amrex::Real,256> new_p;
539  amrex::GpuArray<amrex::Real,256> new_data;
540  amrex::Real* new_z_p = new_z.data();
541  amrex::Real* new_p_p = new_p.data();
542  amrex::Real* new_data_p = new_data.data();
543  for (int k=0; k < kmax_new; k++) {
544  if (stag == 'X') {
545  new_z_p[k] = fourth*(new_z_full(i,j,k)+new_z_full(i,j+1,k)+new_z_full(i,j,k+1)+new_z_full(i,j+1,k+1));
546  } else if (stag == 'Y') {
547  new_z_p[k] = fourth*(new_z_full(i,j,k)+new_z_full(i+1,j,k)+new_z_full(i,j,k+1)+new_z_full(i+1,j,k+1));
548  } else if (stag == 'M') {
549  new_z_p[k] = amrex::Real(0.125)*(new_z_full(i,j,k )+new_z_full(i,j+1,k )+new_z_full(i+1,j,k )+new_z_full(i+1,j+1,k )+
550  new_z_full(i,j,k+1)+new_z_full(i,j+1,k+1)+new_z_full(i+1,j,k+1)+new_z_full(i+1,j+1,k+1));
551  }
552  calc_p_isothermal(new_z_p[k], new_p_p[k]);
553  }
554 
555  amrex::GpuArray<amrex::Real,256> orig_z;
556  amrex::Real* orig_z_p = orig_z.data();
557  for (int k=0; k < kmax_orig; k++) {
558  if (stag == 'M') {
559  orig_z_p[k] = orig_z_full(i,j,k);
560  } else if (stag == 'X') {
561  if (i <= imin_orig) {
562  orig_z_p[k] = orig_z_full(i,j,k);
563  } else if (i >= imax_orig) {
564  orig_z_p[k] = orig_z_full(imax_orig-1,j,k);
565  } else {
566  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k));
567  }
568  } else if (stag == 'Y') {
569  if (j <= jmin_orig) {
570  orig_z_p[k] = orig_z_full(i,j,k);
571  } else if (j >= jmax_orig) {
572  orig_z_p[k] = orig_z_full(i,jmax_orig-1,k);
573  } else {
574  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k));
575  }
576  }
577  }
578 
579  // Check if the data is top-down instead of bottom-up.
580  bool flip_data_required = false;
581  if (orig_z[1] > orig_z[kmax_orig-1]) flip_data_required = true;
582  if (flip_data_required) amrex::Abort("metgrid initialization flip_data_required. Not yet implemented.");
583 
584  // Search for the first level above the surface in the origin data.
585  // This is needed since the origin model topography will be
586  // different than the topography processed by WPS.
587  int k_above_sfc = 0;
588  for (int k=1; k < kmax_orig; k++) {
589  if (orig_z_p[k] > orig_z_p[0]) {
590  k_above_sfc = k;
591  break;
592  }
593  }
594 
595  int kend_order;
596  amrex::GpuArray<amrex::Real,256> ordered_z;
597  amrex::GpuArray<amrex::Real,256> ordered_data;
598  amrex::Real* ordered_z_p = ordered_z.data();
599  amrex::Real* ordered_data_p = ordered_data.data();
600  if (k_above_sfc > 1) {
601  // The levels are not monotonically increasing in height, so
602  // we sort and then make "artistic" quality control choices.
603  int count = 0;
604 
605  // Insert levels that are below the surface.
606  for (int k=1; k < k_above_sfc; k++) {
607  ordered_z_p[count] = orig_z_p[k];
608  ordered_data_p[count] = orig_data(i,j,k,src_comp);
609  count++;
610  }
611 
612  // Check if the level that is nearest to and below the surface
613  // is "too close". If so, we'll ignore the upper level and keep
614  // the lower. Origin data is likely to be on pressure levels
615  // with higher spatial resolution near-surface, which supports
616  // the choice of eliminating levels that are "too close" in
617  // pressure-space. For simplicity, calculate delta P assuming a
618  // baroclinic atmosphere.
619  amrex::Real Pu, Pl;
620  calc_p_isothermal(orig_z_p[0], Pu);
621  calc_p_isothermal(ordered_z_p[count-1], Pl);
622  if (Pl-Pu < metgrid_proximity) {
623  count--;
624  }
625 
626  // Insert the surface level.
627  ordered_z_p[count] = orig_z_p[0];
628  ordered_data_p[count] = orig_data(i,j,0,src_comp);
629  count++;
630 
631  // Quoting WRF's comments, the next level to use is at,
632  // "... ta da, the first level above the surface. I know, wow."
633  int knext = k_above_sfc;
634  // Conditionally more strongly use the surface data by removing
635  // levels between the surface and the height corresponding to a
636  // set number of ERF levels from the surface. This forces the
637  // interpolator to use the surface data up through a number of
638  // ERF levels from the surface.
639  if (metgrid_force_sfc_k > 0) {
640  for (int k=k_above_sfc; k < kmax_orig; k++) {
641  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
642  knext = k;
643  break;
644  }
645  }
646  }
647 
648  // Check if the level that is nearest to and above the surface
649  // is "too close". If so, we'll ignore that level.
650  calc_p_isothermal(orig_z_p[knext], Pu);
651  calc_p_isothermal(ordered_z_p[count-1], Pl);
652  if (Pl-Pu < metgrid_proximity) {
653  knext++;
654  }
655 
656  // Insert levels that are above the surface.
657  for (int k=knext; k < kmax_orig; k++) {
658  ordered_z_p[count] = orig_z_p[k];
659  ordered_data_p[count] = orig_data(i,j,k,src_comp);
660  count++;
661  }
662 
663  kend_order = count;
664  } else {
665  // The surface is the lowest level in the origin data.
666 
667  // Insert the surface.
668  ordered_z_p[0] = orig_z[0];
669  ordered_data_p[0] = orig_data(i,j,0,src_comp);
670 
671  // Similar to above, conditionally more strongly use the
672  // surface data.
673  int count = 1;
674  int knext = count;
675  if (metgrid_force_sfc_k > 0) {
676  for (int k=knext; k < kmax_orig; k++) {
677  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k]) {
678  knext = k;
679  break;
680  }
681  }
682  }
683 
684  // Insert the remaining levels, again ignoring levels that are
685  // "too close" to the prior valid level.
686  for (int k=knext; k < kmax_orig; k++) {
687  amrex::Real Pu, Pl;
688  calc_p_isothermal(orig_z_p[k], Pu);
689  calc_p_isothermal(ordered_z_p[count-1], Pl);
690  if (Pl-Pu < metgrid_proximity) {
691  continue;
692  }
693  ordered_z_p[count] = orig_z_p[k];
694  ordered_data_p[count] = orig_data(i,j,k,src_comp);
695  count++;
696  }
697  kend_order = count;
698  }
699 
700  int ksta(0), kend(0);
701  if (metgrid_use_below_sfc && metgrid_use_sfc) {
702  // Use all levels.
703  ksta = 0;
704  kend = kend_order-1;
705  } else if (metgrid_use_below_sfc && !metgrid_use_sfc) {
706  // Use all levels except for the surface.
707  int ksfc = 0;
708  for (int k=0; k < kmax_orig; k++) {
709  if (ordered_z_p[k] == orig_z_p[0]) {
710  ksfc = k;
711  break;
712  }
713  }
714  for (int k=ksfc; k < kmax_orig-1; k++) {
715  ordered_z_p[k] = ordered_z_p[k+1];
716  ordered_data_p[k] = ordered_data_p[k+1];
717  }
718  ksta = 0;
719  kend = kend_order-2;
720  } else if (!metgrid_use_below_sfc && metgrid_use_sfc) {
721  // Use all levels above and including the surface.
722  int ksfc = 0;
723  for (int k=0; k < kend_order; k++) {
724  if (ordered_z_p[k] == orig_z_p[0]) {
725  ksfc = k;
726  break;
727  }
728  }
729  ksta = ksfc;
730  kend = kend_order-1;
731  } else {
732  // We shouldn't be in here!
733  amrex::Abort("metgrid initialization, !use_levels below_ground && !metgrid_use_sfc");
734  }
735 
736  // Insert the level of maximum winds.
737 // amrex::Real maxw_above_this_level = amrex::Real(30000.0);
738 // amrex::Real maxw_horiz_pres_diff = amrex::Real(5000.0);
739 // if ((flag_maxw == 1) && (use_maxw)) {
740 // amrex::Abort("metgrid initialization, use_maxw not yet implemented");
741 // }
742 
743  // Insert the level of the tropopause.
744 // amrex::Real trop_horiz_pres_diff = amrex::Real(5000.0);
745 // if ((flag_trop == 1) && (use_trop)) {
746 // amrex::Abort("metgrid initialization, use_trop not yet implemented");
747 // }
748 
749  amrex::GpuArray<amrex::Real,256> ordered_p;
750  amrex::Real* ordered_p_p = ordered_p.data();
751  for (int k=0; k < kend_order; k++) {
752  calc_p_isothermal(ordered_z_p[k], ordered_p_p[k]);
753  }
754 
755  // kfinal is the index of the topmost level retained in the final_*
756  // arrays, so those arrays hold kfinal+1 points.
757  int kfinal = 0;
758  amrex::GpuArray<amrex::Real,256> final_z;
759  amrex::GpuArray<amrex::Real,256> final_p;
760  amrex::GpuArray<amrex::Real,256> final_data;
761  amrex::Real* final_z_p = final_z.data();
762  amrex::Real* final_p_p = final_p.data();
763  amrex::Real* final_data_p = final_data.data();
764  final_z_p[0] = ordered_z[ksta];
765  final_p_p[0] = ordered_p[ksta];
766  final_data_p[0] = ordered_data[ksta];
767  for (int k=ksta+1; k <= kend; k++) {
768  // Skip levels that are "too close" to the prior retained level.
769  if ((final_p_p[kfinal]-ordered_p_p[k]) < metgrid_proximity) continue;
770  kfinal++;
771  final_z_p[kfinal] = ordered_z_p[k];
772  final_p_p[kfinal] = ordered_p_p[k];
773  final_data_p[kfinal] = ordered_data_p[k];
774  }
775 
776  // Call the interpolator.
777  lagrange_setup(var_type,
778  exp_interp,
779  kfinal+1,
780  kmax_new,
781  metgrid_order,
782  i,
783  j,
784  final_z_p,
785  final_p_p,
786  final_data_p,
787  new_z_p,
788  new_p_p,
789  new_data_p);
790 
791  // Optionally replace the lowest level of data with the surface
792  // field from the origin data.
793  if (metgrid_retain_sfc) new_data[0] = ordered_data[0];
794 
795  // Save the interpolated data.
796  for (int k=0; k < kmax_new; k++) {
797  new_data_full(i,j,k,dst_comp) = new_data[k];
798  }
799 
800 }
801 
802 /**
803  * Linearly interpolate one metgrid column value onto an ERF point.
804  */
805 AMREX_FORCE_INLINE
806 AMREX_GPU_DEVICE
807 amrex::Real
808 interpolate_column_metgrid_linear (const int& i,
809  const int& j,
810  const int& k,
811  char stag,
812  int src_comp,
813  const amrex::Array4<amrex::Real const>& orig_z,
814  const amrex::Array4<amrex::Real const>& orig_data,
815  const amrex::Array4<amrex::Real const>& new_z)
816 {
817  // This subroutine is a bit ham-handed and can be cleaned up later.
818  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
819  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
820  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
821 
822  amrex::Real z;
823  if (stag == 'X') {
824  z = fourth*(new_z(i,j,k)+new_z(i,j+1,k)+new_z(i,j,k+1)+new_z(i,j+1,k+1));
825  }
826  else if (stag == 'Y') {
827  z = fourth*(new_z(i,j,k)+new_z(i+1,j,k)+new_z(i,j,k+1)+new_z(i+1,j,k+1));
828  }
829  else if (stag == 'M') {
830  z = amrex::Real(0.125)*(new_z(i,j,k )+new_z(i,j+1,k )+new_z(i+1,j,k )+new_z(i+1,j+1,k )+
831  new_z(i,j,k+1)+new_z(i,j+1,k+1)+new_z(i+1,j,k+1)+new_z(i+1,j+1,k+1));
832  }
833 
834  // Initialized to NaN so that any path which fails to set them fails loudly
835  // (with amrex.fpe_trap_invalid=1) rather than returning plausible garbage.
836  amrex::Real z0 = std::numeric_limits<amrex::Real>::quiet_NaN();
837  amrex::Real z1 = std::numeric_limits<amrex::Real>::quiet_NaN();
838  int klow = -1;
839  int khi0 = -1;
840  amrex::Real dzlow = amrex::Real(1.0e12);
841  amrex::Real dzhi0 = -amrex::Real(1.0e12);
842  for (int kk = 0; kk < kmax_orig; kk++) {
843  amrex::Real orig_z_stag = zero;
844  if (stag == 'M') {
845  orig_z_stag = orig_z(i,j,kk);
846  }
847  if (stag == 'X') {
848  if (i == 0) {
849  orig_z_stag = orig_z(i,j,kk);
850  }
851  else if (i == imax_orig) {
852  orig_z_stag = orig_z(imax_orig-1,j,kk);
853  }
854  else {
855  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
856  }
857  }
858  else if (stag == 'Y') {
859  if (j == 0) {
860  orig_z_stag = orig_z(i,j,kk);
861  }
862  else if (j == jmax_orig) {
863  orig_z_stag = orig_z(i,jmax_orig-1,kk);
864  }
865  else {
866  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
867  }
868  }
869 
870  amrex::Real dz = z - orig_z_stag;
871  if ((dz < zero) && (dz > dzhi0)) {
872  dzhi0 = dz;
873  khi0 = kk;
874  z1 = orig_z_stag;
875  }
876  if ((dz >= zero) && (dz < dzlow)) {
877  dzlow = dz;
878  klow = kk;
879  z0 = orig_z_stag;
880  }
881  } // kk
882 
883  // extrapolate below the bottom surface
884  if (klow == -1) {
885  // klow was never found, so z0 was never set above. The two levels used
886  // here are khi0 and khi1, both above z, so z0 is the height at khi0 --
887  // which is exactly what the search loop left in z1 before we clobber it.
888  z0 = z1;
889  int khi1 = -1;
890  amrex::Real dzhi1 = -amrex::Real(1.0e12);
891  for (int kk = 0; kk < kmax_orig; kk++) {
892  amrex::Real orig_z_stag = zero;
893  if (stag == 'M') {
894  orig_z_stag = orig_z(i,j,kk);
895  }
896  else if (stag == 'X') {
897  if (i == 0) {
898  orig_z_stag = orig_z(i,j,kk);
899  }
900  else if (i == imax_orig) {
901  orig_z_stag = orig_z(imax_orig-1,j,kk);
902  }
903  else {
904  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
905  }
906  }
907  else if (stag == 'Y') {
908  if (j == 0) {
909  orig_z_stag = orig_z(i,j,kk);
910  }
911  else if (j == jmax_orig) {
912  orig_z_stag = orig_z(i,jmax_orig-1,kk);
913  }
914  else {
915  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
916  }
917  }
918  amrex::Real dz = z - orig_z_stag;
919  if ((dz < zero) && (dz > dzhi1) && (kk != khi0)) {
920  dzhi1 = dz;
921  khi1 = kk;
922  z1 = orig_z_stag;
923  }
924  }
925  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
926  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
927  return ( y0-(y1-y0)/(z1-z0)*(z0-z) );
928 
929  // Extrapolate above the top surface
930  } else if (khi0 == -1) {
931  khi0 = klow - 1;
932  int khi1 = klow;
933  // khi0 was never found, so z1 was never set above. The upper level used
934  // here is khi1 = klow, whose height is currently held in z0 -- capture
935  // it before z0 is recomputed at the new (lower) khi0.
936  z1 = z0;
937  if (stag == 'M') {
938  z0 = orig_z(i,j,khi0);
939  }
940  else if (stag == 'X') {
941  if (i == 0) {
942  z0 = orig_z(i,j,khi0);
943  }
944  else if (i == imax_orig) {
945  z0 = orig_z(imax_orig-1,j,khi0);
946  }
947  else {
948  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
949  }
950  }
951  else if (stag == 'Y') {
952  if (j == 0) {
953  z0 = orig_z(i,j,khi0);
954  }
955  else if (j == jmax_orig) {
956  z0 = orig_z(i,jmax_orig-1,khi0);
957  }
958  else {
959  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
960  }
961  }
962  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
963  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
964  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
965  } else {
966  // interpolate
967  amrex::Real y0 = orig_data(i,j,klow,src_comp);
968  amrex::Real y1 = orig_data(i,j,khi0,src_comp);
969  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
970 
971  }
972 }
973 
974 /**
975  * Convert relative humidity to vapor mixing ratio using WRF-compatible bounds.
976  */
977 AMREX_FORCE_INLINE
978 AMREX_GPU_DEVICE
979 void
980 rh_to_mxrat (int i,
981  int j,
982  int k,
983  const amrex::Array4<amrex::Real const>& rhum,
984  const amrex::Array4<amrex::Real const>& temp,
985  const amrex::Array4<amrex::Real const>& pres,
986  int src_indx,
987  const amrex::Array4<amrex::Real>& mxrat)
988 {
989  amrex::Real qv_max_p_safe = amrex::Real(10000.0); // WRF default value
990  amrex::Real qv_max_flag = amrex::Real(1.0e-5); // WRF default value
991  amrex::Real qv_max_value = amrex::Real(3.0e-6); // WRF default value
992  amrex::Real qv_min_p_safe = amrex::Real(110000.0); // WRF default value
993  amrex::Real qv_min_flag = amrex::Real(1.0e-6); // WRF default value
994  amrex::Real qv_min_value = amrex::Real(1.0e-6); // WRF default value
995  amrex::Real eps = RdoRv;
996  amrex::Real svp1 = amrex::Real(0.6112);
997  amrex::Real svp2 = amrex::Real(17.67);
998  amrex::Real svp3 = amrex::Real(29.65);
999  amrex::Real svpt0 = amrex::Real(273.15);
1000  // WRF's method when model_config_rec%rh2qv_wrt_liquid=.true. (default behavior)
1001  if (temp(i,j,k) != zero) {
1002  amrex::Real es=amrex::Real(0.01)*rhum(i,j,k)*svp1*amrex::Real(10.0)*std::exp(svp2*(temp(i,j,k)-svpt0)/(temp(i,j,k)-svp3));
1003  if (es >= pres(i,j,k)/amrex::Real(100.0)) {
1004  // vapor pressure exceeds total pressure
1005  mxrat(i,j,k,src_indx) = amrex::Math::powi<-6>(amrex::Real(10.0));
1006  }
1007  else {
1008  mxrat(i,j,k,src_indx) = amrex::max(eps*es/(pres(i,j,k)/amrex::Real(100.0)-es), amrex::Real(1.0e-6));
1009  }
1010  }
1011  else {
1012  // I don't know why there's a fringe case handled in WRF where T is absolute zero...
1013  // Let's just deal with it here in case we also end up needing it.
1014  mxrat(i,j,k,src_indx) = amrex::Real(1.0e-6);
1015  }
1016  // See the below comment from WRF dyn_em/module_initialize_real.F rh_to_mxrat1.
1017  // For pressures above a defined level, reasonable Qv values should be
1018  // a certain value or smaller. If they are larger than this, the input data
1019  // probably had "missing" RH, and we filled in some values. This is an
1020  // attempt to catch those. Also, set the minimum value for the entire
1021  // domain that is above the selected pressure level.
1022  if (pres(i,j,k) < qv_max_p_safe) {
1023  if (mxrat(i,j,k,src_indx) > qv_max_flag) {
1024  mxrat(i,j,k,src_indx) = qv_max_value;
1025  }
1026  }
1027  if (pres(i,j,k) < qv_min_p_safe) {
1028  if (mxrat(i,j,k,src_indx) < qv_min_flag) {
1029  mxrat(i,j,k,src_indx) = qv_min_value;
1030  }
1031  }
1032 }
1033 #endif
ERF_BaseStateParams.H
one
constexpr amrex::Real one
Definition: ERF_Constants.H:9
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
p_0
constexpr amrex::Real p_0
Definition: ERF_Constants.H:61
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:64
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:47
RdoCp
constexpr amrex::Real RdoCp
Definition: ERF_Constants.H:54
RdoRv
constexpr amrex::Real RdoRv
Definition: ERF_Constants.H:57
Coord::y
@ y
Coord::x
@ x
Coord::z
@ z
ERF_EOS.H
ERF_HSEUtils.H
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
use_moisture
const bool use_moisture
Definition: ERF_InitCustomPert_ABL.H:71
init_terrain_from_metgrid
void init_terrain_from_metgrid(const bool &avg_grid_faces_to_nodes, amrex::Geometry &geom, amrex::MultiFab *z_phys_nd, amrex::FArrayBox &NC_hgt_fab)
init_base_state_from_metgrid
void init_base_state_from_metgrid(const bool use_moisture, const bool metgrid_debug_psfc, const amrex::Real l_rdOcp, const amrex::Box &domain, const amrex::Box &valid_bx, const int &flag_psfc, amrex::FArrayBox &state_fab, amrex::FArrayBox &r_hse_fab, amrex::FArrayBox &p_hse_fab, amrex::FArrayBox &pi_hse_fab, amrex::FArrayBox &th_hse_fab, amrex::FArrayBox &qv_hse_fab, amrex::FArrayBox &z_phys_nd_fab, amrex::FArrayBox &z_phys_cc_fab, const amrex::FArrayBox &NC_psfc_fab, const BaseStateParams &bsp)
lagrange_interp
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_interp(const int &order, amrex::Real *x, amrex::Real *y, amrex::Real &new_x, amrex::Real &new_y)
Definition: ERF_MetgridUtils.H:148
interpolate_column_metgrid_linear
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real interpolate_column_metgrid_linear(const int &i, const int &j, const int &k, char stag, int src_comp, const amrex::Array4< amrex::Real const > &orig_z, const amrex::Array4< amrex::Real const > &orig_data, const amrex::Array4< amrex::Real const > &new_z)
Definition: ERF_MetgridUtils.H:808
interpolate_column_metgrid
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void interpolate_column_metgrid(const bool &metgrid_use_below_sfc, const bool &metgrid_use_sfc, const bool &exp_interp, const bool &metgrid_retain_sfc, const amrex::Real &metgrid_proximity, const int &metgrid_order, const int &metgrid_force_sfc_k, const int &i, const int &j, const int &kmax, const int &src_comp, const int &dst_comp, char var_type, char stag, const amrex::Array4< amrex::Real const > &orig_z_full, const amrex::Array4< amrex::Real const > &orig_data, const amrex::Array4< amrex::Real const > &new_z_full, const amrex::Array4< amrex::Real > &new_data_full)
Definition: ERF_MetgridUtils.H:504
read_from_metgrid
void read_from_metgrid(int lev, int itime, const amrex::Box &domain, const std::string &fname, std::string &NC_dateTime, double &NC_epochTime, int &flag_psfc, int &flag_msf, int &flag_sst, int &flag_tsk, int &flag_lmask, int &NC_nx, int &NC_ny, amrex::Real &NC_dx, amrex::Real &NC_dy, amrex::FArrayBox &NC_xvel_fab, amrex::FArrayBox &NC_yvel_fab, amrex::FArrayBox &NC_temp_fab, amrex::FArrayBox &NC_rhum_fab, amrex::FArrayBox &NC_pres_fab, amrex::FArrayBox &NC_ght_fab, amrex::FArrayBox &NC_hgt_fab, amrex::FArrayBox &NC_psfc_fab, amrex::FArrayBox &NC_msfu_fab, amrex::FArrayBox &NC_msfv_fab, amrex::FArrayBox &NC_msfm_fab, amrex::FArrayBox &NC_sst_fab, amrex::FArrayBox &NC_tsk_fab, amrex::FArrayBox &NC_LAT_fab, amrex::FArrayBox &NC_LON_fab, amrex::IArrayBox &NC_lmask_iab, amrex::Geometry &geom)
calc_p_isothermal
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal(const amrex::Real &z, amrex::Real &p)
Definition: ERF_MetgridUtils.H:492
init_msfs_from_metgrid
void init_msfs_from_metgrid(const bool metgrid_debug_msf, amrex::FArrayBox &msfu_fab, amrex::FArrayBox &msfv_fab, amrex::FArrayBox &msfm_fab, const int &flag_msf, amrex::FArrayBox &NC_MSFU_fab, amrex::FArrayBox &NC_MSFV_fab, amrex::FArrayBox &NC_MSFM_fab)
lagrange_setup
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_setup(char var_type, const bool &exp_interp, const int &orig_npts, const int &new_npts, const int &order, const int &i, const int &j, amrex::Real *orig_x_z, amrex::Real *orig_x_p, amrex::Real *orig_y, amrex::Real *new_x_z, amrex::Real *new_x_p, amrex::Real *new_y)
Definition: ERF_MetgridUtils.H:185
init_state_from_metgrid
void init_state_from_metgrid(const int lev, const int itime, const bool use_moisture, const bool interp_theta, const bool metgrid_debug_quiescent, const bool metgrid_debug_isothermal, const bool metgrid_debug_dry, const bool metgrid_basic_linear, const bool metgrid_use_below_sfc, const bool metgrid_use_sfc, const bool metgrid_retain_sfc, const amrex::Real metgrid_proximity, const int metgrid_order, const int metgrid_metgrid_force_sfc_k, const amrex::Real l_rdOcp, amrex::Box &tbxc, amrex::Box &tbxu, amrex::Box &tbxv, amrex::Box &tbxw, amrex::FArrayBox &state_fab, amrex::FArrayBox &x_vel_fab, amrex::FArrayBox &y_vel_fab, amrex::FArrayBox &z_vel_fab, amrex::FArrayBox &z_phys_nd_fab, const amrex::FArrayBox &NC_ght_fab, const amrex::FArrayBox &NC_xvel_fab, const amrex::FArrayBox &NC_yvel_fab, const amrex::FArrayBox &NC_temp_fab, const amrex::FArrayBox &NC_rhum_fab, const amrex::FArrayBox &NC_pres_fab, amrex::FArrayBox &tmp_src_fab, amrex::FArrayBox &tmp_dst_fab, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_xlo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_xhi, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_ylo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_yhi, const amrex::Array4< const int > &mask_c_arr, const amrex::Array4< const int > &mask_u_arr, const amrex::Array4< const int > &mask_v_arr)
rh_to_mxrat
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void rh_to_mxrat(int i, int j, int k, const amrex::Array4< amrex::Real const > &rhum, const amrex::Array4< amrex::Real const > &temp, const amrex::Array4< amrex::Real const > &pres, int src_indx, const amrex::Array4< amrex::Real > &mxrat)
Definition: ERF_MetgridUtils.H:980
ERF_ProbCommon.H
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
ERF_Utils.H
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:27
WDM6SedCellScratch::dz
@ dz
Definition: ERF_AdvanceWDM6.cpp:270
WSM6Ind::p
@ p
Definition: ERF_WSM6.H:191
module_model_constants::svp1
real(c_double), parameter svp1
Definition: ERF_module_model_constants.F90:78
module_model_constants::svp3
real(c_double), parameter svp3
Definition: ERF_module_model_constants.F90:80
module_model_constants::svp2
real(c_double), parameter svp2
Definition: ERF_module_model_constants.F90:79
module_model_constants::svpt0
real(c_double), parameter svpt0
Definition: ERF_module_model_constants.F90:81
BaseStateParams
Definition: ERF_BaseStateParams.H:35