From 40b1355679a805bc44a162db945db67a1f1546e7 Mon Sep 17 00:00:00 2001 From: wiersema1 Date: Fri, 24 May 2024 15:40:56 -0700 Subject: [PATCH] Cleanup and port metgrid-related code from dirty old branches into a single clean branch. These changes include a new vertical interpolation and quality control algorithm for init_type="metgrid", similar to that in WRF's real.exe. Also includes updates to docs and cleanup of some code in touched files. --- Docs/sphinx_doc/Inputs.rst | 75 ++- Exec/DevTests/MetGrid/inputs | 2 +- Source/ERF.H | 9 + Source/ERF.cpp | 9 + Source/IO/NCWpsFile.H | 2 +- Source/IO/ReadFromMetgrid.cpp | 63 +- .../Initialization/ERF_init_from_metgrid.cpp | 449 ++++++++----- Source/Initialization/Metgrid_utils.H | 613 +++++++++++++++++- 8 files changed, 1011 insertions(+), 211 deletions(-) diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index e84316607..9ccf09366 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -954,7 +954,8 @@ List of Parameters | **erf.init_type** | Initialization | “custom”, | “*custom*” | | | type | “ideal”, | | | | | "real", | | -| | |"input_sounding" | | +| | | "input_sounding", | | +| | | "metgrid", | | +----------------------------------+-------------------+--------------------+------------------+ | **erf.input_sounding_file** | Path to WRF-style | String | "input_sounding" | | | input sounding | | | @@ -979,9 +980,70 @@ List of Parameters | | mesoscale data at | | | | | lateral boundaries| | | +----------------------------------+-------------------+--------------------+------------------+ -| **erf.project_initial_velocity** | project initial | Integer | 1 | +| **erf.project_initial_velocity** | Project initial | Integer | 1 | | | velocity? | | | +----------------------------------+-------------------+--------------------+------------------+ +| **erf.real_width** | Lateral boundary | Integer | 0 | +| | total width if | | | +| | use_real_bcs is | | | +| | true | | | ++----------------------------------+-------------------+--------------------+------------------+ +| **erf.real_set_width** | Lateral boundary | Integer | 0 | +| | specified zone | | | +| | width if | | | +| | use_real_bcs is | | | +| | true | | | ++----------------------------------+-------------------+--------------------+------------------+ +| **erf.metgrid_basic_linear** | If init_type is | true or false | false | +| | metgrid, use | | | +| | linear vertical | | | +| | interpolation and | | | +| | no quality | | | +| | control? | | | ++----------------------------------+-------------------+--------------------+------------------+ +| **erf.metgrid_use_below_sfc** | If init_type is | true or false | true | +| | metgrid, use the | | | +| | origin data levels| | | +| | below the surface?| | | ++----------------------------------+-------------------+--------------------+------------------+ +| **erf.metgrid_use_sfc** | If init_type is | true or false | true | +| | metgrid, use the | | | +| | origin data level | | | +| | at the surface? | | | ++----------------------------------+-------------------+--------------------+------------------+ +| **erf.metgrid_retain_sfc** | If init_type is | true or false | false | +| | metgrid, assign | | | +| | the lowest level | | | +| | directly using the| | | +| | surface value from| | | +| | the origin data? | | | ++----------------------------------+-------------------+--------------------+------------------+ +| **erf.metgrid_proximity** | If init_type is | Real | 1000. | +| | metgrid, pressure | | | +| | differential for | | | +| | detecting origin | | | +| | levels that are | | | +| | problematically | | | +| | close together | | | ++----------------------------------+-------------------+--------------------+------------------+ +| **erf.metgrid_order** | If init_type is | Integer | 2 | +| | metgrid, order of | | | +| | the Lagrange | | | +| | polynomial | | | +| | interpolation | | | +| | scheme for | | | +| | vertical | | | +| | interpolation | | | ++----------------------------------+-------------------+--------------------+------------------+ +| **erf.metgrid_force_sfc_k** | If init_type is | Integer | 0 | +| | metgrid, force the| | | +| | origin data | | | +| | surface level to | | | +| | be included in the| | | +| | interpolation for | | | +| | this many ERF | | | +| | vertical levels | | | ++----------------------------------+-------------------+--------------------+------------------+ Notes ----------------- @@ -992,6 +1054,15 @@ If **erf.init_type = real**, the problem is initialized with mesoscale data cont provided via ``erf.nc_init_file``. The mesoscale data are realistic with variation in all three directions. In addition, the lateral boundary conditions must be supplied in a NetCDF files specified by **erf.nc_bdy_file = wrfbdy_d01** +If **erf.init_type = metgrid**, the problem is initialized with data +contained in the first NetCDF file provided via ``erf.nc_init_file_0``. +Lateral boundary conditions are derived from the sequence of NetCDF +files provided via ``erf.nc_init_file_0``. The sequence of +``erf.nc_init_file_0`` should be output from the WRF Preprocessing +System (WPS) listed chronologically starting with the earliest +timestamp. A minimum of two files are required to derive lateral +boundary conditions. + If **erf.init_type = input_sounding**, a WRF-style input sounding is read from ``erf.input_sounding_file``. This text file includes any set of levels that goes at least up to the model top height. The first line includes the surface diff --git a/Exec/DevTests/MetGrid/inputs b/Exec/DevTests/MetGrid/inputs index ccf0ee7dc..6219d7d40 100644 --- a/Exec/DevTests/MetGrid/inputs +++ b/Exec/DevTests/MetGrid/inputs @@ -44,7 +44,7 @@ erf.check_int = 100 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # prefix of plotfile name erf.plot_int_1 = 10 # number of timesteps between plotfiles -erf.plot_vars_1 = density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres +erf.plot_vars_1 = lat_m lon_m density dens_hse rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta z_phys mapfac pres_hse pert_pres # SOLVER CHOICE erf.alpha_T = 1.0 diff --git a/Source/ERF.H b/Source/ERF.H index e299cb61c..89181ae3c 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -882,6 +882,15 @@ private: // used with init_type == "input_sounding" static bool init_sounding_ideal; + // Options for vertical interpolation of met_em*.nc data. + bool metgrid_basic_linear{false}; + bool metgrid_use_below_sfc{true}; + bool metgrid_use_sfc{true}; + bool metgrid_retain_sfc{false}; + amrex::Real metgrid_proximity{1000.0}; + int metgrid_order{2}; + int metgrid_force_sfc_k{0}; + // 1D CDF output (for ingestion in AMR-Wind) static int output_1d_column; static int column_interval; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index c0f1effb5..8b22a1310 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1358,6 +1358,15 @@ ERF::ReadParameters () // Flag to trigger initialization from input_sounding like WRF's ideal.exe pp.query("init_sounding_ideal", init_sounding_ideal); + // Options for vertical interpolation of met_em*.nc data. + pp.query("metgrid_basic_linear", metgrid_basic_linear); + pp.query("metgrid_use_below_sfc", metgrid_use_below_sfc); + pp.query("metgrid_use_sfc", metgrid_use_sfc); + pp.query("metgrid_retain_sfc", metgrid_retain_sfc); + pp.query("metgrid_proximity", metgrid_proximity); + pp.query("metgrid_order", metgrid_order); + pp.query("metgrid_force_sfc_k", metgrid_force_sfc_k); + // Output format pp.query("plotfile_type", plotfile_type); if (plotfile_type != "amrex" && diff --git a/Source/IO/NCWpsFile.H b/Source/IO/NCWpsFile.H index be8fc9577..00d0c41b4 100644 --- a/Source/IO/NCWpsFile.H +++ b/Source/IO/NCWpsFile.H @@ -17,7 +17,7 @@ using PlaneVector = amrex::Vector; NetCDF variables of dimensions Time_BT_SN_WE: "UU", "VV", "TT", "RH", "PRES", "GHT" NetCDF variables of dimensions Time_SN_WE : "HGT", "MAPFAC_U", "MAPFAC_V", "MAPFAC_M", "PSFC" NetCDF global attributes of type int : "FLAG_PSFC", "FLAG_MAPFAC_U", "FLAG_MAPFAC_V", "FLAG_MAPFAC_M", - "FLAG_HGT_M", "WEST-EAST_GRID_DIMENSION", "SOUTH-NORTH_GRID_DIMENSION" + "WEST-EAST_GRID_DIMENSION", "SOUTH-NORTH_GRID_DIMENSION" NetCDF global attributes of type string : "SIMULATION_START_DATE" NetCDF global attributes of type real : "DX", "DY" diff --git a/Source/IO/ReadFromMetgrid.cpp b/Source/IO/ReadFromMetgrid.cpp index 7c66b5b9f..4404b6480 100644 --- a/Source/IO/ReadFromMetgrid.cpp +++ b/Source/IO/ReadFromMetgrid.cpp @@ -10,8 +10,8 @@ using namespace amrex; void read_from_metgrid (int lev, const Box& domain, const std::string& fname, std::string& NC_dateTime, Real& NC_epochTime, - int& flag_psfc, int& flag_msfu, int& flag_msfv, int& flag_msfm, - int& flag_hgt, int& flag_sst, int& flag_lmask, + int& flag_psfc, int& flag_msf, + int& flag_sst, int& flag_lmask, int& NC_nx, int& NC_ny, Real& NC_dx, Real& NC_dy, FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab, @@ -32,24 +32,32 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); { // Global Attributes (int) std::vector attr; - ncf.get_attr("FLAG_PSFC", attr); flag_psfc = attr[0]; - /* UNCOMMENT FOR FLOWMAS - flag_msfu = 0; //ncf.get_attr("FLAG_MAPFAC_U", attr); flag_msfu = attr[0]; - flag_msfv = 0; //ncf.get_attr("FLAG_MAPFAC_V", attr); flag_msfv = attr[0]; - flag_msfm = 0; //ncf.get_attr("FLAG_MAPFAC_M", attr); flag_msfm = attr[0]; - flag_hgt = 1; //ncf.get_attr("FLAG_HGT_M", attr); flag_hgt = attr[0]; - flag_sst = 0; //ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0]; + flag_psfc = 0; //ncf.get_attr("FLAG_PSFC", attr); flag_psfc = attr[0]; + flag_msf = 0; //ncf.get_attr("FLAG_MF_XY", attr); flag_msf = attr[0]; + flag_sst = 0; //ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0]; flag_lmask = 0; //ncf.get_attr("FLAG_LANDMASK", attr); flag_lmask = attr[0]; */ - - ncf.get_attr("FLAG_MAPFAC_U", attr); flag_msfu = attr[0]; - ncf.get_attr("FLAG_MAPFAC_V", attr); flag_msfv = attr[0]; - ncf.get_attr("FLAG_MAPFAC_M", attr); flag_msfm = attr[0]; - ncf.get_attr("FLAG_HGT_M", attr); flag_hgt = attr[0]; - flag_sst = 0; //ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0]; - ncf.get_attr("FLAG_LANDMASK", attr); flag_lmask = attr[0]; - + if (ncf.has_attr("FLAG_PSFC")) { + ncf.get_attr("FLAG_PSFC", attr); flag_psfc = attr[0]; + } else { + flag_psfc = 0; + } + if (ncf.has_attr("FLAG_MF_XY")) { + ncf.get_attr("FLAG_MF_XY", attr); flag_msf = attr[0]; + } else { + flag_msf = 0; + } + if (ncf.has_attr("FLAG_SST")) { + ncf.get_attr("FLAG_SST", attr); flag_sst = attr[0]; + } else { + flag_sst = 0; + } + if (ncf.has_attr("FLAG_LANDMASK")) { + ncf.get_attr("FLAG_LANDMASK", attr); flag_lmask = attr[0]; + } else { + flag_lmask = 0; + } ncf.get_attr("WEST-EAST_GRID_DIMENSION", attr); NC_nx = attr[0]; ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); NC_ny = attr[0]; @@ -85,10 +93,7 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, } int ioproc = ParallelDescriptor::IOProcessorNumber(); // I/O rank ParallelDescriptor::Bcast(&flag_psfc, 1, ioproc); - ParallelDescriptor::Bcast(&flag_msfu, 1, ioproc); - ParallelDescriptor::Bcast(&flag_msfv, 1, ioproc); - ParallelDescriptor::Bcast(&flag_msfm, 1, ioproc); - ParallelDescriptor::Bcast(&flag_hgt, 1, ioproc); + ParallelDescriptor::Bcast(&flag_msf, 1, ioproc); ParallelDescriptor::Bcast(&flag_sst, 1, ioproc); ParallelDescriptor::Bcast(&flag_lmask, 1, ioproc); ParallelDescriptor::Bcast(&NC_nx, 1, ioproc); @@ -112,21 +117,21 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, NC_fabs.push_back(&NC_rhum_fab); NC_fnames.push_back("RH"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); NC_fabs.push_back(&NC_pres_fab); NC_fnames.push_back("PRES"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); NC_fabs.push_back(&NC_ght_fab); NC_fnames.push_back("GHT"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); - NC_fabs.push_back(&NC_LAT_fab); NC_fnames.push_back("XLAT_V"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); - NC_fabs.push_back(&NC_LON_fab); NC_fnames.push_back("XLONG_U"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); + NC_fabs.push_back(&NC_LAT_fab); NC_fnames.push_back("XLAT_M"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); + NC_fabs.push_back(&NC_LON_fab); NC_fnames.push_back("XLONG_M"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); + NC_fabs.push_back(&NC_hgt_fab); NC_fnames.push_back("HGT_M"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); if (flag_psfc) { NC_fabs.push_back(&NC_psfc_fab); NC_fnames.push_back("PSFC"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } - if (flag_msfu) { NC_fabs.push_back(&NC_msfu_fab); NC_fnames.push_back("MAPFAC_U"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } - if (flag_msfv) { NC_fabs.push_back(&NC_msfv_fab); NC_fnames.push_back("MAPFAC_V"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } - if (flag_msfm) { NC_fabs.push_back(&NC_msfm_fab); NC_fnames.push_back("MAPFAC_M"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } - if (flag_hgt) { NC_fabs.push_back(&NC_hgt_fab); NC_fnames.push_back("HGT_M"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } + if (flag_msf) { NC_fabs.push_back(&NC_msfu_fab); NC_fnames.push_back("MAPFAC_U"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } + if (flag_msf) { NC_fabs.push_back(&NC_msfv_fab); NC_fnames.push_back("MAPFAC_V"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } + if (flag_msf) { NC_fabs.push_back(&NC_msfm_fab); NC_fnames.push_back("MAPFAC_M"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } if (flag_sst) { NC_fabs.push_back(&NC_sst_fab); NC_fnames.push_back("SST"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } if (flag_lmask) { NC_iabs.push_back(&NC_lmask_iab); NC_inames.push_back("LANDMASK"); NC_idim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); } // Read the netcdf file and fill these FABs - std::string Lat_var_name = "XLAT_V"; - std::string Lon_var_name = "XLONG_U"; + std::string Lat_var_name = "XLAT_M"; + std::string Lon_var_name = "XLONG_M"; Print() << "Building initial FABS from file " << fname << std::endl; BuildFABsFromNetCDFFile(domain, Latitude, Longitude, Lat_var_name, Lon_var_name, diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 71dcc1748..9c8e7e68d 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -16,19 +16,23 @@ void ERF::init_from_metgrid (int lev) { bool use_moisture = (solverChoice.moisture_type != MoistureType::None); +#ifndef AMREX_USE_GPU if (use_moisture) { Print() << "Init with met_em with valid moisture model." << std::endl; } else { Print() << "Init with met_em without moisture model." << std::endl; } +#endif + + bool interp_theta = false; int ntimes = num_files_at_level[lev]; if (nc_init_file.empty()) - amrex::Error("NetCDF initialization file name must be provided via input"); + Error("NetCDF initialization file name must be provided via input"); if (nc_init_file[lev].empty()) - amrex::Error("NetCDF initialization file name must be provided via input"); + Error("NetCDF initialization file name must be provided via input"); // At least two met_em files are necessary to calculate tendency terms. AMREX_ALWAYS_ASSERT(ntimes >= 2); @@ -58,10 +62,7 @@ ERF::init_from_metgrid (int lev) // *** Variables at this level for holding metgrid file global attributes Vector flag_psfc; flag_psfc.resize( ntimes); - Vector flag_msfu; flag_msfu.resize( ntimes); - Vector flag_msfv; flag_msfv.resize( ntimes); - Vector flag_msfm; flag_msfm.resize( ntimes); - Vector flag_hgt; flag_hgt.resize( ntimes); + Vector flag_msf; flag_msf.resize( ntimes); Vector flag_sst; flag_sst.resize( ntimes); Vector flag_lmask; flag_lmask.resize( ntimes); Vector NC_nx; NC_nx.resize( ntimes); @@ -79,10 +80,13 @@ ERF::init_from_metgrid (int lev) #endif for (int it = 0; it < ntimes; it++) { +#ifndef AMREX_USE_GPU + Print() << " init_from_metgrid: reading nc_init_file[" << lev << "][" << it << "]\t" << nc_init_file[lev][it] << std::endl; +#endif read_from_metgrid(lev, boxes_at_level[lev][0], nc_init_file[lev][it], NC_dateTime[it], NC_epochTime[it], - flag_psfc[it], flag_msfu[it], flag_msfv[it], flag_msfm[it], - flag_hgt[it], flag_sst[it], flag_lmask[it], + flag_psfc[it], flag_msf[it], + flag_sst[it], flag_lmask[it], NC_nx[it], NC_ny[it], NC_dx[it], NC_dy[it], NC_xvel_fab[it], NC_yvel_fab[it], NC_temp_fab[it], NC_rhum_fab[it], NC_pres_fab[it], @@ -106,7 +110,10 @@ ERF::init_from_metgrid (int lev) // Verify that met_em files have even spacing in time. for (int it = 1; it < ntimes; it++) { Real NC_dt = NC_epochTime[it]-NC_epochTime[it-1]; - if (NC_dt != bdy_time_interval) amrex::Error("Time interval between consecutive met_em files must be consistent."); +#ifndef AMREX_USE_GPU + Print() << " " << nc_init_file[lev][it-1] << " / " << nc_init_file[lev][it] << " are " << NC_dt << " seconds apart" << std::endl; +#endif + if (NC_dt != bdy_time_interval) Error("Time interval between consecutive met_em files must be consistent."); } // Set up a FAB for mixing ratio and another for potential temperature. @@ -126,9 +133,6 @@ ERF::init_from_metgrid (int lev) AMREX_ALWAYS_ASSERT(solverChoice.use_terrain); - // Verify that the terrain height (HGT_M) was in each met_em file. - for (int it = 0; it < ntimes; it++) AMREX_ALWAYS_ASSERT(flag_hgt[it] == 1); - z_phys->setVal(0.0); for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { @@ -162,8 +166,8 @@ ERF::init_from_metgrid (int lev) const Array4& src_arr = src.const_array(); ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - int li = amrex::min(amrex::max(i, i_lo), i_hi); - int lj = amrex::min(amrex::max(j, j_lo), j_hi); + int li = min(max(i, i_lo), i_hi); + int lj = min(max(j, j_lo), j_hi); dst_arr(i,j,0) = src_arr(li,lj,0); }); } @@ -183,8 +187,8 @@ ERF::init_from_metgrid (int lev) const Array4& src_arr = src.const_array(); ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - int li = amrex::min(amrex::max(i, i_lo), i_hi); - int lj = amrex::min(amrex::max(j, j_lo), j_hi); + int li = min(max(i, i_lo), i_hi); + int lj = min(max(j, j_lo), j_hi); dst_arr(i,j,0) = src_arr(li,lj,0); }); } @@ -202,8 +206,8 @@ ERF::init_from_metgrid (int lev) const Array4& src_arr = src.const_array(); ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - int li = amrex::min(amrex::max(i, i_lo), i_hi); - int lj = amrex::min(amrex::max(j, j_lo), j_hi); + int li = min(max(i, i_lo), i_hi); + int lj = min(max(j, j_lo), j_hi); dst_arr(i,j,0) = src_arr(li,lj,0); }); } @@ -216,8 +220,8 @@ ERF::init_from_metgrid (int lev) const Array4& src_arr = src.const_array(); ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - int li = amrex::min(amrex::max(i, i_lo), i_hi); - int lj = amrex::min(amrex::max(j, j_lo), j_hi); + int li = min(max(i, i_lo), i_hi); + int lj = min(max(j, j_lo), j_hi); dst_arr(i,j,0) = src_arr(li,lj,0); }); } @@ -244,7 +248,7 @@ ERF::init_from_metgrid (int lev) if (use_moisture) MetGridBdyEnd = MetGridBdyVars::NumTypes; // Zero out fabs_for_bcs on the global domain - amrex::Vector> fabs_for_bcs; + Vector> fabs_for_bcs; fabs_for_bcs.resize(ntimes); for (int it(0); it < ntimes; it++) { fabs_for_bcs[it].resize(MetGridBdyEnd); @@ -256,22 +260,37 @@ ERF::init_from_metgrid (int lev) ldomain = geom[lev].Domain(); ldomain.convert(IntVect(1,0,0)); auto ng = lev_new[Vars::xvel].nGrowVect(); - gdomain = amrex::grow(ldomain,ng); + gdomain = grow(ldomain,ng); } else if (nvar==MetGridBdyVars::V) { ldomain = geom[lev].Domain(); ldomain.convert(IntVect(0,1,0)); auto ng = lev_new[Vars::yvel].nGrowVect(); - gdomain = amrex::grow(ldomain,ng); + gdomain = grow(ldomain,ng); } else { ldomain = geom[lev].Domain(); auto ng = lev_new[Vars::cons].nGrowVect(); - gdomain = amrex::grow(ldomain,ng); + gdomain = grow(ldomain,ng); } fabs_for_bcs[it][nvar].resize(gdomain, 1, Arena_Used); fabs_for_bcs[it][nvar].template setVal(0.0); } } + // Set up FABs to hold interpolated origin variables that are only needed during initialization. + // We will recalculate pressure from z, qv, and theta, but the origin model pressure will be + // necessary if we interpolate temperature (instead of theta) because we'll need pressure to + // calculate theta. + Vector p_interp_fab; + Vector t_interp_fab; + p_interp_fab.resize(ntimes); + t_interp_fab.resize(ntimes); + for (int it(0); it < ntimes; it++) { + Box ldomain = geom[lev].Domain(); + p_interp_fab[it].resize(ldomain, 1, Arena_Used); + p_interp_fab[it].template setVal(0.0); + t_interp_fab[it].resize(ldomain, 1, Arena_Used); + t_interp_fab[it].template setVal(0.0); + } const Real l_rdOcp = solverChoice.rdOcp; std::unique_ptr mask_c = OwnerMask(lev_new[Vars::cons], geom[lev].periodicity());//, lev_new[Vars::cons].nGrowVect()); @@ -302,13 +321,17 @@ ERF::init_from_metgrid (int lev) // z_vel set to 0.0 // theta calculate on origin levels then interpolate // mxrat convert RH -> Q on origin levels then interpolate - init_state_from_metgrid(use_moisture, l_rdOcp, + init_state_from_metgrid(use_moisture, interp_theta, metgrid_basic_linear, + metgrid_use_below_sfc, metgrid_use_sfc, + metgrid_retain_sfc, metgrid_proximity, + metgrid_order, metgrid_force_sfc_k, l_rdOcp, tbxc, tbxu, tbxv, cons_fab, xvel_fab, yvel_fab, zvel_fab, z_phys_cc_fab, NC_hgt_fab, NC_ght_fab, NC_xvel_fab, NC_yvel_fab, NC_temp_fab, NC_rhum_fab, - NC_pres_fab, theta_fab, mxrat_fab, + NC_pres_fab, p_interp_fab, t_interp_fab, + theta_fab, mxrat_fab, fabs_for_bcs, mask_c_arr, mask_u_arr, mask_v_arr); } // mf @@ -323,8 +346,7 @@ ERF::init_from_metgrid (int lev) FArrayBox &msfv_fab = (*mapfac_v[lev])[mfi]; FArrayBox &msfm_fab = (*mapfac_m[lev])[mfi]; - init_msfs_from_metgrid(msfu_fab, msfv_fab, msfm_fab, - flag_msfu[0], flag_msfv[0], flag_msfm[0], + init_msfs_from_metgrid(msfu_fab, msfv_fab, msfm_fab, flag_msf[0], NC_MSFU_fab, NC_MSFV_fab, NC_MSFM_fab); } // mf @@ -342,8 +364,8 @@ ERF::init_from_metgrid (int lev) const Array4& mask_c_arr = mask_c->const_array(mfi); // Fill base state data using origin data (initialization and BC arrays) - // p_hse calculate dry pressure - // r_hse calculate dry density + // p_hse calculate moist hydrostatic pressure + // r_hse calculate moist hydrostatic density // pi_hse calculate Exner term given pressure const Box valid_bx = mfi.validbox(); init_base_state_from_metgrid(use_moisture, l_rdOcp, @@ -367,17 +389,19 @@ ERF::init_from_metgrid (int lev) // complete data set; initialized to 0 above. for (int it(0); it < ntimes; it++) { for (int nvar(0); nvar& R_bcs_arr = fabs_for_bcs[it][MetGridBdyVars::R].const_array(); - for (int ivar(MetGridBdyVars::U); ivar < MetGridBdyEnd; ivar++) { auto xlo_arr = bdy_data_xlo[it][ivar].array(); @@ -473,23 +494,18 @@ ERF::init_from_metgrid (int lev) const Array4& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); if (ivar == MetGridBdyVars::U) { - multiply_rho = false; xlo_plane = xlo_plane_x_stag; xhi_plane = xhi_plane_x_stag; ylo_plane = ylo_plane_x_stag; yhi_plane = yhi_plane_x_stag; } else if (ivar == MetGridBdyVars::V) { - multiply_rho = false; xlo_plane = xlo_plane_y_stag; xhi_plane = xhi_plane_y_stag; ylo_plane = ylo_plane_y_stag; yhi_plane = yhi_plane_y_stag; } else if (ivar == MetGridBdyVars::R) { - multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; } else if (ivar == MetGridBdyVars::T) { - multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; } else if (ivar == MetGridBdyVars::QV) { - multiply_rho = false; xlo_plane = xlo_plane_no_stag; xhi_plane = xhi_plane_no_stag; ylo_plane = ylo_plane_no_stag; yhi_plane = yhi_plane_no_stag; } // MetGridBdyVars::QV @@ -497,26 +513,22 @@ ERF::init_from_metgrid (int lev) // west boundary ParallelFor(xlo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; - xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + xlo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); // xvel at east boundary ParallelFor(xhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; - xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + xhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); // xvel at south boundary ParallelFor(ylo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; - ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + ylo_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); // xvel at north boundary ParallelFor(yhi_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real Factor = (multiply_rho) ? R_bcs_arr(i,j,k) : 1.0; - yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k)*Factor; + yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); }); } // ivar @@ -568,6 +580,15 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, /** * Helper function to initialize state and velocity data read from metgrid data. * + * @param use_moisture bool True if solverChoice.moisture_type != MoistureType::None + * @param interp_theta bool calculate theta on origin levels, then interpolate + * @param metgrid_basic_linear bool linear interpolation without quality control + * @param metgrid_use_below_sfc bool quality control includes points below the surface + * @param metgrid_use_sfc bool quality control includes the point at the surface + * @param metgrid_retain_sfc bool set the lowest level of interpolated field to the surface value + * @param metgrid_proximity Real pressure difference for quality control pruning + * @param metgrid_order int interpolation order + * @param metgrid_force_sfc_k int lower levels pruned by quality control * @param l_rdOcp Real constant specifying Rhydberg constant ($R_d$) divided by specific heat at constant pressure ($c_p$) * @param state_fab FArrayBox holding the state data to initialize * @param x_vel_fab FArrayBox holding the x-velocity data to initialize @@ -582,12 +603,25 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, * @param NC_temp_fab Vector of FArrayBox obects holding metgrid data for temperature * @param NC_rhum_fab Vector of FArrayBox obects holding metgrid data for relative humidity * @param NC_pres_fab Vector of FArrayBox obects holding metgrid data for pressure + * @param p_interp_fab Vector of FArrayBox objects + * @param t_interp_fab Vector of FArrayBox objects * @param theta_fab Vector of FArrayBox obects holding potential temperature calculated from temperature and pressure * @param mxrat_fab Vector of FArrayBox obects holding vapor mixing ratio calculated from relative humidity * @param fabs_for_bcs Vector of Vector of FArrayBox objects holding MetGridBdyVars at each met_em time. + * @param mask_c_arr + * @param mask_u_arr + * @param mask_v_arr */ void init_state_from_metgrid (const bool use_moisture, + const bool interp_theta, + const bool metgrid_basic_linear, + const bool metgrid_use_below_sfc, + const bool metgrid_use_sfc, + const bool metgrid_retain_sfc, + const Real metgrid_proximity, + const int metgrid_order, + const int metgrid_force_sfc_k, const Real l_rdOcp, Box& tbxc, Box& tbxu, @@ -604,20 +638,29 @@ init_state_from_metgrid (const bool use_moisture, const Vector& NC_temp_fab, const Vector& NC_rhum_fab, const Vector& NC_pres_fab, + Vector& p_interp_fab, + Vector& t_interp_fab, Vector& theta_fab, Vector& mxrat_fab, - amrex::Vector>& fabs_for_bcs, - const amrex::Array4& mask_c_arr, - const amrex::Array4& mask_u_arr, - const amrex::Array4& mask_v_arr) + Vector>& fabs_for_bcs, + const Array4& mask_c_arr, + const Array4& mask_u_arr, + const Array4& mask_v_arr) { + bool metgrid_exp_interp = false; // interpolate w.r.t. exp(z) for non-pressure variables. + + // Loop over each time in the origin data. int ntimes = NC_hgt_fab.size(); for (int it = 0; it < ntimes; it++) { + // ******************************************************** // U // ******************************************************** { +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of u-velocity, it = " << it << std::endl; +#endif Box bx2d = NC_xvel_fab[it].box() & tbxu; bx2d.setRange(2,0); auto const orig_data = NC_xvel_fab[it].const_array(); @@ -626,14 +669,22 @@ init_state_from_metgrid (const bool use_moisture, auto bc_data = fabs_for_bcs[it][MetGridBdyVars::U].array(); auto const new_z = z_phys_nd_fab.const_array(); - int kmax = amrex::ubound(tbxu).z; + int kmax = ubound(tbxu).z; ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - for (int k = 0; k<=kmax; k++) { - Real Interp_Val = interpolate_column_metgrid(i,j,k,'X',0,orig_z,orig_data,new_z); - if (mask_u_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; - if (it==0) new_data(i,j,k,0) = Interp_Val; + if (metgrid_basic_linear) { + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid_linear(i,j,k,'X',0,orig_z,orig_data,new_z); + if (mask_u_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; + if (it==0) new_data(i,j,k,0) = Interp_Val; + } + } else { + interpolate_column_metgrid(metgrid_use_below_sfc, metgrid_use_sfc, metgrid_exp_interp, + metgrid_retain_sfc, metgrid_proximity, metgrid_order, + metgrid_force_sfc_k, i, j, 0, it, 'U', 'X', + orig_z, orig_data, new_z, new_data, + true, bc_data, mask_u_arr); } }); } @@ -643,6 +694,9 @@ init_state_from_metgrid (const bool use_moisture, // V // ******************************************************** { +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of v-velocity, it = " << it << std::endl; +#endif Box bx2d = NC_yvel_fab[it].box() & tbxv; bx2d.setRange(2,0); auto const orig_data = NC_yvel_fab[it].const_array(); @@ -651,14 +705,22 @@ init_state_from_metgrid (const bool use_moisture, auto bc_data = fabs_for_bcs[it][MetGridBdyVars::V].array(); auto const new_z = z_phys_nd_fab.const_array(); - int kmax = amrex::ubound(tbxv).z; + int kmax = ubound(tbxv).z; ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - for (int k = 0; k<=kmax; k++) { - Real Interp_Val = interpolate_column_metgrid(i,j,k,'Y',0,orig_z,orig_data,new_z); - if (mask_v_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; - if (it==0) new_data(i,j,k,0) = Interp_Val; + if (metgrid_basic_linear) { + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid_linear(i,j,k,'Y',0,orig_z,orig_data,new_z); + if (mask_v_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; + if (it==0) new_data(i,j,k,0) = Interp_Val; + } + } else { + interpolate_column_metgrid(metgrid_use_below_sfc, metgrid_use_sfc, metgrid_exp_interp, + metgrid_retain_sfc, metgrid_proximity, metgrid_order, + metgrid_force_sfc_k, i, j, 0, it, 'V', 'Y', + orig_z, orig_data, new_z, new_data, + true, bc_data, mask_v_arr); } }); } @@ -682,40 +744,140 @@ init_state_from_metgrid (const bool use_moisture, // ******************************************************** // theta // ******************************************************** - { // calculate potential temperature. - Box bx = NC_rhum_fab[it].box() & tbxc; - auto const temp = NC_temp_fab[it].const_array(); - auto const pres = NC_pres_fab[it].const_array(); - auto theta = theta_fab[it].array(); + if (interp_theta) { + // Calculate potential temperature on the origin model vertical levels + // then interpolate that onto the ERF vertical levels. - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { // calculate potential temperature. + Box bx = NC_rhum_fab[it].box() & tbxc; + auto const temp = NC_temp_fab[it].const_array(); + auto const pres = NC_pres_fab[it].const_array(); + auto theta = theta_fab[it].array(); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + theta(i,j,k) = getThgivenPandT(temp(i,j,k),pres(i,j,k),l_rdOcp); + //theta(i,j,k) = 300.0; // uncomment for an isothermal atmosphere. + }); + } + + { // vertical interpolation of potential temperature. +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of potential temperature, it = " << it << std::endl; +#endif + Box bx2d = NC_temp_fab[it].box() & tbxc; + bx2d.setRange(2,0); + auto const orig_data = theta_fab[it].const_array(); + auto const orig_z = NC_ght_fab[it].const_array(); + auto new_data = state_fab.array(); + auto bc_data = fabs_for_bcs[it][MetGridBdyVars::T].array(); + auto const new_z = z_phys_nd_fab.const_array(); + + int kmax = amrex::ubound(tbxc).z; + + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - theta(i,j,k) = getThgivenPandT(temp(i,j,k),pres(i,j,k),l_rdOcp); - //theta(i,j,k) = 300.0; // TODO: Remove when not needed. Force an isothermal atmosphere for debugging. + if (metgrid_basic_linear) { + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid_linear(i,j,k,'M',0,orig_z,orig_data,new_z); + if (mask_c_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; + if (it==0) new_data(i,j,k,RhoTheta_comp) = Interp_Val; + } + } else { + interpolate_column_metgrid(metgrid_use_below_sfc, metgrid_use_sfc, metgrid_exp_interp, + metgrid_retain_sfc, metgrid_proximity, metgrid_order, + metgrid_force_sfc_k, i, j, RhoTheta_comp, it, 'T', 'M', + orig_z, orig_data, new_z, new_data, + true, bc_data, mask_c_arr); + } }); - } + } - // vertical interpolation of potential temperature. - { - Box bx2d = NC_temp_fab[it].box() & tbxc; - bx2d.setRange(2,0); - auto const orig_data = theta_fab[it].const_array(); - auto const orig_z = NC_ght_fab[it].const_array(); - auto new_data = state_fab.array(); - auto bc_data = fabs_for_bcs[it][MetGridBdyVars::T].array(); - auto const new_z = z_phys_nd_fab.const_array(); + } else { // interp_theta == false - int kmax = amrex::ubound(tbxc).z; + { // vertical interpolation of pressure. +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of pressure, it = " << it << std::endl; +#endif + Box bx2d = p_interp_fab[it].box() & tbxc; + bx2d.setRange(2,0); + auto const orig_data = NC_pres_fab[it].const_array(); + auto const orig_z = NC_ght_fab[it].const_array(); + auto new_data = p_interp_fab[it].array(); + auto const new_z = z_phys_nd_fab.const_array(); + const amrex::Array4 bc_data_unused; - ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept - { - for (int k = 0; k<=kmax; k++) { - Real Interp_Val = interpolate_column_metgrid(i,j,k,'M',0,orig_z,orig_data,new_z); - if (mask_c_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; - if (it==0) new_data(i,j,k,RhoTheta_comp) = Interp_Val; + int kmax = ubound(tbxc).z; + + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + if (metgrid_basic_linear) { + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid_linear(i,j,k,'M',0,orig_z,orig_data,new_z); + new_data(i,j,k) = Interp_Val; + } + } else { + // Interpolate pressure not w.r.t. z but rather p_0*exp(-CONST_GRAV*z/(t_0*R_d)). + // This is akin to interpolating in pressure-space assuming a baroclinic atmosphere. + interpolate_column_metgrid(metgrid_use_below_sfc, metgrid_use_sfc, true, + metgrid_retain_sfc, metgrid_proximity, metgrid_order, + metgrid_force_sfc_k, i, j, 0, 0, 'T', 'M', + orig_z, orig_data, new_z, new_data, + false, bc_data_unused, mask_c_arr); + } + }); } - }); - } + + { // vertical interpolation of temperature. +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of temperature, it = " << it << std::endl; +#endif + Box bx2d = p_interp_fab[it].box() & tbxc; + bx2d.setRange(2,0); + auto const orig_data = NC_temp_fab[it].const_array(); + auto const orig_z = NC_ght_fab[it].const_array(); + auto new_data = t_interp_fab[it].array(); + auto const new_z = z_phys_nd_fab.const_array(); + const amrex::Array4 bc_data_unused; + + int kmax = ubound(tbxc).z; + + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + if (metgrid_basic_linear) { + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid_linear(i,j,k,'M',0,orig_z,orig_data,new_z); + new_data(i,j,k) = Interp_Val; + } + } else { + // According to WRF's code comments, "It is better to + // interpolate temperature and potential temperature + // in LOG(p), regardless of requested default." + interpolate_column_metgrid(metgrid_use_below_sfc, metgrid_use_sfc, false, + metgrid_retain_sfc, metgrid_proximity, metgrid_order, + metgrid_force_sfc_k, i, j, 0, 0, 'T', 'M', + orig_z, orig_data, new_z, new_data, + false, bc_data_unused, mask_c_arr); + } + }); + } + + { // calculate potential temperature on the ERF vertical levels. + auto const temp = t_interp_fab[it].const_array(); + auto const pres = p_interp_fab[it].const_array(); + auto new_data = state_fab.array(); + auto bc_data = fabs_for_bcs[it][MetGridBdyVars::T].array(); + + ParallelFor(tbxc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real Calc_Val = getThgivenPandT(temp(i,j,k),pres(i,j,k),l_rdOcp); + //Calc_Val = 300.0; // uncomment for an isothermal atmosphere. + if (mask_c_arr(i,j,k)) bc_data(i,j,k,0) = Calc_Val; + if (it==0) new_data(i,j,k,RhoTheta_comp) = Calc_Val; + }); + } + + } // interp_theta if (use_moisture) { // ******************************************************** @@ -739,8 +901,10 @@ init_state_from_metgrid (const bool use_moisture, }); } - // vertical interpolation of vapor mixing ratio. - { + { // vertical interpolation of vapor mixing ratio. +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of vapor mixing ratio, it = " << it << std::endl; +#endif Box bx2d = NC_temp_fab[it].box() & tbxc; bx2d.setRange(2,0); auto const orig_data = mxrat_fab[it].const_array(); @@ -749,15 +913,23 @@ init_state_from_metgrid (const bool use_moisture, auto bc_data = fabs_for_bcs[it][MetGridBdyVars::QV].array(); auto const new_z = z_phys_nd_fab.const_array(); - int kmax = amrex::ubound(tbxc).z; + int kmax = ubound(tbxc).z; int state_indx = RhoQ1_comp; ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept { - for (int k = 0; k<=kmax; k++) { - Real Interp_Val = interpolate_column_metgrid(i,j,k,'M',0,orig_z,orig_data,new_z); - if (mask_c_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; - if (it==0) new_data(i,j,k,state_indx) = Interp_Val; + if (metgrid_basic_linear) { + for (int k = 0; k<=kmax; k++) { + Real Interp_Val = interpolate_column_metgrid_linear(i,j,k,'M',0,orig_z,orig_data,new_z); + if (mask_c_arr(i,j,k)) bc_data(i,j,k,0) = Interp_Val; + if (it==0) new_data(i,j,k,state_indx) = Interp_Val; + } + } else { + interpolate_column_metgrid(metgrid_use_below_sfc, metgrid_use_sfc, metgrid_exp_interp, + metgrid_retain_sfc, metgrid_proximity, metgrid_order, + metgrid_force_sfc_k, i, j, state_indx, it, 'Q', 'M', + orig_z, orig_data, new_z, new_data, + true, bc_data, mask_c_arr); } }); } @@ -778,6 +950,7 @@ init_state_from_metgrid (const bool use_moisture, /** * Helper function for initializing hydrostatic base state data from metgrid data * + * @param use_moisture bool True if solverChoice.moisture_type != MoistureType::None * @param l_rdOcp Real constant specifying Rhydberg constant ($R_d$) divided by specific heat at constant pressure ($c_p$) * @param valid_bx Box specifying the index space we are to initialize * @param flag_psfc Vector of Integer 1 if surface pressure is in metgrid data, 0 otherwise @@ -789,6 +962,7 @@ init_state_from_metgrid (const bool use_moisture, * @param NC_ght_fab Vector of FArrayBox objects holding metgrid data for height of cell centers * @param NC_psfc_fab Vector of FArrayBox objects holding metgrid data for surface pressure * @param fabs_for_bcs Vector of Vector of FArrayBox objects holding MetGridBdyVars at each met_em time. + * @param mask_c_arr */ void init_base_state_from_metgrid (const bool use_moisture, @@ -800,13 +974,13 @@ init_base_state_from_metgrid (const bool use_moisture, FArrayBox& p_hse_fab, FArrayBox& pi_hse_fab, FArrayBox& z_phys_cc_fab, - const Vector& /*NC_ght_fab*/, + const Vector& NC_ght_fab, const Vector& NC_psfc_fab, Vector>& fabs_for_bcs, - const amrex::Array4& mask_c_arr) + const Array4& mask_c_arr) { int RhoQ_comp = RhoQ1_comp; - int kmax = amrex::ubound(valid_bx).z; + int kmax = ubound(valid_bx).z; // NOTE: FOEXTRAP is utilized on the validbox but // the FillBoundary call will populate the @@ -828,9 +1002,7 @@ init_base_state_from_metgrid (const bool use_moisture, Gpu::DeviceVector z_vec_d(kmax+2,0); Real* z_vec = z_vec_d.data(); Gpu::DeviceVector Thetad_vec_d(kmax+1,0); Real* Thetad_vec = Thetad_vec_d.data(); Gpu::DeviceVector Thetam_vec_d(kmax+1,0); Real* Thetam_vec = Thetam_vec_d.data(); - Gpu::DeviceVector Rhod_vec_d(kmax+1,0); Real* Rhod_vec = Rhod_vec_d.data(); Gpu::DeviceVector Rhom_vec_d(kmax+1,0); Real* Rhom_vec = Rhom_vec_d.data(); - Gpu::DeviceVector Pd_vec_d(kmax+1,0); Real* Pd_vec = Pd_vec_d.data(); Gpu::DeviceVector Pm_vec_d(kmax+1,0); Real* Pm_vec = Pm_vec_d.data(); Gpu::DeviceVector Q_vec_d(kmax+1,0); Real* Q_vec = Q_vec_d.data(); @@ -852,9 +1024,8 @@ init_base_state_from_metgrid (const bool use_moisture, const Array4& pi_hse_arr = pi_hse_fab.array(); // ******************************************************** - // calculate dry density and dry pressure + // calculate density and pressure for initial conditions. // ******************************************************** - // calculate density and dry pressure on the new grid. Box valid_bx2d = valid_bx; valid_bx2d.setRange(2,0); auto const orig_psfc = NC_psfc_fab[0].const_array(); @@ -866,13 +1037,13 @@ init_base_state_from_metgrid (const bool use_moisture, for (int k=0; k<=kmax; k++) { z_vec[k] = new_z(i,j,k); Thetad_vec[k] = new_data(i,j,k,RhoTheta_comp); - Q_vec[k] = (use_moisture) ? new_data(i,j,k,RhoQ_comp) : 0.0; + Q_vec[k] = (use_moisture) ? new_data(i,j,k,RhoQ_comp) : 0.0; } - z_vec[kmax+1] = new_z(i,j,kmax+1); + z_vec[kmax+1] = new_z(i,j,kmax+1); calc_rho_p(kmax, flag_psfc_vec[0], orig_psfc(i,j,0), Thetad_vec, Thetam_vec, Q_vec, z_vec, - Rhod_vec, Rhom_vec, Pd_vec, Pm_vec); + Rhom_vec, Pm_vec); for (int k=0; k<=kmax; k++) { p_hse_arr(i,j,k) = Pm_vec[k]; @@ -882,7 +1053,7 @@ init_base_state_from_metgrid (const bool use_moisture, ParallelFor(valid_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - // Multiply by Rho to get conserved vars + // Multiply by Rho to get conserved vars Real Qv = 0.0; new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k); new_data(i,j,k,RhoTheta_comp) *= r_hse_arr(i,j,k); @@ -898,20 +1069,23 @@ init_base_state_from_metgrid (const bool use_moisture, pi_hse_arr(i,j,k) = getExnergivenP(p_hse_arr(i,j,k), l_rdOcp); }); - // FOEXTRAP hse arrays + // FOEXTRAP hse arrays to fill ghost cells. FillBoundary is + // called later and will overwrite ghost cell values in the + // interior of the domain, but the FOEXTRAP values will + // remain along the lateral domain boundaries. ParallelFor(gvbx_xlo, gvbx_xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int jj = amrex::max(j ,valid_bx.smallEnd(1)); - jj = amrex::min(jj,valid_bx.bigEnd(1)); + int jj = max(j ,valid_bx.smallEnd(1)); + jj = min(jj,valid_bx.bigEnd(1)); r_hse_arr(i,j,k) = r_hse_arr(i+1,jj,k); p_hse_arr(i,j,k) = p_hse_arr(i+1,jj,k); pi_hse_arr(i,j,k) = pi_hse_arr(i+1,jj,k); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - int jj = amrex::max(j ,valid_bx.smallEnd(1)); - jj = amrex::min(jj,valid_bx.bigEnd(1)); + int jj = max(j ,valid_bx.smallEnd(1)); + jj = min(jj,valid_bx.bigEnd(1)); r_hse_arr(i,j,k) = r_hse_arr(i-1,jj,k); p_hse_arr(i,j,k) = p_hse_arr(i-1,jj,k); pi_hse_arr(i,j,k) = pi_hse_arr(i-1,jj,k); @@ -951,9 +1125,8 @@ init_base_state_from_metgrid (const bool use_moisture, p_hse_bcs_fab.resize(state_fab.box(), 1, Arena_Used); // ******************************************************** - // calculate dry density and dry pressure + // calculate density and pressure for boundary conditions. // ******************************************************** - // calculate density and dry pressure on the new grid. Box valid_bx2d = valid_bx; valid_bx2d.setRange(2,0); auto const orig_psfc = NC_psfc_fab[it].const_array(); @@ -974,8 +1147,9 @@ init_base_state_from_metgrid (const bool use_moisture, calc_rho_p(kmax, flag_psfc_vec[it], orig_psfc(i,j,0), Thetad_vec, Thetam_vec, Q_vec, z_vec, - Rhod_vec, Rhom_vec, Pd_vec, Pm_vec); + Rhom_vec, Pm_vec); + // Multiply by Rho to get conserved vars for (int k=0; k<=kmax; k++) { p_hse_arr(i,j,k) = Pm_vec[k]; if (mask_c_arr(i,j,k)) { @@ -995,9 +1169,7 @@ init_base_state_from_metgrid (const bool use_moisture, * @param msfu_fab FArrayBox specifying x-velocity map factors * @param msfv_fab FArrayBox specifying y-velocity map factors * @param msfm_fab FArrayBox specifying z-velocity map factors - * @param flag_msfu Integer 1 if u-staggered map factor is in metgrid data, 0 otherwise - * @param flag_msfv Integer 1 if v-staggered map factor is in metgrid data, 0 otherwise - * @param flag_msfm Integer 1 if cell center map factor is in metgrid data, 0 otherwise + * @param flag_msf Integer 1 if map factors are in metgrid data, 0 otherwise * @param NC_MSFU_fab Vector of FArrayBox objects holding metgrid data for x-velocity map factors * @param NC_MSFV_fab Vector of FArrayBox objects holding metgrid data for y-velocity map factors * @param NC_MSFM_fab Vector of FArrayBox objects holding metgrid data for z-velocity map factors @@ -1006,9 +1178,7 @@ void init_msfs_from_metgrid (FArrayBox& msfu_fab, FArrayBox& msfv_fab, FArrayBox& msfm_fab, - const int& flag_msfu, - const int& flag_msfv, - const int& flag_msfm, + const int& flag_msf, const Vector& NC_MSFU_fab, const Vector& NC_MSFV_fab, const Vector& NC_MSFM_fab) @@ -1022,35 +1192,18 @@ init_msfs_from_metgrid (FArrayBox& msfu_fab, // // This copies or sets mapfac_m - if (flag_msfm == 1) { + if (flag_msf == 1) { msfm_fab.template copy(NC_MSFM_fab[it]); - } else { -#ifndef AMREX_USE_GPU - Print() << " MAPFAC_M not present in met_em files. Setting to 1.0" << std::endl; -#endif - msfm_fab.template setVal(1.0); - } - - // This copies or sets mapfac_u - if (flag_msfu == 1) { msfu_fab.template copy(NC_MSFU_fab[it]); - } else { -#ifndef AMREX_USE_GPU - Print() << " MAPFAC_U not present in met_em files. Setting to 1.0" << std::endl; -#endif - msfu_fab.template setVal(1.0); - } - - // This copies or sets mapfac_v - if (flag_msfv == 1) { msfv_fab.template copy(NC_MSFV_fab[it]); } else { #ifndef AMREX_USE_GPU - Print() << " MAPFAC_V not present in met_em files. Setting to 1.0" << std::endl; + Print() << " map factors are not present in met_em files. Setting to 1.0" << std::endl; #endif + msfm_fab.template setVal(1.0); + msfu_fab.template setVal(1.0); msfv_fab.template setVal(1.0); } - } // it } #endif // ERF_USE_NETCDF diff --git a/Source/Initialization/Metgrid_utils.H b/Source/Initialization/Metgrid_utils.H index b3c38da7d..a17db4329 100644 --- a/Source/Initialization/Metgrid_utils.H +++ b/Source/Initialization/Metgrid_utils.H @@ -14,10 +14,7 @@ read_from_metgrid (int lev, std::string& NC_dateTime, amrex::Real& NC_epochTime, int& flag_psfc, - int& flag_msfu, - int& flag_msfv, - int& flag_msfm, - int& flag_hgt, + int& flag_msf, int& flag_sst, int& flag_lmask, int& NC_nx, @@ -51,6 +48,14 @@ init_terrain_from_metgrid (amrex::FArrayBox& z_phys_nd_fab, void init_state_from_metgrid (const bool use_moisture, + const bool interp_theta, + 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, @@ -64,9 +69,11 @@ init_state_from_metgrid (const bool use_moisture, const amrex::Vector& NC_ght_fab, const amrex::Vector& NC_xvel_fab, const amrex::Vector& NC_yvel_fab, - const amrex::Vector& NC_zvel_fab, const amrex::Vector& NC_temp_fab, const amrex::Vector& NC_rhum_fab, + const amrex::Vector& NC_pres_fab, + amrex::Vector& p_interp_fab, + amrex::Vector& t_interp_fab, amrex::Vector& theta_fab, amrex::Vector& mxrat_fab, amrex::Vector>& fabs_for_bcs, @@ -78,9 +85,7 @@ void init_msfs_from_metgrid (amrex::FArrayBox& msfu_fab, amrex::FArrayBox& msfv_fab, amrex::FArrayBox& msfm_fab, - const int& flag_msfu, - const int& flag_msfv, - const int& flag_msfm, + const int& flag_msf, const amrex::Vector& NC_MSFU_fab, const amrex::Vector& NC_MSFV_fab, const amrex::Vector& NC_MSFM_fab); @@ -110,12 +115,10 @@ calc_rho_p (const int& kmax, amrex::Real* Thetam_vec, amrex::Real* Q_vec, amrex::Real* z_vec, - amrex::Real* Rhod_vec, amrex::Real* Rhom_vec, - amrex::Real* Pd_vec, amrex::Real* Pm_vec) { - const int maxiter = 10; + const int maxiter = 50; // Calculate or use moist pressure at the surface. amrex::Real Psurf; @@ -132,8 +135,9 @@ calc_rho_p (const int& kmax, amrex::Real half_dz = z_vec[0]; amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[0]; Thetam_vec[0] = Thetad_vec[0]*qvf; + Rhom_vec[0] = 1.0; // an initial guess. for (int it=0; it=0; k--) { - amrex::Real dz = z_vec[k+1]-z_vec[k]; - Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess. - for (int it=0; it < maxiter; it++) { - Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV; - if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0; - Rhod_vec[k] = (p_0/(R_d*Thetam_vec[k]))*std::pow(Pd_vec[k]/p_0, iGamma); - } // it - } // k +AMREX_FORCE_INLINE +AMREX_GPU_DEVICE +void +lagrange_interp (amrex::Vector& x, + amrex::Vector& y, + const int& order, + amrex::Real& new_x, + amrex::Real& new_y) +{ + // Interpolation using Lagrange polynomials. + // P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x) + // where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn) + // --------------------------------------------- + // (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn) + amrex::Real Px = 0.0; + for (int i=0; i <= order; i++) { + amrex::Real numer = 1.0; + amrex::Real denom = 1.0; + for (int k=0; k <= order; k++) { + if (k == i) continue; + numer *= new_x-x[k]; + denom *= x[i]-x[k]; + } + if (denom != 0.0) { + Px += y[i]*numer/denom; + } + } + new_y = Px; } AMREX_FORCE_INLINE AMREX_GPU_DEVICE -amrex::Real -interpolate_column_metgrid (const int& i, +void +lagrange_setup (char var_type, + const bool& exp_interp, + amrex::Vector& orig_x_z, + amrex::Vector& orig_x_p, + amrex::Vector& orig_y, + const int& orig_n, + const int& i, + const int& j, + const int& new_n, + const int& order, + amrex::Vector& new_x_z, + amrex::Vector& new_x_p, + amrex::Vector& new_y) +{ + if (order < 1) amrex::Abort("metgrid initialization, we cannot go lower order than linear"); + + amrex::Real CRC_const1 = 11880.516; // m + amrex::Real CRC_const2 = 0.1902632; + amrex::Real CRC_const3 = 0.0065; // K km-1 + int vboundb = 4; + int vboundt = 0; + + bool debug = false; +// if ((i == 391) && (j == 3)) debug = true; + + for (int new_k=0; new_k < new_n; new_k++) { + if (debug) amrex::Print() << "new_k=" << new_k; + + // Find bounding x values and store the indices. + bool extrapolating = true; + int kl, kr; + for (int ko=0; ko < orig_n-1; ko++) { + amrex::Real a = new_x_z[new_k]-orig_x_z[ko]; + amrex::Real b = new_x_z[new_k]-orig_x_z[ko+1]; + if (a*b <= 0.0) { + kl = ko; + kr = ko+1; + extrapolating = false; + break; + } + } + + if (extrapolating) { + if (var_type == 'T') { + // Assume a standard atmosphere -6.5 K km-1 lapse rate. + // Comparable to the WRF default, t_extrap_type=2. + amrex::Real depth_of_extrap_in_p = new_x_p[new_k]-orig_x_p[1]; + amrex::Real avg_of_extrap_p = 0.5*(new_x_p[new_k]+orig_x_p[1]); + amrex::Real temp_extrap_starting_point = orig_y[1]*std::pow(orig_x_p[1]/100000.0, R_d/Cp_d); + amrex::Real dZdP = CRC_const1*CRC_const2*std::pow(avg_of_extrap_p/100.0, CRC_const2-1.0); + amrex::Real dZ = dZdP*(depth_of_extrap_in_p/100.0); + amrex::Real dT = dZ*CRC_const3; + new_y[new_k] = (temp_extrap_starting_point+dT)*std::pow(100000.0/new_x_p[new_k], R_d/Cp_d); + } else { + // Use a constant value below ground. + // Comparable to the WRF default, extrap_type=2. + new_y[new_k] = orig_y[1]; + } + continue; + } + + if (order%2 != 0) { + if ((kl-((order+1)/2-1) >= 0) && (kr+((order+1)/2-1) <= orig_n-1)) { + int ksta = kl-((order/2)-1); + int kend = ksta+order; + amrex::Real new_x; + amrex::Vector orig_x_sub; + if (exp_interp) { + new_x = new_x_p[new_k]; + orig_x_sub = {orig_x_p.begin()+ksta, orig_x_p.begin()+kend}; + } else { + new_x = new_x_z[new_k]; + orig_x_sub = {orig_x_z.begin()+ksta, orig_x_z.begin()+kend}; + } + amrex::Vector orig_y_sub = {orig_y.begin()+ksta, orig_y.begin()+kend}; + lagrange_interp(orig_x_sub, orig_y_sub, order, new_x, new_y[new_k]); + } else { + amrex::Abort("metgrid initialization, lost in lagrange_setup (odd order)"); + } + } else if ((order%2 == 0) && (new_k >= 1+vboundb) && (new_k < new_n-vboundt)) { + if ((kl-(order/2) >= 0) && (kr+order/2 <= orig_n-1)) { + amrex::Real new_y_l, new_y_r; + { + int ksta = kl-(order/2-1); + int kend = ksta+order; + if (debug) amrex::Print() << " (1a) 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; + amrex::Real new_x; + amrex::Vector orig_x_sub, orig_y_sub; + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_p[k]); + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_z[k]); + } + for (int k=ksta; k <= kend; k++) orig_y_sub.push_back(orig_y[k]); + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < orig_x_sub.size(); k++) amrex::Print() << " " << orig_x_sub[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < orig_y_sub.size(); k++) amrex::Print() << " " << orig_y_sub[k]; + amrex::Print() << "]" << std::endl; + } + lagrange_interp(orig_x_sub, orig_y_sub, order, new_x, new_y_l); + } + { + int ksta = kl-order/2; + int kend = ksta+order; + if (debug) amrex::Print() << "new_k=" << new_k << " (1b) 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; + amrex::Real new_x; + amrex::Vector orig_x_sub, orig_y_sub; + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_p[k]); + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_z[k]); + } + for (int k=ksta; k <= kend; k++) orig_y_sub.push_back(orig_y[k]); + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < orig_x_sub.size(); k++) amrex::Print() << " " << orig_x_sub[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < orig_y_sub.size(); k++) amrex::Print() << " " << orig_y_sub[k]; + amrex::Print() << "]" << std::endl; + } + lagrange_interp(orig_x_sub, orig_y_sub, order, new_x, new_y_r); + } + new_y[new_k] = 0.5*(new_y_l+new_y_r); + if (debug) amrex::Print() << " new_y=" << new_y[new_k] << std::endl; + } else if ((kl-(order/2-1) >= 0) && (kr+order/2 <= orig_n-1)) { + int ksta = kl-(order/2-1); + int kend = ksta+order; + if (debug) amrex::Print() << " (2) 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; + amrex::Real new_x; + amrex::Vector orig_x_sub, orig_y_sub; + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_p[k]); + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_z[k]); + } + for (int k=ksta; k <= kend; k++) orig_y_sub.push_back(orig_y[k]); + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < orig_x_sub.size(); k++) amrex::Print() << " " << orig_x_sub[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < orig_y_sub.size(); k++) amrex::Print() << " " << orig_y_sub[k]; + amrex::Print() << "]" << std::endl; + } + lagrange_interp(orig_x_sub, orig_y_sub, order, new_x, new_y[new_k]); + if (debug) amrex::Print() << " new_y=" << new_y[new_k] << std::endl; + } else if ((kl-order/2 >= 0) && (kr+(order/2-1) <= orig_n-1)) { + int ksta = kl-order/2; + int kend = ksta+order; + if (debug) amrex::Print() << " (3) 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; + amrex::Real new_x; + amrex::Vector orig_x_sub, orig_y_sub; + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_p[k]); + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_z[k]); + } + for (int k=ksta; k <= kend; k++) orig_y_sub.push_back(orig_y[k]); + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < orig_x_sub.size(); k++) amrex::Print() << " " << orig_x_sub[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < orig_y_sub.size(); k++) amrex::Print() << " " << orig_y_sub[k]; + amrex::Print() << "]" << std::endl; + } + lagrange_interp(orig_x_sub, orig_y_sub, order, new_x, new_y[new_k]); + if (debug) amrex::Print() << " new_y=" << new_y[new_k] << std::endl; + } else { + amrex::Abort("metgrid initialization, lost in lagrange_setup (even order)"); + } + } else { + // Linear interpolation. + int ksta = kl; + int kend = kr; + if (debug) amrex::Print() << " (4) 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; + amrex::Real new_x; + amrex::Vector orig_x_sub, orig_y_sub; + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_p[k]); + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) orig_x_sub.push_back(orig_x_z[k]); + } + for (int k=ksta; k <= kend; k++) orig_y_sub.push_back(orig_y[k]); + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < orig_x_sub.size(); k++) amrex::Print() << " " << orig_x_sub[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < orig_y_sub.size(); k++) amrex::Print() << " " << orig_y_sub[k]; + amrex::Print() << "]" << std::endl; + } + lagrange_interp(orig_x_sub, orig_y_sub, 1, new_x, new_y[new_k]); + if (debug) amrex::Print() << " new_y=" << new_y[new_k] << std::endl; + } + } +} + +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& k, + const int& src_comp, + const int& it, + char var_type, char stag, - int src_comp, - const amrex::Array4& orig_z, + const amrex::Array4& orig_z_full, const amrex::Array4& orig_data, - const amrex::Array4& new_z) + const amrex::Array4& new_z_full, + const amrex::Array4& new_data_full, + const bool& update_bc_data, + const amrex::Array4& bc_data_full, + const amrex::Array4& mask) +{ + // Here we closely follow WRF's vert_interp from + // dyn_em/module_initialize_real.F, although changes have been + // made to accomodate interpolation relative to height instead of + // pressure. + int imax_orig = amrex::ubound(amrex::Box(orig_data)).x; + int jmax_orig = amrex::ubound(amrex::Box(orig_data)).y; + int kmax_orig = amrex::ubound(amrex::Box(orig_data)).z; + int kmax_new = amrex::ubound(amrex::Box(new_z_full)).z; + + amrex::Real t_0 = 290.0; + amrex::Vector new_z, new_p, new_data; + new_z.resize(kmax_new); + new_p.resize(kmax_new); + new_data.resize(kmax_new); + for (int k=0; k < kmax_new; k++) { + if (stag == 'X') { + new_z[k] = 0.25*(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)); + } else if (stag == 'Y') { + new_z[k] = 0.25*(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)); + } else if (stag == 'M') { + new_z[k] = 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 )+ + 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)); + } + new_p[k] = p_0*exp(-CONST_GRAV*new_z[k]/(t_0*R_d)); + } + + amrex::Vector orig_z; + orig_z.resize(kmax_orig); + for (int k=0; k < kmax_orig; k++) { + if (stag == 'M') { + orig_z[k] = orig_z_full(i,j,k); + } else if (stag == 'X') { + if (i == 0) { + orig_z[k] = orig_z_full(i,j,k); + } else if (i == imax_orig) { + orig_z[k] = orig_z_full(imax_orig-1,j,k); + } else { + orig_z[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i-1,j,k)); + } + } else if (stag == 'Y') { + if (j == 0) { + orig_z[k] = orig_z_full(i,j,k); + } else if (j == jmax_orig) { + orig_z[k] = orig_z_full(i,jmax_orig-1,k); + } else { + orig_z[k] = 0.5*(orig_z_full(i,j,k)+orig_z_full(i,j-1,k)); + } + } + } + + // Check if the data is top-down instead of bottom-up. + bool flip_data_required = false; + if (orig_z[1] > orig_z[kmax_orig-1]) flip_data_required = true; + if (flip_data_required) amrex::Abort("metgrid initialization flip_data_required. Not yet implemented."); + + // Search for the first level above the surface in the origin data. + // This is needed since the origin model topography will be + // different than the topography processed by WPS. + int k_above_sfc = 0; + for (int k=1; k < kmax_orig; k++) { + if (orig_z[k] > orig_z[0]) { + k_above_sfc = k; + break; + } + } + + int zap = 0; + int zap_below = 0; + int zap_above = 0; + amrex::Vector ordered_z, ordered_data; + + if (k_above_sfc > 1) { + // The levels are not monotonically increasing in height, so + // we sort and then make "artistic" quality control choices. + int count = 0; + + // Insert levels that are below the surface. + for (int k=1; k < k_above_sfc; k++) { + ordered_z.push_back(orig_z[k]); + ordered_data.push_back(orig_data(i,j,k)); + count++; + } + + // Check if the level that is nearest to and below the surface + // is "too close". If so, we'll ignore the upper level and keep + // the lower. Origin data is likely to be on pressure levels + // with higher spatial resolution near-surface, which supports + // the choice of eliminating levels that are "too close" in + // pressure-space. For simplicity, calculate delta P assuming a + // baroclinic atmosphere. + amrex::Real Pu = p_0*exp(-CONST_GRAV*orig_z[0]/(t_0*R_d)); + amrex::Real Pl = p_0*exp(-CONST_GRAV*ordered_z[count-1]/(t_0*R_d)); + if (Pl-Pu < metgrid_proximity) { + count--; + zap = 1; + zap_below = 1; + } + + // Insert the surface level. + ordered_z.push_back(orig_z[0]); + ordered_data.push_back(orig_data(i,j,0)); + count++; + + // Quoting WRF's comments, the next level to use is at, + // "... ta da, the first level above the surface. I know, wow." + int knext = k_above_sfc; + // Conditionally more strongly use the surface data by removing + // levels between the surface and the height corresponding to a + // set number of ERF levels from the surface. This forces the + // interpolator to use the surface data up through a number of + // ERF levels from the surface. + if (metgrid_force_sfc_k > 0) { + for (int k=k_above_sfc; k < kmax_orig; k++) { + if (orig_z[k] > new_z[metgrid_force_sfc_k-1]) { + knext = k; + break; + } else { + zap++; + zap_above++; + } + } + } + + // Check if the level that is nearest to and above the surface + // is "too close". If so, we'll ignore that level. + Pu = p_0*exp(-CONST_GRAV*orig_z[knext]/(t_0*R_d)); + Pl = p_0*exp(-CONST_GRAV*ordered_z[count-1]/(t_0*R_d)); + if (Pl-Pu < metgrid_proximity) { + knext++; + zap++; + zap_above++; + } + + // Insert levels that are above the surface. + for (int k=knext; k < kmax_orig; k++) { + ordered_z.push_back(orig_z[k]); + ordered_data.push_back(orig_data(i,j,k)); + count++; + } + + } else { + // The surface is the lowest level in the origin data. + + // Insert the surface. + ordered_z.push_back(orig_z[0]); + ordered_data.push_back(orig_data(i,j,0)); + + // Similar to above, conditionally more strongly use the + // surface data. + int count = 1; + int knext = count; + if (metgrid_force_sfc_k > 0) { + for (int k=knext; k < kmax_orig; k++) { + if (orig_z[k] > new_z[metgrid_force_sfc_k]) { + knext = k; + break; + } else { + zap++; + zap_above++; + } + } + } + + // Insert the remaining levels, again ignoring levels that are + // "too close" to the prior valid level. + for (int k=knext; k < kmax_orig; k++) { + amrex::Real Pu = p_0*exp(-CONST_GRAV*orig_z[k]/(t_0*R_d)); + amrex::Real Pl = p_0*exp(-CONST_GRAV*ordered_z[count-1]/(t_0*R_d)); + if (Pl-Pu < metgrid_proximity) { + zap++; + zap_above++; + continue; + } + ordered_z.push_back(orig_z[k]); + ordered_data.push_back(orig_data(i,j,k)); + count++; + } + } + + int ksta, kend; + if (metgrid_use_below_sfc && metgrid_use_sfc) { + // Use all levels. + int count = 0; + for (int k=1; k < ordered_z.size(); k++) { + count++; + if (ordered_z[ordered_z.size()-1] == orig_z[k]) break; + } + ksta = 0; + kend = ksta+count-1; + } else if (metgrid_use_below_sfc && !metgrid_use_sfc) { + // Use all levels except for the surface. + int ksfc = 0; + for (int k=0; k < kmax_orig; k++) { + if (ordered_z[k] == orig_z[0]) { + ksfc = k; + break; + } + } + for (int k=ksfc; k < kmax_orig-1; k++) { + ordered_z[k] = ordered_z[k+1]; + ordered_data[k] = ordered_data[k+1]; + } + ordered_z.resize(kmax_orig-1); + ordered_data.resize(kmax_orig-1); + int count = 0; + for (int k=0; k < kmax_orig; k++) { + count++; + if (orig_z[kmax_orig-1] == ordered_z[k]) break; + } + ksta = 0; + kend = ksta+count-1; + + } else if (!metgrid_use_below_sfc && metgrid_use_sfc) { + // Use all levels above and including the surface. + int kcount = k_above_sfc-1-zap_below; + int count = 0; + for (int k=0; k < kmax_orig; k++) { + if (ordered_z[kcount] == orig_z[k]) { + kcount++; + count++; + } + } + ksta = k_above_sfc-1-zap_below; + kend = ksta+count-1; + } else { + // We shouldn't be in here! + amrex::Abort("metgrid initialization, !use_levels below_ground && !metgrid_use_sfc"); + } + + // Insert the level of maximum winds. +// amrex::Real maxw_above_this_level = 30000.0; +// amrex::Real maxw_horiz_pres_diff = 5000.0; +// if ((flag_maxw == 1) && (use_maxw)) { +// amrex::Abort("metgrid initialization, use_maxw not yet implemented"); +// } + + // Insert the level of the tropopause. +// amrex::Real trop_horiz_pres_diff = 5000.0; +// if ((flag_trop == 1) && (use_trop)) { +// amrex::Abort("metgrid initialization, use_trop not yet implemented"); +// } + + amrex::Vector ordered_p; + for (int k=0; k < ordered_z.size(); k++) { + ordered_p.push_back(p_0*exp(-CONST_GRAV*ordered_z[k]/(t_0*R_d))); + } + + int new_n = 1; + int zap_final = 0; + amrex::Vector final_z, final_p, final_data; + final_z.push_back(ordered_z[ksta]); + final_p.push_back(ordered_p[ksta]); + final_data.push_back(ordered_data[ksta]); + for (int k=ksta+1; k <= kend; k++) { + if ((final_p[new_n-1]-ordered_p[k]) < metgrid_proximity) { + zap_final++; + } else { + new_n++; + final_z.push_back(ordered_z[k]); + final_p.push_back(ordered_p[k]); + final_data.push_back(ordered_data[k]); + } + } + kend -= zap_final; + + // Call the interpolator. + lagrange_setup(var_type, + exp_interp, + final_z, + final_p, + final_data, + kend-ksta+1, + i, + j, + kmax_new, + metgrid_order, + new_z, + new_p, + new_data); + + // Optionally replace the lowest level of data with the surface + // field from the origin data. + if (metgrid_retain_sfc) new_data[0] = ordered_data[0]; + + // Save the interpolated data. +// amrex::Print() << "Saving interpolated data i=" << i << " j=" << j << std::endl; + for (int k=0; k < kmax_new; k++) { + if ((mask(i,j,k)) && (update_bc_data)) bc_data_full(i,j,k,0) = new_data[k]; + if (it == 0) new_data_full(i,j,k,src_comp) = new_data[k]; + } + +} + +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& orig_z, + const amrex::Array4& orig_data, + const amrex::Array4& new_z) { // This subroutine is a bit ham-handed and can be cleaned up later. int imax_orig = amrex::ubound(amrex::Box(orig_data)).x;