ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_InputSoundingData.H
Go to the documentation of this file.
1 #ifndef ERF_INPUT_SOUNDING_DATA_H_
2 #define ERF_INPUT_SOUNDING_DATA_H_
3 
4 #include <limits>
5 #include <string>
6 #include <iostream>
7 
8 #include <AMReX_ParmParse.H>
9 #include <AMReX_Print.H>
10 #include <AMReX_Gpu.H>
11 #include <AMReX_Geometry.H>
12 
13 #include <ERF_EOS.H>
14 #include <ERF_Constants.H>
15 #include <ERF_Interpolation_1D.H>
16 #include <ERF_HSEUtils.H>
17 
18 /**
19  * Data structure storing input sounding data. Also
20  * handles reading the input file for sounding data and
21  * hydrostatic column integration.
22  */
24 public:
25  /**
26  * @brief Construct input sounding metadata from the ERF input namespace.
27  */
29  {
30  amrex::ParmParse pp("erf");
31  pp.queryAdd("tau_nudging", tau_nudging);
32 
33  // Read in input_sounding filename
34  n_sounding_files = pp.countval("input_sounding_file");
35  if (n_sounding_files > 0) {
37  pp.queryarr("input_sounding_file", input_sounding_file, 0, n_sounding_files);
38  } else {
39  n_sounding_files = 1;
41  input_sounding_file[0] = "input_sounding";
42  }
43 
44  // Read in input_sounding times
45  n_sounding_times = pp.countval("input_sounding_time");
46 
47  if (n_sounding_times > 0) {
49  pp.queryarr("input_sounding_time", input_sounding_time, 0, n_sounding_times);
50  } else {
51  n_sounding_times = 1;
54  }
55 
56  // If we have more files than times or times than files we just use the minimum
57  int n = std::min(n_sounding_times, n_sounding_files);
58  n_sounding_files = n;
59  n_sounding_times = n;
60  input_sounding_file.resize(n);
61  input_sounding_time.resize(n);
62  }
63 
64  /**
65  * @brief Resize host and device sounding arrays for the configured sounding times.
66  */
67  void resize_arrays ()
68  {
70 
71  z_inp_sound.resize(ntimes);
72  theta_inp_sound.resize(ntimes);
73  qv_inp_sound.resize(ntimes);
74  U_inp_sound.resize(ntimes);
75  V_inp_sound.resize(ntimes);
76 
77  z_inp_sound_d.resize(ntimes);
78  theta_inp_sound_d.resize(ntimes);
79  qv_inp_sound_d.resize(ntimes);
80  U_inp_sound_d.resize(ntimes);
81  V_inp_sound_d.resize(ntimes);
82  }
83 
84  /**
85  * @brief Read and interpolate one input sounding profile onto level-0 heights.
86  * @param geom Geometry defining the level-0 vertical domain.
87  * @param zlevels_stag Staggered vertical coordinates used for interpolation.
88  * @param itime Sounding time/file index to read.
89  * @param is_moist Whether the selected physics permits nonzero water vapor.
90  */
91  void read_from_file (const amrex::Geometry &geom,
92  const amrex::Vector<amrex::Real>& zlevels_stag,
93  int itime, bool is_moist)
94  {
95  const int klo = 0;
96  const int khi = geom.Domain().bigEnd()[AMREX_SPACEDIM-1];
97  const int Nz = geom.Domain().size()[AMREX_SPACEDIM-1];
98 
99  const amrex::Real zbot = zlevels_stag[klo];
100  const amrex::Real ztop = zlevels_stag[khi+1];
101 
102  z_inp_sound[itime].resize(Nz+2);
103  theta_inp_sound[itime].resize(Nz+2);
104  qv_inp_sound[itime].resize(Nz+2);
105  U_inp_sound[itime].resize(Nz+2);
106  V_inp_sound[itime].resize(Nz+2);
107 
108  // Read the input_sounding file
109  amrex::Print() << "input_sounding file location : " << input_sounding_file[itime] << std::endl;
110  std::ifstream input_sounding_reader(input_sounding_file[itime]);
111  if(!input_sounding_reader.is_open()) {
112  amrex::Error("Error opening the input_sounding file\n");
113  }
114  else {
115  // Read the contents of the input_sounding file
116  amrex::Print() << "Successfully opened the input_sounding file. Now reading... " << std::endl;
117  std::string line;
118 
119  // First, read the input data into temp vectors; then, interpolate vectors to the
120  // domain lo/hi and cell centers (from level 0)
121  amrex::Vector<amrex::Real> z_inp_sound_tmp, theta_inp_sound_tmp, qv_inp_sound_tmp,
122  U_inp_sound_tmp, V_inp_sound_tmp;
123 
124  // Read surface quantities from the first line
125  std::getline(input_sounding_reader, line);
126  std::istringstream iss(line);
128  press_ref_inp_sound *= 100; // convert from hPa to Pa
129  qv_ref_inp_sound *= amrex::Real(0.001); // convert from g/kg to kg/kg
130 
131  // Add surface
132  z_inp_sound_tmp.push_back(zbot); // height above sea level [m]
133  theta_inp_sound_tmp.push_back(theta_ref_inp_sound);
134  qv_inp_sound_tmp.push_back(qv_ref_inp_sound);
135  U_inp_sound_tmp.push_back(0);
136  V_inp_sound_tmp.push_back(0);
137 
138  // *************************************************************************************************
139  // Make sure that if we specify surface values for temp in the inputs file and in the input sounding file,
140  // those values must match
141  // *************************************************************************************************
142  amrex::ParmParse pp_erf("erf");
143 
144  // Read with query, not queryAdd: ERF_SurfaceLayer.H reads the same keys and
145  // treats a key that is present as given, so adding an absent key here would
146  // make it look given there (and be rejected as non-positive). ERF_SurfaceLayer.H
147  // also rejects a given value that is not positive, so the sentinel test below
148  // is the same as testing that the key was given.
149  // Only a surface potential temperature (most.pot_surf_temp, or the deprecated
150  // most.surf_temp) is comparable to the sounding's surface theta;
151  // most.abs_surf_temp is not checked.
152  amrex::Real surf_temp = amrex::Real(-1.0);
153  amrex::Real old_surf_temp = amrex::Real(-1.0);
154  pp_erf.query("most.pot_surf_temp", surf_temp);
155  pp_erf.query("most.surf_temp", old_surf_temp);
156  if (old_surf_temp > amrex::Real(0.0)) { surf_temp = old_surf_temp; }
157  bool surf_temp_specified = (surf_temp > amrex::Real(0.0));
158  if ( surf_temp_specified && (surf_temp != theta_inp_sound_tmp[0]) ) {
159  amrex::Print() << "Input sounding temp does not match most.pot_surf_temp" << std::endl;
160  amrex::Print() << "Modify one or both to match." << std::endl;
161  amrex::Abort();
162  }
163 
164  // *************************************************************************************************
165  // Make sure that if we specify surface values for qv in the inputs file and in the input sounding file,
166  // those values must match
167  // *************************************************************************************************
168  // Negative sentinel rather than the queryAdd return value: a mixing ratio is
169  // never negative, and ERF_SurfaceLayer.H parses this same key.
170  amrex::Real surf_moist = amrex::Real(-1.0);
171  pp_erf.queryAdd("most.surf_moist", surf_moist);
172  bool surf_moist_specified = (surf_moist >= amrex::Real(0.0));
173  if ( surf_moist_specified && (surf_moist != qv_inp_sound_tmp[0]) ) {
174  amrex::Print() << "Input sounding qv does not match most.surf_moist" << std::endl;
175  amrex::Print() << "Modify one or both to match." << std::endl;
176  amrex::Abort();
177  }
178 
179  // Read the vertical profile at each given height
180  amrex::Real z, theta, qv, U, V;
181 
182  while(std::getline(input_sounding_reader, line)) {
183  std::istringstream iss_z(line);
184  iss_z >> z >> theta >> qv >> U >> V;
185 
186  // Dont read in non-zero qv if not using a moisture model
187  AMREX_ALWAYS_ASSERT(qv == zero || is_moist);
188 
189  if (z == zbot) {
190  AMREX_ALWAYS_ASSERT(theta == theta_inp_sound_tmp[0]);
191  AMREX_ALWAYS_ASSERT(qv*amrex::Real(0.001) == qv_inp_sound_tmp[0]); // convert from g/kg to kg/kg
192  U_inp_sound_tmp[0] = U;
193  V_inp_sound_tmp[0] = V;
194  } else {
195  AMREX_ALWAYS_ASSERT(z > z_inp_sound_tmp[z_inp_sound_tmp.size()-1]); // sounding is increasing in height
196  z_inp_sound_tmp.push_back(z);
197 
198  theta_inp_sound_tmp.push_back(theta);
199  qv_inp_sound_tmp.push_back(qv*amrex::Real(0.001)); // convert from g/kg to kg/kg
200  U_inp_sound_tmp.push_back(U);
201  V_inp_sound_tmp.push_back(V);
202  if (z >= ztop) break;
203  }
204  }
205 
206  // At this point, we have an input_sounding from zbot up to
207  // z_inp_sound_tmp[N-1] >= ztop. Now, interpolate to grid level 0 heights
208  const int Ninp = z_inp_sound_tmp.size();
209  z_inp_sound[itime][0] = zbot;
210  theta_inp_sound[itime][0] = theta_inp_sound_tmp[0];
211  qv_inp_sound[itime][0] = qv_inp_sound_tmp[0];
212  U_inp_sound[itime][0] = U_inp_sound_tmp[0];
213  V_inp_sound[itime][0] = V_inp_sound_tmp[0];
214  for (int k=0; k < Nz; ++k) {
215  z_inp_sound[itime][k+1] = myhalf * (zlevels_stag[k] + zlevels_stag[k+1]);
216  theta_inp_sound[itime][k+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), theta_inp_sound_tmp.dataPtr(), z_inp_sound[itime][k+1], Ninp);
217  qv_inp_sound[itime][k+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), qv_inp_sound_tmp.dataPtr(), z_inp_sound[itime][k+1], Ninp);
218  U_inp_sound[itime][k+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), U_inp_sound_tmp.dataPtr(), z_inp_sound[itime][k+1], Ninp);
219  V_inp_sound[itime][k+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), V_inp_sound_tmp.dataPtr(), z_inp_sound[itime][k+1], Ninp);
220  }
221  z_inp_sound[itime][Nz+1] = ztop;
222  theta_inp_sound[itime][Nz+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), theta_inp_sound_tmp.dataPtr(), ztop, Ninp);
223  qv_inp_sound[itime][Nz+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), qv_inp_sound_tmp.dataPtr(), ztop, Ninp);
224  U_inp_sound[itime][Nz+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), U_inp_sound_tmp.dataPtr(), ztop, Ninp);
225  V_inp_sound[itime][Nz+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), V_inp_sound_tmp.dataPtr(), ztop, Ninp);
226  }
227 
228  amrex::Print() << "Successfully read the " << itime << "th input_sounding file..." << std::endl;
229  input_sounding_reader.close();
230 
231  host_to_device(itime);
232  }
233 
234  /**
235  * @brief Hydrostatically integrate density and pressure for one sounding.
236  * @param itime Sounding time/file index to integrate.
237  */
238  void calc_rho_p (int itime)
239  {
240  /* Calculate density and pressure, roughly following the procedure in
241  * WRF dyn_em/module_initialize_ideal.F. We integrate hydrostatically
242  * from the surface up through the air column to get the dry density
243  * and moist pressure.
244  */
245  const int Ninp = size(itime);
246  pm_integ.resize(Ninp);
247  rhod_integ.resize(Ninp);
248 
249  // evaluate surface quantities (k=0): total pressure and dry air
253  RdoCp,
255 
256  amrex::Print() << "ideal sounding init: surface dry air density = "
257  << std::setprecision(15)
258  << rhod_integ[0] << " kg/m^3" << std::endl;
259 
260  // Note:
261  // p_dry = rho_d R_d T
262  // p_tot = rho_m R_d T_v
263  // = rho_d(1 + q_v) R_d T_v
264 
265 #if 0 // Printing
266  // In this absence of moisture, this moist profile will match the
267  // following dry profile
268  amrex::Print() << "z p_m rho_d theta qv U V" << std::endl;
269  amrex::Print() << z_inp_sound[itime][0]
270  << " " << pm_integ[0]
271  << " " << rhod_integ[0]
272  << " " << theta_inp_sound[itime][0]
273  << " " << qv_inp_sound[itime][0]
274  << " " << U_inp_sound[itime][0]
275  << " " << V_inp_sound[itime][0]
276  << std::endl;
277 #endif
278 
279  // integrate from surface to domain top
280  amrex::Real dz, F, C;
281  amrex::Real T_hi;
282  amrex::Real rho_tot_hi, rho_tot_lo;
283  for (int k=1; k < size(itime); ++k)
284  {
285  // Vertical grid spacing
286  dz = z_inp_sound[itime][k] - z_inp_sound[itime][k-1];
287 
288  // The Newton residual F below is a pressure, so the tolerance it is
289  // tested against has to be scaled to the pressure being solved for.
290  // The iterate cannot do better than settle into a one-ulp limit
291  // cycle about the true root, which puts a hard floor of roughly
292  // ulp(p) on |F| -- 3.6e-12 for a 1.6e4 Pa pressure in double, and
293  // ~1e-2 Pa at 1e5 Pa in single. The fixed 1e-12 Pa tolerance used
294  // here previously sat below that floor, so the iteration burned all
295  // of its allowed passes and reported a spurious non-convergence
296  // warning for a column that had in fact converged to roundoff.
297 #ifdef AMREX_USE_FLOAT
298  const amrex::Real tol = amrex::Real(4)
300  * pm_integ[k-1];
301 #else
302  const amrex::Real tol = amrex::Real(1.0e-12);
303 #endif
304 
305  // Establish known constant
306  rho_tot_lo = rhod_integ[k-1] * (one + qv_inp_sound[itime][k-1]);
307  C = -pm_integ[k-1] + myhalf*rho_tot_lo*CONST_GRAV*dz;
308 
309  // Initial guess and residual
310  pm_integ[k] = pm_integ[k-1];
311  T_hi = getTgivenPandTh(pm_integ[k], theta_inp_sound[itime][k], RdoCp);
313  pm_integ[k],
314  RdoCp,
315  qv_inp_sound[itime][k]);
316  rho_tot_hi = rhod_integ[k] * (one + qv_inp_sound[itime][k]);
317  F = pm_integ[k] + myhalf*rho_tot_hi*CONST_GRAV*dz + C;
318 
319  // Do iterations
320  if (std::abs(F)>tol) {
321  bool maintain_Th = true;
323  CONST_GRAV, C, theta_inp_sound[itime][k], T_hi,
324  qv_inp_sound[itime][k], qv_inp_sound[itime][k],
325  pm_integ[k], rhod_integ[k], F, maintain_Th);
326  }
327 
328 #if 0 // Printing
329  amrex::Print() << z_inp_sound[itime][k]
330  << " " << pm_integ[k]
331  << " " << rhod_integ[k]
332  << " " << theta_inp_sound[itime][k]
333  << " " << qv_inp_sound[itime][k]
334  << " " << U_inp_sound[itime][k]
335  << " " << V_inp_sound[itime][k]
336  << std::endl;
337 #endif
338  }
339  // Note: at this point, the surface pressure, density of the dry air
340  // column is stored in pm_integ[0], rhod_integ[0]
341 
342  // update
343  host_to_device(itime);
344  }
345 
346  /**
347  * @brief Compute density and pressure for an isentropic sounding profile.
348  * @param itime Sounding time/file index to integrate.
349  */
350  void calc_rho_p_isentropic (int itime)
351  {
352  /* Calculate density and pressure assuming isentropic (constant theta)
353  background conditions. This does not use Newton-Raphson iterations
354  to calculate rho and p.
355  */
356  const int Ninp = size(itime);
357  pm_integ.resize(Ninp);
358  rhod_integ.resize(Ninp);
359 
360  // evaluate surface quantities (k=0): total pressure and dry air
364  RdoCp,
366 
367  const amrex::Real th0 = theta_ref_inp_sound;
368  const amrex::Real p0_pow = std::pow(p_0, (Gamma-1)/Gamma);
369 
370  amrex::Print() << "isentropic sounding init: surface dry air density = "
371  << std::setprecision(15)
372  << rhod_integ[0] << " kg/m^3" << std::endl;
373 
374 #if 0 // Printing
375  // In this absence of moisture, this moist profile will match the
376  // following dry profile
377  amrex::Print() << "z p_m rho_d theta qv U V" << std::endl;
378  amrex::Print() << z_inp_sound[itime][0]
379  << " " << pm_integ[0]
380  << " " << rhod_integ[0]
381  << " " << theta_inp_sound[itime][0]
382  << " " << qv_inp_sound[itime][0]
383  << " " << U_inp_sound[itime][0]
384  << " " << V_inp_sound[itime][0]
385  << std::endl;
386 #endif
387 
388  // integrate from surface to domain top
389  amrex::Real dz;
390  for (int k=1; k < size(itime); ++k)
391  {
392  // Vertical grid spacing
393  dz = z_inp_sound[itime][k] - z_inp_sound[itime][k-1];
394 
395  // Isentropic pressure at next level
396  amrex::Real qvmean = (assume_dry) ? zero : myhalf*(qv_inp_sound[itime][k-1] + qv_inp_sound[itime][k]);
397  amrex::Real thm0 = th0 * (one + RvoRd * qvmean);
398  amrex::Real plo_pow = std::pow(pm_integ[k-1], (Gamma-1)/Gamma);
399  pm_integ[k] = (Gamma-1)/Gamma *
400  (-CONST_GRAV * p0_pow / (R_d * thm0) * (1 + qvmean) * dz
401  + Gamma/(Gamma-1) * plo_pow);
402  pm_integ[k] = std::pow(pm_integ[k], Gamma / (Gamma-1));
403 
404  // Note: The pressure calculated here is copied to `p_inp_sound_d`
405  // and currently not used.
406 
407  // Get corresponding dry air density
409  pm_integ[k],
410  RdoCp,
411  qv_inp_sound[itime][k]);
412 
413 #if 0 // Printing
414  amrex::Print() << z_inp_sound[itime][k]
415  << " " << pm_integ[k]
416  << " " << rhod_integ[k]
417  << " " << theta_inp_sound[itime][k]
418  << " " << qv_inp_sound[itime][k]
419  << " " << U_inp_sound[itime][k]
420  << " " << V_inp_sound[itime][k]
421  << std::endl;
422 #endif
423  }
424  // Note: at this point, the surface pressure, density of the dry air
425  // column is stored in pm_integ[0], rhod_integ[0]
426 
427  // update
428  host_to_device(itime);
429  }
430 
431  /**
432  * @brief Copy host sounding arrays for one time index to device storage.
433  * @param itime Sounding time/file index to copy.
434  */
435  void host_to_device (int itime)
436  {
437  const int Ninp = size(itime);
438  z_inp_sound_d[itime].resize(Ninp);
439  theta_inp_sound_d[itime].resize(Ninp);
440  qv_inp_sound_d[itime].resize(Ninp);
441  U_inp_sound_d[itime].resize(Ninp);
442  V_inp_sound_d[itime].resize(Ninp);
443 
444  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
445  z_inp_sound[itime].begin(), z_inp_sound[itime].end(),
446  z_inp_sound_d[itime].begin());
447  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
448  theta_inp_sound[itime].begin(), theta_inp_sound[itime].end(),
449  theta_inp_sound_d[itime].begin());
450  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
451  qv_inp_sound[itime].begin(), qv_inp_sound[itime].end(),
452  qv_inp_sound_d[itime].begin());
453  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
454  U_inp_sound[itime].begin(), U_inp_sound[itime].end(),
455  U_inp_sound_d[itime].begin());
456  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
457  V_inp_sound[itime].begin(), V_inp_sound[itime].end(),
458  V_inp_sound_d[itime].begin());
459 
460  if (rhod_integ.size() > 0)
461  {
462  //amrex::Print() << "Copying rho_d, p_d to device" << std::endl;
463  rho_inp_sound_d.resize(size(itime)+2);
464  p_inp_sound_d.resize(size(itime)+2);
465  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
466  rhod_integ.begin(), rhod_integ.end(),
467  rho_inp_sound_d.begin());
468  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
469  pm_integ.begin(), pm_integ.end(),
470  p_inp_sound_d.begin());
471  }
472  }
473 
474  /**
475  * @brief Return the number of interpolated samples for one sounding time.
476  * @param itime Sounding time/file index to query.
477  * @return Number of stored sounding profile levels.
478  */
479  int size (int itime) const
480  {
482  AMREX_ALWAYS_ASSERT(z_inp_sound[itime].size() == qv_inp_sound[itime].size());
483  AMREX_ALWAYS_ASSERT(z_inp_sound[itime].size() == U_inp_sound[itime].size());
484  AMREX_ALWAYS_ASSERT(z_inp_sound[itime].size() == V_inp_sound[itime].size());
485  return z_inp_sound[itime].size();
486  }
487 
488  void set_start_time (amrex::Real start_time)
489  {
490  for (int i = 0; i < n_sounding_times; i++)
491  {
492  input_sounding_time[i] += start_time;
493  }
494  }
495 
496  // Members
497  int ntimes;
498 
499  amrex::Real tau_nudging = amrex::Real(5.0); // time scale used for nudging
500 
501  amrex::Vector<std::string> input_sounding_file = {};
502  amrex::Vector<amrex::Real> input_sounding_time = {};
505 
506  bool assume_dry{false};
507 
508  // - read from file
510 
511  // This is a vector (over time) of Vectors
512  amrex::Vector<amrex::Vector<amrex::Real>> z_inp_sound, theta_inp_sound, qv_inp_sound, U_inp_sound, V_inp_sound;
513 
514  // This is a vector (over time) of DeviceVectors
515  amrex::Vector<amrex::Gpu::DeviceVector<amrex::Real>> z_inp_sound_d, theta_inp_sound_d, qv_inp_sound_d, U_inp_sound_d, V_inp_sound_d;
516 
517  // - moist profiles
518  amrex::Vector<amrex::Real> pm_integ; // from integrating up air column
519  // - dry profiles
520  amrex::Vector<amrex::Real> rhod_integ; // from integrating down air column
521  // - to set solution fields
522  amrex::Gpu::DeviceVector<amrex::Real> p_inp_sound_d, rho_inp_sound_d;
523 };
524 
525 /**
526  * @var InputSoundingData::ntimes
527  * @brief Number of sounding times represented in the arrays.
528  * @var InputSoundingData::tau_nudging
529  * @brief Nudging time scale for relaxing toward input sounding data.
530  * @var InputSoundingData::input_sounding_file
531  * @brief Sounding file path for each configured sounding time.
532  * @var InputSoundingData::input_sounding_time
533  * @brief Physical time associated with each input sounding file.
534  * @var InputSoundingData::n_sounding_files
535  * @brief Number of input sounding files.
536  * @var InputSoundingData::n_sounding_times
537  * @brief Number of input sounding times.
538  * @var InputSoundingData::assume_dry
539  * @brief Whether to ignore moisture in isentropic pressure integration.
540  * @var InputSoundingData::press_ref_inp_sound
541  * @brief Surface reference pressure read from the input sounding [Pa].
542  * @var InputSoundingData::theta_ref_inp_sound
543  * @brief Surface reference potential temperature read from the input sounding [K].
544  * @var InputSoundingData::qv_ref_inp_sound
545  * @brief Surface reference water vapor mixing ratio read from the input sounding [kg/kg].
546  * @var InputSoundingData::z_inp_sound
547  * @brief Interpolated host sounding heights for each sounding time [m].
548  * @var InputSoundingData::theta_inp_sound
549  * @brief Interpolated host potential temperature profiles [K].
550  * @var InputSoundingData::qv_inp_sound
551  * @brief Interpolated host water vapor mixing ratio profiles [kg/kg].
552  * @var InputSoundingData::U_inp_sound
553  * @brief Interpolated host x velocity profiles.
554  * @var InputSoundingData::V_inp_sound
555  * @brief Interpolated host y velocity profiles.
556  * @var InputSoundingData::z_inp_sound_d
557  * @brief Device sounding heights for each sounding time [m].
558  * @var InputSoundingData::theta_inp_sound_d
559  * @brief Device potential temperature profiles [K].
560  * @var InputSoundingData::qv_inp_sound_d
561  * @brief Device water vapor mixing ratio profiles [kg/kg].
562  * @var InputSoundingData::U_inp_sound_d
563  * @brief Device x velocity profiles.
564  * @var InputSoundingData::V_inp_sound_d
565  * @brief Device y velocity profiles.
566  * @var InputSoundingData::pm_integ
567  * @brief Hydrostatically integrated moist pressure profile.
568  * @var InputSoundingData::rhod_integ
569  * @brief Hydrostatically integrated dry density profile.
570  * @var InputSoundingData::p_inp_sound_d
571  * @brief Device pressure profile used to initialize solution fields.
572  * @var InputSoundingData::rho_inp_sound_d
573  * @brief Device density profile used to initialize solution fields.
574  */
575 #endif
constexpr amrex::Real p_0
Definition: ERF_Constants.H:69
constexpr amrex::Real RvoRd
Definition: ERF_Constants.H:59
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:72
constexpr amrex::Real R_d
Definition: ERF_Constants.H:49
constexpr amrex::Real RdoCp
Definition: ERF_Constants.H:57
constexpr amrex::Real Gamma
Definition: ERF_Constants.H:70
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhogivenThetaPress(const amrex::Real th, const amrex::Real p, const amrex::Real rdOcp, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:96
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenPandTh(const amrex::Real P, const amrex::Real th, const amrex::Real rdOcp)
Definition: ERF_EOS.H:32
const Real ztop
Definition: ERF_InitCustomPertVels_ParticleTests.H:4
ParmParse pp("prob")
const int khi
Definition: ERF_InitCustomPert_Bubble.H:25
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
const int klo
Definition: ERF_InitCustomPert_KE.H:9
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real interpolate_1d(const amrex::Real *alpha, const amrex::Real *beta, const amrex::Real alpha_interp, const int alpha_size)
Definition: ERF_Interpolation_1D.H:14
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
amrex::Real Real
Definition: ERF_ShocInterface.H:19
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Newton_Raphson_hse(const Real &m_tol, const Real &RdoCp, const Real &dz, const Real &g, const Real &C, const Real &Th, const Real &T, const Real &qt, const Real &qv, Real &P, Real &rd, Real &F, const bool &maintain_Th)
Definition: ERF_HSEUtils.H:46
@ theta
Definition: ERF_Kessler.H:26
@ qv
Definition: ERF_Kessler.H:31
@ U
Definition: ERF_IndexDefines.H:166
@ V
Definition: ERF_IndexDefines.H:167
@ dz
Definition: ERF_AdvanceWDM6.cpp:273
real(c_double), parameter epsilon
Definition: ERF_module_model_constants.F90:12
Definition: ERF_InputSoundingData.H:23
void host_to_device(int itime)
Copy host sounding arrays for one time index to device storage.
Definition: ERF_InputSoundingData.H:435
amrex::Gpu::DeviceVector< amrex::Real > rho_inp_sound_d
Device density profile used to initialize solution fields.
Definition: ERF_InputSoundingData.H:522
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > qv_inp_sound_d
Device water vapor mixing ratio profiles [kg/kg].
Definition: ERF_InputSoundingData.H:515
amrex::Vector< amrex::Vector< amrex::Real > > V_inp_sound
Interpolated host y velocity profiles.
Definition: ERF_InputSoundingData.H:512
amrex::Vector< amrex::Real > pm_integ
Hydrostatically integrated moist pressure profile.
Definition: ERF_InputSoundingData.H:518
amrex::Real qv_ref_inp_sound
Surface reference water vapor mixing ratio read from the input sounding [kg/kg].
Definition: ERF_InputSoundingData.H:509
amrex::Vector< amrex::Real > input_sounding_time
Physical time associated with each input sounding file.
Definition: ERF_InputSoundingData.H:502
void resize_arrays()
Resize host and device sounding arrays for the configured sounding times.
Definition: ERF_InputSoundingData.H:67
bool assume_dry
Whether to ignore moisture in isentropic pressure integration.
Definition: ERF_InputSoundingData.H:506
int n_sounding_times
Number of input sounding times.
Definition: ERF_InputSoundingData.H:504
amrex::Vector< amrex::Vector< amrex::Real > > theta_inp_sound
Interpolated host potential temperature profiles [K].
Definition: ERF_InputSoundingData.H:512
amrex::Vector< amrex::Real > rhod_integ
Hydrostatically integrated dry density profile.
Definition: ERF_InputSoundingData.H:520
void set_start_time(amrex::Real start_time)
Definition: ERF_InputSoundingData.H:488
void calc_rho_p(int itime)
Hydrostatically integrate density and pressure for one sounding.
Definition: ERF_InputSoundingData.H:238
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > V_inp_sound_d
Device y velocity profiles.
Definition: ERF_InputSoundingData.H:515
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > z_inp_sound_d
Device sounding heights for each sounding time [m].
Definition: ERF_InputSoundingData.H:515
void calc_rho_p_isentropic(int itime)
Compute density and pressure for an isentropic sounding profile.
Definition: ERF_InputSoundingData.H:350
int ntimes
Number of sounding times represented in the arrays.
Definition: ERF_InputSoundingData.H:497
void read_from_file(const amrex::Geometry &geom, const amrex::Vector< amrex::Real > &zlevels_stag, int itime, bool is_moist)
Read and interpolate one input sounding profile onto level-0 heights.
Definition: ERF_InputSoundingData.H:91
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > theta_inp_sound_d
Device potential temperature profiles [K].
Definition: ERF_InputSoundingData.H:515
amrex::Vector< amrex::Vector< amrex::Real > > U_inp_sound
Interpolated host x velocity profiles.
Definition: ERF_InputSoundingData.H:512
amrex::Real theta_ref_inp_sound
Surface reference potential temperature read from the input sounding [K].
Definition: ERF_InputSoundingData.H:509
amrex::Vector< std::string > input_sounding_file
Sounding file path for each configured sounding time.
Definition: ERF_InputSoundingData.H:501
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > U_inp_sound_d
Device x velocity profiles.
Definition: ERF_InputSoundingData.H:515
amrex::Gpu::DeviceVector< amrex::Real > p_inp_sound_d
Device pressure profile used to initialize solution fields.
Definition: ERF_InputSoundingData.H:522
amrex::Vector< amrex::Vector< amrex::Real > > z_inp_sound
Interpolated host sounding heights for each sounding time [m].
Definition: ERF_InputSoundingData.H:512
amrex::Vector< amrex::Vector< amrex::Real > > qv_inp_sound
Interpolated host water vapor mixing ratio profiles [kg/kg].
Definition: ERF_InputSoundingData.H:512
amrex::Real tau_nudging
Nudging time scale for relaxing toward input sounding data.
Definition: ERF_InputSoundingData.H:499
InputSoundingData()
Construct input sounding metadata from the ERF input namespace.
Definition: ERF_InputSoundingData.H:28
amrex::Real press_ref_inp_sound
Surface reference pressure read from the input sounding [Pa].
Definition: ERF_InputSoundingData.H:509
int size(int itime) const
Return the number of interpolated samples for one sounding time.
Definition: ERF_InputSoundingData.H:479
int n_sounding_files
Number of input sounding files.
Definition: ERF_InputSoundingData.H:503