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
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
227 {
228  const Real* z_inp_sound = inputSoundingData.z_inp_sound_d[0].dataPtr();
229  const Real* theta_inp_sound = inputSoundingData.theta_inp_sound_d[0].dataPtr();
230  const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d[0].dataPtr();
231  const int inp_sound_size = inputSoundingData.size(0);
232 
233  // Geometry
234  const Real* prob_lo = geomdata.ProbLo();
235  const Real* dx = geomdata.CellSize();
236  const Real z_lo = prob_lo[2];
237  const Real dz = dx[2];
238 
239  // We want to set the lateral BC values, too
240  Box gbx = bx; // Copy constructor
241  gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
242 
243  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
244  const Real z = (z_cc_arr) ? z_cc_arr(i,j,k) : z_lo + (k + myhalf) * dz;
245 
246  Real rho_0 = one;
247 
248  // Set the density
249  state(i, j, k, Rho_comp) = rho_0;
250 
251  // Initial Rho0*Theta0
252  state(i, j, k, RhoTheta_comp) = rho_0 * interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size);
253 
254  // total nonprecipitating water (Q1) == water vapor (Qv), i.e., there is no cloud water or cloud ice
255  if (l_moist) {
256  state(i, j, k, RhoQ1_comp) = rho_0 * interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size);
257  }
258  });
259 }
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:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
Real rho_0
Definition: ERF_InitCustomPert_ABL.H:4
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
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_AdvanceWSM6.cpp:104

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 ,
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
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
297 {
298  const Real* z_inp_sound = inputSoundingData.z_inp_sound_d[0].dataPtr();
299  const Real* rho_inp_sound = inputSoundingData.rho_inp_sound_d.dataPtr();
300  const Real* theta_inp_sound = inputSoundingData.theta_inp_sound_d[0].dataPtr();
301  const Real* qv_inp_sound = inputSoundingData.qv_inp_sound_d[0].dataPtr();
302  const int inp_sound_size = inputSoundingData.size(0);
303  const bool anel_assume_dry = inputSoundingData.assume_dry;
304 
305  // Geometry
306  const Real* prob_lo = geomdata.ProbLo();
307  const Real* dx = geomdata.CellSize();
308  const Real z_lo = prob_lo[2];
309  const Real dz = dx[2];
310 
311  int kbot = geomdata.Domain().smallEnd(2);
312  int ktop = geomdata.Domain().bigEnd(2);
313 
314  // We want to set the lateral BC values, too
315  Box gbx = bx; // Copy constructor
316  gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
317 
318  ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
319  const Real z = (z_cc_arr) ? z_cc_arr(i,j,k)
320  : z_lo + (k + myhalf) * dz;
321 
322  Real rho_k = interpolate_1d(z_inp_sound, rho_inp_sound, z, inp_sound_size);
323  Real rhoTh_k = rho_k * interpolate_1d(z_inp_sound, theta_inp_sound, z, inp_sound_size);
324 
325  Real rho_k_base = rho_k;
326  if (l_isentropic) {
327  // `rho_inp_sound` previously calculated in calc_rho_p_isentropic()
328  // is in HSE and, when multiplied by the specified
329  // `theta_input_sound`, give p, T, and theta that are consistent
330  // with each other.
331  //
332  // Here, we do not require thermodynamic consistency between the
333  // initial base state rho and the prognostic variables. Instead,
334  // we calculate a `rho_hse` that is consistent with the isentropic
335  // (constant theta) assumption.
336  rho_k_base = rhoTh_k / theta_inp_sound[0];
337  }
338 
339  // Set the density
340  state(i, j, k, Rho_comp) = rho_k;
341 
342  // Initial Rho0*Theta0
343  state(i, j, k, RhoTheta_comp) = rhoTh_k;
344 
345  // Total nonprecipitating water (Q1) == water vapor (Qv), i.e., there
346  // is no cloud water or cloud ice
347  Real qv_k = zero;
348  if (l_moist) {
349  qv_k = interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size);
350  state(i, j, k, RhoQ1_comp) = rho_k * qv_k;
351  }
352 
353  // Update hse quantities with values calculated from InputSoundingData.calc_rho_p()
354  if (anel_assume_dry) qv_k = 0;
355  r_hse_arr (i,j,k) = rho_k_base;
356  p_hse_arr (i,j,k) = getPgivenRTh(rhoTh_k, qv_k);
357  pi_hse_arr(i,j,k) = getExnergivenRTh(rhoTh_k, l_rdOcp, qv_k);
358  th_hse_arr(i,j,k) = getRhoThetagivenP(p_hse_arr(i,j,k), qv_k) / rho_k_base;
359  qv_hse_arr(i,j,k) = qv_k;
360 
361  if (l_isentropic) {
362 #if 0
363  if (i==0 && j==0) {
364  Print() << "HSE rho,p,T=pi*th,th at " << IntVect(i,j,k) << " : "
365  << r_hse_arr(i,j,k) << " "
366  << p_hse_arr(i,j,k) << " "
367  << pi_hse_arr(i,j,k)*th_hse_arr(i,j,k) << " "
368  << th_hse_arr(i,j,k)
369  << " with rho,rhotheta=" << rho_k << " " << rhoTh_k
370  << std::endl;
371  }
372 #endif
373  // If everything above is thermodynamically consistent, this should be constant
374  AMREX_ALWAYS_ASSERT(std::abs(th_hse_arr(i,j,k) - theta_inp_sound[0]) < 1e-12);
375  }
376 
377  // FOEXTRAP hse arrays
378  if (k==kbot)
379  {
380  for (int kk = 1; kk <= ngz; kk++) {
381  r_hse_arr(i, j, k-kk) = r_hse_arr(i,j,k);
382  p_hse_arr(i, j, k-kk) = p_hse_arr(i,j,k);
383  pi_hse_arr(i, j, k-kk) = pi_hse_arr(i,j,k);
384  th_hse_arr(i, j, k-kk) = th_hse_arr(i,j,k);
385  qv_hse_arr(i, j, k-kk) = qv_hse_arr(i,j,k);
386  }
387  }
388  else if (k==ktop)
389  {
390  for (int kk = 1; kk <= ngz; kk++) {
391  r_hse_arr(i, j, k+kk) = r_hse_arr(i,j,k);
392  p_hse_arr(i, j, k+kk) = p_hse_arr(i,j,k);
393  pi_hse_arr(i, j, k+kk) = pi_hse_arr(i,j,k);
394  th_hse_arr(i, j, k+kk) = th_hse_arr(i,j,k);
395  qv_hse_arr(i, j, k+kk) = qv_hse_arr(i,j,k);
396  }
397  }
398  });
399 }
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
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
420 {
421  const Real* z_inp_sound = inputSoundingData.z_inp_sound_d[0].dataPtr();
422  const Real* U_inp_sound = inputSoundingData.U_inp_sound_d[0].dataPtr();
423  const Real* V_inp_sound = inputSoundingData.V_inp_sound_d[0].dataPtr();
424  const int inp_sound_size = inputSoundingData.size(0);
425 
426  // Geometry
427  const Real* prob_lo = geomdata.ProbLo();
428  const Real* dx = geomdata.CellSize();
429  const Real z_lo = prob_lo[2];
430  const Real dz = dx[2];
431 
432  // We want to set the lateral BC values, too
433  Box gbx = bx; // Copy constructor
434  gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions
435 
436  // Construct a box that is on x-faces
437  const Box& xbx = surroundingNodes(gbx,0);
438  // Construct a box that is on y-faces
439  const Box& ybx = surroundingNodes(gbx,1);
440  // Construct a box that is on z-faces
441  const Box& zbx = surroundingNodes(gbx,2);
442 
443  // Set the x,y,z-velocities
444  ParallelFor(xbx, ybx, zbx,
445  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
446  // Note that this is called on a box of x-faces
447  const Real z = (z_nd_arr) ? fourth*( z_nd_arr(i,j ,k )
448  + z_nd_arr(i,j+1,k )
449  + z_nd_arr(i,j ,k+1)
450  + z_nd_arr(i,j+1,k+1))
451  : z_lo + (k + myhalf) * dz;
452 
453  // Set the x-velocity
454  x_vel(i, j, k) = interpolate_1d(z_inp_sound, U_inp_sound, z, inp_sound_size);
455  },
456  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
457  // Note that this is called on a box of y-faces
458  const Real z = (z_nd_arr) ? fourth*( z_nd_arr(i ,j,k )
459  + z_nd_arr(i+1,j,k )
460  + z_nd_arr(i ,j,k+1)
461  + z_nd_arr(i+1,j,k+1))
462  : z_lo + (k + myhalf) * dz;
463 
464  // Set the y-velocity
465  y_vel(i, j, k) = interpolate_1d(z_inp_sound, V_inp_sound, z, inp_sound_size);
466  },
467  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
468  // Note that this is called on a box of z-faces
469  // Set the z-velocity
470  z_vel(i, j, k) = zero;
471  });
472 }
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: