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_state_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_state_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_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_state_from_input_sounding()

void init_state_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
[in]bxBox specifying the indices we are initializing
[out]stateArray4 specifying the state data we are to initialize
[in]geomdataGeometryData object specifying the domain geometry
[in]z_cc_arrArray4 specifying cell-centered terrain heights, if present
[in]l_moistWhether moisture variables should be initialized
[in]inputSoundingDataInputSoundingData object we are to initialize from

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
z_cc_arrArray4 specifying cell-centered terrain heights, if present
l_moistWhether moisture variables should be initialized
inputSoundingDataInputSoundingData object we are to initialize from
271 {
272  const Real* z_inp_sound = inputSoundingData.z_inp_sound_d[0].dataPtr();
273  const Real* theta_inp_sound = inputSoundingData.theta_inp_sound_d[0].dataPtr();
274  const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d[0].dataPtr();
275  const int inp_sound_size = inputSoundingData.size(0);
276 
277  // Geometry
278  const Real* prob_lo = geomdata.ProbLo();
279  const Real* dx = geomdata.CellSize();
280  const Real z_lo = prob_lo[2];
281  const Real dz = dx[2];
282 
283  // We want to set the lateral BC values, too
284  Box gbx = bx; // Copy constructor
285  gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
286 
287  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
288  const Real z = (z_cc_arr) ? z_cc_arr(i,j,k) : z_lo + (k + myhalf) * dz;
289 
290  Real rho_0 = one;
291 
292  // Set the density
293  state(i, j, k, Rho_comp) = rho_0;
294 
295  // Initial Rho0*Theta0
296  state(i, j, k, RhoTheta_comp) = rho_0 * interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size);
297 
298  // total nonprecipitating water (Q1) == water vapor (Qv), i.e., there is no cloud water or cloud ice
299  if (l_moist) {
300  state(i, j, k, RhoQ1_comp) = rho_0 * interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size);
301  }
302  });
303 }
constexpr amrex::Real one
Definition: ERF_Constants.H:9
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
#define Rho_comp
Definition: ERF_IndexDefines.H:39
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:40
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:45
Real rho_0
Definition: ERF_InitCustomPert_ABL.H:18
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
const GpuArray< Real, AMREX_SPACEDIM > prob_lo
Definition: ERF_InitCustomPert_CloudChamber.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:14
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
amrex::Real Real
Definition: ERF_ShocInterface.H:19
@ dz
Definition: ERF_AdvanceWDM6.cpp:270

Referenced by ERF::init_from_input_sounding().

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

◆ init_state_from_input_sounding_hse()

void init_state_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 
)

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

Parameters
[in]bxBox specifying the indices we are initializing
[out]stateArray4 specifying the state data we are to initialize
[out]r_hse_arrArray4 specifying the density HSE base state data we are to initialize
[out]p_hse_arrArray4 specifying the pressure HSE base state data we are to initialize
[out]pi_hse_arrArray4 specifying the Exner pressure HSE base state data we are to initialize
[out]th_hse_arrArray4 specifying the base state potential temperature we are to initialize
[out]qv_hse_arrArray4 specifying the base state water vapor mixing ratio we are to initialize
[in]geomdataGeometryData object specifying the domain geometry
[in]z_cc_arrArray4 specifying cell-centered terrain heights, if present
[in]l_gravityGravity constant (retained for interface compatibility).
[in]l_rdOcpReal number specifying the Rhydberg constant ($R_d$) divided by specific heat at constant pressure ($c_p$)
[in]l_moistWhether moisture variables should be initialized
[in]inputSoundingDataInputSoundingData object we are to initialize from
[in]l_isentropicWhether to construct an isentropic HSE base state
[in]ngzNumber of vertical ghost cells to fill in HSE arrays

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
qv_hse_arrArray4 specifying the base state water vapor mixing ratio we are to initialize
geomdataGeometryData object specifying the domain geometry The unused gravity argument is retained for interface compatibility.
z_cc_arrArray4 specifying cell-centered terrain heights, if present
l_rdOcpReal number specifying the Rhydberg constant ($R_d$) divided by specific heat at constant pressure ($c_p$)
l_moistWhether moisture variables should be initialized
inputSoundingDataInputSoundingData object we are to initialize from
l_isentropicWhether to construct an isentropic HSE base state
ngzNumber of vertical ghost cells to fill in HSE arrays
341 {
342  const Real* z_inp_sound = inputSoundingData.z_inp_sound_d[0].dataPtr();
343  const Real* rho_inp_sound = inputSoundingData.rho_inp_sound_d.dataPtr();
344  const Real* theta_inp_sound = inputSoundingData.theta_inp_sound_d[0].dataPtr();
345  const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d[0].dataPtr();
346  const int inp_sound_size = inputSoundingData.size(0);
347  const bool anel_assume_dry = inputSoundingData.assume_dry;
348 
349  // Geometry
350  const Real* prob_lo = geomdata.ProbLo();
351  const Real* dx = geomdata.CellSize();
352  const Real z_lo = prob_lo[2];
353  const Real dz = dx[2];
354 
355  int kbot = geomdata.Domain().smallEnd(2);
356  int ktop = geomdata.Domain().bigEnd(2);
357 
358  // We want to set the lateral BC values, too
359  Box gbx = bx; // Copy constructor
360  gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
361 
362  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
363  const Real z = (z_cc_arr) ? z_cc_arr(i,j,k)
364  : z_lo + (k + myhalf) * dz;
365 
366  Real rho_k = interpolate_1d(z_inp_sound, rho_inp_sound, z, inp_sound_size);
367  Real rhoTh_k = rho_k * interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size);
368 
369  Real rho_k_base = rho_k;
370  if (l_isentropic) {
371  // `rho_inp_sound` previously calculated in calc_rho_p_isentropic()
372  // is in HSE and, when multiplied by the specified
373  // `theta_input_sound`, give p, T, and theta that are consistent
374  // with each other.
375  //
376  // Here, we do not require thermodynamic consistency between the
377  // initial base state rho and the prognostic variables. Instead,
378  // we calculate a `rho_hse` that is consistent with the isentropic
379  // (constant theta) assumption.
380  rho_k_base = rhoTh_k / theta_inp_sound[0];
381  }
382 
383  // Set the density
384  state(i, j, k, Rho_comp) = rho_k;
385 
386  // Initial Rho0*Theta0
387  state(i, j, k, RhoTheta_comp) = rhoTh_k;
388 
389  // Total nonprecipitating water (Q1) == water vapor (Qv), i.e., there
390  // is no cloud water or cloud ice
391  Real qv_k = zero;
392  if (l_moist) {
393  qv_k = interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size);
394  state(i, j, k, RhoQ1_comp) = rho_k * qv_k;
395  }
396 
397  // Update hse quantities with values calculated from InputSoundingData.calc_rho_p()
398  if (anel_assume_dry) qv_k = 0;
399  r_hse_arr (i,j,k) = rho_k_base;
400  p_hse_arr (i,j,k) = getPgivenRTh(rhoTh_k, qv_k);
401  pi_hse_arr(i,j,k) = getExnergivenRTh(rhoTh_k, l_rdOcp, qv_k);
402  th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k), qv_k) / rho_k_base;
403  qv_hse_arr(i,j,k) = qv_k;
404 
405  if (l_isentropic) {
406 #if 0
407  if (i==0 && j==0) {
408  Print() << "HSE rho,p,T=pi*th,th at " << IntVect(i,j,k) << " : "
409  << r_hse_arr(i,j,k) << " "
410  << p_hse_arr(i,j,k) << " "
411  << pi_hse_arr(i,j,k)*th_hse_arr(i,j,k) << " "
412  << th_hse_arr(i,j,k)
413  << " with rho,rhotheta=" << rho_k << " " << rhoTh_k
414  << std::endl;
415  }
416 #endif
417  // If everything above is thermodynamically consistent, this should be constant
418  AMREX_ALWAYS_ASSERT(std::abs(th_hse_arr(i,j,k) - theta_inp_sound[0]) < 1e-12);
419  }
420 
421  // FOEXTRAP hse arrays
422  if (k==kbot)
423  {
424  for (int kk = 1; kk <= ngz; kk++) {
425  r_hse_arr(i, j, k-kk) = r_hse_arr(i,j,k);
426  p_hse_arr(i, j, k-kk) = p_hse_arr(i,j,k);
427  pi_hse_arr(i, j, k-kk) = pi_hse_arr(i,j,k);
428  th_hse_arr(i, j, k-kk) = th_hse_arr(i,j,k);
429  qv_hse_arr(i, j, k-kk) = qv_hse_arr(i,j,k);
430  }
431  }
432  else if (k==ktop)
433  {
434  for (int kk = 1; kk <= ngz; kk++) {
435  r_hse_arr(i, j, k+kk) = r_hse_arr(i,j,k);
436  p_hse_arr(i, j, k+kk) = p_hse_arr(i,j,k);
437  pi_hse_arr(i, j, k+kk) = pi_hse_arr(i,j,k);
438  th_hse_arr(i, j, k+kk) = th_hse_arr(i,j,k);
439  qv_hse_arr(i, j, k+kk) = qv_hse_arr(i,j,k);
440  }
441  }
442  });
443 }
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:172
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:156
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=amrex::Real(0))
Definition: ERF_EOS.H:81
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_velocities_from_input_sounding()

void init_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
[in]bxBox specifying the indices we are initializing
[out]x_velArray4 specifying the x-velocity data we are to initialize
[out]y_velArray4 specifying the y-velocity data we are to initialize
[out]z_velArray4 specifying the z-velocity data we are to initialize
[in]geomdataGeometryData object specifying the domain geometry
[in]z_nd_arrArray4 specifying node-centered terrain heights, if present
[in]inputSoundingDataInputSoundingData object we are to initialize from

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
z_nd_arrArray4 specifying node-centered terrain heights, if present
inputSoundingDataInputSoundingData object we are to initialize from
464 {
465  const Real* z_inp_sound = inputSoundingData.z_inp_sound_d[0].dataPtr();
466  const Real* U_inp_sound = inputSoundingData.U_inp_sound_d[0].dataPtr();
467  const Real* V_inp_sound = inputSoundingData.V_inp_sound_d[0].dataPtr();
468  const int inp_sound_size = inputSoundingData.size(0);
469 
470  // Geometry
471  const Real* prob_lo = geomdata.ProbLo();
472  const Real* dx = geomdata.CellSize();
473  const Real z_lo = prob_lo[2];
474  const Real dz = dx[2];
475 
476  // We want to set the lateral BC values, too
477  Box gbx = bx; // Copy constructor
478  gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
479 
480  // Construct a box that is on x-faces
481  const Box& xbx = surroundingNodes(gbx,0);
482  // Construct a box that is on y-faces
483  const Box& ybx = surroundingNodes(gbx,1);
484  // Construct a box that is on z-faces
485  const Box& zbx = surroundingNodes(gbx,2);
486 
487  // Set the x,y,z-velocities
488  ParallelFor(xbx, ybx, zbx,
489  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
490  // Note that this is called on a box of x-faces
491  const Real z = (z_nd_arr) ? fourth*( z_nd_arr(i,j ,k )
492  + z_nd_arr(i,j+1,k )
493  + z_nd_arr(i,j ,k+1)
494  + z_nd_arr(i,j+1,k+1))
495  : z_lo + (k + myhalf) * dz;
496 
497  // Set the x-velocity
498  x_vel(i, j, k) = interpolate_1d(z_inp_sound, U_inp_sound, z, inp_sound_size);
499  },
500  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
501  // Note that this is called on a box of y-faces
502  const Real z = (z_nd_arr) ? fourth*( z_nd_arr(i ,j,k )
503  + z_nd_arr(i+1,j,k )
504  + z_nd_arr(i ,j,k+1)
505  + z_nd_arr(i+1,j,k+1))
506  : z_lo + (k + myhalf) * dz;
507 
508  // Set the y-velocity
509  y_vel(i, j, k) = interpolate_1d(z_inp_sound, V_inp_sound, z, inp_sound_size);
510  },
511  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
512  // Note that this is called on a box of z-faces
513  // Set the z-velocity
514  z_vel(i, j, k) = zero;
515  });
516 }
constexpr amrex::Real fourth
Definition: ERF_Constants.H:14

Referenced by ERF::init_from_input_sounding().

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