ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_DataStruct.H
Go to the documentation of this file.
1 #ifndef ERF_DATA_STRUCT_H_
2 #define ERF_DATA_STRUCT_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_Constants.H>
13 #include <ERF_IndexDefines.H>
14 #include <ERF_AdvStruct.H>
15 #include <ERF_DiffStruct.H>
16 #include <ERF_SpongeStruct.H>
17 #include <ERF_TurbStruct.H>
18 #include <ERF_TurbPertStruct.H>
19 
20 enum MapFacType {
21 // This version assumes isotropic
22  m_x, u_x, v_x, num,
23  m_y = 0, u_y = 1, v_y = 2
24 // This version allows for non-isotropic
25 // m_x, u_x, v_x,
26 // m_y, u_y, v_y, num
27 };
28 
29 enum TauType {
31 };
32 
33 AMREX_ENUM(InitType,
34  None, Input_Sounding, NCFile, WRFInput, Metgrid, Uniform
35 );
36 
37 AMREX_ENUM(SoundingType,
38  ConstantDensity, Ideal, Isentropic, DryIsentropic
39 );
40 
41 AMREX_ENUM(ABLDriverType,
42  None, PressureGradient, GeostrophicWind
43 );
44 
45 AMREX_ENUM(CouplingType,
46  OneWay, TwoWay
47 );
48 
49 AMREX_ENUM(SubsteppingType,
50  None, Explicit, Implicit
51 );
52 
53 AMREX_ENUM(MeshType,
54  ConstantDz, StretchedDz, VariableDz
55 );
56 
57 AMREX_ENUM(TerrainType,
58  None, StaticFittedMesh, MovingFittedMesh, EB, ImmersedForcing
59 );
60 
61 AMREX_ENUM(MoistureModelType,
62  Eulerian, Lagrangian, Undefined
63 );
64 
65 AMREX_ENUM(MoistureType,
66  SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler, Kessler_NoRain, SatAdj, Morrison, Morrison_NoIce, None
67 );
68 
69 AMREX_ENUM(WindFarmType,
70  Fitch, EWP, SimpleAD, GeneralAD, None
71 );
72 
73 AMREX_ENUM(WindFarmLocType,
74  lat_lon, x_y, None
75 );
76 
77 AMREX_ENUM(LandSurfaceType,
78  SLM, MM5, None, NOAHMP
79 );
80 
81 AMREX_ENUM(RadiationType,
82  None, RRTMGP
83 );
84 
85 enum struct Coord {
86  x, y, z
87 };
88 
89 // These are used as integers so must be enum not enum struct
90 enum Rayleigh {
92 };
93 
94 // These are used as integers so must be enum not enum struct
95 enum Sponge {
97 };
98 
100  int qv = -1; // Water vapor
101  int qc = -1; // Cloud liquid water
102  int qi = -1; // Cloud ice
103  int qr = -1; // Rain
104  int qs = -1; // Snow
105  int qg = -1; // Graupel
106 
107  // Constructor for easy initialization
108  MoistureComponentIndices (int qv_comp, int qc_comp,
109  int qi_comp=-1,
110  int qr_comp=-1,
111  int qs_comp=-1,
112  int qg_comp=-1)
113  : qv(qv_comp), qc(qc_comp), qi(qi_comp), qr(qr_comp), qs(qs_comp), qg(qg_comp) {}
114 
115  // Default constructor
117 };
118 
119 /**
120  * Container holding many of the algorithmic options and parameters
121  */
122 
123 struct SolverChoice {
124  public:
125  void init_params (int max_level, std::string pp_prefix)
126  {
127  amrex::ParmParse pp(pp_prefix);
128 
129  bool bogus;
130  if (pp.query("use_terrain",bogus) > 0) {
131  amrex::Abort("The input use_terrain is deprecated. Set terrain_type instead.");
132  }
133 
134  pp.query("grid_stretching_ratio", grid_stretching_ratio);
135  if (grid_stretching_ratio != 0) {
136  AMREX_ASSERT_WITH_MESSAGE((grid_stretching_ratio >= 1.),
137  "The grid stretching ratio must be greater than 1");
138  }
139  if (grid_stretching_ratio >= 1) {
140  if (mesh_type == MeshType::ConstantDz) {
141  mesh_type = MeshType::StretchedDz;
142  }
143  if (terrain_type != TerrainType::StaticFittedMesh) {
144  amrex::Print() << "Turning terrain on to enable grid stretching" << std::endl;
145  terrain_type = TerrainType::StaticFittedMesh;
146  }
147  pp.query("zsurface", zsurf);
148  if (zsurf != 0.0) {
149  amrex::Print() << "Nominal zsurface height != 0, may result in unexpected behavior"
150  << std::endl;
151  }
152  pp.get("initial_dz", dz0);
153  }
154 
155  // Do we set map scale factors to 0.5 instead of 1 for testing?
156  pp.query("test_mapfactor", test_mapfactor);
157 
158  // What type of moisture model to use?
159  moisture_type = MoistureType::None; // Default
160  if (pp.query("moisture_type",moisture_type) > 0) {
161  amrex::Abort("The input moisture_type is deprecated. Set moisture_model instead.");
162  }
163  pp.query_enum_case_insensitive("moisture_model",moisture_type);
164  if ( (moisture_type == MoistureType::Morrison) ||
165  (moisture_type == MoistureType::SAM) ) {
167  RhoQ1_comp, // water vapor
168  RhoQ2_comp, // cloud water
169  RhoQ3_comp, // cloud ice
170  RhoQ4_comp, // rain
171  RhoQ5_comp, // snow
172  RhoQ6_comp // graupel
173  );
174  } else if ( (moisture_type == MoistureType::Morrison_NoIce) ||
175  (moisture_type == MoistureType::SAM_NoIce) ) {
177  RhoQ1_comp, // water vapor
178  RhoQ2_comp, // cloud water
179  -1, // cloud ice
180  RhoQ4_comp // rain
181  );
182  } else if ( (moisture_type == MoistureType::SAM_NoPrecip_NoIce) ||
183  (moisture_type == MoistureType::Kessler_NoRain) ||
184  (moisture_type == MoistureType::SatAdj) ) {
186  RhoQ1_comp, // water vapor
187  RhoQ2_comp // cloud water
188  );
189  } else if (moisture_type == MoistureType::Kessler) {
191  RhoQ1_comp, // water vapor
192  RhoQ2_comp, // cloud water
193  RhoQ3_comp // rain
194  );
195  }
196 
197  // TODO: should we set default for dry??
198  // Set a different default for moist vs dry
199  if (moisture_type != MoistureType::None) {
200  if (moisture_type == MoistureType::Kessler_NoRain ||
201  moisture_type == MoistureType::SAM ||
202  moisture_type == MoistureType::SAM_NoIce ||
203  moisture_type == MoistureType::SAM_NoPrecip_NoIce ||
204  moisture_type == MoistureType::Morrison ||
205  moisture_type == MoistureType::Morrison_NoIce ||
206  moisture_type == MoistureType::SatAdj)
207  {
208  buoyancy_type = 1; // uses Rhoprime
209  } else {
210  buoyancy_type = 2; // uses Tprime
211  }
212 
213  pp.query("moisture_tight_coupling",moisture_tight_coupling);
214  }
215 
216  // Which expression (1,2/3 or 4) to use for buoyancy
217  pp.query("buoyancy_type", buoyancy_type);
218 
219  // What type of land surface model to use
220  lsm_type = LandSurfaceType::None; // Default
221  pp.query_enum_case_insensitive("land_surface_model",lsm_type);
222 
223  // What type of radiation model to use
224  rad_type = RadiationType::None; // Default
225  pp.query_enum_case_insensitive("radiation_model", rad_type);
226 
227  // Verify that radiation model cannot be RRTMGP if ERF was not compiled with RRTMGP
228 #ifndef ERF_USE_RRTMGP
229  if (rad_type == RadiationType::RRTMGP)
230  {
231  amrex::Abort("ERF was not compiled with RRTMGP enabled!");
232  }
233 #endif
234 
235  // Is the terrain none, static or moving?
236  std::string terrain_type_temp = "";
237  pp.query("terrain_type", terrain_type_temp);
238  if (terrain_type_temp == "Moving") {
239  amrex::Warning("erf.terrain_type = Moving is deprecated; please replace Moving by MovingFittedMesh");
240  terrain_type = TerrainType::MovingFittedMesh;
241  } else if (terrain_type_temp == "Static") {
242  amrex::Warning("erf.terrain_type = Static is deprecated; please replace Static by StaticFittedMesh");
243  terrain_type = TerrainType::StaticFittedMesh;
244  } else {
245  pp.query_enum_case_insensitive("terrain_type",terrain_type);
246  }
247 
248  //
249  // Read the init_type here to make sure we correctly set the mesh and terrain types
250  //
251  std::string init_type_temp_string;
252  pp.query("init_type",init_type_temp_string);
253  if ( (init_type_temp_string == "Real") || (init_type_temp_string == "real") ) {
254  amrex::Warning("erf.init_type = Real is deprecated; please replace Real by WRFInput");
255  init_type = InitType::WRFInput;
256  use_real_bcs = true;
257  pp.query("upwind_real_bcs",upwind_real_bcs);
258  } else if ( (init_type_temp_string == "Ideal") || (init_type_temp_string == "ideal") ) {
259  amrex::Warning("erf.init_type = Ideal is deprecated; please replace Ideal by WRFInput");
260  init_type = InitType::WRFInput;
261  use_real_bcs = false;
262  } else {
263  pp.query_enum_case_insensitive("init_type",init_type);
264  use_real_bcs = ( (init_type == InitType::WRFInput) || (init_type == InitType::Metgrid) );
265  if (use_real_bcs) {
266  pp.query("upwind_real_bcs",upwind_real_bcs);
267  }
268  }
269 
270  if ( (init_type == InitType::WRFInput) || (init_type == InitType::Metgrid) ) {
271  if (terrain_type != TerrainType::StaticFittedMesh) {
272  amrex::Abort("Only terrain_type = StaticFittedMesh are allowed with init_type = WRFInput or Metgrid");
273  }
274  }
275 
276  if (init_type == InitType::WRFInput) {
277  if (moisture_type == MoistureType::None) {
278  amrex::Abort("Can't have moisture_type = None with init_type = WRFInput");
279  }
280 
281  // NetCDF wrfbdy lateral boundary file
282  std::string nc_bdy_file_temp_string;
283  bool has_bdy = pp.query("nc_bdy_file", nc_bdy_file_temp_string);
284  if (!has_bdy) use_real_bcs = false;
285 
286  bool use_real_bcs_temp = use_real_bcs;
287  pp.query("use_real_bcs", use_real_bcs_temp);
288  if (use_real_bcs && !use_real_bcs_temp) {
289  use_real_bcs = false;
290  }
291  if (use_real_bcs) {
292  pp.query("upwind_real_bcs",upwind_real_bcs);
293  }
294  }
295 
296  // Check for rebalancing with wrfinput
297  if (init_type == InitType::WRFInput) {
298  pp.query("rebalance_wrfinput",rebalance_wrfinput);
299  }
300 
301  // How to interpret input_sounding
302  if (init_type == InitType::Input_Sounding) {
303  pp.query_enum_case_insensitive("sounding_type",sounding_type);
304  }
305 
306  if (terrain_type == TerrainType::StaticFittedMesh ||
307  terrain_type == TerrainType::MovingFittedMesh) {
308  mesh_type = MeshType::VariableDz;
309  }
310 
311  int n_zlevels = pp.countval("terrain_z_levels");
312  if (n_zlevels > 0)
313  {
314  if (terrain_type == TerrainType::None) {
315  terrain_type = TerrainType::StaticFittedMesh;
316  }
317  if (mesh_type == MeshType::ConstantDz) {
318  mesh_type = MeshType::StretchedDz;
319  }
320  }
321 
322  // Use lagged_delta_rt in the fast integrator?
323  pp.query("use_lagged_delta_rt", use_lagged_delta_rt);
324 
325  // These default to true but are used for unit testing
326  pp.query("use_gravity", use_gravity);
328 
329  pp.query("c_p", c_p);
330  rdOcp = R_d / c_p;
331 
332  read_int_string(max_level, "anelastic", anelastic, 0);
333 
334  // *******************************************************************************
335  // Read substepping_type and allow for different values at each level
336  // *******************************************************************************
337  substepping_type.resize(max_level+1);
338 
339  for (int i = 0; i <= max_level; i++) {
340  substepping_type[i] = SubsteppingType::Implicit;
341  }
342 
343  int nvals = pp.countval("substepping_type");
344  AMREX_ALWAYS_ASSERT(nvals == 0 || nvals == 1 || nvals >= max_level+1);
345 
346  if (nvals == 1) {
347  pp.query_enum_case_insensitive("substepping_type",substepping_type[0]);
348  for (int i = 1; i <= max_level; i++) {
350  }
351  } else if (nvals > 1) { // in this case we have asserted nvals >= max_level+1
352  for (int i = 0; i <= max_level; i++) {
353  pp.query_enum_case_insensitive("substepping_type",substepping_type[i],i);
354  }
355  }
356 
357  // *******************************************************************************
358  // Error check on deprecated input
359  // *******************************************************************************
360  int nvals_old = pp.countval("no_substepping");
361  if (nvals_old > 0) {
362  amrex::Abort("The no_substepping flag is deprecated -- set substepping_type instead");
363  }
364  // *******************************************************************************
365 
366  bool any_anelastic = false;
367  for (int i = 0; i <= max_level; ++i) {
368  if (anelastic[i] == 1) any_anelastic = true;
369  }
370 
371  if (any_anelastic == 1) {
373  fixed_density = true; // We default to true but are allowed to override below
374  buoyancy_type = 3; // (This isn't actually used when anelastic is set)
375  } else {
376  pp.query("project_initial_velocity", project_initial_velocity);
377  }
378 
379  pp.query("fixed_density", fixed_density);
380 
381  // *******************************************************************************
382 
383  pp.query("ncorr", ncorr);
384  pp.query("poisson_abstol", poisson_abstol);
385  pp.query("poisson_reltol", poisson_reltol);
386 
387  for (int lev = 0; lev <= max_level; lev++) {
388  if (anelastic[lev] != 0)
389  {
390  substepping_type[lev] = SubsteppingType::None;
391  }
392  }
393 
394  pp.query("force_stage1_single_substep", force_stage1_single_substep);
395 
396  // Include Coriolis forcing?
397  pp.query("use_coriolis", use_coriolis);
398  pp.query("has_lat_lon", has_lat_lon);
399  pp.query("variable_coriolis", variable_coriolis);
400 
401  // Include Rayleigh damping (separate flags for each variable)
402  pp.query("rayleigh_damp_U", rayleigh_damp_U);
403  pp.query("rayleigh_damp_V", rayleigh_damp_V);
404  pp.query("rayleigh_damp_W", rayleigh_damp_W);
405  pp.query("rayleigh_damp_T", rayleigh_damp_T);
406  pp.query("rayleigh_dampcoef", rayleigh_dampcoef);
407  pp.query("rayleigh_zdamp", rayleigh_zdamp);
408  pp.query("rayleigh_damp_substep", rayleigh_damp_substep); // apply Rayleigh damping source terms in substep only
409 
410  // flags for whether to apply other source terms in substep only
411  pp.query("immersed_forcing_substep", immersed_forcing_substep); // apply Rayleigh damping source terms in substep only
412  pp.query("forest_substep", forest_substep); // apply Rayleigh damping source terms in substep only
413 
414  // Flag to do MOST rotations with terrain
415  pp.query("use_rotate_surface_flux",use_rotate_surface_flux);
417  AMREX_ASSERT_WITH_MESSAGE(terrain_type != TerrainType::None,"MOST stress rotations are only valid with terrain!");
418  }
419 
420  // Which external forcings?
421  abl_driver_type = ABLDriverType::None; // Default: no ABL driver for simulating classical fluid dynamics problems
422  pp.query_enum_case_insensitive("abl_driver_type",abl_driver_type);
423  pp.query("const_massflux_u", const_massflux_u);
424  pp.query("const_massflux_v", const_massflux_v);
425  pp.query("const_massflux_tau", const_massflux_tau);
426  pp.query("const_massflux_layer_lo", const_massflux_layer_lo);
427  pp.query("const_massflux_layer_hi", const_massflux_layer_hi);
428 
429  // Which type of inflow turbulent generation
430  pert_type = PerturbationType::None; // Default
431  pp.query_enum_case_insensitive("perturbation_type",pert_type);
432 
433  amrex::Vector<amrex::Real> abl_pressure_grad_in = {0.0, 0.0, 0.0};
434  pp.queryarr("abl_pressure_grad",abl_pressure_grad_in);
435  for(int i = 0; i < AMREX_SPACEDIM; ++i) abl_pressure_grad[i] = abl_pressure_grad_in[i];
436 
437  amrex::Vector<amrex::Real> abl_geo_forcing_in = {0.0, 0.0, 0.0};
438  if(pp.queryarr("abl_geo_forcing",abl_geo_forcing_in)) {
439  amrex::Print() << "Specified abl_geo_forcing: (";
440  for (int i = 0; i < AMREX_SPACEDIM; ++i) {
441  abl_geo_forcing[i] = abl_geo_forcing_in[i];
442  amrex::Print() << abl_geo_forcing[i] << " ";
443  }
444  amrex::Print() << ")" << std::endl;
445  }
446 
447  if (use_coriolis)
448  {
450  }
451 
452  pp.query("add_custom_rhotheta_forcing", custom_rhotheta_forcing);
453  pp.query("add_custom_moisture_forcing", custom_moisture_forcing);
454  pp.query("add_custom_w_subsidence", custom_w_subsidence);
455  pp.query("add_custom_geostrophic_profile", custom_geostrophic_profile);
456  pp.query("custom_forcing_uses_primitive_vars", custom_forcing_prim_vars);
457 
458  pp.query("nudging_from_input_sounding", nudging_from_input_sounding);
459 
461  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!(!abl_geo_wind_table.empty() && custom_geostrophic_profile),
462  "Should not have both abl_geo_wind_table and custom_geostrophic_profile set.");
463 
464  pp.query("Ave_Plane", ave_plane);
465 
466  pp.query("use_moist_background", use_moist_background);
467 
468  // Use numerical diffusion?
469  pp.query("num_diff_coeff",num_diff_coeff);
470  AMREX_ASSERT_WITH_MESSAGE(( (num_diff_coeff >= 0.) && (num_diff_coeff <= 1.) ),
471  "Numerical diffusion coefficient must be between 0 & 1.");
473  if (use_num_diff) {
474  amrex::Print() << "6th-order numerical diffusion turned on with coefficient = "
475  << num_diff_coeff << std::endl;
476  num_diff_coeff *= std::pow(2.0,-6);
477  }
478 
479  advChoice.init_params(pp_prefix);
480  diffChoice.init_params(pp_prefix);
481  spongeChoice.init_params(pp_prefix);
482 
483  turbChoice.resize(max_level+1);
484  for (int lev = 0; lev <= max_level; lev++) {
485  turbChoice[lev].init_params(lev,max_level,pp_prefix);
486  }
487 
488  // YSU PBL: use consistent coriolis frequency
489  for (int lev = 0; lev <= max_level; lev++) {
490  if (turbChoice[lev].pbl_ysu_use_consistent_coriolis) {
491  if (use_coriolis) {
492  turbChoice[lev].pbl_ysu_coriolis_freq = coriolis_factor * sinphi;
493  if (lev == 0) {
494  amrex::Print() << "YSU PBL using ERF coriolis frequency: " << turbChoice[lev].pbl_ysu_coriolis_freq << std::endl;
495  }
496  } else {
497  amrex::Abort("YSU cannot use ERF coriolis frequency if not using coriolis");
498  }
499  }
500  }
501  // MRF
502  for (int lev = 0; lev <= max_level; lev++) {
503  if (turbChoice[lev].pbl_ysu_use_consistent_coriolis) {
504  if (use_coriolis) {
505  turbChoice[lev].pbl_ysu_coriolis_freq = coriolis_factor * sinphi;
506  if (lev == 0) {
507  amrex::Print() << "MRF PBL using ERF coriolis frequency: " << turbChoice[lev].pbl_ysu_coriolis_freq << std::endl;
508  }
509  } else {
510  amrex::Abort("MRF cannot use ERF coriolis frequency if not using coriolis");
511  }
512  }
513  }
514 
515  // Are we using SHOC?
516  pp.query("use_shoc", use_shoc);
517 #ifndef ERF_USE_SHOC
518  if (use_shoc) {
519  amrex::Abort("You set use_shoc to true but didn't build with SHOC; you must rebuild the executable");
520  }
521 #endif
522 
523  // Which type of multilevel coupling
524  coupling_type = CouplingType::TwoWay; // Default
525  pp.query_enum_case_insensitive("coupling_type",coupling_type);
526 
527  // Which type of windfarm model
528  windfarm_type = WindFarmType::None; // Default
529  pp.query_enum_case_insensitive("windfarm_type",windfarm_type);
530 
531  static std::string windfarm_loc_type_string = "None";
532  windfarm_loc_type = WindFarmLocType::None;
533  pp.query_enum_case_insensitive("windfarm_loc_type",windfarm_loc_type);
534 
535  pp.query("windfarm_loc_table", windfarm_loc_table);
536  pp.query("windfarm_spec_table", windfarm_spec_table);
537  pp.query("windfarm_blade_table", windfarm_blade_table);
538  pp.query("windfarm_airfoil_tables", windfarm_airfoil_tables);
539  pp.query("windfarm_spec_table_extra", windfarm_spec_table_extra);
540 
541  // Sampling distance upstream of the turbine to find the
542  // incoming free stream velocity as a factor of the diameter of the
543  // turbine. ie. the sampling distance will be this number multiplied
544  // by the diameter of the turbine
545  pp.query("sampling_distance_by_D", sampling_distance_by_D);
546  pp.query("turb_disk_angle_from_x", turb_disk_angle);
547 
548  pp.query("windfarm_x_shift",windfarm_x_shift);
549  pp.query("windfarm_y_shift",windfarm_y_shift);
550  // Test if time averaged data is to be output
551  pp.query("time_avg_vel",time_avg_vel);
552 
553  pp.query("do_hurricane_simulation", do_hurricane_simulation);
554 
555  check_params(max_level);
556  }
557 
558  void check_params (int max_level)
559  {
560  // Warn for PBL models and moisture - these may not yet be compatible
561  for (int lev = 0; lev <= max_level; lev++) {
562  if ((moisture_type != MoistureType::None) && (turbChoice[lev].pbl_type != PBLType::None)) {
563  amrex::Warning("\n*** WARNING: Moisture may not yet be compatible with PBL models, \n proceed with caution ***");
564  }
565  }
566  //
567  // Buoyancy type check
568  //
569  if (buoyancy_type != 1 && buoyancy_type != 2 && buoyancy_type != 3 && buoyancy_type != 4) {
570  amrex::Abort("buoyancy_type must be 1, 2, 3 or 4");
571  }
572 
573  if (!use_lagged_delta_rt && !(terrain_type == TerrainType::MovingFittedMesh)) {
574  amrex::Error("Can't turn off lagged_delta_rt when terrain not moving");
575  }
576 
577  //
578  // Wind farm checks
579  //
580  if (windfarm_type==WindFarmType::SimpleAD and sampling_distance_by_D < 0.0) {
581  amrex::Abort("To use simplified actuator disks, you need to provide a variable"
582  " erf.sampling_distance_by_D in the inputs which specifies the upstream"
583  " distance as a factor of the turbine diameter at which the incoming free stream"
584  " velocity will be computed at.");
585  }
586  if ( (windfarm_type==WindFarmType::SimpleAD ||
587  windfarm_type==WindFarmType::GeneralAD ) && turb_disk_angle < 0.0) {
588  amrex::Abort("To use simplified actuator disks, you need to provide a variable"
589  " erf.turb_disk_angle_from_x in the inputs which is the angle of the face of the"
590  " turbine disk from the x-axis. A turbine facing an oncoming flow in the x-direction"
591  " will have turb_disk_angle value of 90 deg.");
592  }
593  if (windfarm_loc_type == WindFarmLocType::lat_lon and (windfarm_x_shift < 0.0 or windfarm_y_shift < 0.0)) {
594  amrex::Abort("You are using windfarms with latitude-logitude option to position the turbines."
595  " For this you should provide the inputs erf.windfarm_x_shift and"
596  " erf.windfarm_y_shift which are the values by which the bounding box of the"
597  " windfarm is shifted from the x and the y axes.");
598  }
599  }
600 
601  void display (int max_level, std::string pp_prefix)
602  {
603  amrex::Print() << "SOLVER CHOICE: " << std::endl;
604  amrex::Print() << "force_stage1_single_substep : " << force_stage1_single_substep << std::endl;
605  for (int lev = 0; lev <= max_level; lev++) {
606  if (anelastic[lev]) {
607  amrex::Print() << "Level " << lev << " is anelastic" << std::endl;
608  }
609  if (substepping_type[lev] == SubsteppingType::None) {
610  amrex::Print() << "No substepping at level " << lev << std::endl;
611  } else if (substepping_type[lev] == SubsteppingType::Explicit) {
612  amrex::Print() << "Explicit substepping at level " << lev << std::endl;
613  } else if (substepping_type[lev] == SubsteppingType::Implicit) {
614  amrex::Print() << "Implicit substepping at level " << lev << std::endl;
615  }
616  }
617  amrex::Print() << "fixed_density : " << fixed_density << std::endl;
618  amrex::Print() << "use_coriolis : " << use_coriolis << std::endl;
619  amrex::Print() << "use_gravity : " << use_gravity << std::endl;
620 
621  if (moisture_type == MoistureType::SAM) {
622  amrex::Print() << "Moisture Model: SAM" << std::endl;
623  } else if (moisture_type == MoistureType::SAM_NoIce) {
624  amrex::Print() << "Moisture Model: SAM No Ice" << std::endl;
625  } else if (moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
626  amrex::Print() << "Moisture Model: SAM No Precip No Ice" << std::endl;
627  } else if (moisture_type == MoistureType::Morrison) {
628  amrex::Print() << "Moisture Model: Morrison" << std::endl;
629  } else if (moisture_type == MoistureType::Morrison_NoIce) {
630  amrex::Print() << "Moisture Model: Morrison_NoIce" << std::endl;
631  } else if (moisture_type == MoistureType::Kessler) {
632  amrex::Print() << "Moisture Model: Kessler" << std::endl;
633  } else if (moisture_type == MoistureType::Kessler_NoRain) {
634  amrex::Print() << "Moisture Model: Kessler No Rain" << std::endl;
635  } else if (moisture_type == MoistureType::SatAdj) {
636  amrex::Print() << "Moisture Model: Saturation Adjustment" << std::endl;
637  } else {
638  amrex::Print() << "Moisture Model: None" << std::endl;
639  }
640 
641  if (terrain_type == TerrainType::StaticFittedMesh) {
642  amrex::Print() << "Terrain Type: StaticFittedMesh" << std::endl;
643  } else if (terrain_type == TerrainType::MovingFittedMesh) {
644  amrex::Print() << "Terrain Type: MovingFittedMesh" << std::endl;
645  } else if (terrain_type == TerrainType::EB) {
646  amrex::Print() << "Terrain Type: EB" << std::endl;
647  } else if (terrain_type == TerrainType::ImmersedForcing) {
648  amrex::Print() << "Terrain Type: ImmersedForcing" << std::endl;
649  } else {
650  amrex::Print() << "Terrain Type: None" << std::endl;
651  }
652 
653  if (mesh_type == MeshType::ConstantDz) {
654  amrex::Print() << " Mesh Type: ConstantDz" << std::endl;
655  } else if (mesh_type == MeshType::StretchedDz) {
656  amrex::Print() << " Mesh Type: StretchedDz" << std::endl;
657  } else if (mesh_type == MeshType::VariableDz) {
658  amrex::Print() << " Mesh Type: VariableDz" << std::endl;
659  } else {
660  amrex::Abort("No mesh_type set!");
661  }
662 
663  amrex::Print() << "ABL Driver Type: " << std::endl;
664  if (abl_driver_type == ABLDriverType::None) {
665  amrex::Print() << " None" << std::endl;
666  } else if (abl_driver_type == ABLDriverType::PressureGradient) {
667  amrex::Print() << " Pressure Gradient "
668  << amrex::RealVect(abl_pressure_grad[0],abl_pressure_grad[1],abl_pressure_grad[2])
669  << std::endl;
670  } else if (abl_driver_type == ABLDriverType::GeostrophicWind) {
671  amrex::Print() << " Geostrophic Wind "
672  << amrex::RealVect(abl_geo_forcing[0],abl_geo_forcing[1],abl_geo_forcing[2])
673  << std::endl;
674  }
675 
676  if (max_level > 0) {
677  amrex::Print() << "Coupling Type: " << std::endl;
678  if (coupling_type == CouplingType::TwoWay) {
679  amrex::Print() << " Two-way" << std::endl;
680  } else if (coupling_type == CouplingType::OneWay) {
681  amrex::Print() << " One-way" << std::endl;
682  }
683  }
684 
685  if (rad_type == RadiationType::RRTMGP) {
686  amrex::Print() << "Radiation Model: RRTMGP" << std::endl;
687  } else {
688  amrex::Print() << "Radiation Model: None" << std::endl;
689  }
690 
691  amrex::Print() << "Buoyancy_type : " << buoyancy_type << std::endl;
692 
693  advChoice.display(pp_prefix);
696 
697  for (int lev = 0; lev <= max_level; lev++) {
698  turbChoice[lev].display(lev);
699  }
700  }
701 
702  void build_coriolis_forcings_const_lat (std::string pp_prefix)
703  {
704  amrex::ParmParse pp(pp_prefix);
705 
706  // Read the rotational time period (in seconds)
707  amrex::Real rot_time_period = 86400.0;
708  pp.query("rotational_time_period", rot_time_period);
709 
710  coriolis_factor = 2.0 * 2.0 * PI / rot_time_period;
711 
712  amrex::Real latitude = 90.0;
713  pp.query("latitude", latitude);
714 
715  pp.query("coriolis_3d", coriolis_3d);
716 
717  // Convert to radians
718  latitude *= (PI/180.);
719  sinphi = std::sin(latitude);
720  if (coriolis_3d) {
721  cosphi = std::cos(latitude);
722  }
723 
724  amrex::Print() << "Coriolis frequency, f = " << coriolis_factor * sinphi << " 1/s" << std::endl;
725 
726  if (abl_driver_type == ABLDriverType::GeostrophicWind) {
727  // Read in the geostrophic wind -- we only use this to construct
728  // the forcing term so no need to keep it
729  amrex::Vector<amrex::Real> abl_geo_wind(3);
730  pp.queryarr("abl_geo_wind",abl_geo_wind);
731 
732  if(!pp.query("abl_geo_wind_table",abl_geo_wind_table)) {
733  abl_geo_forcing = {
734  -coriolis_factor * (abl_geo_wind[1]*sinphi - abl_geo_wind[2]*cosphi),
735  coriolis_factor * abl_geo_wind[0]*sinphi,
736  -coriolis_factor * abl_geo_wind[0]*cosphi
737  };
738  } else {
739  amrex::Print() << "NOTE: abl_geo_wind_table provided, ignoring input abl_geo_wind" << std::endl;
740  }
741  }
742  }
743 
744  void read_int_string (int max_level, const char* string_to_read,
745  amrex::Vector<int>& vec_to_fill, int default_int)
746  {
747  amrex::ParmParse pp("erf");
748  int nvals = pp.countval(string_to_read);
749  AMREX_ALWAYS_ASSERT(nvals == 0 || nvals == 1 || nvals >= max_level+1);
750  amrex::Vector<int> temp; temp.resize(nvals);
751  pp.queryarr(string_to_read,temp);
752  if (nvals == 0) {
753  for (int i = 0; i <= max_level; ++i) vec_to_fill.push_back(default_int);
754  } else if (nvals == 1) {
755  for (int i = 0; i <= max_level; ++i) vec_to_fill.push_back(temp[0]);
756  } else {
757  for (int i = 0; i <= max_level; ++i) vec_to_fill.push_back(temp[i]);
758  }
759  }
760 
761  inline static
762  InitType init_type = InitType::None;
763 
764  inline static
765  SoundingType sounding_type = SoundingType::Ideal;
766 
767  inline static
768  TerrainType terrain_type = TerrainType::None;
769 
770  inline static
771  bool use_real_bcs = false;
772 
773  inline static
774  bool upwind_real_bcs = false;
775 
776  inline static
777  MeshType mesh_type = MeshType::ConstantDz;
778 
779  static
780  void set_mesh_type (MeshType new_mesh_type)
781  {
782  mesh_type = new_mesh_type;
783  }
784 
788  amrex::Vector<TurbChoice> turbChoice;
789 
791 
792  amrex::Vector<SubsteppingType> substepping_type;
793  amrex::Vector<int> anelastic;
794 
795  bool fixed_density = false;
796  int ncorr = 1;
799 
800  bool test_mapfactor = false;
801 
802  int buoyancy_type = 1; // uses rhoprime directly
803 
804  // Specify what additional physics/forcing modules we use
805  bool use_gravity = false;
806  bool use_coriolis = false;
807  bool coriolis_3d = true;
808 
809  bool rayleigh_damp_U = false;
810  bool rayleigh_damp_V = false;
811  bool rayleigh_damp_W = false;
812  bool rayleigh_damp_T = false;
813  amrex::Real rayleigh_dampcoef = 0.2; // inverse time scale [1/s]
814  amrex::Real rayleigh_zdamp = 500.0; // damping layer depth [m]
816  bool rayleigh_damp_substep = false; // use Rayleigh on substep only?
817 
818  // Specify whether to apply other various source terms on substep only
820  bool forest_substep = false;
821 
822  // This defaults to true but can be set to false for moving terrain cases only
823  bool use_lagged_delta_rt = true;
824 
825  // Constants
827  amrex::Real c_p = Cp_d; // specific heat at constant pressure for dry air [J/(kg-K)]
829 
830  // Staggered z levels for vertical grid stretching
834 
836 
837  // Coriolis forcing
841 
842  // User-specified forcings in problem definition
845  bool custom_w_subsidence = false;
848 
849  // Do we use source terms to nudge the solution towards
850  // the time-varying data provided in input sounding files?
852 
853  // MOST stress rotations
855 
856  // Should we use SHOC?
857  bool use_shoc = false;
858 
859  // User wishes to output time averaged velocity fields
860  bool time_avg_vel = false;
861 
862  // Type of perturbation
863  PerturbationType pert_type;
864 
865  // Numerical diffusion
866  bool use_num_diff{false};
868 
869  // Rebalance wrfinput
870  bool rebalance_wrfinput{false};
871 
872  CouplingType coupling_type;
873  MoistureType moisture_type;
874  WindFarmType windfarm_type;
875  WindFarmLocType windfarm_loc_type;
876  LandSurfaceType lsm_type;
877  RadiationType rad_type;
878 
879  ABLDriverType abl_driver_type;
880  amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> abl_pressure_grad;
881  amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> abl_geo_forcing;
882  std::string abl_geo_wind_table;
883  bool have_geo_wind_profile {false};
884 
885  bool has_lat_lon{false};
886  bool variable_coriolis{false};
887 
888  int ave_plane {2};
889 
890  // Microphysics params
891  bool use_moist_background {false};
893 
895 
902 
903  // Use forest canopy model?
904  bool do_forest_drag {false};
905 
906  // Enforce constant mass flux?
912  int massflux_klo {0}; // these are updated in ERF.cpp
913  int massflux_khi {0};
914 
916 
917 };
918 #endif
constexpr amrex::Real Cp_d
Definition: ERF_Constants.H:12
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
constexpr amrex::Real R_d
Definition: ERF_Constants.H:10
TauType
Definition: ERF_DataStruct.H:29
@ tau12
Definition: ERF_DataStruct.H:30
@ tau23
Definition: ERF_DataStruct.H:30
@ tau33
Definition: ERF_DataStruct.H:30
@ tau22
Definition: ERF_DataStruct.H:30
@ tau11
Definition: ERF_DataStruct.H:30
@ tau32
Definition: ERF_DataStruct.H:30
@ tau31
Definition: ERF_DataStruct.H:30
@ tau21
Definition: ERF_DataStruct.H:30
@ tau13
Definition: ERF_DataStruct.H:30
Rayleigh
Definition: ERF_DataStruct.H:90
@ ubar
Definition: ERF_DataStruct.H:91
@ wbar
Definition: ERF_DataStruct.H:91
@ nvars
Definition: ERF_DataStruct.H:91
@ vbar
Definition: ERF_DataStruct.H:91
@ thetabar
Definition: ERF_DataStruct.H:91
Sponge
Definition: ERF_DataStruct.H:95
@ nvars_sponge
Definition: ERF_DataStruct.H:96
@ vbar_sponge
Definition: ERF_DataStruct.H:96
@ ubar_sponge
Definition: ERF_DataStruct.H:96
MapFacType
Definition: ERF_DataStruct.H:20
@ v_x
Definition: ERF_DataStruct.H:22
@ num
Definition: ERF_DataStruct.H:22
@ u_y
Definition: ERF_DataStruct.H:23
@ v_y
Definition: ERF_DataStruct.H:23
@ m_y
Definition: ERF_DataStruct.H:23
@ u_x
Definition: ERF_DataStruct.H:22
@ m_x
Definition: ERF_DataStruct.H:22
Coord
Definition: ERF_DataStruct.H:85
AMREX_ENUM(InitType, None, Input_Sounding, NCFile, WRFInput, Metgrid, Uniform)
#define RhoQ4_comp
Definition: ERF_IndexDefines.H:45
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
#define RhoQ3_comp
Definition: ERF_IndexDefines.H:44
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
#define RhoQ6_comp
Definition: ERF_IndexDefines.H:47
#define RhoQ5_comp
Definition: ERF_IndexDefines.H:46
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real pp(amrex::Real y)
Definition: ERF_MicrophysicsUtils.H:230
amrex::Real Real
Definition: ERF_ShocInterface.H:16
Definition: ERF_EWP.H:9
Definition: ERF_Fitch.H:9
Definition: ERF_GeneralAD.H:8
Definition: ERF_Kessler.H:38
Definition: ERF_MM5.H:26
Definition: ERF_Morrison.H:58
Definition: ERF_NOAHMP.H:37
Definition: ERF_SAM.H:53
Definition: ERF_SLM.H:26
Definition: ERF_SatAdj.H:41
Definition: ERF_SimpleAD.H:8
@ bogus
Definition: ERF_IndexDefines.H:204
Definition: ERF_AdvStruct.H:19
void display(std::string &pp_prefix)
Definition: ERF_AdvStruct.H:211
void init_params(std::string pp_prefix)
Definition: ERF_AdvStruct.H:21
Definition: ERF_DiffStruct.H:19
void init_params(std::string pp_prefix)
Definition: ERF_DiffStruct.H:21
void display()
Definition: ERF_DiffStruct.H:67
Definition: ERF_DataStruct.H:99
int qs
Definition: ERF_DataStruct.H:104
int qr
Definition: ERF_DataStruct.H:103
MoistureComponentIndices()=default
int qi
Definition: ERF_DataStruct.H:102
int qv
Definition: ERF_DataStruct.H:100
int qc
Definition: ERF_DataStruct.H:101
int qg
Definition: ERF_DataStruct.H:105
MoistureComponentIndices(int qv_comp, int qc_comp, int qi_comp=-1, int qr_comp=-1, int qs_comp=-1, int qg_comp=-1)
Definition: ERF_DataStruct.H:108
Definition: ERF_DataStruct.H:123
bool rayleigh_damp_T
Definition: ERF_DataStruct.H:812
amrex::Real dz0
Definition: ERF_DataStruct.H:833
amrex::Real const_massflux_layer_lo
Definition: ERF_DataStruct.H:910
bool use_lagged_delta_rt
Definition: ERF_DataStruct.H:823
amrex::Real coriolis_factor
Definition: ERF_DataStruct.H:838
static MeshType mesh_type
Definition: ERF_DataStruct.H:777
amrex::Real windfarm_x_shift
Definition: ERF_DataStruct.H:900
bool rayleigh_damp_V
Definition: ERF_DataStruct.H:810
void display(int max_level, std::string pp_prefix)
Definition: ERF_DataStruct.H:601
bool rebalance_wrfinput
Definition: ERF_DataStruct.H:870
amrex::Real poisson_reltol
Definition: ERF_DataStruct.H:798
bool rayleigh_damp_substep
Definition: ERF_DataStruct.H:816
void build_coriolis_forcings_const_lat(std::string pp_prefix)
Definition: ERF_DataStruct.H:702
amrex::Real rayleigh_zdamp
Definition: ERF_DataStruct.H:814
amrex::Real rdOcp
Definition: ERF_DataStruct.H:828
RadiationType rad_type
Definition: ERF_DataStruct.H:877
void read_int_string(int max_level, const char *string_to_read, amrex::Vector< int > &vec_to_fill, int default_int)
Definition: ERF_DataStruct.H:744
std::string windfarm_spec_table
Definition: ERF_DataStruct.H:896
DiffChoice diffChoice
Definition: ERF_DataStruct.H:786
amrex::Real const_massflux_v
Definition: ERF_DataStruct.H:908
bool use_gravity
Definition: ERF_DataStruct.H:805
int ncorr
Definition: ERF_DataStruct.H:796
int force_stage1_single_substep
Definition: ERF_DataStruct.H:790
std::string windfarm_spec_table_extra
Definition: ERF_DataStruct.H:896
amrex::Real cosphi
Definition: ERF_DataStruct.H:839
LandSurfaceType lsm_type
Definition: ERF_DataStruct.H:876
amrex::Real c_p
Definition: ERF_DataStruct.H:827
std::string windfarm_loc_table
Definition: ERF_DataStruct.H:896
amrex::Real gravity
Definition: ERF_DataStruct.H:826
void check_params(int max_level)
Definition: ERF_DataStruct.H:558
bool custom_rhotheta_forcing
Definition: ERF_DataStruct.H:843
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_geo_forcing
Definition: ERF_DataStruct.H:881
bool use_shoc
Definition: ERF_DataStruct.H:857
WindFarmLocType windfarm_loc_type
Definition: ERF_DataStruct.H:875
int massflux_klo
Definition: ERF_DataStruct.H:912
bool moisture_tight_coupling
Definition: ERF_DataStruct.H:894
bool custom_w_subsidence
Definition: ERF_DataStruct.H:845
bool nudging_from_input_sounding
Definition: ERF_DataStruct.H:851
bool rayleigh_damp_U
Definition: ERF_DataStruct.H:809
bool custom_geostrophic_profile
Definition: ERF_DataStruct.H:846
amrex::Real rayleigh_ztop
Definition: ERF_DataStruct.H:815
bool immersed_forcing_substep
Definition: ERF_DataStruct.H:819
amrex::Real grid_stretching_ratio
Definition: ERF_DataStruct.H:831
amrex::Real sinphi
Definition: ERF_DataStruct.H:840
bool have_geo_wind_profile
Definition: ERF_DataStruct.H:883
amrex::Real const_massflux_u
Definition: ERF_DataStruct.H:907
amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > abl_pressure_grad
Definition: ERF_DataStruct.H:880
void init_params(int max_level, std::string pp_prefix)
Definition: ERF_DataStruct.H:125
amrex::Vector< SubsteppingType > substepping_type
Definition: ERF_DataStruct.H:792
bool coriolis_3d
Definition: ERF_DataStruct.H:807
bool use_num_diff
Definition: ERF_DataStruct.H:866
amrex::Real sampling_distance_by_D
Definition: ERF_DataStruct.H:898
bool test_mapfactor
Definition: ERF_DataStruct.H:800
bool use_coriolis
Definition: ERF_DataStruct.H:806
static SoundingType sounding_type
Definition: ERF_DataStruct.H:765
bool custom_moisture_forcing
Definition: ERF_DataStruct.H:844
amrex::Real num_diff_coeff
Definition: ERF_DataStruct.H:867
std::string windfarm_blade_table
Definition: ERF_DataStruct.H:897
bool fixed_density
Definition: ERF_DataStruct.H:795
amrex::Real zsurf
Definition: ERF_DataStruct.H:832
amrex::Vector< TurbChoice > turbChoice
Definition: ERF_DataStruct.H:788
bool variable_coriolis
Definition: ERF_DataStruct.H:886
bool project_initial_velocity
Definition: ERF_DataStruct.H:835
amrex::Vector< int > anelastic
Definition: ERF_DataStruct.H:793
static bool upwind_real_bcs
Definition: ERF_DataStruct.H:774
AdvChoice advChoice
Definition: ERF_DataStruct.H:785
bool use_moist_background
Definition: ERF_DataStruct.H:891
MoistureType moisture_type
Definition: ERF_DataStruct.H:873
bool custom_forcing_prim_vars
Definition: ERF_DataStruct.H:847
std::string abl_geo_wind_table
Definition: ERF_DataStruct.H:882
static TerrainType terrain_type
Definition: ERF_DataStruct.H:768
ABLDriverType abl_driver_type
Definition: ERF_DataStruct.H:879
bool rayleigh_damp_W
Definition: ERF_DataStruct.H:811
PerturbationType pert_type
Definition: ERF_DataStruct.H:863
SpongeChoice spongeChoice
Definition: ERF_DataStruct.H:787
WindFarmType windfarm_type
Definition: ERF_DataStruct.H:874
static InitType init_type
Definition: ERF_DataStruct.H:762
amrex::Real const_massflux_layer_hi
Definition: ERF_DataStruct.H:911
static bool use_real_bcs
Definition: ERF_DataStruct.H:771
int buoyancy_type
Definition: ERF_DataStruct.H:802
amrex::Real poisson_abstol
Definition: ERF_DataStruct.H:797
MoistureComponentIndices moisture_indices
Definition: ERF_DataStruct.H:892
amrex::Real turb_disk_angle
Definition: ERF_DataStruct.H:899
amrex::Real windfarm_y_shift
Definition: ERF_DataStruct.H:901
bool do_hurricane_simulation
Definition: ERF_DataStruct.H:915
bool has_lat_lon
Definition: ERF_DataStruct.H:885
bool use_rotate_surface_flux
Definition: ERF_DataStruct.H:854
bool do_forest_drag
Definition: ERF_DataStruct.H:904
amrex::Real const_massflux_tau
Definition: ERF_DataStruct.H:909
int massflux_khi
Definition: ERF_DataStruct.H:913
bool time_avg_vel
Definition: ERF_DataStruct.H:860
bool forest_substep
Definition: ERF_DataStruct.H:820
amrex::Real rayleigh_dampcoef
Definition: ERF_DataStruct.H:813
CouplingType coupling_type
Definition: ERF_DataStruct.H:872
std::string windfarm_airfoil_tables
Definition: ERF_DataStruct.H:897
static void set_mesh_type(MeshType new_mesh_type)
Definition: ERF_DataStruct.H:780
int ave_plane
Definition: ERF_DataStruct.H:888
Definition: ERF_SpongeStruct.H:15
void display()
Definition: ERF_SpongeStruct.H:45
void init_params(std::string pp_prefix)
Definition: ERF_SpongeStruct.H:17