diff --git a/Source/ERF.H b/Source/ERF.H index 89e16009a..a1aee5aa2 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -493,6 +493,8 @@ private: amrex::Vector> bdy_data_yhi; amrex::Real bdy_time_interval; + amrex::Real Latitude; + amrex::Real Longitude; #endif // ERF_USE_NETCDF // Struct for working with the sounding data we take as an input diff --git a/Source/IO/NCWpsFile.H b/Source/IO/NCWpsFile.H index fd949ca9e..be8fc9577 100644 --- a/Source/IO/NCWpsFile.H +++ b/Source/IO/NCWpsFile.H @@ -223,6 +223,10 @@ void ReadNetCDFFile (const std::string& fname, amrex::Vector names, template void fill_fab_from_arrays (int iv, + amrex::Real& Latitude, + amrex::Real& Longitude, + std::string& Lat_var_name, + std::string& Lon_var_name, amrex::Vector>& nc_arrays, const std::string& var_name, NC_Data_Dims_Type& NC_dim_type, @@ -280,6 +284,9 @@ fill_fab_from_arrays (int iv, int i = n - k*(ns2*ns3) - (j-joff) * ns3 + ioff; fab_arr(i,j,k,0) = static_cast(*(nc_arrays[iv].get_data()+n)); } + + if (var_name == Lat_var_name) Latitude = fab_arr(0,0,0); + if (var_name == Lon_var_name) Longitude = fab_arr(0,0,0); } /** @@ -293,6 +300,10 @@ fill_fab_from_arrays (int iv, template void BuildFABsFromNetCDFFile (const amrex::Box& domain, + amrex::Real& Latitude, + amrex::Real& Longitude, + std::string& Lat_var_name, + std::string& Lon_var_name, const std::string &fname, amrex::Vector nc_var_names, amrex::Vector NC_dim_types, @@ -311,7 +322,10 @@ BuildFABsFromNetCDFFile (const amrex::Box& domain, { FAB tmp; if (amrex::ParallelDescriptor::IOProcessor()) { - fill_fab_from_arrays(iv, nc_arrays, nc_var_names[iv], NC_dim_types[iv], tmp); + fill_fab_from_arrays(iv, Latitude, Longitude, + Lat_var_name, Lon_var_name, + nc_arrays, nc_var_names[iv], + NC_dim_types[iv], tmp); } int ncomp = tmp.nComp(); diff --git a/Source/IO/ReadFromMetgrid.cpp b/Source/IO/ReadFromMetgrid.cpp index 0cd0939da..08cceae6d 100644 --- a/Source/IO/ReadFromMetgrid.cpp +++ b/Source/IO/ReadFromMetgrid.cpp @@ -1,6 +1,7 @@ #include #include #include +#include using namespace amrex; @@ -19,9 +20,13 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, FArrayBox& NC_hgt_fab, FArrayBox& NC_psfc_fab, FArrayBox& NC_msfu_fab, FArrayBox& NC_msfv_fab, FArrayBox& NC_msfm_fab, FArrayBox& NC_sst_fab, - IArrayBox& NC_lmask_iab) + FArrayBox& NC_LAT_fab, FArrayBox& NC_LON_fab, + IArrayBox& NC_lmask_iab, + Real& Latitude, + Real& Longitude, + Geometry& geom) { - Print() << "Loading initial data from NetCDF file at level " << lev << std::endl; + Print() << "Loading header data from NetCDF file at level " << lev << std::endl; if (ParallelDescriptor::IOProcessor()) { auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); @@ -60,6 +65,23 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, ncf.get_attr("DY", attr); NC_dy = attr[0]; } ncf.close(); + + // Verify the inputs geometry matches what the NETCDF file has + Real tol = 1.0e-3; + Real Len_x = NC_dx * Real(NC_nx-1); + Real Len_y = NC_dy * Real(NC_ny-1); + if (std::fabs(Len_x - (geom.ProbHi(0) - geom.ProbLo(0))) > tol) { + Print() << "X problem extent " << (geom.ProbHi(0) - geom.ProbLo(0)) << " does not match NETCDF file " + << Len_x << "!\n"; + Print() << "dx: " << NC_dx << ' ' << "Nx: " << NC_nx-1 << "\n"; + Abort("Domain specification error"); + } + if (std::fabs(Len_y - (geom.ProbHi(1) - geom.ProbLo(1))) > tol) { + Print() << "Y problem extent " << (geom.ProbHi(1) - geom.ProbLo(1)) << " does not match NETCDF file " + << Len_y << "!\n"; + Print() << "dy: " << NC_dy << ' ' << "Ny: " << NC_ny-1 << "\n"; + Abort("Domain specification error"); + } } int ioproc = ParallelDescriptor::IOProcessorNumber(); // I/O rank ParallelDescriptor::Bcast(&flag_psfc, 1, ioproc); @@ -75,6 +97,8 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, ParallelDescriptor::Bcast(&NC_dx, 1, ioproc); ParallelDescriptor::Bcast(&NC_dy, 1, ioproc); + Print() << "Loading initial data from NetCDF file at level " << lev << std::endl; + Vector NC_fabs; Vector NC_iabs; Vector NC_fnames; @@ -88,6 +112,8 @@ 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); 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); } @@ -99,12 +125,18 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname, 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"; Print() << "Building initial FABS from file " << fname << std::endl; - BuildFABsFromNetCDFFile(domain, fname, NC_fnames, NC_fdim_types, NC_fabs); + BuildFABsFromNetCDFFile(domain, Latitude, Longitude, + Lat_var_name, Lon_var_name, + fname, NC_fnames, NC_fdim_types, NC_fabs); // Read the netcdf file and fill these IABs Print() << "Building initial IABS from file " << fname << std::endl; - BuildFABsFromNetCDFFile(domain, fname, NC_inames, NC_idim_types, NC_iabs); + BuildFABsFromNetCDFFile(domain, Latitude, Longitude, + Lat_var_name, Lon_var_name, + fname, NC_inames, NC_idim_types, NC_iabs); // TODO: FIND OUT IF WE NEED TO DIVIDE VELS BY MAPFAC // diff --git a/Source/IO/ReadFromWRFInput.cpp b/Source/IO/ReadFromWRFInput.cpp index 28ff84e34..30dd1fc63 100644 --- a/Source/IO/ReadFromWRFInput.cpp +++ b/Source/IO/ReadFromWRFInput.cpp @@ -25,43 +25,100 @@ read_from_wrfinput (int lev, FArrayBox& NC_PHB_fab, FArrayBox& NC_ALB_fab, FArrayBox& NC_PB_fab, - MoistureType moisture_type) + FArrayBox& NC_LAT_fab, + FArrayBox& NC_LON_fab, + MoistureType moisture_type, + Real& Latitude, + Real& Longitude, + Geometry& geom) { + Print() << "Loading header data from NetCDF file at level " << lev << std::endl; + + if (ParallelDescriptor::IOProcessor()) { + int NC_nx, NC_ny; + Real NC_dx, NC_dy; + Real NC_epochTime; + std::string NC_dateTime; + auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4); + { // Global Attributes (int) + std::vector attr; + ncf.get_attr("WEST-EAST_GRID_DIMENSION", attr); NC_nx = attr[0]; + ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); NC_ny = attr[0]; + } + { // Global Attributes (Real) + std::vector attr; + ncf.get_attr("DX", attr); NC_dx = attr[0]; + ncf.get_attr("DY", attr); NC_dy = attr[0]; + } + { // Global Attributes (string) + NC_dateTime = ncf.get_attr("SIMULATION_START_DATE")+"UTC"; + const std::string dateTimeFormat = "%Y-%m-%d_%H:%M:%S%Z"; + NC_epochTime = getEpochTime(NC_dateTime, dateTimeFormat); + } + ncf.close(); + + amrex::ignore_unused(NC_dateTime); amrex::ignore_unused(NC_epochTime); + + // Verify the inputs geometry matches what the NETCDF file has + Real tol = 1.0e-3; + Real Len_x = NC_dx * Real(NC_nx-1); + Real Len_y = NC_dy * Real(NC_ny-1); + if (std::fabs(Len_x - (geom.ProbHi(0) - geom.ProbLo(0))) > tol) { + Print() << "X problem extent " << (geom.ProbHi(0) - geom.ProbLo(0)) << " does not match NETCDF file " + << Len_x << "!\n"; + Print() << "dx: " << NC_dx << ' ' << "Nx: " << NC_nx-1 << "\n"; + Abort("Domain specification error"); + } + if (std::fabs(Len_y - (geom.ProbHi(1) - geom.ProbLo(1))) > tol) { + Print() << "Y problem extent " << (geom.ProbHi(1) - geom.ProbLo(1)) << " does not match NETCDF file " + << Len_y << "!\n"; + Print() << "dy: " << NC_dy << ' ' << "Ny: " << NC_ny-1 << "\n"; + Abort("Domain specification error"); + } + } + Print() << "Loading initial data from NetCDF file at level " << lev << std::endl; Vector NC_fabs; Vector NC_names; Vector NC_dim_types; - NC_fabs.push_back(&NC_xvel_fab); NC_names.push_back("U"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 0 - NC_fabs.push_back(&NC_yvel_fab); NC_names.push_back("V"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 1 - NC_fabs.push_back(&NC_zvel_fab); NC_names.push_back("W"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 2 - NC_fabs.push_back(&NC_rho_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 3 - NC_fabs.push_back(&NC_rhop_fab), NC_names.push_back("AL"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 4 - NC_fabs.push_back(&NC_rhotheta_fab); NC_names.push_back("T"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 5 - NC_fabs.push_back(&NC_PH_fab); NC_names.push_back("PH"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 6 - NC_fabs.push_back(&NC_PHB_fab); NC_names.push_back("PHB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 7 - NC_fabs.push_back(&NC_PB_fab); NC_names.push_back("PB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 9 - NC_fabs.push_back(&NC_P_fab); NC_names.push_back("P"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 10 - NC_fabs.push_back(&NC_ALB_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 11 - NC_fabs.push_back(&NC_MUB_fab); NC_names.push_back("MUB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 12 - NC_fabs.push_back(&NC_MSFU_fab); NC_names.push_back("MAPFAC_UY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 13 - NC_fabs.push_back(&NC_MSFV_fab); NC_names.push_back("MAPFAC_VY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 14 - NC_fabs.push_back(&NC_MSFM_fab); NC_names.push_back("MAPFAC_MY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 15 - NC_fabs.push_back(&NC_SST_fab); NC_names.push_back("SST"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 16 - NC_fabs.push_back(&NC_C1H_fab); NC_names.push_back("C1H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 17 - NC_fabs.push_back(&NC_C2H_fab); NC_names.push_back("C2H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 18 - NC_fabs.push_back(&NC_RDNW_fab); NC_names.push_back("RDNW"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 19 + NC_fabs.push_back(&NC_xvel_fab); NC_names.push_back("U"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 0 + NC_fabs.push_back(&NC_yvel_fab); NC_names.push_back("V"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 1 + NC_fabs.push_back(&NC_zvel_fab); NC_names.push_back("W"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 2 + NC_fabs.push_back(&NC_rho_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 3 + NC_fabs.push_back(&NC_rhop_fab), NC_names.push_back("AL"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 4 + NC_fabs.push_back(&NC_rhotheta_fab); NC_names.push_back("T"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 5 + NC_fabs.push_back(&NC_PH_fab); NC_names.push_back("PH"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 6 + NC_fabs.push_back(&NC_PHB_fab); NC_names.push_back("PHB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 7 + NC_fabs.push_back(&NC_PB_fab); NC_names.push_back("PB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 9 + NC_fabs.push_back(&NC_P_fab); NC_names.push_back("P"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 10 + NC_fabs.push_back(&NC_ALB_fab); NC_names.push_back("ALB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 11 + NC_fabs.push_back(&NC_MUB_fab); NC_names.push_back("MUB"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 12 + NC_fabs.push_back(&NC_MSFU_fab); NC_names.push_back("MAPFAC_UY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 13 + NC_fabs.push_back(&NC_MSFV_fab); NC_names.push_back("MAPFAC_VY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 14 + NC_fabs.push_back(&NC_MSFM_fab); NC_names.push_back("MAPFAC_MY"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 15 + NC_fabs.push_back(&NC_SST_fab); NC_names.push_back("SST"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 16 + NC_fabs.push_back(&NC_C1H_fab); NC_names.push_back("C1H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 17 + NC_fabs.push_back(&NC_C2H_fab); NC_names.push_back("C2H"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 18 + NC_fabs.push_back(&NC_RDNW_fab); NC_names.push_back("RDNW"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT); // 19 + NC_fabs.push_back(&NC_LAT_fab); NC_names.push_back("XLAT_V"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 20 + NC_fabs.push_back(&NC_LON_fab); NC_names.push_back("XLONG_U"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_SN_WE); // 21 if (moisture_type != MoistureType::None) { - NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 20 - NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 21 - NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 22 + NC_fabs.push_back(&NC_QVAPOR_fab); NC_names.push_back("QVAPOR"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 22 + NC_fabs.push_back(&NC_QCLOUD_fab); NC_names.push_back("QCLOUD"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 23 + NC_fabs.push_back(&NC_QRAIN_fab); NC_names.push_back("QRAIN"); NC_dim_types.push_back(NC_Data_Dims_Type::Time_BT_SN_WE); // 24 } // Read the netcdf file and fill these FABs + std::string Lat_var_name = "XLAT_V"; + std::string Lon_var_name = "XLONG_U"; Print() << "Building initial FABS from file " << fname << std::endl; - BuildFABsFromNetCDFFile(domain, fname, NC_names, NC_dim_types, NC_fabs); + BuildFABsFromNetCDFFile(domain, Latitude, Longitude, + Lat_var_name, Lon_var_name, + fname, NC_names, NC_dim_types, NC_fabs); + // // Convert the velocities using the map factors diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 4eea8f9c3..51bfe8042 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -51,6 +51,8 @@ ERF::init_from_metgrid (int lev) Vector NC_MSFV_fab; NC_MSFV_fab.resize(ntimes); Vector NC_MSFM_fab; NC_MSFM_fab.resize(ntimes); Vector NC_sst_fab; NC_sst_fab.resize (ntimes); + Vector NC_LAT_fab; NC_LAT_fab.resize (ntimes); + Vector NC_LON_fab; NC_LON_fab.resize (ntimes); // *** IArrayBox's at this level for holding mask data Vector NC_lmask_iab; NC_lmask_iab.resize(ntimes); @@ -87,7 +89,8 @@ ERF::init_from_metgrid (int lev) 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_lmask_iab[it]); + NC_sst_fab[it], NC_LAT_fab[it], NC_LON_fab[it], + NC_lmask_iab[it], Latitude, Longitude, geom[lev]); } // it // Verify that files in nc_init_file[lev] are ordered from earliest to latest. diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp index 26c5a47ac..cd6eb626e 100644 --- a/Source/Initialization/ERF_init_from_wrfinput.cpp +++ b/Source/Initialization/ERF_init_from_wrfinput.cpp @@ -31,7 +31,12 @@ read_from_wrfinput (int lev, const Box& domain, const std::string& fname, FArrayBox& NC_PHB_fab, FArrayBox& NC_ALB_fab, FArrayBox& NC_PB_fab, - MoistureType moisture_type); + FArrayBox& NC_LAT_fab, + FArrayBox& NC_LON_fab, + MoistureType moisture_type, + Real& Latitude, + Real& Longitude, + Geometry& geom); Real read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain, @@ -79,6 +84,12 @@ init_msfs_from_wrfinput (int lev, FArrayBox& msfu_fab, const Vector& NC_MSFU_fab, const Vector& NC_MSFV_fab, const Vector& NC_MSFM_fab); + +void +verify_terrain_top_boundary (const Real& z_top, + const Vector& NC_PH_fab, + const Vector& NC_PHB_fab); + void init_terrain_from_wrfinput (int lev, const Real& z_top, const Box& domain, FArrayBox& z_phys, @@ -109,28 +120,30 @@ void ERF::init_from_wrfinput (int lev) { // *** FArrayBox's at this level for holding the INITIAL data - Vector NC_xvel_fab ; NC_xvel_fab.resize(num_boxes_at_level[lev]); - Vector NC_yvel_fab ; NC_yvel_fab.resize(num_boxes_at_level[lev]); - Vector NC_zvel_fab ; NC_zvel_fab.resize(num_boxes_at_level[lev]); - Vector NC_rho_fab ; NC_rho_fab.resize(num_boxes_at_level[lev]); - Vector NC_rhop_fab ; NC_rhop_fab.resize(num_boxes_at_level[lev]); - Vector NC_rhoth_fab; NC_rhoth_fab.resize(num_boxes_at_level[lev]); - Vector NC_MUB_fab ; NC_MUB_fab.resize(num_boxes_at_level[lev]); - Vector NC_MSFU_fab ; NC_MSFU_fab.resize(num_boxes_at_level[lev]); - Vector NC_MSFV_fab ; NC_MSFV_fab.resize(num_boxes_at_level[lev]); - Vector NC_MSFM_fab ; NC_MSFM_fab.resize(num_boxes_at_level[lev]); - Vector NC_SST_fab ; NC_SST_fab.resize(num_boxes_at_level[lev]); - Vector NC_C1H_fab ; NC_C1H_fab.resize(num_boxes_at_level[lev]); - Vector NC_C2H_fab ; NC_C2H_fab.resize(num_boxes_at_level[lev]); - Vector NC_RDNW_fab ; NC_RDNW_fab.resize(num_boxes_at_level[lev]); - Vector NC_PH_fab ; NC_PH_fab.resize(num_boxes_at_level[lev]); - Vector NC_P_fab ; NC_P_fab.resize(num_boxes_at_level[lev]); - Vector NC_PHB_fab ; NC_PHB_fab.resize(num_boxes_at_level[lev]); - Vector NC_ALB_fab ; NC_ALB_fab.resize(num_boxes_at_level[lev]); - Vector NC_PB_fab ; NC_PB_fab.resize(num_boxes_at_level[lev]); + Vector NC_xvel_fab; NC_xvel_fab.resize(num_boxes_at_level[lev]); + Vector NC_yvel_fab; NC_yvel_fab.resize(num_boxes_at_level[lev]); + Vector NC_zvel_fab; NC_zvel_fab.resize(num_boxes_at_level[lev]); + Vector NC_rho_fab; NC_rho_fab.resize(num_boxes_at_level[lev]); + Vector NC_rhop_fab; NC_rhop_fab.resize(num_boxes_at_level[lev]); + Vector NC_rhoth_fab; NC_rhoth_fab.resize(num_boxes_at_level[lev]); + Vector NC_MUB_fab; NC_MUB_fab.resize(num_boxes_at_level[lev]); + Vector NC_MSFU_fab; NC_MSFU_fab.resize(num_boxes_at_level[lev]); + Vector NC_MSFV_fab; NC_MSFV_fab.resize(num_boxes_at_level[lev]); + Vector NC_MSFM_fab; NC_MSFM_fab.resize(num_boxes_at_level[lev]); + Vector NC_SST_fab; NC_SST_fab.resize(num_boxes_at_level[lev]); + Vector NC_C1H_fab; NC_C1H_fab.resize(num_boxes_at_level[lev]); + Vector NC_C2H_fab; NC_C2H_fab.resize(num_boxes_at_level[lev]); + Vector NC_RDNW_fab; NC_RDNW_fab.resize(num_boxes_at_level[lev]); + Vector NC_PH_fab; NC_PH_fab.resize(num_boxes_at_level[lev]); + Vector NC_P_fab; NC_P_fab.resize(num_boxes_at_level[lev]); + Vector NC_PHB_fab; NC_PHB_fab.resize(num_boxes_at_level[lev]); + Vector NC_ALB_fab; NC_ALB_fab.resize(num_boxes_at_level[lev]); + Vector NC_PB_fab; NC_PB_fab.resize(num_boxes_at_level[lev]); Vector NC_QVAPOR_fab; NC_QVAPOR_fab.resize(num_boxes_at_level[lev]); Vector NC_QCLOUD_fab; NC_QCLOUD_fab.resize(num_boxes_at_level[lev]); Vector NC_QRAIN_fab ; NC_QRAIN_fab.resize(num_boxes_at_level[lev]); + Vector NC_LAT_fab; NC_LAT_fab.resize(num_boxes_at_level[lev]); + Vector NC_LON_fab; NC_LON_fab.resize(num_boxes_at_level[lev]); // Print() << "Building initial FABS from file " << nc_init_file[lev][idx] << std::endl; if (nc_init_file.empty()) @@ -146,7 +159,8 @@ ERF::init_from_wrfinput (int lev) NC_QVAPOR_fab[idx], NC_QCLOUD_fab[idx], NC_QRAIN_fab[idx], NC_PH_fab[idx] , NC_P_fab[idx] , NC_PHB_fab[idx] , NC_ALB_fab[idx] , NC_PB_fab[idx] , - solverChoice.moisture_type); + NC_LAT_fab[idx] , NC_LON_fab[idx], + solverChoice.moisture_type, Latitude, Longitude, geom[lev]); } auto& lev_new = vars_new[lev]; @@ -190,6 +204,10 @@ ERF::init_from_wrfinput (int lev) const Real& z_top = geom[lev].ProbHi(2); if (solverChoice.use_terrain) { + if (ParallelDescriptor::IOProcessor()) { + verify_terrain_top_boundary(z_top, NC_PH_fab, NC_PHB_fab); + } + std::unique_ptr& z_phys = z_phys_nd[lev]; for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { @@ -451,6 +469,56 @@ init_base_state_from_wrfinput (int lev, } // idx } +/** + * Helper function for verifying the top boundary is valid. + * + * @param z_top Real user specified top boundary + * @param NC_PH_fab Vector of FArrayBox objects storing WRF terrain coordinate data (PH) + * @param NC_PHB_fab Vector of FArrayBox objects storing WRF terrain coordinate data (PHB) + */ +void +verify_terrain_top_boundary (const Real& z_top, + const Vector& NC_PH_fab, + const Vector& NC_PHB_fab) +{ + int nboxes = NC_PH_fab.size(); + for (int idx = 0; idx < nboxes; idx++) { + Real zmin = 1.0e16; + Real zmax = -1.0e16; + + Box Fab2dBox (NC_PHB_fab[idx].box()); Fab2dBox.setSmall(2,Fab2dBox.bigEnd(2)); + + Box nodal_box = amrex::surroundingNodes(NC_PHB_fab[idx].box()); + int ilo = nodal_box.smallEnd()[0]; + int ihi = nodal_box.bigEnd()[0]; + int jlo = nodal_box.smallEnd()[1]; + int jhi = nodal_box.bigEnd()[1]; + + auto const& phb = NC_PHB_fab[idx].const_array(); + auto const& ph = NC_PH_fab[idx].const_array(); + + ParallelFor(Fab2dBox, [=,&zmin,&zmax] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + int ii = std::max(std::min(i,ihi-1),ilo+1); + int jj = std::max(std::min(j,jhi-1),jlo+1); + Real z_calc = 0.25 * ( ph (ii,jj ,k) + ph (ii-1,jj ,k) + + ph (ii,jj-1,k) + ph (ii-1,jj-1,k) + + phb(ii,jj ,k) + phb(ii-1,jj ,k) + + phb(ii,jj-1,k) + phb(ii-1,jj-1,k) ) / CONST_GRAV; + zmin = amrex::min(zmin,z_calc); + zmax = amrex::max(zmax,z_calc); + }); + + Real tol = 1.0e-2; + if (std::fabs((z_top-zmin)/zmin) > tol && + std::fabs((z_top-zmax)/zmax) > tol ) { + Print() << "Z problem extent " << z_top << " does not match NETCDF file min " + << zmin << " and max " << zmax << "!\n"; + Abort("Domain specification error"); + } + } +} + /** * Helper function for initializing terrain coordinates from a WRF dataset. * @@ -466,8 +534,7 @@ init_terrain_from_wrfinput (int lev, const Real& z_top, const Vector& NC_PHB_fab) { int nboxes = NC_PH_fab.size(); - for (int idx = 0; idx < nboxes; idx++) - { + for (int idx = 0; idx < nboxes; idx++) { // This copies from NC_zphys on z-faces to z_phys_nd on nodes const Array4& z_arr = z_phys.array(); const Array4& nc_phb_arr = NC_PHB_fab[idx].const_array(); diff --git a/Source/Initialization/Metgrid_utils.H b/Source/Initialization/Metgrid_utils.H index 82d95dbe3..b3c38da7d 100644 --- a/Source/Initialization/Metgrid_utils.H +++ b/Source/Initialization/Metgrid_utils.H @@ -36,7 +36,12 @@ read_from_metgrid (int lev, amrex::FArrayBox& NC_msfv_fab, amrex::FArrayBox& NC_msfm_fab, amrex::FArrayBox& NC_sst_fab, - amrex::IArrayBox& NC_lmask_iab); + amrex::FArrayBox& NC_LAT_fab, + amrex::FArrayBox& NC_LON_fab, + amrex::IArrayBox& NC_lmask_iab, + amrex::Real& Latitude, + amrex::Real& Longitude, + amrex::Geometry& geom);