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 <string>
5 #include <iostream>
6 
7 #include <AMReX_ParmParse.H>
8 #include <AMReX_Print.H>
9 #include <AMReX_Gpu.H>
10 #include <AMReX_Geometry.H>
11 
12 #include <ERF_EOS.H>
13 #include <ERF_Constants.H>
14 #include <ERF_Interpolation_1D.H>
15 #include <ERF_HSE_utils.H>
16 
17 /**
18  * Data structure storing input sounding data. Also
19  * handles reading the input file for sounding data and
20  * hydrostatic column integration.
21  */
23 public:
25  {
26  amrex::ParmParse pp("erf");
27  pp.query("tau_nudging", tau_nudging);
28 
29  // Read in input_sounding filename
30  n_sounding_files = pp.countval("input_sounding_file");
31  if (n_sounding_files > 0) {
33  pp.queryarr("input_sounding_file", input_sounding_file, 0, n_sounding_files);
34  } else {
35  n_sounding_files = 1;
37  input_sounding_file[0] = "input_sounding";
38  }
39 
40  // Read in input_sounding times
41  n_sounding_times = pp.countval("input_sounding_time");
42 
43  if (n_sounding_times > 0) {
45  pp.queryarr("input_sounding_time", input_sounding_time, 0, n_sounding_times);
46  } else {
47  n_sounding_times = 1;
49  input_sounding_time[0] = 0.0;
50  }
51 
52  // If we have more files than times or times than files we just use the minimum
53  int n = std::min(n_sounding_times, n_sounding_files);
54  n_sounding_files = n;
55  n_sounding_times = n;
56  input_sounding_file.resize(n);
57  input_sounding_time.resize(n);
58  }
59 
60  void resize_arrays ()
61  {
63 
64  z_inp_sound.resize(ntimes);
65  theta_inp_sound.resize(ntimes);
66  qv_inp_sound.resize(ntimes);
67  U_inp_sound.resize(ntimes);
68  V_inp_sound.resize(ntimes);
69 
70  z_inp_sound_d.resize(ntimes);
71  theta_inp_sound_d.resize(ntimes);
72  qv_inp_sound_d.resize(ntimes);
73  U_inp_sound_d.resize(ntimes);
74  V_inp_sound_d.resize(ntimes);
75  }
76 
77  void read_from_file (const amrex::Geometry &geom,
78  const amrex::Vector<amrex::Real>& zlevels_stag,
79  int itime)
80  {
81  const int klo = 0;
82  const int khi = geom.Domain().bigEnd()[AMREX_SPACEDIM-1];
83  const int Nz = geom.Domain().size()[AMREX_SPACEDIM-1];
84  const amrex::Real dz = geom.CellSize()[AMREX_SPACEDIM-1];
85 
86  const bool use_terrain = (zlevels_stag.size() > 0);
87  const amrex::Real zbot = (use_terrain) ? zlevels_stag[klo] : geom.ProbLo(AMREX_SPACEDIM-1);
88  const amrex::Real ztop = (use_terrain) ? zlevels_stag[khi+1] : geom.ProbHi(AMREX_SPACEDIM-1);
89 
90  z_inp_sound[itime].resize(Nz+2);
91  theta_inp_sound[itime].resize(Nz+2);
92  qv_inp_sound[itime].resize(Nz+2);
93  U_inp_sound[itime].resize(Nz+2);
94  V_inp_sound[itime].resize(Nz+2);
95 
96  // Read the input_sounding file
97  amrex::Print() << "input_sounding file location : " << input_sounding_file[itime] << std::endl;
98  std::ifstream input_sounding_reader(input_sounding_file[itime]);
99  if(!input_sounding_reader.is_open()) {
100  amrex::Error("Error opening the input_sounding file\n");
101  }
102  else {
103  // Read the contents of the input_sounding file
104  amrex::Print() << "Successfully opened the input_sounding file. Now reading... " << std::endl;
105  std::string line;
106 
107  // First, read the input data into temp vectors; then, interpolate vectors to the
108  // domain lo/hi and cell centers (from level 0)
109  amrex::Vector<amrex::Real> z_inp_sound_tmp, theta_inp_sound_tmp, qv_inp_sound_tmp,
110  U_inp_sound_tmp, V_inp_sound_tmp;
111 
112  // Read surface quantities from the first line
113  std::getline(input_sounding_reader, line);
114  std::istringstream iss(line);
116  press_ref_inp_sound *= 100; // convert from hPa to Pa
117  qv_ref_inp_sound *= 0.001; // convert from g/kg to kg/kg
118 
119  // Add surface
120  z_inp_sound_tmp.push_back(zbot); // height above sea level [m]
121  theta_inp_sound_tmp.push_back(theta_ref_inp_sound);
122  qv_inp_sound_tmp.push_back(qv_ref_inp_sound);
123  U_inp_sound_tmp.push_back(0);
124  V_inp_sound_tmp.push_back(0);
125 
126  // Read the vertical profile at each given height
127  amrex::Real z, theta, qv, U, V;
128  while(std::getline(input_sounding_reader, line)) {
129  std::istringstream iss_z(line);
130  iss_z >> z >> theta >> qv >> U >> V;
131  if (z == zbot) {
132  AMREX_ALWAYS_ASSERT(theta == theta_inp_sound_tmp[0]);
133  AMREX_ALWAYS_ASSERT(qv*0.001 == qv_inp_sound_tmp[0]); // convert from g/kg to kg/kg
134  U_inp_sound_tmp[0] = U;
135  V_inp_sound_tmp[0] = V;
136  } else {
137  AMREX_ALWAYS_ASSERT(z > z_inp_sound_tmp[z_inp_sound_tmp.size()-1]); // sounding is increasing in height
138  z_inp_sound_tmp.push_back(z);
139  theta_inp_sound_tmp.push_back(theta);
140  qv_inp_sound_tmp.push_back(qv*0.001); // convert from g/kg to kg/kg
141  U_inp_sound_tmp.push_back(U);
142  V_inp_sound_tmp.push_back(V);
143  if (z >= ztop) break;
144  }
145  }
146 
147  // At this point, we have an input_sounding from zbot up to
148  // z_inp_sound_tmp[N-1] >= ztop. Now, interpolate to grid level 0 heights
149  const int Ninp = z_inp_sound_tmp.size();
150  z_inp_sound[itime][0] = zbot;
151  theta_inp_sound[itime][0] = theta_inp_sound_tmp[0];
152  qv_inp_sound[itime][0] = qv_inp_sound_tmp[0];
153  U_inp_sound[itime][0] = U_inp_sound_tmp[0];
154  V_inp_sound[itime][0] = V_inp_sound_tmp[0];
155  for (int k=0; k < Nz; ++k) {
156  z_inp_sound[itime][k+1] = (use_terrain) ? 0.5 * (zlevels_stag[k] + zlevels_stag[k+1])
157  : zbot + (k + 0.5) * dz;
158  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);
159  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);
160  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);
161  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);
162  }
163  z_inp_sound[itime][Nz+1] = ztop;
164  theta_inp_sound[itime][Nz+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), theta_inp_sound_tmp.dataPtr(), ztop, Ninp);
165  qv_inp_sound[itime][Nz+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), qv_inp_sound_tmp.dataPtr(), ztop, Ninp);
166  U_inp_sound[itime][Nz+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), U_inp_sound_tmp.dataPtr(), ztop, Ninp);
167  V_inp_sound[itime][Nz+1] = interpolate_1d(z_inp_sound_tmp.dataPtr(), V_inp_sound_tmp.dataPtr(), ztop, Ninp);
168  }
169 
170  amrex::Print() << "Successfully read the " << itime << "th input_sounding file..." << std::endl;
171  input_sounding_reader.close();
172 
173  host_to_device(itime);
174  }
175 
176  void calc_rho_p (int itime)
177  {
178  /* Calculate density and pressure, roughly following the procedure in
179  * WRF dyn_em/module_initialize_ideal.F. We integrate hydrostatically
180  * from the surface up through the air column to get the dry density
181  * and moist pressure.
182  */
183  const amrex::Real tol = 1.0e-12;
184  const int Ninp = size(itime);
185  pm_integ.resize(Ninp);
186  rhod_integ.resize(Ninp);
187 
188  // evaluate surface quantities (k=0): total pressure and dry air
192  R_d/Cp_d,
194 
195  amrex::Print() << "ideal sounding init: surface density of moist air = "
196  << rhod_integ[0] << " kg/m^3" << std::endl;
197 
198  // Note:
199  // p_dry = rho_d R_d T
200  // p_tot = rho_m R_d T_v
201  // = rho_d(1 + q_v) R_d T_v
202 
203 #if 0 // Printing
204  // In this absence of moisture, this moist profile will match the
205  // following dry profile
206  amrex::Print() << "z p_m rho_d theta qv U V" << std::endl;
207  amrex::Print() << z_inp_sound[itime][0]
208  << " " << pm_integ[0]
209  << " " << rhod_integ[0]
210  << " " << theta_inp_sound[itime][0]
211  << " " << qv_inp_sound[itime][0]
212  << " " << U_inp_sound[itime][0]
213  << " " << V_inp_sound[itime][0]
214  << std::endl;
215 #endif
216 
217  // integrate from surface to domain top
218  amrex::Real dz, F, C;
219  amrex::Real rho_tot_hi, rho_tot_lo;
220  for (int k=1; k < size(itime); ++k)
221  {
222  // Vertical grid spacing
223  dz = z_inp_sound[itime][k] - z_inp_sound[itime][k-1];
224 
225  // Establish known constant
226  rho_tot_lo = rhod_integ[k-1] * (1. + qv_inp_sound[itime][k-1]);
227  C = -pm_integ[k-1] + 0.5*rho_tot_lo*CONST_GRAV*dz;
228 
229  // Initial guess and residual
230  pm_integ[k] = pm_integ[k-1];
232  pm_integ[k],
233  R_d/Cp_d,
234  qv_inp_sound[itime][k]);
235  rho_tot_hi = rhod_integ[k] * (1. + qv_inp_sound[itime][k]);
236  F = pm_integ[k] + 0.5*rho_tot_hi*CONST_GRAV*dz + C;
237 
238  // Do iterations
239  if (std::abs(F)>tol) HSEutils::Newton_Raphson_hse(tol, R_d/Cp_d, dz,
240  CONST_GRAV, C, theta_inp_sound[itime][k],
241  qv_inp_sound[itime][k], qv_inp_sound[itime][k],
242  pm_integ[k], rhod_integ[k], F);
243 #if 0 // Printing
244  amrex::Print() << z_inp_sound[itime][k]
245  << " " << pm_integ[k]
246  << " " << rhod_integ[k]
247  << " " << theta_inp_sound[itime][k]
248  << " " << qv_inp_sound[itime][k]
249  << " " << U_inp_sound[itime][k]
250  << " " << V_inp_sound[itime][k]
251  << std::endl;
252 #endif
253  }
254  // Note: at this point, the surface pressure, density of the dry air
255  // column is stored in pm_integ[0], rhod_integ[0]
256 
257  // update
258  host_to_device(itime);
259  }
260 
261  void host_to_device (int itime)
262  {
263  const int Ninp = size(itime);
264  z_inp_sound_d[itime].resize(Ninp);
265  theta_inp_sound_d[itime].resize(Ninp);
266  qv_inp_sound_d[itime].resize(Ninp);
267  U_inp_sound_d[itime].resize(Ninp);
268  V_inp_sound_d[itime].resize(Ninp);
269 
270  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
271  z_inp_sound[itime].begin(), z_inp_sound[itime].end(),
272  z_inp_sound_d[itime].begin());
273  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
274  theta_inp_sound[itime].begin(), theta_inp_sound[itime].end(),
275  theta_inp_sound_d[itime].begin());
276  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
277  qv_inp_sound[itime].begin(), qv_inp_sound[itime].end(),
278  qv_inp_sound_d[itime].begin());
279  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
280  U_inp_sound[itime].begin(), U_inp_sound[itime].end(),
281  U_inp_sound_d[itime].begin());
282  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
283  V_inp_sound[itime].begin(), V_inp_sound[itime].end(),
284  V_inp_sound_d[itime].begin());
285 
286  if (rhod_integ.size() > 0)
287  {
288  //amrex::Print() << "Copying rho_d, p_d to device" << std::endl;
289  rho_inp_sound_d.resize(size(itime)+2);
290  p_inp_sound_d.resize(size(itime)+2);
291  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
292  rhod_integ.begin(), rhod_integ.end(),
293  rho_inp_sound_d.begin());
294  amrex::Gpu::copy(amrex::Gpu::hostToDevice,
295  pm_integ.begin(), pm_integ.end(),
296  p_inp_sound_d.begin());
297  }
298  }
299 
300  int size (int itime) const
301  {
302  AMREX_ALWAYS_ASSERT(z_inp_sound[itime].size() == theta_inp_sound[itime].size());
303  AMREX_ALWAYS_ASSERT(z_inp_sound[itime].size() == qv_inp_sound[itime].size());
304  AMREX_ALWAYS_ASSERT(z_inp_sound[itime].size() == U_inp_sound[itime].size());
305  AMREX_ALWAYS_ASSERT(z_inp_sound[itime].size() == V_inp_sound[itime].size());
306  return z_inp_sound[itime].size();
307  }
308 
309  // Members
310  int ntimes;
311 
312  amrex::Real tau_nudging = 5.0; // time scale used for nudging
313 
314  amrex::Vector<std::string> input_sounding_file = {};
315  amrex::Vector<amrex::Real> input_sounding_time = {};
318 
319  // - read from file
321 
322  // This is a vector (over time) of Vectors
323  amrex::Vector<amrex::Vector<amrex::Real>> z_inp_sound, theta_inp_sound, qv_inp_sound, U_inp_sound, V_inp_sound;
324 
325  // This is a vector (over time) of DeviceVectors
326  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;
327 
328  // - moist profiles
329  amrex::Vector<amrex::Real> pm_integ; // from integrating up air column
330  // - dry profiles
331  amrex::Vector<amrex::Real> rhod_integ; // from integrating down air column
332  // - to set solution fields
333  amrex::Gpu::DeviceVector<amrex::Real> p_inp_sound_d, rho_inp_sound_d;
334 };
335 #endif
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:12
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
constexpr amrex::Real R_d
Definition: ERF_Constants.H:10
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhogivenThetaPress(const amrex::Real theta, const amrex::Real p, const amrex::Real rdOcp, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:89
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:12
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_Microphysics_Utils.H:203
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 &T, const Real &qt, const Real &qv, Real &P, Real &rd, Real &F)
Definition: ERF_HSE_utils.H:35
@ theta
Definition: ERF_MM5.H:20
@ qv
Definition: ERF_Kessler.H:36
@ U
Definition: ERF_IndexDefines.H:97
@ V
Definition: ERF_IndexDefines.H:98
Definition: ERF_InputSoundingData.H:22
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > theta_inp_sound_d
Definition: ERF_InputSoundingData.H:326
amrex::Real press_ref_inp_sound
Definition: ERF_InputSoundingData.H:320
amrex::Gpu::DeviceVector< amrex::Real > p_inp_sound_d
Definition: ERF_InputSoundingData.H:333
void host_to_device(int itime)
Definition: ERF_InputSoundingData.H:261
amrex::Real theta_ref_inp_sound
Definition: ERF_InputSoundingData.H:320
amrex::Vector< amrex::Real > pm_integ
Definition: ERF_InputSoundingData.H:329
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > V_inp_sound_d
Definition: ERF_InputSoundingData.H:326
amrex::Vector< amrex::Real > rhod_integ
Definition: ERF_InputSoundingData.H:331
void resize_arrays()
Definition: ERF_InputSoundingData.H:60
int n_sounding_times
Definition: ERF_InputSoundingData.H:317
amrex::Vector< amrex::Real > input_sounding_time
Definition: ERF_InputSoundingData.H:315
int n_sounding_files
Definition: ERF_InputSoundingData.H:316
amrex::Real tau_nudging
Definition: ERF_InputSoundingData.H:312
void read_from_file(const amrex::Geometry &geom, const amrex::Vector< amrex::Real > &zlevels_stag, int itime)
Definition: ERF_InputSoundingData.H:77
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > qv_inp_sound_d
Definition: ERF_InputSoundingData.H:326
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > z_inp_sound_d
Definition: ERF_InputSoundingData.H:326
amrex::Vector< std::string > input_sounding_file
Definition: ERF_InputSoundingData.H:314
amrex::Gpu::DeviceVector< amrex::Real > rho_inp_sound_d
Definition: ERF_InputSoundingData.H:333
amrex::Vector< amrex::Vector< amrex::Real > > theta_inp_sound
Definition: ERF_InputSoundingData.H:323
void calc_rho_p(int itime)
Definition: ERF_InputSoundingData.H:176
amrex::Real qv_ref_inp_sound
Definition: ERF_InputSoundingData.H:320
amrex::Vector< amrex::Vector< amrex::Real > > z_inp_sound
Definition: ERF_InputSoundingData.H:323
amrex::Vector< amrex::Vector< amrex::Real > > qv_inp_sound
Definition: ERF_InputSoundingData.H:323
amrex::Vector< amrex::Gpu::DeviceVector< amrex::Real > > U_inp_sound_d
Definition: ERF_InputSoundingData.H:326
amrex::Vector< amrex::Vector< amrex::Real > > U_inp_sound
Definition: ERF_InputSoundingData.H:323
amrex::Vector< amrex::Vector< amrex::Real > > V_inp_sound
Definition: ERF_InputSoundingData.H:323
InputSoundingData()
Definition: ERF_InputSoundingData.H:24
int size(int itime) const
Definition: ERF_InputSoundingData.H:300
int ntimes
Definition: ERF_InputSoundingData.H:310