379 amrex::Vector<amrex::Real> m_xterrain,m_yterrain,m_zterrain;
381 int nx = 0;
int ny = 0;
383 if (amrex::ParallelDescriptor::IOProcessor()) {
385 amrex::Print()<<
"Reading terrain file: "<< fname<< std::endl;
386 std::ifstream file(fname);
388 if (!file.is_open()) {
389 amrex::Abort(
"Error: Could not open the file " + fname+
"\n");
393 if (file.peek() == std::ifstream::traits_type::eof()) {
394 amrex::Abort(
"Error: The file " + fname+
" is empty.\n");
402 file >> lon_min >> lat_min;
404 amrex::Error(
"The value of longitude for entry in the first line in " + fname
405 +
" should not exceed amrex::Real(180.) It is " + std::to_string(lon_min));
408 amrex::Error(
"The value of latitude for entry in the first line in " + fname
409 +
" should not exceed amrex::Real(90.) It is " + std::to_string(lat_min));
415 while (file >> value1 >> value2 >> value3) {
416 m_xterrain.push_back(value1);
418 m_yterrain.push_back(value2);
420 m_zterrain.push_back(value3);
423 AMREX_ASSERT(m_xterrain.size() ==
static_cast<long int>(nx*ny));
424 AMREX_ASSERT(m_yterrain.size() ==
static_cast<long int>(ny));
425 AMREX_ASSERT(m_zterrain.size() ==
static_cast<long int>(nx*ny));
429 nx = erf_get_single_value<int>(file,cnt); cnt++;
430 ny = erf_get_single_value<int>(file,cnt); cnt++;
431 amrex::Print()<<
"Expecting " << nx <<
" values of x, " <<
432 ny <<
" values of y, and " <<
433 nx*ny <<
" values of z" << std::endl;
436 m_xterrain.resize(nx);
437 m_yterrain.resize(ny);
438 m_zterrain.resize(nx * ny);
439 for (
int n = 0; n < nx; n++) {
440 m_xterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
443 for (
int n = 0; n < ny; n++) {
444 m_yterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
447 for (
int n = 0; n < nx * ny; n++) {
448 m_zterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
457 amrex::ParallelDescriptor::Bcast(&nx,1,amrex::ParallelDescriptor::IOProcessorNumber());
458 amrex::ParallelDescriptor::Bcast(&ny,1,amrex::ParallelDescriptor::IOProcessorNumber());
461 int nx_vals = is_usgs ? (nx * ny) : nx;
463 m_xterrain.resize(nx_vals);
464 m_yterrain.resize(ny);
465 m_zterrain.resize(nz);
467 amrex::ParallelDescriptor::Bcast(m_xterrain.data(),nx_vals,amrex::ParallelDescriptor::IOProcessorNumber());
468 amrex::ParallelDescriptor::Bcast(m_yterrain.data(),ny,amrex::ParallelDescriptor::IOProcessorNumber());
469 amrex::ParallelDescriptor::Bcast(m_zterrain.data(),nz,amrex::ParallelDescriptor::IOProcessorNumber());
472 amrex::Gpu::DeviceVector<amrex::Real> d_xterrain(nx_vals),d_yterrain(ny),d_zterrain(nz);
473 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin());
474 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin());
475 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin());
481 auto dx = geom.CellSizeArray();
482 auto ProbLoArr = geom.ProbLoArray();
484 int ilo = geom.Domain().smallEnd(0);
485 int jlo = geom.Domain().smallEnd(1);
486 int klo = geom.Domain().smallEnd(2);
487 int ihi = geom.Domain().bigEnd(0) + 1;
488 int jhi = geom.Domain().bigEnd(1) + 1;
490 amrex::Box zbx = terrain_fab.box();
491 amrex::Array4<amrex::Real>
const& z_arr = terrain_fab.array();
496 int ii = amrex::min(amrex::max(i,ilo),ihi);
497 int jj = amrex::min(amrex::max(j,jlo),jhi);
503 int ind11, ind12, ind21, ind22;
506 int iindex_terrain=-1;
507 int jindex_terrain=-1;
513 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
515 jindex_terrain = it-1;
518 if (jindex_terrain == -1) {
519 jindex_terrain = ny-1;
521 jindex_terrain = amrex::min(amrex::max(jindex_terrain,0), ny-2);
523 int gstart = (jindex_terrain )*nx;
524 int gend = (jindex_terrain+1)*nx-1;
525 for (
int it = gstart; it <= gend && iindex_terrain == -1; it++) {
527 iindex_terrain = it-gstart-1;
531 iindex_terrain = amrex::min(amrex::max(iindex_terrain,0), nx-2);
534 ind11 = jindex_terrain*nx + iindex_terrain;
541 y1 = d_yt[jindex_terrain];
542 y2 = d_yt[jindex_terrain+1];
551 z_arr(i,j,klo) = w_11*d_zt[ind11] + w_12*d_zt[ind12] + w_21*d_zt[ind21] + w_22*d_zt[ind22];
555 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
557 jindex_terrain = it-1;
560 if (jindex_terrain == -1) {
561 jindex_terrain = ny-1;
564 for (
int it = 0; it < nx && iindex_terrain == -1; it++) {
566 iindex_terrain = it-1;
569 if (iindex_terrain == -1) {
570 iindex_terrain = nx-1;
574 x1 = d_xt[iindex_terrain];
575 x2 = d_xt[iindex_terrain+1];
576 y1 = d_yt[jindex_terrain];
577 y2 = d_yt[jindex_terrain+1];
583 ind11 = iindex_terrain * ny + jindex_terrain;
584 ind21 = std::min(iindex_terrain+1,nx-1) * ny + jindex_terrain;
586 ind12 = iindex_terrain * ny + std::min(jindex_terrain+1,ny-1);
587 ind22 = std::min(iindex_terrain+1,nx-1) * ny + std::min(jindex_terrain+1,ny-1);
593 ind11 = jindex_terrain * nx + iindex_terrain;
594 ind12 = std::min(jindex_terrain+1,ny-1) * nx + iindex_terrain;
596 ind21 = jindex_terrain * nx + std::min(iindex_terrain+1,nx-1);
597 ind22 = std::min(jindex_terrain+1,ny-1) * nx + std::min(iindex_terrain+1,nx-1);
600 if (iindex_terrain == nx-1 && jindex_terrain == ny-1)
602 z_arr(i,j,klo) = d_zt[ind11];
604 else if (iindex_terrain != nx-1 && jindex_terrain == ny-1)
609 z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_21*d_zt[ind21])/denom;
611 else if (iindex_terrain == nx-1 && jindex_terrain != ny-1)
616 z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12])/denom;
625 z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12] + w_21*d_zt[ind21] + w_22*d_zt[ind22]) / denom;
const Real dx
Definition: ERF_InitCustomPert_ABL.H:23
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
amrex::Real Real
Definition: ERF_ShocInterface.H:19