Skip to content

Updates to the MetGrid initialization and forcing pathway #2806

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wiersema1 wants to merge 1 commit into
base: development
Choose a base branch
from
Draft

Updates to the MetGrid initialization and forcing pathway #2806

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions Source/ERF_IndexDefines.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,22 @@ namespace MetGridBdyVars {
};
}

namespace MetGridTmpSrcVars {
enum {
Theta = 0,
QV = 1,
NumTypes
};
}

namespace MetGridTmpDstVars {
enum {
P = 0,
T = 1,
NumTypes
};
}

namespace Vars {
enum {
cons = 0,
Expand Down
45 changes: 36 additions & 9 deletions Source/ERF_Tagging.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ using namespace amrex;
#ifdef ERF_USE_NETCDF
Box read_subdomain_from_wrfinput (int lev, const std::string& fname, int& ratio);
Real read_start_time_from_wrfinput (int lev, const std::string& fname);
Box read_subdomain_from_metgrid (int lev, const std::string& fname, int& ratio, int& klo, int& khi);
#endif

void
Expand All @@ -23,41 +24,62 @@ ERF::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)
const int tagval = TagBox::SET;

#ifdef ERF_USE_NETCDF
if (solverChoice.init_type == InitType::WRFInput) {
if ((solverChoice.init_type == InitType::WRFInput) || (solverChoice.init_type == InitType::Metgrid)) {
int ratio;
Box subdomain;

if (!nc_init_file[levc+1].empty())
{
Real levc_start_time = read_start_time_from_wrfinput(levc , nc_init_file[levc ][0]);
amrex::Print() << " WRFInput time at level " << levc << " is " << levc_start_time << std::endl;
if (solverChoice.init_type == InitType::WRFInput) {
amrex::Print() << " WRFInput time at level " << levc << " is " << levc_start_time << std::endl;
} else if (solverChoice.init_type == InitType::Metgrid) {
amrex::Print() << " met_em time at level " << levc << " is " << levc_start_time << std::endl;
}

for (int isub = 0; isub < nc_init_file[levc+1].size(); isub++) {
if (!have_read_nc_init_file[levc+1][isub])
{
Real levf_start_time = read_start_time_from_wrfinput(levc+1, nc_init_file[levc+1][isub]);
amrex::Print() << " WRFInput start_time at level " << levc+1 << " is " << levf_start_time << std::endl;
if (solverChoice.init_type == InitType::WRFInput) {
amrex::Print() << " WRFInput start_time at level " << levc+1 << " is " << levf_start_time << std::endl;
} else if (solverChoice.init_type == InitType::Metgrid) {
amrex::Print() << " met_em start time at level " << levc+1 << " is " << levf_start_time << std::endl;
}

// We assume there is only one subdomain at levc; otherwise we don't know
// which one is the parent of the fine region we are trying to create
AMREX_ALWAYS_ASSERT(subdomains[levc].size() == 1);

if ( (ref_ratio[levc][2]) != 1) {
if ((solverChoice.init_type == InitType::WRFInput) && ((ref_ratio[levc][2]) != 1)) {
amrex::Abort("The ref_ratio specified in the inputs file must have 1 in the z direction; please use ref_ratio_vect rather than ref_ratio");
}

if ( levf_start_time <= (levc_start_time + t_new[levc]) ) {
amrex::Print() << " WRFInput file to read: " << nc_init_file[levc+1][isub] << std::endl;
subdomain = read_subdomain_from_wrfinput(levc, nc_init_file[levc+1][isub], ratio);
amrex::Print() << " WRFInput subdomain " << isub << " at level " << levc+1 << " is " << subdomain << std::endl;
if (solverChoice.init_type == InitType::WRFInput) {
amrex::Print() << " WRFInput file to read: " << nc_init_file[levc+1][isub] << std::endl;
subdomain = read_subdomain_from_wrfinput(levc, nc_init_file[levc+1][isub], ratio);
amrex::Print() << " WRFInput subdomain " << isub << " at level " << levc+1 << " is " << subdomain << std::endl;
} else if (solverChoice.init_type == InitType::Metgrid) {
amrex::Print() << "met_em file to read: " << nc_init_file[levc+1][0] << std::endl;
const Box& domain = geom[levc].Domain();
int klo = domain.smallEnd(2);
int khi = domain.bigEnd(2);
subdomain = read_subdomain_from_metgrid(levc, nc_init_file[levc+1][0], ratio, klo, khi);
amrex::Print() << " met_em subdomain at level " << levc+1 << " is " << subdomain << std::endl;
}

if ( (ratio != ref_ratio[levc][0]) || (ratio != ref_ratio[levc][1]) ) {
amrex::Print() << "File " << nc_init_file[levc+1][0] << " has refinement ratio = " << ratio << std::endl;
amrex::Print() << "The inputs file has refinement ratio = " << ref_ratio[levc] << std::endl;
amrex::Abort("These must be the same -- please edit your inputs file and try again.");
}

subdomain.coarsen(IntVect(ratio,ratio,1));
if (solverChoice.init_type == InitType::WRFInput) {
subdomain.coarsen(IntVect(ratio,ratio,1));
} else if (solverChoice.init_type == InitType::Metgrid) {
subdomain.coarsen(ref_ratio[levc]);
}

Box coarser_level(subdomains[levc][isub].minimalBox());
subdomain.shift(coarser_level.smallEnd());
Expand All @@ -66,7 +88,12 @@ ERF::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)
amrex::Print() << " Crse subdomain to be tagged is" << subdomain << std::endl;
}

Box new_fine(subdomain); new_fine.refine(IntVect(ratio,ratio,1));
Box new_fine(subdomain);
if (solverChoice.init_type == InitType::WRFInput) {
new_fine.refine(IntVect(ratio,ratio,1));
} else if (solverChoice.init_type == InitType::Metgrid) {
new_fine.refine(ref_ratio[levc]);
}
num_boxes_at_level[levc+1] = 1;
boxes_at_level[levc+1].push_back(new_fine);

Expand Down
106 changes: 83 additions & 23 deletions Source/IO/ERF_ReadFromMetgrid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,43 @@ using namespace amrex;

#ifdef ERF_USE_NETCDF

Box
read_subdomain_from_metgrid(int /*lev*/, const std::string& fname, int& ratio, int& klo, int& khi)
{
int is, js;
int nx, ny, nz;
if (ParallelDescriptor::IOProcessor()) {
auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
{ // Global Attributes (int)
std::vector<int> attr;
ncf.get_attr("WEST-EAST_GRID_DIMENSION" , attr); nx = attr[0]-1;
ncf.get_attr("SOUTH-NORTH_GRID_DIMENSION", attr); ny = attr[0]-1;
ncf.get_attr("BOTTOM-TOP_GRID_DIMENSION" , attr); nz = attr[0]-1;
#ifndef AMREX_USE_GPU
Print() << "Have read (nx,ny,nz) = " << nx << " " << ny << " " << nz << std::endl;
#endif
ncf.get_attr("i_parent_start", attr) ; is = attr[0]-1;
ncf.get_attr("j_parent_start", attr) ; js = attr[0]-1;
ncf.get_attr("parent_grid_ratio", attr); ratio = attr[0];
}
ncf.close();
#ifndef AMREX_USE_GPU
amrex::Print() << "Have read (parent_ilo,parent_jlo) = " << is << " " << js << std::endl;
amrex::Print() << "Have read refinement ratio = " << ratio << std::endl;
#endif
}
amrex::ParallelDescriptor::Bcast(&is ,1,amrex::ParallelDescriptor::IOProcessorNumber());
amrex::ParallelDescriptor::Bcast(&js ,1,amrex::ParallelDescriptor::IOProcessorNumber());
amrex::ParallelDescriptor::Bcast(&nx ,1,amrex::ParallelDescriptor::IOProcessorNumber());
amrex::ParallelDescriptor::Bcast(&ny ,1,amrex::ParallelDescriptor::IOProcessorNumber());
amrex::ParallelDescriptor::Bcast(&nz ,1,amrex::ParallelDescriptor::IOProcessorNumber());
amrex::ParallelDescriptor::Bcast(&ratio,1,amrex::ParallelDescriptor::IOProcessorNumber());
return Box( IntVect(ratio*is, ratio*js, klo), IntVect(ratio*is+nx-1, ratio*js+ny-1, khi) );
}

void
read_from_metgrid (int lev, const Box& domain, const std::string& fname,
read_from_metgrid (int lev, int itime,
const Box& domain, const std::string& fname,
std::string& NC_dateTime, Real& NC_epochTime,
int& flag_psfc, int& flag_msf,
int& flag_sst, int& flag_tsk, int& flag_lmask,
Expand All @@ -25,7 +60,9 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
IArrayBox& NC_lmask_iab,
Geometry& geom)
{
#ifndef AMREX_USE_GPU
Print() << "Loading header data from NetCDF file at level " << lev << std::endl;
#endif

if (ParallelDescriptor::IOProcessor()) {
auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
Expand All @@ -50,30 +87,39 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
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");
}
}
if (lev == 0) {
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) {
#ifndef AMREX_USE_GPU
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";
#endif
Abort("Domain specification error");
}
if (std::fabs(Len_y - (geom.ProbHi(1) - geom.ProbLo(1))) > tol) {
#ifndef AMREX_USE_GPU
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";
#endif
Abort("Domain specification error");
}
} // lev == 0
} // IOProc

int ioproc = ParallelDescriptor::IOProcessorNumber(); // I/O rank
ParallelDescriptor::Bcast(&NC_nx, 1, ioproc);
ParallelDescriptor::Bcast(&NC_ny, 1, ioproc);
ParallelDescriptor::Bcast(&NC_epochTime, 1, ioproc);
ParallelDescriptor::Bcast(&NC_dx, 1, ioproc);
ParallelDescriptor::Bcast(&NC_dy, 1, ioproc);

#ifndef AMREX_USE_GPU
Print() << "Loading initial data from NetCDF file at level " << lev << std::endl;
#endif

Vector<FArrayBox*> NC_fabs;
Vector<IArrayBox*> NC_iabs;
Expand All @@ -92,7 +138,9 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
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);

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 (itime == 0) {
NC_fabs.push_back(&NC_psfc_fab); NC_fnames.push_back("PSFC"); NC_fdim_types.push_back(NC_Data_Dims_Type::Time_SN_WE);
}
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);
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);
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);
Expand All @@ -102,24 +150,36 @@ read_from_metgrid (int lev, const Box& domain, const std::string& fname,
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
#ifndef AMREX_USE_GPU
Print() << "Building initial FABS from file " << fname << std::endl;
#endif
Vector<int> success; success.resize(NC_fabs.size());
BuildFABsFromNetCDFFile<FArrayBox,Real>(domain, fname, NC_fnames, NC_fdim_types, NC_fabs, success);
for (int i = 0; i < success.size(); i++) {
if (NC_fnames[i] == "PSFC" && success[i] == 1) {flag_psfc = 1;}
if (NC_fnames[i] == "SST" && success[i] == 1) {flag_sst = 1;}
if (NC_fnames[i] == "SKINTEMP" && success[i] == 1) {flag_tsk = 1;}
if (NC_fnames[i] == "MAPFAC_M" && success[i] == 1) {flag_msf = 1;}
flag_psfc = (NC_fnames[i] == "PSFC" && success[i] == 1) ? 1 : 0;
flag_sst = (NC_fnames[i] == "SST" && success[i] == 1) ? 1 : 0;
flag_tsk = (NC_fnames[i] == "SKINTEMP" && success[i] == 1) ? 1 : 0;
flag_msf = (NC_fnames[i] == "MAPFAC_M" && success[i] == 1) ? 1 : 0;
}

// Read the netcdf file and fill these IABs
#ifndef AMREX_USE_GPU
Print() << "Building initial IABS from file " << fname << std::endl;
#endif
Vector<int> success_i; success_i.resize(NC_iabs.size());
BuildFABsFromNetCDFFile<IArrayBox,int>(domain, fname, NC_inames, NC_idim_types, NC_iabs, success_i);
for (int i = 0; i < success_i.size(); i++) {
if (NC_inames[i] == "LANDMASK" && success_i[i] == 1) {flag_lmask = 1;}
flag_lmask = (NC_inames[i] == "LANDMASK" && success_i[i] == 1) ? 1 : 0;
}

#ifndef AMREX_USE_GPU
Print() << " flag_psfc: " << flag_psfc << std::endl;
Print() << " flag_sst: " << flag_sst << std::endl;
Print() << " flag_tsk: " << flag_tsk << std::endl;
Print() << " flag_msf: " << flag_msf << std::endl;
Print() << " flag_lmask: " << flag_lmask << std::endl;
#endif

// TODO: FIND OUT IF WE NEED TO DIVIDE VELS BY MAPFAC
//
// Convert the velocities using the map factors
Expand Down
Loading
Loading