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