ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_InitFromInputSounding.cpp File Reference
#include <ERF.H>
#include <ERF_EOS.H>
#include <ERF_Constants.H>
#include <ERF_Utils.H>
#include <ERF_ProbCommon.H>
Include dependency graph for ERF_InitFromInputSounding.cpp:

Functions

void init_bx_scalars_from_input_sounding (const Box &bx, Array4< Real > const &state, GeometryData const &geomdata, Array4< const Real > const &z_cc_arr, const bool &l_moist, InputSoundingData const &inputSoundingData)
 
void init_bx_scalars_from_input_sounding_hse (const Box &bx, Array4< Real > const &state, Array4< Real > const &r_hse_arr, Array4< Real > const &p_hse_arr, Array4< Real > const &pi_hse_arr, Array4< Real > const &th_hse_arr, Array4< Real > const &qv_hse_arr, GeometryData const &geomdata, Array4< const Real > const &z_cc_arr, const Real &l_gravity, const Real &l_rdOcp, const bool &l_moist, InputSoundingData const &inputSoundingData, const bool &l_isentropic, const int &ngz)
 
void init_bx_velocities_from_input_sounding (const Box &bx, Array4< Real > const &x_vel, Array4< Real > const &y_vel, Array4< Real > const &z_vel, GeometryData const &geomdata, Array4< const Real > const &z_nd_arr, InputSoundingData const &inputSoundingData)
 

Function Documentation

◆ init_bx_scalars_from_input_sounding()

void init_bx_scalars_from_input_sounding ( const Box &  bx,
Array4< Real > const &  state,
GeometryData const &  geomdata,
Array4< const Real > const &  z_cc_arr,
const bool &  l_moist,
InputSoundingData const &  inputSoundingData 
)

Box level wrapper for initializing scalar data from input sounding data.

Parameters
bxBox specifying the indices we are initializing
stateArray4 specifying the state data we are to initialize
geomdataGeometryData object specifying the domain geometry
inputSoundingDataInputSoundingData object we are to initialize from
185 {
186  const Real* z_inp_sound = inputSoundingData.z_inp_sound_d[0].dataPtr();
187  const Real* theta_inp_sound = inputSoundingData.theta_inp_sound_d[0].dataPtr();
188  const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d[0].dataPtr();
189  const int inp_sound_size = inputSoundingData.size(0);
190 
191  // Geometry
192  const Real* prob_lo = geomdata.ProbLo();
193  const Real* dx = geomdata.CellSize();
194  const Real z_lo = prob_lo[2];
195  const Real dz = dx[2];
196 
197  // We want to set the lateral BC values, too
198  Box gbx = bx; // Copy constructor
199  gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
200 
201  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
202  const Real z = (z_cc_arr) ? z_cc_arr(i,j,k)
203  : z_lo + (k + 0.5) * dz;
204 
205  Real rho_0 = 1.0;
206 
207  // Set the density
208  state(i, j, k, Rho_comp) = rho_0;
209 
210  // Initial Rho0*Theta0
211  state(i, j, k, RhoTheta_comp) = rho_0 * interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size);
212 
213  // Initialize all scalars to 0.
214  for (int n = 0; n < NSCALARS; n++) {
215  state(i, j, k, RhoScalar_comp+n) = 0;
216  }
217 
218  // total nonprecipitating water (Q1) == water vapor (Qv), i.e., there is no cloud water or cloud ice
219  if (l_moist) {
220  state(i, j, k, RhoQ1_comp) = rho_0 * interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size);
221  }
222  });
223 }
#define RhoScalar_comp
Definition: ERF_IndexDefines.H:40
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define NSCALARS
Definition: ERF_IndexDefines.H:16
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
const auto dx
Definition: ERF_InitCustomPertVels_ParticleTests.H:15
ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const auto prob_lo=geomdata.ProbLo();const auto dx=geomdata.CellSize();const Real x=(prob_lo[0]+(i+0.5) *dx[0])/mf_m(i, j, 0);const Real z=z_cc(i, j, k);Real L=std::sqrt(std::pow((x - x_c)/x_r, 2)+std::pow((z - z_c)/z_r, 2));if(L<=1.0) { Real dT=T_pert *(std::cos(PI *L)+1.0)/2.0;Real Tbar_hse=p_hse(i, j, k)/(R_d *r_hse(i, j, k));Real theta_perturbed=(Tbar_hse+dT) *std::pow(p_0/p_hse(i, j, k), rdOcp);Real theta_0=(Tbar_hse) *std::pow(p_0/p_hse(i, j, k), rdOcp);if(const_rho) { state_pert(i, j, k, RhoTheta_comp)=r_hse(i, j, k) *(theta_perturbed - theta_0);} else { state_pert(i, j, k, Rho_comp)=getRhoThetagivenP(p_hse(i, j, k))/theta_perturbed - r_hse(i, j, k);} } })
const amrex::Real * prob_lo
Definition: ERF_InitCustomPert_IsentropicVortex.H:16
amrex::Real rho_0
Definition: ERF_InitCustomPert_IsentropicVortex.H:33
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::Real Real
Definition: ERF_ShocInterface.H:19

Referenced by ERF::init_from_input_sounding().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ init_bx_scalars_from_input_sounding_hse()

void init_bx_scalars_from_input_sounding_hse ( const Box &  bx,
Array4< Real > const &  state,
Array4< Real > const &  r_hse_arr,
Array4< Real > const &  p_hse_arr,
Array4< Real > const &  pi_hse_arr,
Array4< Real > const &  th_hse_arr,
Array4< Real > const &  qv_hse_arr,
GeometryData const &  geomdata,
Array4< const Real > const &  z_cc_arr,
const Real ,
const Real l_rdOcp,
const bool &  l_moist,
InputSoundingData const &  inputSoundingData,
const bool &  l_isentropic,
const int &  ngz 
)

Box level wrapper for initializing scalar and hydrostatic base state data from input sounding data.

Parameters
bxBox specifying the indices we are initializing
stateArray4 specifying the state data we are to initialize
r_hse_arrArray4 specifying the density HSE base state data we are to initialize
p_hse_arrArray4 specifying the pressure HSE base state data we are to initialize
pi_hse_arrArray4 specifying the Exner pressure HSE base state data we are to initialize
th_hse_arrArray4 specifying the base state potential temperature we are to initialize
geomdataGeometryData object specifying the domain geometry
l_gravityReal number specifying the gravitational acceleration constant
l_rdOcpReal number specifying the Rhydberg constant ($R_d$) divided by specific heat at constant pressure ($c_p$)
inputSoundingDataInputSoundingData object we are to initialize from
256 {
257  const Real* z_inp_sound = inputSoundingData.z_inp_sound_d[0].dataPtr();
258  const Real* rho_inp_sound = inputSoundingData.rho_inp_sound_d.dataPtr();
259  const Real* theta_inp_sound = inputSoundingData.theta_inp_sound_d[0].dataPtr();
260  const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d[0].dataPtr();
261  const int inp_sound_size = inputSoundingData.size(0);
262  const bool anel_assume_dry = inputSoundingData.assume_dry;
263 
264  // Geometry
265  const Real* prob_lo = geomdata.ProbLo();
266  const Real* dx = geomdata.CellSize();
267  const Real z_lo = prob_lo[2];
268  const Real dz = dx[2];
269 
270  int kbot = geomdata.Domain().smallEnd(2);
271  int ktop = geomdata.Domain().bigEnd(2);
272 
273  // We want to set the lateral BC values, too
274  Box gbx = bx; // Copy constructor
275  gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
276 
277  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
278  const Real z = (z_cc_arr) ? z_cc_arr(i,j,k)
279  : z_lo + (k + 0.5) * dz;
280 
281  Real rho_k = interpolate_1d(z_inp_sound, rho_inp_sound, z, inp_sound_size);
282  Real rhoTh_k = rho_k * interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size);
283 
284  Real rho_k_base = rho_k;
285  if (l_isentropic) {
286  // `rho_inp_sound` previously calculated in calc_rho_p_isentropic()
287  // is in HSE and, when multiplied by the specified
288  // `theta_input_sound`, give p, T, and theta that are consistent
289  // with each other.
290  //
291  // Here, we do not require thermodynamic consistency between the
292  // initial base state rho and the prognostic variables. Instead,
293  // we calculate a `rho_hse` that is consistent with the isentropic
294  // (constant theta) assumption.
295  rho_k_base = rhoTh_k / theta_inp_sound[0];
296  }
297 
298  // Set the density
299  state(i, j, k, Rho_comp) = rho_k;
300 
301  // Initial Rho0*Theta0
302  state(i, j, k, RhoTheta_comp) = rhoTh_k;
303 
304  // Initialize all scalars to 0.
305  for (int n = 0; n < NSCALARS; n++) {
306  state(i, j, k, RhoScalar_comp+n) = 0;
307  }
308 
309  // total nonprecipitating water (Q1) == water vapor (Qv), i.e., there
310  // is no cloud water or cloud ice
311  Real qv_k = 0.0;
312  if (l_moist) {
313  qv_k = interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size);
314  state(i, j, k, RhoQ1_comp) = rho_k * qv_k;
315  }
316 
317  // Update hse quantities with values calculated from InputSoundingData.calc_rho_p()
318  if (anel_assume_dry) qv_k = 0;
319  r_hse_arr (i,j,k) = rho_k_base * (1.0 + qv_k);
320  p_hse_arr (i,j,k) = getPgivenRTh(rhoTh_k, qv_k);
321  pi_hse_arr(i,j,k) = getExnergivenRTh(rhoTh_k, l_rdOcp, qv_k);
322  th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k), qv_k) / rho_k_base;
323  qv_hse_arr(i,j,k) = qv_k;
324  if (l_isentropic) {
325 #if 0
326  if (i==0 && j==0) {
327  Print() << "HSE rho,p,T=pi*th,th at " << IntVect(i,j,k) << " : "
328  << r_hse_arr(i,j,k) << " "
329  << p_hse_arr(i,j,k) << " "
330  << pi_hse_arr(i,j,k)*th_hse_arr(i,j,k) << " "
331  << th_hse_arr(i,j,k)
332  << " with rho,rhotheta=" << rho_k << " " << rhoTh_k
333  << std::endl;
334  }
335 #endif
336  // If everything above is thermodynamically consistent, this should be constant
337  AMREX_ALWAYS_ASSERT(std::abs(th_hse_arr(i,j,k) - theta_inp_sound[0]) < 1e-12);
338  }
339 
340  // FOEXTRAP hse arrays
341  if (k==kbot)
342  {
343  for (int kk = 1; kk <= ngz; kk++) {
344  r_hse_arr(i, j, k-kk) = r_hse_arr(i,j,k);
345  p_hse_arr(i, j, k-kk) = p_hse_arr(i,j,k);
346  pi_hse_arr(i, j, k-kk) = pi_hse_arr(i,j,k);
347  th_hse_arr(i, j, k-kk) = th_hse_arr(i,j,k);
348  qv_hse_arr(i, j, k-kk) = qv_hse_arr(i,j,k);
349  }
350  }
351  else if (k==ktop)
352  {
353  for (int kk = 1; kk <= ngz; kk++) {
354  r_hse_arr(i, j, k+kk) = r_hse_arr(i,j,k);
355  p_hse_arr(i, j, k+kk) = p_hse_arr(i,j,k);
356  pi_hse_arr(i, j, k+kk) = pi_hse_arr(i,j,k);
357  th_hse_arr(i, j, k+kk) = th_hse_arr(i,j,k);
358  qv_hse_arr(i, j, k+kk) = qv_hse_arr(i,j,k);
359  }
360  }
361  });
362 }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=0.)
Definition: ERF_EOS.H:81
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:156
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:172
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)

Referenced by ERF::init_from_input_sounding().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ init_bx_velocities_from_input_sounding()

void init_bx_velocities_from_input_sounding ( const Box &  bx,
Array4< Real > const &  x_vel,
Array4< Real > const &  y_vel,
Array4< Real > const &  z_vel,
GeometryData const &  geomdata,
Array4< const Real > const &  z_nd_arr,
InputSoundingData const &  inputSoundingData 
)

Box level wrapper for initializing velocities from input sounding data.

Parameters
bxBox specifying the indices we are initializing
x_velArray4 specifying the x-velocity data we are to initialize
y_velArray4 specifying the y-velocity data we are to initialize
z_velArray4 specifying the z-velocity data we are to initialize
geomdataGeometryData object specifying the domain geometry
inputSoundingDataInputSoundingData object we are to initialize from
382 {
383  const Real* z_inp_sound = inputSoundingData.z_inp_sound_d[0].dataPtr();
384  const Real* U_inp_sound = inputSoundingData.U_inp_sound_d[0].dataPtr();
385  const Real* V_inp_sound = inputSoundingData.V_inp_sound_d[0].dataPtr();
386  const int inp_sound_size = inputSoundingData.size(0);
387 
388  // Geometry
389  const Real* prob_lo = geomdata.ProbLo();
390  const Real* dx = geomdata.CellSize();
391  const Real z_lo = prob_lo[2];
392  const Real dz = dx[2];
393 
394  // We want to set the lateral BC values, too
395  Box gbx = bx; // Copy constructor
396  gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
397 
398  // Construct a box that is on x-faces
399  const Box& xbx = surroundingNodes(gbx,0);
400  // Construct a box that is on y-faces
401  const Box& ybx = surroundingNodes(gbx,1);
402  // Construct a box that is on z-faces
403  const Box& zbx = surroundingNodes(gbx,2);
404 
405  // Set the x,y,z-velocities
407  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
408  // Note that this is called on a box of x-faces
409  const Real z = (z_nd_arr) ? 0.25*( z_nd_arr(i,j ,k )
410  + z_nd_arr(i,j+1,k )
411  + z_nd_arr(i,j ,k+1)
412  + z_nd_arr(i,j+1,k+1))
413  : z_lo + (k + 0.5) * dz;
414 
415  // Set the x-velocity
416  x_vel(i, j, k) = interpolate_1d(z_inp_sound, U_inp_sound, z, inp_sound_size);
417  },
418  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
419  // Note that this is called on a box of y-faces
420  const Real z = (z_nd_arr) ? 0.25*( z_nd_arr(i ,j,k )
421  + z_nd_arr(i+1,j,k )
422  + z_nd_arr(i ,j,k+1)
423  + z_nd_arr(i+1,j,k+1))
424  : z_lo + (k + 0.5) * dz;
425 
426  // Set the y-velocity
427  y_vel(i, j, k) = interpolate_1d(z_inp_sound, V_inp_sound, z, inp_sound_size);
428  },
429  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
430  // Note that this is called on a box of z-faces
431  // Set the z-velocity
432  z_vel(i, j, k) = 0.0;
433  });
434 }
const Box zbx
Definition: ERF_SetupDiff.H:9
const Box xbx
Definition: ERF_SetupDiff.H:7
const Box ybx
Definition: ERF_SetupDiff.H:8

Referenced by ERF::init_from_input_sounding().

Here is the call graph for this function:
Here is the caller graph for this function: