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