ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_MetgridUtils.H File Reference
#include <ERF.H>
#include <ERF_EOS.H>
#include <ERF_Utils.H>
#include <ERF_ProbCommon.H>
#include <ERF_HSEUtils.H>
Include dependency graph for ERF_MetgridUtils.H:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

void read_from_metgrid (int lev, int itime, const amrex::Box &domain, const std::string &fname, std::string &NC_dateTime, amrex::Real &NC_epochTime, int &flag_psfc, int &flag_msf, int &flag_sst, int &flag_tsk, int &flag_lmask, int &NC_nx, int &NC_ny, amrex::Real &NC_dx, amrex::Real &NC_dy, amrex::FArrayBox &NC_xvel_fab, amrex::FArrayBox &NC_yvel_fab, amrex::FArrayBox &NC_temp_fab, amrex::FArrayBox &NC_rhum_fab, amrex::FArrayBox &NC_pres_fab, amrex::FArrayBox &NC_ght_fab, amrex::FArrayBox &NC_hgt_fab, amrex::FArrayBox &NC_psfc_fab, amrex::FArrayBox &NC_msfu_fab, amrex::FArrayBox &NC_msfv_fab, amrex::FArrayBox &NC_msfm_fab, amrex::FArrayBox &NC_sst_fab, amrex::FArrayBox &NC_tsk_fab, amrex::FArrayBox &NC_LAT_fab, amrex::FArrayBox &NC_LON_fab, amrex::IArrayBox &NC_lmask_iab, amrex::Geometry &geom)
 
void init_terrain_from_metgrid (amrex::FArrayBox &z_phys_nd_fab, amrex::FArrayBox &NC_hgt_fab)
 
void init_state_from_metgrid (const int lev, const int itime, const bool use_moisture, const bool interp_theta, const bool metgrid_debug_quiescent, const bool metgrid_debug_isothermal, const bool metgrid_debug_dry, const bool metgrid_basic_linear, const bool metgrid_use_below_sfc, const bool metgrid_use_sfc, const bool metgrid_retain_sfc, const amrex::Real metgrid_proximity, const int metgrid_order, const int metgrid_metgrid_force_sfc_k, const amrex::Real l_rdOcp, amrex::Box &tbxc, amrex::Box &tbxu, amrex::Box &tbxv, amrex::FArrayBox &state_fab, amrex::FArrayBox &x_vel_fab, amrex::FArrayBox &y_vel_fab, amrex::FArrayBox &z_vel_fab, amrex::FArrayBox &z_phys_nd_fab, const amrex::FArrayBox &NC_ght_fab, const amrex::FArrayBox &NC_xvel_fab, const amrex::FArrayBox &NC_yvel_fab, const amrex::FArrayBox &NC_temp_fab, const amrex::FArrayBox &NC_rhum_fab, const amrex::FArrayBox &NC_pres_fab, amrex::FArrayBox &tmp_src_fab, amrex::FArrayBox &tmp_dst_fab, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_xlo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_xhi, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_ylo, amrex::Vector< amrex::Vector< amrex::FArrayBox >> &fabs_for_bcs_yhi, const amrex::Array4< const int > &mask_c_arr, const amrex::Array4< const int > &mask_u_arr, const amrex::Array4< const int > &mask_v_arr)
 
void init_msfs_from_metgrid (const bool metgrid_debug_msf, amrex::FArrayBox &msfu_fab, amrex::FArrayBox &msfv_fab, amrex::FArrayBox &msfm_fab, const int &flag_msf, amrex::FArrayBox &NC_MSFU_fab, amrex::FArrayBox &NC_MSFV_fab, amrex::FArrayBox &NC_MSFM_fab)
 
void init_base_state_from_metgrid (const bool use_moisture, const bool metgrid_debug_psfc, const amrex::Real l_rdOcp, const amrex::Box &valid_bx, const int &flag_psfc, amrex::FArrayBox &state_fab, amrex::FArrayBox &r_hse_fab, amrex::FArrayBox &p_hse_fab, amrex::FArrayBox &pi_hse_fab, amrex::FArrayBox &th_hse_fab, amrex::FArrayBox &qv_hse_fab, amrex::FArrayBox &z_phys_nd_fab, amrex::FArrayBox &z_phys_cc_fab, const amrex::FArrayBox &NC_psfc_fab)
 
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_interp (const int &order, amrex::Real *x, amrex::Real *y, amrex::Real &new_x, amrex::Real &new_y)
 
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_setup (char var_type, const bool &exp_interp, const int &orig_n, const int &new_n, const int &order, const int &i, const int &j, amrex::Real *orig_x_z, amrex::Real *orig_x_p, amrex::Real *orig_y, amrex::Real *new_x_z, amrex::Real *new_x_p, amrex::Real *new_y)
 
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal (const amrex::Real &z, amrex::Real &p)
 
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void interpolate_column_metgrid (const bool &metgrid_use_below_sfc, const bool &metgrid_use_sfc, const bool &exp_interp, const bool &metgrid_retain_sfc, const amrex::Real &metgrid_proximity, const int &metgrid_order, const int &metgrid_force_sfc_k, const int &i, const int &j, const int &kmax, const int &src_comp, const int &dst_comp, char var_type, char stag, const amrex::Array4< amrex::Real const > &orig_z_full, const amrex::Array4< amrex::Real const > &orig_data, const amrex::Array4< amrex::Real const > &new_z_full, const amrex::Array4< amrex::Real > &new_data_full)
 
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real interpolate_column_metgrid_linear (const int &i, const int &j, const int &k, char stag, int src_comp, const amrex::Array4< amrex::Real const > &orig_z, const amrex::Array4< amrex::Real const > &orig_data, const amrex::Array4< amrex::Real const > &new_z)
 
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void rh_to_mxrat (int i, int j, int k, const amrex::Array4< amrex::Real const > &rhum, const amrex::Array4< amrex::Real const > &temp, const amrex::Array4< amrex::Real const > &pres, int src_indx, const amrex::Array4< amrex::Real > &mxrat)
 

Function Documentation

◆ calc_p_isothermal()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal ( const amrex::Real &  z,
amrex::Real &  p 
)
453 {
454  p = p_0*std::exp(-CONST_GRAV*z/(amrex::Real(290.0)*R_d));
455 }
p_0
constexpr amrex::Real p_0
Definition: ERF_Constants.H:29
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:32
R_d
constexpr amrex::Real R_d
Definition: ERF_Constants.H:21
p
Real * p
Definition: ERF_InitCustomPert_SquallLine.H:61
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19

Referenced by interpolate_column_metgrid().

Here is the caller graph for this function:

◆ init_base_state_from_metgrid()

void init_base_state_from_metgrid ( const bool  use_moisture,
const bool  metgrid_debug_psfc,
const amrex::Real  l_rdOcp,
const amrex::Box &  valid_bx,
const int &  flag_psfc,
amrex::FArrayBox &  state_fab,
amrex::FArrayBox &  r_hse_fab,
amrex::FArrayBox &  p_hse_fab,
amrex::FArrayBox &  pi_hse_fab,
amrex::FArrayBox &  th_hse_fab,
amrex::FArrayBox &  qv_hse_fab,
amrex::FArrayBox &  z_phys_nd_fab,
amrex::FArrayBox &  z_phys_cc_fab,
const amrex::FArrayBox &  NC_psfc_fab 
)

◆ init_msfs_from_metgrid()

void init_msfs_from_metgrid ( const bool  metgrid_debug_msf,
amrex::FArrayBox &  msfu_fab,
amrex::FArrayBox &  msfv_fab,
amrex::FArrayBox &  msfm_fab,
const int &  flag_msf,
amrex::FArrayBox &  NC_MSFU_fab,
amrex::FArrayBox &  NC_MSFV_fab,
amrex::FArrayBox &  NC_MSFM_fab 
)

◆ init_state_from_metgrid()

void init_state_from_metgrid ( const int  lev,
const int  itime,
const bool  use_moisture,
const bool  interp_theta,
const bool  metgrid_debug_quiescent,
const bool  metgrid_debug_isothermal,
const bool  metgrid_debug_dry,
const bool  metgrid_basic_linear,
const bool  metgrid_use_below_sfc,
const bool  metgrid_use_sfc,
const bool  metgrid_retain_sfc,
const amrex::Real  metgrid_proximity,
const int  metgrid_order,
const int  metgrid_metgrid_force_sfc_k,
const amrex::Real  l_rdOcp,
amrex::Box &  tbxc,
amrex::Box &  tbxu,
amrex::Box &  tbxv,
amrex::FArrayBox &  state_fab,
amrex::FArrayBox &  x_vel_fab,
amrex::FArrayBox &  y_vel_fab,
amrex::FArrayBox &  z_vel_fab,
amrex::FArrayBox &  z_phys_nd_fab,
const amrex::FArrayBox &  NC_ght_fab,
const amrex::FArrayBox &  NC_xvel_fab,
const amrex::FArrayBox &  NC_yvel_fab,
const amrex::FArrayBox &  NC_temp_fab,
const amrex::FArrayBox &  NC_rhum_fab,
const amrex::FArrayBox &  NC_pres_fab,
amrex::FArrayBox &  tmp_src_fab,
amrex::FArrayBox &  tmp_dst_fab,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs_xlo,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs_xhi,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs_ylo,
amrex::Vector< amrex::Vector< amrex::FArrayBox >> &  fabs_for_bcs_yhi,
const amrex::Array4< const int > &  mask_c_arr,
const amrex::Array4< const int > &  mask_u_arr,
const amrex::Array4< const int > &  mask_v_arr 
)

◆ init_terrain_from_metgrid()

void init_terrain_from_metgrid ( amrex::FArrayBox &  z_phys_nd_fab,
amrex::FArrayBox &  NC_hgt_fab 
)

◆ interpolate_column_metgrid()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE void interpolate_column_metgrid ( const bool &  metgrid_use_below_sfc,
const bool &  metgrid_use_sfc,
const bool &  exp_interp,
const bool &  metgrid_retain_sfc,
const amrex::Real &  metgrid_proximity,
const int &  metgrid_order,
const int &  metgrid_force_sfc_k,
const int &  i,
const int &  j,
const int &  kmax,
const int &  src_comp,
const int &  dst_comp,
char  var_type,
char  stag,
const amrex::Array4< amrex::Real const > &  orig_z_full,
const amrex::Array4< amrex::Real const > &  orig_data,
const amrex::Array4< amrex::Real const > &  new_z_full,
const amrex::Array4< amrex::Real > &  new_data_full 
)
478 {
479  // Here we closely follow WRF's vert_interp from
480  // dyn_em/module_initialize_real.F, although changes have been
481  // made to accommodate interpolation relative to height instead of
482  // pressure.
483  int imin_orig = amrex::lbound(amrex::Box(orig_data)).x;
484  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
485  int jmin_orig = amrex::lbound(amrex::Box(orig_data)).y;
486  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
487  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
488  int kmax_new = kmax + 1;
489 
490  AMREX_ASSERT(kmax_orig < 256);
491  AMREX_ASSERT(kmax_new < 256);
492 
493  amrex::GpuArray<amrex::Real,256> new_z;
494  amrex::GpuArray<amrex::Real,256> new_p;
495  amrex::GpuArray<amrex::Real,256> new_data;
496  amrex::Real* new_z_p = new_z.data();
497  amrex::Real* new_p_p = new_p.data();
498  amrex::Real* new_data_p = new_data.data();
499  for (int k=0; k < kmax_new; k++) {
500  if (stag == 'X') {
501  new_z_p[k] = fourth*(new_z_full(i,j,k)+new_z_full(i,j+1,k)+new_z_full(i,j,k+1)+new_z_full(i,j+1,k+1));
502  } else if (stag == 'Y') {
503  new_z_p[k] = fourth*(new_z_full(i,j,k)+new_z_full(i+1,j,k)+new_z_full(i,j,k+1)+new_z_full(i+1,j,k+1));
504  } else if (stag == 'M') {
505  new_z_p[k] = amrex::Real(0.125)*(new_z_full(i,j,k )+new_z_full(i,j+1,k )+new_z_full(i+1,j,k )+new_z_full(i+1,j+1,k )+
506  new_z_full(i,j,k+1)+new_z_full(i,j+1,k+1)+new_z_full(i+1,j,k+1)+new_z_full(i+1,j+1,k+1));
507  }
508  calc_p_isothermal(new_z_p[k], new_p_p[k]);
509  }
510 
511  amrex::GpuArray<amrex::Real,256> orig_z;
512  amrex::Real* orig_z_p = orig_z.data();
513  for (int k=0; k < kmax_orig; k++) {
514  if (stag == 'M') {
515  orig_z_p[k] = orig_z_full(i,j,k);
516  } else if (stag == 'X') {
517  if (i <= imin_orig) {
518  orig_z_p[k] = orig_z_full(i,j,k);
519  } else if (i >= imax_orig) {
520  orig_z_p[k] = orig_z_full(imax_orig-1,j,k);
521  } else {
522  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k));
523  }
524  } else if (stag == 'Y') {
525  if (j <= jmin_orig) {
526  orig_z_p[k] = orig_z_full(i,j,k);
527  } else if (j >= jmax_orig) {
528  orig_z_p[k] = orig_z_full(i,jmax_orig-1,k);
529  } else {
530  orig_z_p[k] = myhalf*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k));
531  }
532  }
533  }
534 
535  // Check if the data is top-down instead of bottom-up.
536  bool flip_data_required = false;
537  if (orig_z[1] > orig_z[kmax_orig-1]) flip_data_required = true;
538  if (flip_data_required) amrex::Abort("metgrid initialization flip_data_required. Not yet implemented.");
539 
540  // Search for the first level above the surface in the origin data.
541  // This is needed since the origin model topography will be
542  // different than the topography processed by WPS.
543  int k_above_sfc = 0;
544  for (int k=1; k < kmax_orig; k++) {
545  if (orig_z_p[k] > orig_z_p[0]) {
546  k_above_sfc = k;
547  break;
548  }
549  }
550 
551  int kend_order;
552  amrex::GpuArray<amrex::Real,256> ordered_z;
553  amrex::GpuArray<amrex::Real,256> ordered_data;
554  amrex::Real* ordered_z_p = ordered_z.data();
555  amrex::Real* ordered_data_p = ordered_data.data();
556  if (k_above_sfc > 1) {
557  // The levels are not monotonically increasing in height, so
558  // we sort and then make "artistic" quality control choices.
559  int count = 0;
560 
561  // Insert levels that are below the surface.
562  for (int k=1; k < k_above_sfc; k++) {
563  ordered_z_p[count] = orig_z_p[k];
564  ordered_data_p[count] = orig_data(i,j,k,src_comp);
565  count++;
566  }
567 
568  // Check if the level that is nearest to and below the surface
569  // is "too close". If so, we'll ignore the upper level and keep
570  // the lower. Origin data is likely to be on pressure levels
571  // with higher spatial resolution near-surface, which supports
572  // the choice of eliminating levels that are "too close" in
573  // pressure-space. For simplicity, calculate delta P assuming a
574  // baroclinic atmosphere.
575  amrex::Real Pu, Pl;
576  calc_p_isothermal(orig_z_p[0], Pu);
577  calc_p_isothermal(ordered_z_p[count-1], Pl);
578  if (Pl-Pu < metgrid_proximity) {
579  count--;
580  }
581 
582  // Insert the surface level.
583  ordered_z_p[count] = orig_z_p[0];
584  ordered_data_p[count] = orig_data(i,j,0,src_comp);
585  count++;
586 
587  // Quoting WRF's comments, the next level to use is at,
588  // "... ta da, the first level above the surface. I know, wow."
589  int knext = k_above_sfc;
590  // Conditionally more strongly use the surface data by removing
591  // levels between the surface and the height corresponding to a
592  // set number of ERF levels from the surface. This forces the
593  // interpolator to use the surface data up through a number of
594  // ERF levels from the surface.
595  if (metgrid_force_sfc_k > 0) {
596  for (int k=k_above_sfc; k < kmax_orig; k++) {
597  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k-1]) {
598  knext = k;
599  break;
600  }
601  }
602  }
603 
604  // Check if the level that is nearest to and above the surface
605  // is "too close". If so, we'll ignore that level.
606  calc_p_isothermal(orig_z_p[knext], Pu);
607  calc_p_isothermal(ordered_z_p[count-1], Pl);
608  if (Pl-Pu < metgrid_proximity) {
609  knext++;
610  }
611 
612  // Insert levels that are above the surface.
613  for (int k=knext; k < kmax_orig; k++) {
614  ordered_z_p[count] = orig_z_p[k];
615  ordered_data_p[count] = orig_data(i,j,k,src_comp);
616  count++;
617  }
618 
619  kend_order = count;
620  } else {
621  // The surface is the lowest level in the origin data.
622 
623  // Insert the surface.
624  ordered_z_p[0] = orig_z[0];
625  ordered_data_p[0] = orig_data(i,j,0,src_comp);
626 
627  // Similar to above, conditionally more strongly use the
628  // surface data.
629  int count = 1;
630  int knext = count;
631  if (metgrid_force_sfc_k > 0) {
632  for (int k=knext; k < kmax_orig; k++) {
633  if (orig_z_p[k] > new_z_p[metgrid_force_sfc_k]) {
634  knext = k;
635  break;
636  }
637  }
638  }
639 
640  // Insert the remaining levels, again ignoring levels that are
641  // "too close" to the prior valid level.
642  for (int k=knext; k < kmax_orig; k++) {
643  amrex::Real Pu, Pl;
644  calc_p_isothermal(orig_z_p[k], Pu);
645  calc_p_isothermal(ordered_z_p[count-1], Pl);
646  if (Pl-Pu < metgrid_proximity) {
647  continue;
648  }
649  ordered_z_p[count] = orig_z_p[k];
650  ordered_data_p[count] = orig_data(i,j,k,src_comp);
651  count++;
652  }
653  kend_order = count;
654  }
655 
656  int ksta(0), kend(0);
657  if (metgrid_use_below_sfc && metgrid_use_sfc) {
658  // Use all levels.
659  ksta = 0;
660  kend = kend_order-1;
661  } else if (metgrid_use_below_sfc && !metgrid_use_sfc) {
662  // Use all levels except for the surface.
663  int ksfc = 0;
664  for (int k=0; k < kmax_orig; k++) {
665  if (ordered_z_p[k] == orig_z_p[0]) {
666  ksfc = k;
667  break;
668  }
669  }
670  for (int k=ksfc; k < kmax_orig-1; k++) {
671  ordered_z_p[k] = ordered_z_p[k+1];
672  ordered_data_p[k] = ordered_data_p[k+1];
673  }
674  ksta = 0;
675  kend = kend_order-2;
676  } else if (!metgrid_use_below_sfc && metgrid_use_sfc) {
677  // Use all levels above and including the surface.
678  int ksfc = 0;
679  for (int k=0; k < kend_order; k++) {
680  if (ordered_z_p[k] == orig_z_p[0]) {
681  ksfc = k;
682  break;
683  }
684  }
685  ksta = ksfc;
686  kend = kend_order-1;
687  } else {
688  // We shouldn't be in here!
689  amrex::Abort("metgrid initialization, !use_levels below_ground && !metgrid_use_sfc");
690  }
691 
692  // Insert the level of maximum winds.
693 // amrex::Real maxw_above_this_level = amrex::Real(30000.0);
694 // amrex::Real maxw_horiz_pres_diff = amrex::Real(5000.0);
695 // if ((flag_maxw == 1) && (use_maxw)) {
696 // amrex::Abort("metgrid initialization, use_maxw not yet implemented");
697 // }
698 
699  // Insert the level of the tropopause.
700 // amrex::Real trop_horiz_pres_diff = amrex::Real(5000.0);
701 // if ((flag_trop == 1) && (use_trop)) {
702 // amrex::Abort("metgrid initialization, use_trop not yet implemented");
703 // }
704 
705  amrex::GpuArray<amrex::Real,256> ordered_p;
706  amrex::Real* ordered_p_p = ordered_p.data();
707  for (int k=0; k < kend_order; k++) {
708  calc_p_isothermal(ordered_z_p[k], ordered_p_p[k]);
709  }
710 
711  int new_n = 0;
712  int zap_final = 0;
713  amrex::GpuArray<amrex::Real,256> final_z;
714  amrex::GpuArray<amrex::Real,256> final_p;
715  amrex::GpuArray<amrex::Real,256> final_data;
716  amrex::Real* final_z_p = final_z.data();
717  amrex::Real* final_p_p = final_p.data();
718  amrex::Real* final_data_p = final_data.data();
719  final_z_p[0] = ordered_z[ksta];
720  final_p_p[0] = ordered_p[ksta];
721  final_data_p[0] = ordered_data[ksta];
722  for (int k=ksta+1; k <= kend; k++) {
723  if ((final_p_p[new_n]-ordered_p_p[k]) < metgrid_proximity) {
724  zap_final++;
725  } else {
726  new_n++;
727  final_z_p[new_n] = ordered_z_p[k];
728  final_p_p[new_n] = ordered_p_p[k];
729  final_data_p[new_n] = ordered_data_p[k];
730  }
731  }
732  kend -= zap_final;
733 
734  // Call the interpolator.
735  lagrange_setup(var_type,
736  exp_interp,
737  kend-ksta,
738  kmax_new,
739  metgrid_order,
740  i,
741  j,
742  final_z_p,
743  final_p_p,
744  final_data_p,
745  new_z_p,
746  new_p_p,
747  new_data_p);
748 
749  // Optionally replace the lowest level of data with the surface
750  // field from the origin data.
751  if (metgrid_retain_sfc) new_data[0] = ordered_data[0];
752 
753  // Save the interpolated data.
754  for (int k=0; k < kmax_new; k++) {
755  new_data_full(i,j,k,dst_comp) = new_data[k];
756  }
757 
758 }
fourth
constexpr amrex::Real fourth
Definition: ERF_Constants.H:12
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:11
calc_p_isothermal
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void calc_p_isothermal(const amrex::Real &z, amrex::Real &p)
Definition: ERF_MetgridUtils.H:451
lagrange_setup
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_setup(char var_type, const bool &exp_interp, const int &orig_n, const int &new_n, const int &order, const int &i, const int &j, amrex::Real *orig_x_z, amrex::Real *orig_x_p, amrex::Real *orig_y, amrex::Real *new_x_z, amrex::Real *new_x_p, amrex::Real *new_y)
Definition: ERF_MetgridUtils.H:147
Here is the call graph for this function:

◆ interpolate_column_metgrid_linear()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real interpolate_column_metgrid_linear ( const int &  i,
const int &  j,
const int &  k,
char  stag,
int  src_comp,
const amrex::Array4< amrex::Real const > &  orig_z,
const amrex::Array4< amrex::Real const > &  orig_data,
const amrex::Array4< amrex::Real const > &  new_z 
)
771 {
772  // This subroutine is a bit ham-handed and can be cleaned up later.
773  int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;
774  int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y;
775  int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z;
776 
777  amrex::Real z;
778  if (stag == 'X') {
779  z = fourth*(new_z(i,j,k)+new_z(i,j+1,k)+new_z(i,j,k+1)+new_z(i,j+1,k+1));
780  }
781  else if (stag == 'Y') {
782  z = fourth*(new_z(i,j,k)+new_z(i+1,j,k)+new_z(i,j,k+1)+new_z(i+1,j,k+1));
783  }
784  else if (stag == 'M') {
785  z = amrex::Real(0.125)*(new_z(i,j,k )+new_z(i,j+1,k )+new_z(i+1,j,k )+new_z(i+1,j+1,k )+
786  new_z(i,j,k+1)+new_z(i,j+1,k+1)+new_z(i+1,j,k+1)+new_z(i+1,j+1,k+1));
787  }
788 
789  amrex::Real z0, z1;
790  int klow = -1;
791  int khi0 = -1;
792  amrex::Real dzlow = amrex::Real(1.0e12);
793  amrex::Real dzhi0 = -amrex::Real(1.0e12);
794  for (int kk = 0; kk < kmax_orig; kk++) {
795  amrex::Real orig_z_stag = zero;
796  if (stag == 'M') {
797  orig_z_stag = orig_z(i,j,kk);
798  }
799  if (stag == 'X') {
800  if (i == 0) {
801  orig_z_stag = orig_z(i,j,kk);
802  }
803  else if (i == imax_orig) {
804  orig_z_stag = orig_z(imax_orig-1,j,kk);
805  }
806  else {
807  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
808  }
809  }
810  else if (stag == 'Y') {
811  if (j == 0) {
812  orig_z_stag = orig_z(i,j,kk);
813  }
814  else if (j == jmax_orig) {
815  orig_z_stag = orig_z(i,jmax_orig-1,kk);
816  }
817  else {
818  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
819  }
820  }
821 
822  amrex::Real dz = z - orig_z_stag;
823  if ((dz < zero) && (dz > dzhi0)) {
824  dzhi0 = dz;
825  khi0 = kk;
826  z1 = orig_z_stag;
827  }
828  if ((dz >= zero) && (dz < dzlow)) {
829  dzlow = dz;
830  klow = kk;
831  z0 = orig_z_stag;
832  }
833  } // kk
834 
835  // extrapolate below the bottom surface
836  if (klow == -1) {
837  int khi1 = -1;
838  amrex::Real dzhi1 = -amrex::Real(1.0e12);
839  for (int kk = 0; kk < kmax_orig; kk++) {
840  amrex::Real orig_z_stag = zero;
841  if (stag == 'M') {
842  orig_z_stag = orig_z(i,j,kk);
843  }
844  else if (stag == 'X') {
845  if (i == 0) {
846  orig_z_stag = orig_z(i,j,kk);
847  }
848  else if (i == imax_orig) {
849  orig_z_stag = orig_z(imax_orig-1,j,kk);
850  }
851  else {
852  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i-1,j,kk));
853  }
854  }
855  else if (stag == 'Y') {
856  if (j == 0) {
857  orig_z_stag = orig_z(i,j,kk);
858  }
859  else if (j == jmax_orig) {
860  orig_z_stag = orig_z(i,jmax_orig-1,kk);
861  }
862  else {
863  orig_z_stag = myhalf*(orig_z(i,j,kk)+orig_z(i,j-1,kk));
864  }
865  }
866  amrex::Real dz = z - orig_z_stag;
867  if ((dz < zero) && (dz > dzhi1) && (kk != khi0)) {
868  dzhi1 = dz;
869  khi1 = kk;
870  z1 = orig_z_stag;
871  }
872  }
873  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
874  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
875  return ( y0-(y1-y0)/(z1-z0)*(z0-z) );
876 
877  // Extrapolate above the top surface
878  } else if (khi0 == -1) {
879  khi0 = klow - 1;
880  int khi1 = klow;
881  if (stag == 'M') {
882  z0 = orig_z(i,j,khi0);
883  }
884  else if (stag == 'X') {
885  if (i == 0) {
886  z0 = orig_z(i,j,khi0);
887  }
888  else if (i == imax_orig) {
889  z0 = orig_z(imax_orig-1,j,khi0);
890  }
891  else {
892  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i-1,j,khi0));
893  }
894  }
895  else if (stag == 'Y') {
896  if (j == 0) {
897  z0 = orig_z(i,j,khi0);
898  }
899  else if (j == jmax_orig) {
900  z0 = orig_z(i,jmax_orig-1,khi0);
901  }
902  else {
903  z0 = myhalf*(orig_z(i,j,khi0)+orig_z(i,j-1,khi0));
904  }
905  }
906  amrex::Real y0 = orig_data(i,j,khi0,src_comp);
907  amrex::Real y1 = orig_data(i,j,khi1,src_comp);
908  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
909  } else {
910  // interpolate
911  amrex::Real y0 = orig_data(i,j,klow,src_comp);
912  amrex::Real y1 = orig_data(i,j,khi0,src_comp);
913  return ( y0+(y1-y0)/(z1-z0)*(z-z0) );
914 
915  }
916 }
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:6
Coord::z
@ z
z0
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
dz
const Real dz
Definition: ERF_InitCustomPert_Bubble.H:25

◆ lagrange_interp()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_interp ( const int &  order,
amrex::Real *  x,
amrex::Real *  y,
amrex::Real &  new_x,
amrex::Real &  new_y 
)
122 {
123  // Interpolation using Lagrange polynomials.
124  // P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x)
125  // where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn)
126  // ---------------------------------------------
127  // (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn)
128  amrex::Real Px = zero;
129  for (int i=0; i <= order; i++) {
130  amrex::Real n = one;
131  amrex::Real d = one;
132  for (int k=0; k <= order; k++) {
133  if (k == i) continue;
134  n *= new_x-x[k];
135  d *= x[i]-x[k];
136  }
137  if (d != zero) {
138  Px += y[i]*n/d;
139  }
140  }
141  new_y = Px;
142 }
one
constexpr amrex::Real one
Definition: ERF_Constants.H:7
Coord::y
@ y
Coord::x
@ x

Referenced by lagrange_setup().

Here is the caller graph for this function:

◆ lagrange_setup()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_setup ( char  var_type,
const bool &  exp_interp,
const int &  orig_n,
const int &  new_n,
const int &  order,
const int &  i,
const int &  j,
amrex::Real *  orig_x_z,
amrex::Real *  orig_x_p,
amrex::Real *  orig_y,
amrex::Real *  new_x_z,
amrex::Real *  new_x_p,
amrex::Real *  new_y 
)
160 {
161 
162  amrex::Real CRC_const1 = amrex::Real(11880.516); // m
163  amrex::Real CRC_const2 = amrex::Real(0.1902632);
164  amrex::Real CRC_const3 = amrex::Real(0.0065); // K km-1
165 
166  amrex::ignore_unused(i,j);
167 #ifndef AMREX_USE_GPU
168  bool debug = false;
169 #endif
170 
171  for (int new_k=0; new_k < new_n; new_k++) {
172 #ifndef AMREX_USE_GPU
173  if (debug) amrex::Print() << "new_k=" << new_k;
174 #endif
175  // Determine if interpolating or extrapolating.
176  // If interpolating, find bounding x values and store the indices.
177  bool extrapolating = true;
178  int kl, kr;
179  for (int ko=0; ko < orig_n-1; ko++) {
180  amrex::Real a = new_x_z[new_k]-orig_x_z[ko];
181  amrex::Real b = new_x_z[new_k]-orig_x_z[ko+1];
182  if (a*b <= zero) {
183  kl = ko;
184  kr = ko+1;
185  extrapolating = false;
186  break;
187  }
188  }
189 
190  if (extrapolating) {
191  if (var_type == 'T') {
192  // Assume a standard atmosphere -amrex::Real(6.5) K km-1 lapse rate.
193  // Comparable to the WRF default, t_extrap_type=two
194  amrex::Real depth_of_extrap_in_p = new_x_p[new_k]-orig_x_p[0];
195  amrex::Real avg_of_extrap_p = myhalf*(new_x_p[new_k]+orig_x_p[0]);
196  amrex::Real temp_extrap_starting_point = orig_y[0]*std::pow(orig_x_p[0]/amrex::Real(100000.0), R_d/Cp_d);
197  amrex::Real dZdP = CRC_const1*CRC_const2*std::pow(avg_of_extrap_p/amrex::Real(100.0), CRC_const2-one);
198  amrex::Real dZ = dZdP*(depth_of_extrap_in_p/amrex::Real(100.0));
199  amrex::Real dT = dZ*CRC_const3;
200  new_y[new_k] = (temp_extrap_starting_point+dT)*std::pow(amrex::Real(100000.0)/new_x_p[new_k], R_d/Cp_d);
201  } else {
202  // Use a constant value below ground.
203  // Comparable to the WRF default, extrap_type=two
204  new_y[new_k] = orig_y[0];
205  }
206  continue;
207  }
208 
209  if (order%2 != 0) {
210  if ((kl-((order+1)/2-1) >= 0) && (kr+((order+1)/2-1) <= orig_n-1)) {
211  // Odd order interpolation.
212  int ksta = kl-(((order+1)/2)-1);
213  int kend = ksta+order;
214 #ifndef AMREX_USE_GPU
215  int ksize = kend-ksta;
216  if (debug) amrex::Print() << " (1a) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
217 #endif
218  amrex::Real new_x;
219  amrex::GpuArray<amrex::Real,9> orig_x_sub;
220  amrex::Real* orig_x_sub_p = orig_x_sub.data();
221  if (exp_interp) {
222  new_x = new_x_p[new_k];
223  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
224  } else {
225  new_x = new_x_z[new_k];
226  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
227  }
228  amrex::GpuArray<amrex::Real,9> orig_y_sub;
229  amrex::Real* orig_y_sub_p = orig_y_sub.data();
230  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
231 #ifndef AMREX_USE_GPU
232  if (debug) {
233  amrex::Print() << " orig_x_sub_p = [";
234  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
235  amrex::Print() << "]" << std::endl;
236  amrex::Print() << " orig_y_sub_p = [";
237  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
238  amrex::Print() << "]" << std::endl;
239  }
240 #endif
241  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
242  } else {
243  // Interpolation stencil is too big due to proximity to the surface or model top.
244  // We have bounding points so resort to a simple linear interpolation.
245  int ksta = kl;
246  int kend = kr;
247 #ifndef AMREX_USE_GPU
248  int ksize = kend-ksta+1;
249  if (debug) amrex::Print() << " (1b) order=" << order << " new_x=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
250 #endif
251  amrex::Real new_x;
252  amrex::GpuArray<amrex::Real,2> orig_x_sub;
253  amrex::GpuArray<amrex::Real,2> orig_y_sub;
254  amrex::Real* orig_x_sub_p = orig_x_sub.data();
255  amrex::Real* orig_y_sub_p = orig_y_sub.data();
256  if (exp_interp) {
257  new_x = new_x_p[new_k];
258  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
259  } else {
260  new_x = new_x_z[new_k];
261  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
262  }
263  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
264 #ifndef AMREX_USE_GPU
265  if (debug) {
266  amrex::Print() << " orig_x_sub_p = [";
267  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
268  amrex::Print() << "]" << std::endl;
269  amrex::Print() << " orig_y_sub_p = [";
270  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
271  amrex::Print() << "]" << std::endl;
272  }
273 #endif
274  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
275  }
276  } else if (order%2 == 0) {
277  if ((kl-(order/2) >= 0) && (kr+order/2 <= orig_n-1)) {
278  // Even order interpolation.
279  // Interpolate twice. Offset the range by 1 the 2nd time. Average the results.
280  amrex::Real new_y_l, new_y_r;
281  {
282  int ksta = kl-(order/2-1);
283  int kend = ksta+order;
284 #ifndef AMREX_USE_GPU
285  int ksize = kend-ksta;
286  if (debug) amrex::Print() << " (2a) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
287 #endif
288  amrex::Real new_x;
289  amrex::GpuArray<amrex::Real,10> orig_x_sub;
290  amrex::GpuArray<amrex::Real,10> orig_y_sub;
291  amrex::Real* orig_x_sub_p = orig_x_sub.data();
292  amrex::Real* orig_y_sub_p = orig_y_sub.data();
293  if (exp_interp) {
294  new_x = new_x_p[new_k];
295  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
296  } else {
297  new_x = new_x_z[new_k];
298  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
299  }
300  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
301 #ifndef AMREX_USE_GPU
302  if (debug) {
303  amrex::Print() << " orig_x_sub_p = [";
304  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
305  amrex::Print() << "]" << std::endl;
306  amrex::Print() << " orig_y_sub_p = [";
307  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
308  amrex::Print() << "]" << std::endl;
309  }
310 #endif
311  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_l);
312  }
313  {
314  int ksta = kl-order/2;
315  int kend = ksta+order;
316 #ifndef AMREX_USE_GPU
317  int ksize = kend-ksta;
318  if (debug) amrex::Print() << "new_k=" << new_k << " (2b) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
319 #endif
320  amrex::Real new_x;
321  amrex::GpuArray<amrex::Real,10> orig_x_sub;
322  amrex::GpuArray<amrex::Real,10> orig_y_sub;
323  amrex::Real* orig_x_sub_p = orig_x_sub.data();
324  amrex::Real* orig_y_sub_p = orig_y_sub.data();
325  if (exp_interp) {
326  new_x = new_x_p[new_k];
327  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
328  } else {
329  new_x = new_x_z[new_k];
330  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
331  }
332  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
333 #ifndef AMREX_USE_GPU
334  if (debug) {
335  amrex::Print() << " orig_x_sub_p = [";
336  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
337  amrex::Print() << "]" << std::endl;
338  amrex::Print() << " orig_y_sub_p = [";
339  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
340  amrex::Print() << "]" << std::endl;
341  }
342 #endif
343  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_r);
344  }
345  new_y[new_k] = myhalf*(new_y_l+new_y_r);
346  } else if ((kl-(order/2-1) >= 0) && (kr+order/2 <= orig_n-1)) {
347  int ksta = kl-(order/2-1);
348  int kend = ksta+order;
349 #ifndef AMREX_USE_GPU
350  int ksize = kend-ksta;
351  if (debug) amrex::Print() << " (3) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
352 #endif
353  amrex::Real new_x;
354  amrex::GpuArray<amrex::Real,10> orig_x_sub;
355  amrex::GpuArray<amrex::Real,10> orig_y_sub;
356  amrex::Real* orig_x_sub_p = orig_x_sub.data();
357  amrex::Real* orig_y_sub_p = orig_y_sub.data();
358  if (exp_interp) {
359  new_x = new_x_p[new_k];
360  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
361  } else {
362  new_x = new_x_z[new_k];
363  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
364  }
365  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
366 #ifndef AMREX_USE_GPU
367  if (debug) {
368  amrex::Print() << " orig_x_sub_p = [";
369  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
370  amrex::Print() << "]" << std::endl;
371  amrex::Print() << " orig_y_sub_p = [";
372  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
373  amrex::Print() << "]" << std::endl;
374  }
375 #endif
376  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
377  } else if ((kl-order/2 >= 0) && (kr+order/2-1 <= orig_n-1)) {
378  int ksta = kl-order/2;
379  int kend = ksta+order;
380 #ifndef AMREX_USE_GPU
381  int ksize = kend-ksta;
382  if (debug) amrex::Print() << " (4) order=" << order << " new_x_z=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
383 #endif
384  amrex::Real new_x;
385  amrex::GpuArray<amrex::Real,10> orig_x_sub;
386  amrex::GpuArray<amrex::Real,10> orig_y_sub;
387  amrex::Real* orig_x_sub_p = orig_x_sub.data();
388  amrex::Real* orig_y_sub_p = orig_y_sub.data();
389  if (exp_interp) {
390  new_x = new_x_p[new_k];
391  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
392  } else {
393  new_x = new_x_z[new_k];
394  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
395  }
396  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
397 #ifndef AMREX_USE_GPU
398  if (debug) {
399  amrex::Print() << " orig_x_sub_p = [";
400  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
401  amrex::Print() << "]" << std::endl;
402  amrex::Print() << " orig_y_sub_p = [";
403  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
404  amrex::Print() << "]" << std::endl;
405  }
406 #endif
407  lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
408  } else {
409  // Linear interpolation.
410  int ksta = kl;
411  int kend = kr;
412 #ifndef AMREX_USE_GPU
413  int ksize = kend-ksta+1;
414  if (debug) amrex::Print() << " (5) order=" << order << " new_x=" << new_x_z[new_k] << " new_x_p=" << new_x_p[new_k] << " kl=" << kl << " kr=" << kr << " ksta=" << ksta << " kend=" << kend << std::endl;
415 #endif
416  amrex::Real new_x;
417  amrex::GpuArray<amrex::Real,2> orig_x_sub;
418  amrex::GpuArray<amrex::Real,2> orig_y_sub;
419  amrex::Real* orig_x_sub_p = orig_x_sub.data();
420  amrex::Real* orig_y_sub_p = orig_y_sub.data();
421  if (exp_interp) {
422  new_x = new_x_p[new_k];
423  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; }
424  } else {
425  new_x = new_x_z[new_k];
426  for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; }
427  }
428  for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; }
429 #ifndef AMREX_USE_GPU
430  if (debug) {
431  amrex::Print() << " orig_x_sub_p = [";
432  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k];
433  amrex::Print() << "]" << std::endl;
434  amrex::Print() << " orig_y_sub_p = [";
435  for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k];
436  amrex::Print() << "]" << std::endl;
437  }
438 #endif
439  lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]);
440  }
441  }
442 #ifndef AMREX_USE_GPU
443  if (debug) amrex::Print() << " new_y[" << new_k << "]=" << new_y[new_k] << std::endl;
444 #endif
445  }
446 }
Cp_d
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:23
lagrange_interp
AMREX_FORCE_INLINE AMREX_GPU_DEVICE void lagrange_interp(const int &order, amrex::Real *x, amrex::Real *y, amrex::Real &new_x, amrex::Real &new_y)
Definition: ERF_MetgridUtils.H:117

Referenced by interpolate_column_metgrid().

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

◆ read_from_metgrid()

void read_from_metgrid ( int  lev,
int  itime,
const amrex::Box &  domain,
const std::string &  fname,
std::string &  NC_dateTime,
amrex::Real &  NC_epochTime,
int &  flag_psfc,
int &  flag_msf,
int &  flag_sst,
int &  flag_tsk,
int &  flag_lmask,
int &  NC_nx,
int &  NC_ny,
amrex::Real &  NC_dx,
amrex::Real &  NC_dy,
amrex::FArrayBox &  NC_xvel_fab,
amrex::FArrayBox &  NC_yvel_fab,
amrex::FArrayBox &  NC_temp_fab,
amrex::FArrayBox &  NC_rhum_fab,
amrex::FArrayBox &  NC_pres_fab,
amrex::FArrayBox &  NC_ght_fab,
amrex::FArrayBox &  NC_hgt_fab,
amrex::FArrayBox &  NC_psfc_fab,
amrex::FArrayBox &  NC_msfu_fab,
amrex::FArrayBox &  NC_msfv_fab,
amrex::FArrayBox &  NC_msfm_fab,
amrex::FArrayBox &  NC_sst_fab,
amrex::FArrayBox &  NC_tsk_fab,
amrex::FArrayBox &  NC_LAT_fab,
amrex::FArrayBox &  NC_LON_fab,
amrex::IArrayBox &  NC_lmask_iab,
amrex::Geometry &  geom 
)

◆ rh_to_mxrat()

AMREX_FORCE_INLINE AMREX_GPU_DEVICE void rh_to_mxrat ( int  i,
int  j,
int  k,
const amrex::Array4< amrex::Real const > &  rhum,
const amrex::Array4< amrex::Real const > &  temp,
const amrex::Array4< amrex::Real const > &  pres,
int  src_indx,
const amrex::Array4< amrex::Real > &  mxrat 
)
929 {
930  amrex::Real qv_max_p_safe = amrex::Real(10000.0); // WRF default value
931  amrex::Real qv_max_flag = amrex::Real(1.0e-5); // WRF default value
932  amrex::Real qv_max_value = amrex::Real(3.0e-6); // WRF default value
933  amrex::Real qv_min_p_safe = amrex::Real(110000.0); // WRF default value
934  amrex::Real qv_min_flag = amrex::Real(1.0e-6); // WRF default value
935  amrex::Real qv_min_value = amrex::Real(1.0e-6); // WRF default value
936  amrex::Real eps = amrex::Real(0.622);
937  amrex::Real svp1 = amrex::Real(0.6112);
938  amrex::Real svp2 = amrex::Real(17.67);
939  amrex::Real svp3 = amrex::Real(29.65);
940  amrex::Real svpt0 = amrex::Real(273.15);
941  // WRF's method when model_config_rec%rh2qv_wrt_liquid=.true. (default behavior)
942  if (temp(i,j,k) != zero) {
943  amrex::Real es=amrex::Real(0.01)*rhum(i,j,k)*svp1*amrex::Real(10.0)*std::exp(svp2*(temp(i,j,k)-svpt0)/(temp(i,j,k)-svp3));
944  if (es >= pres(i,j,k)/amrex::Real(100.0)) {
945  // vapor pressure exceeds total pressure
946  mxrat(i,j,k,src_indx) = amrex::Math::powi<-6>(amrex::Real(10.0));
947  }
948  else {
949  mxrat(i,j,k,src_indx) = amrex::max(eps*es/(pres(i,j,k)/amrex::Real(100.0)-es), amrex::Real(1.0e-6));
950  }
951  }
952  else {
953  // I don't know why there's a fringe case handled in WRF where T is absolute zero...
954  // Let's just deal with it here in case we also end up needing it.
955  mxrat(i,j,k,src_indx) = amrex::Real(1.0e-6);
956  }
957  // See the below comment from WRF dyn_em/module_initialize_real.F rh_to_mxrat1.
958  // For pressures above a defined level, reasonable Qv values should be
959  // a certain value or smaller. If they are larger than this, the input data
960  // probably had "missing" RH, and we filled in some values. This is an
961  // attempt to catch those. Also, set the minimum value for the entire
962  // domain that is above the selected pressure level.
963  if (pres(i,j,k) < qv_max_p_safe) {
964  if (mxrat(i,j,k,src_indx) > qv_max_flag) {
965  mxrat(i,j,k,src_indx) = qv_max_value;
966  }
967  }
968  if (pres(i,j,k) < qv_min_p_safe) {
969  if (mxrat(i,j,k,src_indx) < qv_min_flag) {
970  mxrat(i,j,k,src_indx) = qv_min_value;
971  }
972  }
973 }
MicVar_Kess::pres
@ pres
Definition: ERF_Kessler.H:25
module_model_constants::svp1
real(c_double), parameter svp1
Definition: ERF_module_model_constants.F90:78
module_model_constants::svp3
real(c_double), parameter svp3
Definition: ERF_module_model_constants.F90:80
module_model_constants::svp2
real(c_double), parameter svp2
Definition: ERF_module_model_constants.F90:79
module_model_constants::svpt0
real(c_double), parameter svpt0
Definition: ERF_module_model_constants.F90:81