From 41f6533f518aefc580a96678db8885cda5f6e400 Mon Sep 17 00:00:00 2001 From: wiersema1 <46633497+wiersema1@users.noreply.github.com> Date: Fri, 6 Dec 2024 11:53:49 -0800 Subject: [PATCH] Metgrid interpolator (#1630) * Porting init_type="metgrid" updates into a clean branch of development. This includes,\n * an updated vertical interpolator used when interpolating met_em files onto the ERF vertical grid.\n * option to either interpolate potential temperature or interpolate pressure and temperature then calculate potential temperature.\n * debugging runtime options to overwrite incoming met_em data with idealized hard-coded values. * Updated documentation for init_type="metgrid" options. Made interp_theta an input parameter. * Changed temporary FABs used in init_type="metgrid" from a vector of FABs for each met_em timestamp to a single FAB since we only ever need access to these FABs for whatever met_em time is currently being worked on. Minor clean up of the metgrid initialization code. * Fixed a bug introduced when cleaning up code for init_type="metgrid" * Undoing the recent switch from dynamically sized vectors to Array1D in the vertical interpolation scheme used for met_em data * Make this GPU compatible. * Fix remaining issues with bounds and unused vars. * Remove size var since we went to fixed length GpuArrays. * Reduce memory overhead and fix GPU warning. * Added assertions during initialization for init_type of "metgrid" and "real" that there are sufficient points in a relaxation zone for computing the laplacian. --------- Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com> Co-authored-by: AMLattanzi Co-authored-by: Aaron Lattanzi --- Docs/sphinx_doc/Inputs.rst | 110 ++ Source/ERF.H | 15 + Source/ERF.cpp | 15 + Source/IO/ERF_ReadFromMetgrid.cpp | 63 +- Source/Initialization/ERF_InitFromMetgrid.cpp | 1044 +++++++++-------- .../Initialization/ERF_InitFromWRFInput.cpp | 5 + Source/Initialization/ERF_MetgridUtils.H | 690 ++++++++++- Submodules/AMReX | 2 +- 8 files changed, 1418 insertions(+), 526 deletions(-) diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 29136d252..46ec4148a 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -1219,6 +1219,107 @@ List of Parameters | **erf.project_initial_velocity** | project initial | true or false | true if anelastic; | | | velocity? | | false if compressible | +----------------------------------+-------------------+--------------------+-----------------------+ +| **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_debug_quiescent** | If init_type is | true or false | false | +| | metgrid, overwrite| | | +| | initial conditions| | | +| | and boundary | | | +| | conditions to be | | | +| | quiescent. | | | ++----------------------------------+-------------------+--------------------+-----------------------+ +| **erf.metgrid_debug_isothermal** | If init_type is | true or false | false | +| | metgrid, overwrite| | | +| | theta to be 300 in| | | +| | initial conditions| | | +| | and boundary | | | +| | conditions. | | | ++----------------------------------+-------------------+--------------------+-----------------------+ +| **erf.metgrid_debug_dry** | If init_type is | true or false | false | +| | metgrid, overwrite| | | +| | qv to be dry in | | | +| | initial conditions| | | +| | and boundary | | | +| | conditions. | | | ++----------------------------------+-------------------+--------------------+-----------------------+ +| **erf.metgrid_debug_msf** | If init_type is | true or false | false | +| | metgrid, overwrite| | | +| | map scale factors | | | +| | to be 1. | | | ++----------------------------------+-------------------+--------------------+-----------------------+ +| **erf.metgrid_debug_psfc** | If init_type is | true or false | false | +| | metgrid, overwrite| | | +| | surface pressure | | | +| | to be 10**5. | | | ++----------------------------------+-------------------+--------------------+-----------------------+ +| **erf.metgrid_interp_theta** | If init_type is | true or false | false | +| | metgrid, calculate| | | +| | theta on origin | | | +| | model vertical | | | +| | levels and then | | | +| | interpolate onto | | | +| | the ERF vertical | | | +| | levels. | | | ++----------------------------------+-------------------+--------------------+-----------------------+ +| **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 ----------------- @@ -1233,6 +1334,15 @@ The extent of the relaxation zone may be controlled with ``erf.real_width`` (cor and ``erf.real_set_width`` (corresponding to WRF's **spec_zone**, typically set to 1), which corresponds to a relaxation zone with a width of **real_width - real_set_width**. +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/Source/ERF.H b/Source/ERF.H index 9556d5316..d98d2be63 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -1050,6 +1050,21 @@ private: // used with init_type == InitType::Input_Sounding static bool init_sounding_ideal; + // Options for vertical interpolation of met_em*.nc data. + bool metgrid_debug_quiescent{false}; + bool metgrid_debug_isothermal{false}; + bool metgrid_debug_dry{false}; + bool metgrid_debug_psfc{false}; + bool metgrid_debug_msf{false}; + bool metgrid_interp_theta{false}; + 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 7b427e1ad..306eaf40f 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -1503,6 +1503,21 @@ 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_debug_quiescent", metgrid_debug_quiescent); + pp.query("metgrid_debug_isothermal", metgrid_debug_isothermal); + pp.query("metgrid_debug_dry", metgrid_debug_dry); + pp.query("metgrid_debug_psfc", metgrid_debug_psfc); + pp.query("metgrid_debug_msf", metgrid_debug_msf); + pp.query("metgrid_interp_theta", metgrid_interp_theta); + 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); + // Set default to FullState for now ... later we will try Perturbation interpolation_type = StateInterpType::FullState; pp.query_enum_case_insensitive("interpolation_type" ,interpolation_type); diff --git a/Source/IO/ERF_ReadFromMetgrid.cpp b/Source/IO/ERF_ReadFromMetgrid.cpp index b2fa08971..63681cd27 100644 --- a/Source/IO/ERF_ReadFromMetgrid.cpp +++ b/Source/IO/ERF_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_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp index 2abfaceab..020e87b71 100644 --- a/Source/Initialization/ERF_InitFromMetgrid.cpp +++ b/Source/Initialization/ERF_InitFromMetgrid.cpp @@ -16,23 +16,29 @@ 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 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); + // At least three points are necessary if there is a relaxation zone. + if (real_width > real_set_width) + AMREX_ALWAYS_ASSERT(real_width-real_set_width >= 3); + // Size the SST and LANDMASK sst_lev[lev].resize(ntimes); lmask_lev[lev].resize(ntimes); @@ -58,10 +64,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); @@ -78,22 +81,35 @@ ERF::init_from_metgrid (int lev) Arena_Used = The_Pinned_Arena(); #endif - for (int it = 0; it < ntimes; it++) { - 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], - 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], - NC_ght_fab[it], NC_hgt_fab[it], NC_psfc_fab[it], - NC_MSFU_fab[it], NC_MSFV_fab[it], NC_MSFM_fab[it], - NC_sst_fab[it], NC_LAT_fab[it], NC_LON_fab[it], - NC_lmask_iab[it], Latitude, Longitude, geom[lev]); - } // it +#ifndef AMREX_USE_GPU + if (metgrid_debug_quiescent) Print() << "metgrid_debug_quiescent = true" << std::endl; + if (metgrid_debug_isothermal) Print() << "metgrid_debug_isothermal = true" << std::endl; + if (metgrid_debug_dry) Print() << "metgrid_debug_dry = true" << std::endl; + if (metgrid_debug_psfc) Print() << "metgrid_debug_psfc = true" << std::endl; + if (metgrid_debug_msf) Print() << "metgrid_debug_msf = true" << std::endl; + if (metgrid_interp_theta) Print() << "metgrid_interp_theta = true" << std::endl; + if (metgrid_basic_linear) Print() << "metgrid_basic_linear = true" << std::endl; +#endif + + for (int itime(0); itime < ntimes; itime++) { +#ifndef AMREX_USE_GPU + Print() << " init_from_metgrid: reading nc_init_file[" << lev << "][" << itime << "]\t" << nc_init_file[lev][itime] << std::endl; +#endif + read_from_metgrid(lev, boxes_at_level[lev][0], nc_init_file[lev][itime], + NC_dateTime[itime], NC_epochTime[itime], + flag_psfc[itime], flag_msf[itime], + flag_sst[itime], flag_lmask[itime], + NC_nx[itime], NC_ny[itime], NC_dx[itime], NC_dy[itime], + NC_xvel_fab[itime], NC_yvel_fab[itime], + NC_temp_fab[itime], NC_rhum_fab[itime], NC_pres_fab[itime], + NC_ght_fab[itime], NC_hgt_fab[itime], NC_psfc_fab[itime], + NC_MSFU_fab[itime], NC_MSFV_fab[itime], NC_MSFM_fab[itime], + NC_sst_fab[itime], NC_LAT_fab[itime], NC_LON_fab[itime], + NC_lmask_iab[itime], Latitude, Longitude, geom[lev]); + } // itime // Verify that files in nc_init_file[lev] are ordered from earliest to latest. - for (int it = 1; it < ntimes; it++) AMREX_ALWAYS_ASSERT(NC_epochTime[it] > NC_epochTime[it-1]); + for (int itime(1); itime < ntimes; itime++) AMREX_ALWAYS_ASSERT(NC_epochTime[itime] > NC_epochTime[itime-1]); // Start at the earliest time in nc_init_file[lev]. start_bdy_time = NC_epochTime[0]; @@ -104,20 +120,12 @@ ERF::init_from_metgrid (int lev) bdy_time_interval = NC_epochTime[1]-NC_epochTime[0]; // 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."); - } - - // Set up a FAB for mixing ratio and another for potential temperature. - // Necessary because the input data has relative humidity and temperature, not mixing ratio and potential temperature. - // TODO: add alternate pathways for other origin models where different combinations of variables may be present. - Vector mxrat_fab; mxrat_fab.resize(ntimes); - Vector theta_fab; theta_fab.resize(ntimes); - for (int it = 0; it < ntimes; it++) { - Box NC_box_unstag = NC_rhum_fab[it].box(); - mxrat_fab[it].resize(NC_box_unstag, 1, Arena_Used); - theta_fab[it].resize(NC_box_unstag, 1, Arena_Used); + for (int itime(1); itime < ntimes; itime++) { + Real NC_dt = NC_epochTime[itime]-NC_epochTime[itime-1]; +#ifndef AMREX_USE_GPU + Print() << " " << nc_init_file[lev][itime-1] << " / " << nc_init_file[lev][itime] << " 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."); } auto& lev_new = vars_new[lev]; @@ -126,9 +134,6 @@ ERF::init_from_metgrid (int lev) AMREX_ALWAYS_ASSERT(SolverChoice::terrain_type != TerrainType::None); - // 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 ) { @@ -152,44 +157,44 @@ ERF::init_from_metgrid (int lev) int i_lo = geom[lev].Domain().smallEnd(0); int i_hi = geom[lev].Domain().bigEnd(0); int j_lo = geom[lev].Domain().smallEnd(1); int j_hi = geom[lev].Domain().bigEnd(1); if (flag_sst[0]) { - for (int it = 0; it < ntimes; ++it) { - sst_lev[lev][it] = std::make_unique(ba2d,dm,1,ngv); - for ( MFIter mfi(*(sst_lev[lev][it]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + for (int itime(0); itime < ntimes; ++itime) { + sst_lev[lev][itime] = std::make_unique(ba2d,dm,1,ngv); + for ( MFIter mfi(*(sst_lev[lev][itime]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Box gtbx = mfi.growntilebox(); - FArrayBox& dst = (*(sst_lev[lev][it]))[mfi]; - FArrayBox& src = NC_sst_fab[it]; + FArrayBox& dst = (*(sst_lev[lev][itime]))[mfi]; + FArrayBox& src = NC_sst_fab[itime]; const Array4< Real>& dst_arr = dst.array(); 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); }); } - sst_lev[lev][it]->FillBoundary(geom[lev].periodicity()); + sst_lev[lev][itime]->FillBoundary(geom[lev].periodicity()); } } else { - for (int it = 0; it < ntimes; ++it) sst_lev[lev][it] = nullptr; + for (int itime(0); itime < ntimes; ++itime) sst_lev[lev][itime] = nullptr; } if (flag_lmask[0]) { - for (int it = 0; it < ntimes; ++it) { - lmask_lev[lev][it] = std::make_unique(ba2d,dm,1,ngv); - for ( MFIter mfi(*(lmask_lev[lev][it]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + for (int itime(0); itime < ntimes; ++itime) { + lmask_lev[lev][itime] = std::make_unique(ba2d,dm,1,ngv); + for ( MFIter mfi(*(lmask_lev[lev][itime]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Box gtbx = mfi.growntilebox(); - IArrayBox& dst = (*(lmask_lev[lev][it]))[mfi]; - IArrayBox& src = NC_lmask_iab[it]; + IArrayBox& dst = (*(lmask_lev[lev][itime]))[mfi]; + IArrayBox& src = NC_lmask_iab[itime]; const Array4< int>& dst_arr = dst.array(); 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); }); } - lmask_lev[lev][it]->FillBoundary(geom[lev].periodicity()); + lmask_lev[lev][itime]->FillBoundary(geom[lev].periodicity()); } } @@ -202,8 +207,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); }); } @@ -217,21 +222,21 @@ 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); }); } - for (int it = 0; it < ntimes; it++) { + for (int itime(0); itime < ntimes; itime++) { // Verify that the grid size and resolution from met_em file matches that in geom (from ERF inputs file). - AMREX_ALWAYS_ASSERT(geom[lev].CellSizeArray()[0] == NC_dx[it]); - AMREX_ALWAYS_ASSERT(geom[lev].CellSizeArray()[1] == NC_dy[it]); + AMREX_ALWAYS_ASSERT(geom[lev].CellSizeArray()[0] == NC_dx[itime]); + AMREX_ALWAYS_ASSERT(geom[lev].CellSizeArray()[1] == NC_dy[itime]); // NC_nx-2 because NC_nx is the number of staggered grid points indexed from 1. - AMREX_ALWAYS_ASSERT(geom[lev].Domain().bigEnd(0) == NC_nx[it]-2); + AMREX_ALWAYS_ASSERT(geom[lev].Domain().bigEnd(0) == NC_nx[itime]-2); // NC_ny-2 because NC_ny is the number of staggered grid points indexed from 1. - AMREX_ALWAYS_ASSERT(geom[lev].Domain().bigEnd(1) == NC_ny[it]-2); - } // it + AMREX_ALWAYS_ASSERT(geom[lev].Domain().bigEnd(1) == NC_ny[itime]-2); + } // itime // This makes the Jacobian. make_J(geom[lev],*z_phys, *detJ_cc[lev]); @@ -244,35 +249,96 @@ ERF::init_from_metgrid (int lev) int MetGridBdyEnd = MetGridBdyVars::NumTypes-1; if (use_moisture) MetGridBdyEnd = MetGridBdyVars::NumTypes; - // Zero out fabs_for_bcs on the global domain - amrex::Vector> fabs_for_bcs; - fabs_for_bcs.resize(ntimes); - for (int it(0); it < ntimes; it++) { - fabs_for_bcs[it].resize(MetGridBdyEnd); + // Zero out the bdy data to start off + bdy_data_xlo.resize(ntimes); + bdy_data_xhi.resize(ntimes); + bdy_data_ylo.resize(ntimes); + bdy_data_yhi.resize(ntimes); + + for (int itime(0); itime < ntimes; itime++) { + bdy_data_xlo[itime].resize(MetGridBdyEnd); + bdy_data_xhi[itime].resize(MetGridBdyEnd); + bdy_data_ylo[itime].resize(MetGridBdyEnd); + bdy_data_yhi[itime].resize(MetGridBdyEnd); + + // Build the boxes for each BDY FAB + const auto& lo = geom[lev].Domain().loVect(); + const auto& hi = geom[lev].Domain().hiVect(); + IntVect plo(lo); + IntVect phi(hi); + + plo[0] = lo[0]; plo[1] = lo[1]; plo[2] = lo[2]; + phi[0] = lo[0]+real_width-1; phi[1] = hi[1]; phi[2] = hi[2]; + const Box pbx_xlo(plo, phi); + Box xlo_plane_no_stag(pbx_xlo); + Box xlo_plane_x_stag = pbx_xlo; xlo_plane_x_stag.shiftHalf(0,-1); + Box xlo_plane_y_stag = convert(pbx_xlo, {0, 1, 0}); + + plo[0] = hi[0]-real_width+1; plo[1] = lo[1]; plo[2] = lo[2]; + phi[0] = hi[0]; phi[1] = hi[1]; phi[2] = hi[2]; + const Box pbx_xhi(plo, phi); + Box xhi_plane_no_stag(pbx_xhi); + Box xhi_plane_x_stag = pbx_xhi; xhi_plane_x_stag.shiftHalf(0,1); + Box xhi_plane_y_stag = convert(pbx_xhi, {0, 1, 0}); + + plo[1] = lo[1]; plo[0] = lo[0]; plo[2] = lo[2]; + phi[1] = lo[1]+real_width-1; phi[0] = hi[0]; phi[2] = hi[2]; + const Box pbx_ylo(plo, phi); + Box ylo_plane_no_stag(pbx_ylo); + Box ylo_plane_x_stag = convert(pbx_ylo, {1, 0, 0}); + Box ylo_plane_y_stag = pbx_ylo; ylo_plane_y_stag.shiftHalf(1,-1); + + plo[1] = hi[1]-real_width+1; plo[0] = lo[0]; plo[2] = lo[2]; + phi[1] = hi[1]; phi[0] = hi[0]; phi[2] = hi[2]; + const Box pbx_yhi(plo, phi); + Box yhi_plane_no_stag(pbx_yhi); + Box yhi_plane_x_stag = convert(pbx_yhi, {1, 0, 0}); + Box yhi_plane_y_stag = pbx_yhi; yhi_plane_y_stag.shiftHalf(1,1); - Box gdomain; - Box ldomain; for (int nvar(0); nvar(0.0); + bdy_data_xlo[itime][nvar].template setVal(0.0); + bdy_data_xhi[itime][nvar].template setVal(0.0); + bdy_data_ylo[itime][nvar].template setVal(0.0); + bdy_data_yhi[itime][nvar].template setVal(0.0); } } + // Set up FABs for mixing ratio and potential temperature on the origin model vertical grid. + // This is necessary for input data that has relative humidity and temperature, but does not + // include mixing ratio and potential temperature. + // TODO: add alternate pathways for other origin models where different combinations of variables may be present. + Box NC_box_unstag = NC_rhum_fab[0].box(); + FArrayBox mxrat_fab, theta_fab; + mxrat_fab.resize(NC_box_unstag, 1, Arena_Used); + theta_fab.resize(NC_box_unstag, 1, Arena_Used); + + // Set up FABs to hold interpolated origin variables that are only needed during initialization. + // Pressure will be iteratively calculated from z, qv, and theta, but the origin model pressure + // is necessary if metgrid_interp_theta=True and we interpolate temperature instead of theta. + FArrayBox p_interp_fab, t_interp_fab; + if (!metgrid_interp_theta) { + Box ldomain = geom[lev].Domain(); + p_interp_fab.resize(ldomain, 1, Arena_Used); + p_interp_fab.template setVal(0.0); + t_interp_fab.resize(ldomain, 1, Arena_Used); + t_interp_fab.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()); @@ -303,14 +369,22 @@ 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, metgrid_interp_theta, metgrid_debug_quiescent, + metgrid_debug_isothermal, metgrid_debug_dry, + 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, - fabs_for_bcs, mask_c_arr, mask_u_arr, mask_v_arr); + NC_pres_fab, p_interp_fab, t_interp_fab, + theta_fab, mxrat_fab, + bdy_data_xlo, bdy_data_xhi, + bdy_data_ylo, bdy_data_yhi, + mask_c_arr, mask_u_arr, mask_v_arr); } // mf @@ -324,8 +398,8 @@ 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(metgrid_debug_msf, + msfu_fab, msfv_fab, msfm_fab, flag_msf[0], NC_MSFU_fab, NC_MSFV_fab, NC_MSFM_fab); } // mf @@ -342,25 +416,23 @@ ERF::init_from_metgrid (int lev) FArrayBox& cons_fab = lev_new[Vars::cons][mfi]; FArrayBox& z_phys_nd_fab = (*z_phys)[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 // th_hse calculate potential temperature const Box valid_bx = mfi.validbox(); - init_base_state_from_metgrid(use_moisture, l_rdOcp, - valid_bx, - flag_psfc, + init_base_state_from_metgrid(use_moisture, metgrid_debug_psfc, l_rdOcp, + valid_bx, flag_psfc, cons_fab, r_hse_fab, p_hse_fab, pi_hse_fab, th_hse_fab, - z_phys_nd_fab, NC_ght_fab, NC_psfc_fab, - fabs_for_bcs); + z_phys_nd_fab, NC_psfc_fab); } // mf // FillBoundary to populate the internal halo cells r_hse.FillBoundary(geom[lev].periodicity()); p_hse.FillBoundary(geom[lev].periodicity()); pi_hse.FillBoundary(geom[lev].periodicity()); - th_hse.FillBoundary(geom[lev].periodicity()); // NOTE: fabs_for_bcs is defined over the whole domain on each rank. // However, the operations needed to define the data on the ERF @@ -368,141 +440,30 @@ ERF::init_from_metgrid (int lev) // So when we save the data in fabs_for_bc, only regions owned // by the rank are populated. Use an allreduce sum to make the // complete data set; initialized to 0 above. - for (int it(0); it < ntimes; it++) { + for (int itime(0); itime < ntimes; itime++) { for (int nvar(0); nvar& fabs_for_bcs_arr = fabs_for_bcs[it][ivar].const_array(); - - if (ivar == MetGridBdyVars::U) { - 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) { - 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::T) { - 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) { - 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 - - // west boundary - ParallelFor(xlo_plane, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - 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 - { - 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 - { - 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 - { - yhi_arr(i,j,k,0) = fabs_for_bcs_arr(i,j,k); - }); - - } // ivar - } // it } /** @@ -517,12 +478,12 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, { int ntimes = 1; // Use terrain from the first met_em file. - for (int it = 0; it < ntimes; it++) { + for (int itime(0); itime < ntimes; itime++) { // This copies from NC_zphys on z-faces to z_phys_nd on nodes const Array4& z_arr = z_phys_nd_fab.array(); - const Array4& nc_hgt_arr = NC_hgt_fab[it].const_array(); + const Array4& nc_hgt_arr = NC_hgt_fab[itime].const_array(); - const Box z_hgt_box = NC_hgt_fab[it].box(); + const Box z_hgt_box = NC_hgt_fab[itime].box(); int ilo = z_hgt_box.smallEnd()[0]; int ihi = z_hgt_box.bigEnd()[0]; @@ -530,7 +491,7 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, int jhi = z_hgt_box.bigEnd()[1]; Box z_phys_box = z_phys_nd_fab.box(); - Box from_box = surroundingNodes(NC_hgt_fab[it].box()); + Box from_box = surroundingNodes(NC_hgt_fab[itime].box()); from_box.growHi(2,-1); Box bx = z_phys_box & from_box; @@ -544,12 +505,24 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, z_arr(i,j,k) = 0.25 * ( nc_hgt_arr (ii,jj ,k) + nc_hgt_arr(ii-1,jj ,k) + nc_hgt_arr (ii,jj-1,k) + nc_hgt_arr(ii-1,jj-1,k) ); }); - } // it + } // itime } /** * Helper function to initialize state and velocity data read from metgrid data. * + * @param use_moisture bool True if solverChoice.moisture_type != MoistureType::None + * @param metgrid_interp_theta bool calculate theta on origin levels, then interpolate + * @param metgrid_debug_quiescent bool overwrite u and v with 0.0 + * @param metgrid_debug_isothermal bool overwrite theta with 300.0 + * @param metgrid_debug_dry bool overwrite qv with 0.0 + * @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 @@ -564,12 +537,28 @@ init_terrain_from_metgrid (FArrayBox& z_phys_nd_fab, * @param NC_temp_fab Vector of FArrayBox objects holding metgrid data for temperature * @param NC_rhum_fab Vector of FArrayBox objects holding metgrid data for relative humidity * @param NC_pres_fab Vector of FArrayBox objects holding metgrid data for pressure - * @param theta_fab Vector of FArrayBox objects holding potential temperature calculated from temperature and pressure - * @param mxrat_fab Vector of FArrayBox objects holding vapor mixing ratio calculated from relative humidity + * @param p_interp_fab FArrayBox object + * @param t_interp_fab FArrayBox object + * @param theta_fab FArrayBox object holding potential temperature calculated from temperature and pressure + * @param mxrat_fab FArrayBox object 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 metgrid_interp_theta, + const bool metgrid_debug_quiescent, + const bool metgrid_debug_isothermal, + const bool metgrid_debug_dry, + const bool metgrid_basic_linear, + const bool metgrid_use_below_sfc, + const bool metgrid_use_sfc, + const bool metgrid_retain_sfc, + const Real metgrid_proximity, + const int metgrid_order, + const int metgrid_force_sfc_k, const Real l_rdOcp, Box& tbxc, Box& tbxu, @@ -586,36 +575,76 @@ init_state_from_metgrid (const bool use_moisture, const Vector& NC_temp_fab, const Vector& NC_rhum_fab, const Vector& NC_pres_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) + FArrayBox& p_interp_fab, + FArrayBox& t_interp_fab, + FArrayBox& theta_fab, + FArrayBox& mxrat_fab, + Vector>& fabs_for_bcs_xlo, + Vector>& fabs_for_bcs_xhi, + Vector>& fabs_for_bcs_ylo, + Vector>& fabs_for_bcs_yhi, + 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++) + for (int itime(0); itime < ntimes; itime++) { + // ******************************************************** // U // ******************************************************** { - Box bx2d = NC_xvel_fab[it].box() & tbxu; +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of u-velocity, itime = " << itime << std::endl; +#endif + Box bx2d = NC_xvel_fab[itime].box() & tbxu; bx2d.setRange(2,0); - auto const orig_data = NC_xvel_fab[it].const_array(); - auto const orig_z = NC_ght_fab[it].const_array(); + auto const orig_data = NC_xvel_fab[itime].const_array(); + auto const orig_z = NC_ght_fab[itime].const_array(); auto new_data = x_vel_fab.array(); - 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; + Box bx_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::U].box(); + Box bx_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::U].box(); + Box bx_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::U].box(); + Box bx_yhi = fabs_for_bcs_yhi[itime][MetGridBdyVars::U].box(); + auto bc_data_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::U].array(); + auto bc_data_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::U].array(); + auto bc_data_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::U].array(); + auto bc_data_yhi = fabs_for_bcs_yhi[itime][MetGridBdyVars::U].array(); + + 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_debug_quiescent) { // Debugging option to run quiescent. + for (int k(0); k<=kmax; k++) { + if (mask_u_arr(i,j,k) && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = 0.0; + if (mask_u_arr(i,j,k) && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = 0.0; + if (mask_u_arr(i,j,k) && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = 0.0; + if (mask_u_arr(i,j,k) && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = 0.0; + if (itime == 0) new_data(i,j,k,0) = 0.0; + } + } else if (metgrid_basic_linear) { // Linear interpolation with no quality control. + 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) && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = Interp_Val; + if (mask_u_arr(i,j,k) && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = Interp_Val; + if (mask_u_arr(i,j,k) && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = Interp_Val; + if (mask_u_arr(i,j,k) && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = Interp_Val; + if (itime == 0) new_data(i,j,k,0) = Interp_Val; + } + } else { // Vertical interpolation and quality control similar to that from WRF. + 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, itime, 'U', 'X', + orig_z, orig_data, new_z, new_data, true, + bc_data_xlo, bc_data_xhi, bc_data_ylo, bc_data_yhi, + bx_xlo, bx_xhi, bx_ylo, bx_yhi, mask_u_arr); } }); } @@ -625,22 +654,53 @@ init_state_from_metgrid (const bool use_moisture, // V // ******************************************************** { - Box bx2d = NC_yvel_fab[it].box() & tbxv; +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of v-velocity, itime = " << itime << std::endl; +#endif + Box bx2d = NC_yvel_fab[itime].box() & tbxv; bx2d.setRange(2,0); - auto const orig_data = NC_yvel_fab[it].const_array(); - auto const orig_z = NC_ght_fab[it].const_array(); + auto const orig_data = NC_yvel_fab[itime].const_array(); + auto const orig_z = NC_ght_fab[itime].const_array(); auto new_data = y_vel_fab.array(); - 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; + Box bx_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::V].box(); + Box bx_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::V].box(); + Box bx_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::V].box(); + Box bx_yhi = fabs_for_bcs_yhi[itime][MetGridBdyVars::V].box(); + auto bc_data_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::V].array(); + auto bc_data_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::V].array(); + auto bc_data_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::V].array(); + auto bc_data_yhi = fabs_for_bcs_yhi[itime][MetGridBdyVars::V].array(); + + 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_debug_quiescent) { // Debugging option to run quiescent. + for (int k(0); k<=kmax; k++) { + if (mask_v_arr(i,j,k) && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = 0.0; + if (mask_v_arr(i,j,k) && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = 0.0; + if (mask_v_arr(i,j,k) && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = 0.0; + if (mask_v_arr(i,j,k) && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = 0.0; + if (itime == 0) new_data(i,j,k,0) = 0.0; + } + } else if (metgrid_basic_linear) { // Linear interpolation with no quality control. + 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) && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = Interp_Val; + if (mask_v_arr(i,j,k) && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = Interp_Val; + if (mask_v_arr(i,j,k) && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = Interp_Val; + if (mask_v_arr(i,j,k) && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = Interp_Val; + if (itime == 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, itime, 'V', 'Y', + orig_z, orig_data, new_z, new_data, true, + bc_data_xlo, bc_data_xhi, bc_data_ylo, bc_data_yhi, + bx_xlo, bx_xhi, bx_ylo, bx_yhi, mask_v_arr); } }); } @@ -649,14 +709,15 @@ init_state_from_metgrid (const bool use_moisture, // ******************************************************** // W // ******************************************************** - if (it == 0) { // update at initialization + if (itime == 0) { // update at initialization z_vel_fab.template setVal(0.0); } + // ******************************************************** // Initialize all state_fab variables to zero // ******************************************************** - if (it == 0) { // update at initialization + if (itime == 0) { // update at initialization state_fab.template setVal(0.0); } @@ -664,40 +725,176 @@ 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 (metgrid_interp_theta) { + // Calculate potential temperature on the origin model vertical levels + // then interpolate that onto the ERF vertical levels. + + { // calculate potential temperature. + Box bx = NC_rhum_fab[itime].box() & tbxc; + auto const temp = NC_temp_fab[itime].const_array(); + auto const pres = NC_pres_fab[itime].const_array(); + auto theta = theta_fab.array(); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + 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); + }); + } + + { // vertical interpolation of potential temperature. +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of potential temperature, itime = " << itime << std::endl; +#endif + Box bx2d = NC_temp_fab[itime].box() & tbxc; + bx2d.setRange(2,0); + auto const orig_data = theta_fab.const_array(); + auto const orig_z = NC_ght_fab[itime].const_array(); + auto new_data = state_fab.array(); + auto const new_z = z_phys_nd_fab.const_array(); + + Box bx_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::T].box(); + Box bx_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::T].box(); + Box bx_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::T].box(); + Box bx_yhi = fabs_for_bcs_yhi[itime][MetGridBdyVars::T].box(); + auto bc_data_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::T].array(); + auto bc_data_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::T].array(); + auto bc_data_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::T].array(); + auto bc_data_yhi = fabs_for_bcs_yhi[itime][MetGridBdyVars::T].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_debug_isothermal) { // Debugging option to run isothermal. + for (int k(0); k<=kmax; k++) { + if (mask_c_arr(i,j,k) && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = 300.0; + if (mask_c_arr(i,j,k) && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = 300.0; + if (mask_c_arr(i,j,k) && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = 300.0; + if (mask_c_arr(i,j,k) && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = 300.0; + if (itime == 0) new_data(i,j,k,RhoTheta_comp) = 300.0; + } + } else if (metgrid_basic_linear) { // Linear interpolation with no quality control. + 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) && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = Interp_Val; + if (mask_c_arr(i,j,k) && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = Interp_Val; + if (mask_c_arr(i,j,k) && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = Interp_Val; + if (mask_c_arr(i,j,k) && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = Interp_Val; + if (itime == 0) new_data(i,j,k,RhoTheta_comp) = Interp_Val; + } + } else { // Vertical interpolation and quality control similar to that from WRF. + 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, itime, 'T', 'M', + orig_z, orig_data, new_z, new_data, true, + bc_data_xlo, bc_data_xhi, bc_data_ylo, bc_data_yhi, + bx_xlo, bx_xhi, bx_ylo, bx_yhi, 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 { // metgrid_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, itime = " << itime << std::endl; +#endif + Box bx2d = p_interp_fab.box() & tbxc; + bx2d.setRange(2,0); + auto const orig_data = NC_pres_fab[itime].const_array(); + auto const orig_z = NC_ght_fab[itime].const_array(); + auto new_data = p_interp_fab.array(); + auto const new_z = z_phys_nd_fab.const_array(); - 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; + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + const amrex::Array4 bc_data_xlo, bc_data_xhi, bc_data_ylo, bc_data_yhi; + + int kmax = ubound(tbxc).z; + + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + if (metgrid_basic_linear) { // Linear interpolation with no quality control. + 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 { // Vertical interpolation and quality control similar to that from WRF. + // 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_xlo, bc_data_xhi, bc_data_ylo, bc_data_yhi, + bx_xlo, bx_xhi, bx_ylo, bx_yhi, mask_c_arr); + } + }); } - }); - } + + { // vertical interpolation of temperature. +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of temperature, itime = " << itime << std::endl; +#endif + Box bx2d = p_interp_fab.box() & tbxc; + bx2d.setRange(2,0); + auto const orig_data = NC_temp_fab[itime].const_array(); + auto const orig_z = NC_ght_fab[itime].const_array(); + auto new_data = t_interp_fab.array(); + auto const new_z = z_phys_nd_fab.const_array(); + + Box bx_xlo, bx_xhi, bx_ylo, bx_yhi; + const amrex::Array4 bc_data_xlo, bc_data_xhi, bc_data_ylo, bc_data_yhi; + + int kmax = ubound(tbxc).z; + + ParallelFor(bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept + { + if (metgrid_basic_linear) { // Linear interpolation with no quality control. + 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 { // Vertical interpolation and quality control similar to that from WRF. + // 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_xlo, bc_data_xhi, bc_data_ylo, bc_data_yhi, + bx_xlo, bx_xhi, bx_ylo, bx_yhi, mask_c_arr); + } + }); + } + + { // calculate potential temperature on the ERF vertical levels. + auto const temp = t_interp_fab.const_array(); + auto const pres = p_interp_fab.const_array(); + auto new_data = state_fab.array(); + + Box bx_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::T].box(); + Box bx_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::T].box(); + Box bx_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::T].box(); + Box bx_yhi = fabs_for_bcs_yhi[itime][MetGridBdyVars::T].box(); + auto bc_data_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::T].array(); + auto bc_data_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::T].array(); + auto bc_data_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::T].array(); + auto bc_data_yhi = fabs_for_bcs_yhi[itime][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); + if (metgrid_debug_isothermal) Calc_Val = 300.0; // Debugging option to run isothermal. + if (mask_c_arr(i,j,k) && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = Calc_Val; + if (mask_c_arr(i,j,k) && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = Calc_Val; + if (mask_c_arr(i,j,k) && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = Calc_Val; + if (mask_c_arr(i,j,k) && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = Calc_Val; + if (itime == 0) new_data(i,j,k,RhoTheta_comp) = Calc_Val; + }); + } + + } // metgrid_interp_theta if (use_moisture) { // ******************************************************** @@ -709,11 +906,11 @@ init_state_from_metgrid (const bool use_moisture, // could be specific humidity or a mixing ratio. // { // calculate vapor mixing ratio from relative humidity. - Box bx = NC_temp_fab[it].box() & tbxc; - auto const rhum = NC_rhum_fab[it].const_array(); - auto const temp = NC_temp_fab[it].const_array(); - auto const pres = NC_pres_fab[it].const_array(); - auto mxrat = mxrat_fab[it].array(); + Box bx = NC_temp_fab[itime].box() & tbxc; + auto const rhum = NC_rhum_fab[itime].const_array(); + auto const temp = NC_temp_fab[itime].const_array(); + auto const pres = NC_pres_fab[itime].const_array(); + auto mxrat = mxrat_fab.array(); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -721,45 +918,68 @@ init_state_from_metgrid (const bool use_moisture, }); } - // vertical interpolation of vapor mixing ratio. - { - Box bx2d = NC_temp_fab[it].box() & tbxc; + { // vertical interpolation of vapor mixing ratio. +#ifndef AMREX_USE_GPU + Print() << "[init_state_from_metgrid] vertical interpolation of vapor mixing ratio, itime = " << itime << std::endl; +#endif + Box bx2d = NC_temp_fab[itime].box() & tbxc; bx2d.setRange(2,0); - auto const orig_data = mxrat_fab[it].const_array(); - auto const orig_z = NC_ght_fab[it].const_array(); + auto const orig_data = mxrat_fab.const_array(); + auto const orig_z = NC_ght_fab[itime].const_array(); auto new_data = state_fab.array(); - 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; + Box bx_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::QV].box(); + Box bx_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::QV].box(); + Box bx_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::QV].box(); + Box bx_yhi = fabs_for_bcs_yhi[itime][MetGridBdyVars::QV].box(); + auto bc_data_xlo = fabs_for_bcs_xlo[itime][MetGridBdyVars::QV].array(); + auto bc_data_xhi = fabs_for_bcs_xhi[itime][MetGridBdyVars::QV].array(); + auto bc_data_ylo = fabs_for_bcs_ylo[itime][MetGridBdyVars::QV].array(); + auto bc_data_yhi = fabs_for_bcs_yhi[itime][MetGridBdyVars::QV].array(); + + 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_debug_dry) { // Debugging option to run dry. + for (int k(0); k<=kmax; k++) { + if (mask_c_arr(i,j,k) && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = 0.0; + if (mask_c_arr(i,j,k) && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = 0.0; + if (mask_c_arr(i,j,k) && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = 0.0; + if (mask_c_arr(i,j,k) && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = 0.0; + if (itime == 0) new_data(i,j,k,state_indx) = 0.0; + } + } else if (metgrid_basic_linear) { // Linear interpolation with no quality control. + 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) && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = Interp_Val; + if (mask_c_arr(i,j,k) && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = Interp_Val; + if (mask_c_arr(i,j,k) && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = Interp_Val; + if (mask_c_arr(i,j,k) && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = Interp_Val; + if (itime == 0) new_data(i,j,k,state_indx) = Interp_Val; + } + } else { // Vertical interpolation and quality control similar to that from WRF. + 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, itime, 'Q', 'M', + orig_z, orig_data, new_z, new_data, true, + bc_data_xlo, bc_data_xhi, bc_data_ylo, bc_data_yhi, + bx_xlo, bx_xhi, bx_ylo, bx_yhi, mask_c_arr); } }); } - } - - // TODO: TEMPORARY CODE TO RUN QUIESCENT, REMOVE WHEN NOT NEEDED. -// if (it == 0) { -// x_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. -// y_vel_fab.template setVal(0.0); // TODO: temporary code to initialize with quiescent atmosphere. -// } -// fabs_for_bcs[it][MetGridBdyVars::U].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. -// fabs_for_bcs[it][MetGridBdyVars::V].template setVal(0.0); // TODO: temporary code to force with quiescent atmosphere. + } // use_moisture - } // it + } // itime } /** * 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 @@ -769,12 +989,13 @@ init_state_from_metgrid (const bool use_moisture, * @param pi_hse_fab FArrayBox holding the hydrostatic base Exner pressure we are initializing * @param th_hse_fab FArrayBox holding the base state potential temperature we are initializing * @param z_phys_nd_fab FArrayBox holding nodal z coordinate data for terrain - * @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, + const bool metgrid_debug_psfc, const Real l_rdOcp, const Box& valid_bx, const Vector& flag_psfc, @@ -784,12 +1005,10 @@ init_base_state_from_metgrid (const bool use_moisture, FArrayBox& pi_hse_fab, FArrayBox& th_hse_fab, FArrayBox& z_phys_cc_fab, - const Vector& /*NC_ght_fab*/, - const Vector& NC_psfc_fab, - Vector>& fabs_for_bcs) + const Vector& NC_psfc_fab) { 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 @@ -812,12 +1031,6 @@ init_base_state_from_metgrid (const bool use_moisture, Gpu::copy(Gpu::hostToDevice, flag_psfc.begin(), flag_psfc.end(), flag_psfc_d.begin()); int* flag_psfc_vec = flag_psfc_d.data(); - // Define the arena to be used for data allocation - Arena* Arena_Used = The_Arena(); -#ifdef AMREX_USE_GPU - // Inside MFiter use async arena - Arena_Used = The_Async_Arena(); -#endif // Expose for copy to GPU Real grav = CONST_GRAV; @@ -829,9 +1042,8 @@ init_base_state_from_metgrid (const bool use_moisture, auto psfc_flag = flag_psfc_vec[0]; // ******************************************************** - // 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(); @@ -852,7 +1064,9 @@ init_base_state_from_metgrid (const bool use_moisture, Real thetad_lo, thetad_hi; // Calculate or use pressure at the surface. - if (psfc_flag == 1) { + if (metgrid_debug_psfc) { + psurf = std::pow(10, 5); + } else if (psfc_flag == 1) { psurf = orig_psfc(i,j,0); } else { z_lo = new_z(i,j,0); @@ -870,7 +1084,7 @@ init_base_state_from_metgrid (const bool use_moisture, Real half_dz = z_lo; Real qvf = 1.0+(R_v/R_d)*qv_lo; Real thetam = thetad_lo*qvf; - for (int it=0; it{}; - auto p_hse_arr = p_hse_bcs_fab.array(); - - ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept - { - const int maxiter = 10; - const amrex::Real tol = 1.0e-10; - - // Low and Hi column variables - Real psurf; - Real z_lo, z_hi; - Real p_lo, p_hi; - Real qv_lo, qv_hi; - Real rd_lo, rd_hi; - Real thetad_lo, thetad_hi; - - // Calculate or use pressure at the surface. - if (psfc_flag == 1) { - psurf = orig_psfc(i,j,0); - } else { - z_lo = new_z(i,j,0); - Real t_0 = 290.0; // WRF's model_config_rec%base_temp - Real a = 50.0; // WRF's model_config_rec%base_lapse - psurf = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2.)-2.0*grav*z_lo/(a*R_d)), 0.5)); - } - - // Iterations for the first CC point that is 1/2 dz off the surface - { - z_lo = new_z(i,j,0); - qv_lo = (use_moisture) ? Q_arr(i,j,0) : 0.0; - rd_lo = 0.0; // initial guess - thetad_lo = Theta_arr(i,j,0); - Real half_dz = z_lo; - Real qvf = 1.0+(R_v/R_d)*qv_lo; - Real thetam = thetad_lo*qvf; - for (int it=0; ittol) HSEutils::Newton_Raphson_hse(tol, R_d/Cp_d, dz, - grav, C, thetad_hi, - qv_hi, qv_hi, p_hi, - rd_hi, F); - - // Copy solution to base state - p_hse_arr(i,j,k) = p_hi; - - // Copy hi to lo - z_lo = z_hi; - p_lo = p_hi; - qv_lo = qv_hi; - rd_lo = rd_hi; - thetad_lo = thetad_hi; - } - - }); - } // itime } @@ -1107,62 +1221,42 @@ 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 */ void -init_msfs_from_metgrid (FArrayBox& msfu_fab, +init_msfs_from_metgrid (const bool metgrid_debug_msf, + 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) { // int ntimes = NC_MSFU_fab.size(); int ntimes = 1; - for (int it = 0; it < ntimes; it++) { + for (int itime(0); itime < ntimes; itime++) { // // FArrayBox to FArrayBox copy does "copy on intersection" // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks // // This copies or sets mapfac_m - if (flag_msfm == 1) { - msfm_fab.template copy(NC_MSFM_fab[it]); + if ((flag_msf == 1) and (!metgrid_debug_msf)) { + msfm_fab.template copy(NC_MSFM_fab[itime]); + msfu_fab.template copy(NC_MSFU_fab[itime]); + msfv_fab.template copy(NC_MSFV_fab[itime]); } else { #ifndef AMREX_USE_GPU - Print() << " MAPFAC_M 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); - } - - // 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; -#endif msfv_fab.template setVal(1.0); } - - } // it + } // itime } #endif // ERF_USE_NETCDF diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp index 325712670..da74b3fea 100644 --- a/Source/Initialization/ERF_InitFromWRFInput.cpp +++ b/Source/Initialization/ERF_InitFromWRFInput.cpp @@ -305,6 +305,11 @@ ERF::init_from_wrfinput (int lev) if (nc_bdy_file.empty()) { amrex::Error("NetCDF boundary file name must be provided via input"); } + + // At least three points are necessary if there is a relaxation zone. + if (real_width > real_set_width) + AMREX_ALWAYS_ASSERT(real_width-real_set_width >= 3); + bdy_time_interval = read_from_wrfbdy(nc_bdy_file,geom[0].Domain(), bdy_data_xlo,bdy_data_xhi,bdy_data_ylo,bdy_data_yhi, real_width, start_bdy_time); diff --git a/Source/Initialization/ERF_MetgridUtils.H b/Source/Initialization/ERF_MetgridUtils.H index 0b6a0a6d3..e17bf2f81 100644 --- a/Source/Initialization/ERF_MetgridUtils.H +++ b/Source/Initialization/ERF_MetgridUtils.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,17 @@ 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_debug_quiescent, + const bool metgrid_debug_isothermal, + const bool metgrid_debug_dry, + const bool metgrid_basic_linear, + const bool metgrid_use_below_sfc, + const bool metgrid_use_sfc, + const bool metgrid_retain_sfc, + const amrex::Real metgrid_proximity, + const int metgrid_order, + const int metgrid_metgrid_force_sfc_k, const amrex::Real l_rdOcp, amrex::Box& tbxc, amrex::Box& tbxu, @@ -64,29 +72,34 @@ 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, - amrex::Vector& theta_fab, - amrex::Vector& mxrat_fab, - amrex::Vector>& fabs_for_bcs, + const amrex::Vector& NC_pres_fab, + amrex::FArrayBox& p_interp_fab, + amrex::FArrayBox& t_interp_fab, + amrex::FArrayBox& theta_fab, + amrex::FArrayBox& mxrat_fab, + amrex::Vector>& fabs_for_bcs_xlo, + amrex::Vector>& fabs_for_bcs_xhi, + amrex::Vector>& fabs_for_bcs_ylo, + amrex::Vector>& fabs_for_bcs_yhi, const amrex::Array4& mask_c_arr, const amrex::Array4& mask_u_arr, const amrex::Array4& mask_v_arr); void -init_msfs_from_metgrid (amrex::FArrayBox& msfu_fab, +init_msfs_from_metgrid (const bool metgrid_debug_msf, + 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); void init_base_state_from_metgrid (const bool use_moisture, + const bool metgrid_debug_psfc, const amrex::Real l_rdOcp, const amrex::Box& valid_bx, const amrex::Vector& flag_psfc, @@ -96,21 +109,656 @@ init_base_state_from_metgrid (const bool use_moisture, amrex::FArrayBox& pi_hse_fab, amrex::FArrayBox& th_hse_fab, amrex::FArrayBox& z_phys_cc_fab, - const amrex::Vector& NC_ght_fab, - const amrex::Vector& NC_psfc_fab, - amrex::Vector>& fabs_for_bcs); + const amrex::Vector& NC_psfc_fab); AMREX_FORCE_INLINE AMREX_GPU_DEVICE -amrex::Real -interpolate_column_metgrid (const int& i, +void +lagrange_interp (const int& order, + amrex::Real* x, + amrex::Real* y, + 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 n = 1.0; + amrex::Real d = 1.0; + for (int k=0; k <= order; k++) { + if (k == i) continue; + n *= new_x-x[k]; + d *= x[i]-x[k]; + } + if (d != 0.0) { + Px += y[i]*n/d; + } + } + new_y = Px; +} + +AMREX_FORCE_INLINE +AMREX_GPU_DEVICE +void +lagrange_setup (char var_type, + const bool& exp_interp, + const int& orig_n, + const int& new_n, + const int& order, + amrex::Real* orig_x_z, + amrex::Real* orig_x_p, + amrex::Real* orig_y, + amrex::Real* new_x_z, + amrex::Real* new_x_p, + amrex::Real* new_y) +{ + 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; + +#ifndef AMREX_USE_GPU + bool debug = false; +#endif + + for (int new_k=0; new_k < new_n; new_k++) { +#ifndef AMREX_USE_GPU + if (debug) amrex::Print() << "new_k=" << new_k; +#endif + // 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::GpuArray orig_x_sub; + amrex::Real* orig_x_sub_p = orig_x_sub.data(); + if (exp_interp) { + new_x = new_x_p[new_k]; + orig_x_sub_p[0] = orig_x_p[ksta]; + orig_x_sub_p[1] = orig_x_p[kend]; + } else { + new_x = new_x_z[new_k]; + orig_x_sub_p[0] = orig_x_z[ksta]; + orig_x_sub_p[1] = orig_x_z[kend]; + } + amrex::GpuArray orig_y_sub; + amrex::Real* orig_y_sub_p = orig_y_sub.data(); + orig_y_sub_p[0] = orig_y[ksta]; + orig_y_sub_p[1] = orig_y[kend]; + lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, 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; +#ifndef AMREX_USE_GPU + int ksize = kend-ksta; + 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; +#endif + amrex::Real new_x; + amrex::GpuArray orig_x_sub; + amrex::GpuArray orig_y_sub; + amrex::Real* orig_x_sub_p = orig_x_sub.data(); + amrex::Real* orig_y_sub_p = orig_y_sub.data(); + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; } + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; } + } + for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; } +#ifndef AMREX_USE_GPU + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k]; + amrex::Print() << "]" << std::endl; + } +#endif + lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_l); + } + { + int ksta = kl-order/2; + int kend = ksta+order; +#ifndef AMREX_USE_GPU + int ksize = kend-ksta; + 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; +#endif + amrex::Real new_x; + amrex::GpuArray orig_x_sub; + amrex::GpuArray orig_y_sub; + amrex::Real* orig_x_sub_p = orig_x_sub.data(); + amrex::Real* orig_y_sub_p = orig_y_sub.data(); + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; } + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; } + } + for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; } +#ifndef AMREX_USE_GPU + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k]; + amrex::Print() << "]" << std::endl; + } +#endif + lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y_r); + } + new_y[new_k] = 0.5*(new_y_l+new_y_r); +#ifndef AMREX_USE_GPU + if (debug) amrex::Print() << " new_y=" << new_y[new_k] << std::endl; +#endif + } else if ((kl-(order/2-1) >= 0) && (kr+order/2 <= orig_n-1)) { + int ksta = kl-(order/2-1); + int kend = ksta+order; +#ifndef AMREX_USE_GPU + int ksize = kend-ksta; + 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; +#endif + amrex::Real new_x; + amrex::GpuArray orig_x_sub; + amrex::GpuArray orig_y_sub; + amrex::Real* orig_x_sub_p = orig_x_sub.data(); + amrex::Real* orig_y_sub_p = orig_y_sub.data(); + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; } + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; } + } + for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; } +#ifndef AMREX_USE_GPU + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k]; + amrex::Print() << "]" << std::endl; + } +#endif + lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]); +#ifndef AMREX_USE_GPU + if (debug) amrex::Print() << " new_y=" << new_y[new_k] << std::endl; +#endif + } else if ((kl-order/2 >= 0) && (kr+(order/2-1) <= orig_n-1)) { + int ksta = kl-order/2; + int kend = ksta+order; +#ifndef AMREX_USE_GPU + int ksize = kend-ksta; + 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; +#endif + amrex::Real new_x; + amrex::GpuArray orig_x_sub; + amrex::GpuArray orig_y_sub; + amrex::Real* orig_x_sub_p = orig_x_sub.data(); + amrex::Real* orig_y_sub_p = orig_y_sub.data(); + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; } + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; } + } + for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; } +#ifndef AMREX_USE_GPU + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k]; + amrex::Print() << "]" << std::endl; + } +#endif + lagrange_interp(order, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]); +#ifndef AMREX_USE_GPU + if (debug) amrex::Print() << " new_y=" << new_y[new_k] << std::endl; +#endif + } else { + amrex::Abort("metgrid initialization, lost in lagrange_setup (even order)"); + } + } else { + // Linear interpolation. + int ksta = kl; + int kend = kr; +#ifndef AMREX_USE_GPU + int ksize = kend-ksta; + 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; +#endif + amrex::Real new_x; + amrex::GpuArray orig_x_sub; + amrex::GpuArray orig_y_sub; + amrex::Real* orig_x_sub_p = orig_x_sub.data(); + amrex::Real* orig_y_sub_p = orig_y_sub.data(); + if (exp_interp) { + new_x = new_x_p[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_p[k]; } + } else { + new_x = new_x_z[new_k]; + for (int k=ksta; k <= kend; k++) { orig_x_sub_p[k-ksta] = orig_x_z[k]; } + } + for (int k=ksta; k <= kend; k++) { orig_y_sub_p[k-ksta] = orig_y[k]; } +#ifndef AMREX_USE_GPU + if (debug) { + amrex::Print() << " orig_x_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_x_sub_p[k]; + amrex::Print() << "]" << std::endl; + amrex::Print() << " orig_y_sub = ["; + for (int k=0; k < ksize; k++) amrex::Print() << " " << orig_y_sub_p[k]; + amrex::Print() << "]" << std::endl; + } +#endif + lagrange_interp(1, orig_x_sub_p, orig_y_sub_p, new_x, new_y[new_k]); +#ifndef AMREX_USE_GPU + if (debug) amrex::Print() << " new_y=" << new_y[new_k] << std::endl; +#endif + } + } +} + +AMREX_FORCE_INLINE +AMREX_GPU_DEVICE +void +calc_p_isothermal (const amrex::Real& z, + amrex::Real& p) +{ + p = p_0*exp(-CONST_GRAV*z/(290.0*R_d)); +} + +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& itime, + 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_xlo, + const amrex::Array4& bc_data_xhi, + const amrex::Array4& bc_data_ylo, + const amrex::Array4& bc_data_yhi, + const amrex::Box& bx_xlo, + const amrex::Box& bx_xhi, + const amrex::Box& bx_ylo, + const amrex::Box& bx_yhi, + 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 accommodate 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_ASSERT(kmax_orig < 256); + AMREX_ASSERT(kmax_new < 256); + + amrex::GpuArray new_z; + amrex::GpuArray new_p; + amrex::GpuArray new_data; + amrex::Real* new_z_p = new_z.data(); + amrex::Real* new_p_p = new_p.data(); + amrex::Real* new_data_p = new_data.data(); + for (int k=0; k < kmax_new; k++) { + if (stag == 'X') { + new_z_p[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_p[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_p[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)); + } + calc_p_isothermal(new_z_p[k], new_p_p[k]); + } + + amrex::GpuArray orig_z; + amrex::Real* orig_z_p = orig_z.data(); + for (int k=0; k < kmax_orig; k++) { + if (stag == 'M') { + orig_z_p[k] = orig_z_full(i,j,k); + } else if (stag == 'X') { + if (i == 0) { + orig_z_p[k] = orig_z_full(i,j,k); + } else if (i == imax_orig) { + orig_z_p[k] = orig_z_full(imax_orig-1,j,k); + } else { + orig_z_p[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_p[k] = orig_z_full(i,j,k); + } else if (j == jmax_orig) { + orig_z_p[k] = orig_z_full(i,jmax_orig-1,k); + } else { + orig_z_p[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_p[k] > orig_z_p[0]) { + k_above_sfc = k; + break; + } + } + + int zap = 0; + int zap_below = 0; + int zap_above = 0; + int kend_order; + amrex::ignore_unused(zap,zap_above); + amrex::GpuArray ordered_z; + amrex::GpuArray ordered_data; + amrex::Real* ordered_z_p = ordered_z.data(); + amrex::Real* ordered_data_p = ordered_data.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_p[count] = orig_z_p[k]; + ordered_data_p[count] = 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, Pl; + calc_p_isothermal(orig_z_p[0], Pu); + calc_p_isothermal(ordered_z_p[count-1], Pl); + if (Pl-Pu < metgrid_proximity) { + count--; + zap = 1; + zap_below = 1; + } + + // Insert the surface level. + ordered_z_p[count] = orig_z_p[0]; + ordered_data_p[count] = 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_p[k] > new_z_p[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. + calc_p_isothermal(orig_z_p[knext], Pu); + calc_p_isothermal(ordered_z_p[count-1], Pl); + 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_p[count] = orig_z_p[k]; + ordered_data_p[count] = orig_data(i,j,k); + count++; + } + + kend_order = count; + } else { + // The surface is the lowest level in the origin data. + + // Insert the surface. + ordered_z_p[0] = orig_z[0]; + ordered_data_p[0] = 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_p[k] > new_z_p[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, Pl; + calc_p_isothermal(orig_z_p[k], Pu); + calc_p_isothermal(ordered_z_p[count-1], Pl); + if (Pl-Pu < metgrid_proximity) { + zap++; + zap_above++; + continue; + } + ordered_z_p[count] = orig_z_p[k]; + ordered_data_p[count] = orig_data(i,j,k); + count++; + } + kend_order = count; + } + + int ksta(0), kend(0); + if (metgrid_use_below_sfc && metgrid_use_sfc) { + // Use all levels. + ksta = 0; + kend = ksta+kend_order-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_p[k] == orig_z_p[0]) { + ksfc = k; + break; + } + } + for (int k=ksfc; k < kmax_orig-1; k++) { + ordered_z_p[k] = ordered_z_p[k+1]; + ordered_data_p[k] = ordered_data_p[k+1]; + } + ksta = 0; + kend = ksta+kend_order-2; + } 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_p[kcount] == orig_z_p[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::GpuArray ordered_p; + amrex::Real* ordered_p_p = ordered_p.data(); + for (int k=0; k < kend_order; k++) { + calc_p_isothermal(ordered_z_p[k], ordered_p_p[k]); + } + + int new_n = 0; + int zap_final = 0; + amrex::GpuArray final_z; + amrex::GpuArray final_p; + amrex::GpuArray final_data; + amrex::Real* final_z_p = final_z.data(); + amrex::Real* final_p_p = final_p.data(); + amrex::Real* final_data_p = final_data.data(); + final_z_p[0] = ordered_z[ksta]; + final_p_p[0] = ordered_p[ksta]; + final_data_p[0] = ordered_data[ksta]; + for (int k=ksta+1; k <= kend; k++) { + if ((final_p_p[new_n]-ordered_p_p[k]) < metgrid_proximity) { + zap_final++; + } else { + new_n++; + final_z_p[new_n] = ordered_z_p[k]; + final_p_p[new_n] = ordered_p_p[k]; + final_data_p[new_n] = ordered_data_p[k]; + } + } + kend -= zap_final; + + // Call the interpolator. + lagrange_setup(var_type, + exp_interp, + kend-ksta, + kmax_new, + metgrid_order, + final_z_p, + final_p_p, + final_data_p, + new_z_p, + new_p_p, + new_data_p); + + // 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. + for (int k=0; k < kmax_new; k++) { + if (mask(i,j,k) && update_bc_data && bx_xlo.contains(i,j,k)) bc_data_xlo(i,j,k,0) = new_data[k]; + if (mask(i,j,k) && update_bc_data && bx_xhi.contains(i,j,k)) bc_data_xhi(i,j,k,0) = new_data[k]; + if (mask(i,j,k) && update_bc_data && bx_ylo.contains(i,j,k)) bc_data_ylo(i,j,k,0) = new_data[k]; + if (mask(i,j,k) && update_bc_data && bx_yhi.contains(i,j,k)) bc_data_yhi(i,j,k,0) = new_data[k]; + if (itime == 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; diff --git a/Submodules/AMReX b/Submodules/AMReX index 2f10b41e3..f7f3ee9f5 160000 --- a/Submodules/AMReX +++ b/Submodules/AMReX @@ -1 +1 @@ -Subproject commit 2f10b41e39747e1c38e01c39e08251e20be11309 +Subproject commit f7f3ee9f57170e227af1fee1c6b95e5b41a45af8