337 amrex::Vector<amrex::Real> m_xterrain,m_yterrain,m_zterrain;
339 int nx = 0;
int ny = 0;
341 if (amrex::ParallelDescriptor::IOProcessor()) {
343 amrex::Print()<<
"Reading terrain file: "<< fname<< std::endl;
344 std::ifstream file(fname);
346 if (!file.is_open()) {
347 amrex::Abort(
"Error: Could not open the file " + fname+
"\n");
351 if (file.peek() == std::ifstream::traits_type::eof()) {
352 amrex::Abort(
"Error: The file " + fname+
" is empty.\n");
360 file >> lon_min >> lat_min;
361 if(std::fabs(lon_min) > 180.0) {
362 amrex::Error(
"The value of longitude for entry in the first line in " + fname
363 +
" should not exceed 180. It is " + std::to_string(lon_min));
365 if(std::fabs(lat_min) > 90.0) {
366 amrex::Error(
"The value of latitude for entry in the first line in " + fname
367 +
" should not exceed 90. It is " + std::to_string(lat_min));
373 while (file >> value1 >> value2 >> value3) {
374 m_xterrain.push_back(value1);
376 m_yterrain.push_back(value2);
378 m_zterrain.push_back(value3);
381 AMREX_ASSERT(m_xterrain.size() ==
static_cast<long int>(nx*ny));
382 AMREX_ASSERT(m_yterrain.size() ==
static_cast<long int>(ny));
383 AMREX_ASSERT(m_zterrain.size() ==
static_cast<long int>(nx*ny));
387 nx = erf_get_single_value<int>(file,cnt); cnt++;
388 ny = erf_get_single_value<int>(file,cnt); cnt++;
389 amrex::Print()<<
"Expecting " << nx <<
" values of x, " <<
390 ny <<
" values of y, and " <<
391 nx*ny <<
" values of z" << std::endl;
392 AMREX_ALWAYS_ASSERT(nx > 0);
393 AMREX_ALWAYS_ASSERT(ny > 0);
394 m_xterrain.resize(nx);
395 m_yterrain.resize(ny);
396 m_zterrain.resize(nx * ny);
397 for (
int n = 0; n < nx; n++) {
398 m_xterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
401 for (
int n = 0; n < ny; n++) {
402 m_yterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
405 for (
int n = 0; n < nx * ny; n++) {
406 m_zterrain[n] = erf_get_single_value<amrex::Real>(file,cnt);
415 amrex::ParallelDescriptor::Bcast(&nx,1,amrex::ParallelDescriptor::IOProcessorNumber());
416 amrex::ParallelDescriptor::Bcast(&ny,1,amrex::ParallelDescriptor::IOProcessorNumber());
420 m_xterrain.resize(nx);
421 m_yterrain.resize(ny);
422 m_zterrain.resize(nz);
424 amrex::ParallelDescriptor::Bcast(m_xterrain.data(),nx,amrex::ParallelDescriptor::IOProcessorNumber());
425 amrex::ParallelDescriptor::Bcast(m_yterrain.data(),ny,amrex::ParallelDescriptor::IOProcessorNumber());
426 amrex::ParallelDescriptor::Bcast(m_zterrain.data(),nz,amrex::ParallelDescriptor::IOProcessorNumber());
429 amrex::Gpu::DeviceVector<amrex::Real> d_xterrain(nx),d_yterrain(ny),d_zterrain(nz);
430 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), d_xterrain.begin());
431 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), d_yterrain.begin());
432 amrex::Gpu::copy(amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), d_zterrain.begin());
438 auto dx = geom.CellSizeArray();
439 auto ProbLoArr = geom.ProbLoArray();
441 int ilo = geom.Domain().smallEnd(0);
442 int jlo = geom.Domain().smallEnd(1);
443 int klo = geom.Domain().smallEnd(2);
444 int ihi = geom.Domain().bigEnd(0) + 1;
445 int jhi = geom.Domain().bigEnd(1) + 1;
447 amrex::Box
zbx = terrain_fab.box();
448 amrex::Array4<amrex::Real>
const& z_arr = terrain_fab.array();
450 amrex::ParallelFor(
zbx, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
453 int ii = amrex::min(amrex::max(i,ilo),ihi);
454 int jj = amrex::min(amrex::max(j,jlo),jhi);
460 int ind11, ind12, ind21, ind22;
463 int iindex_terrain=-1;
464 int jindex_terrain=-1;
470 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
472 jindex_terrain = it-1;
475 if (jindex_terrain == -1) {
476 jindex_terrain = ny-1;
479 int gstart = (jindex_terrain )*nx;
480 int gend = (jindex_terrain+1)*nx-1;
481 for (
int it = gstart; it <= gend && iindex_terrain == -1; it++) {
483 iindex_terrain = it-gstart;
488 ind11 = jindex_terrain*nx + iindex_terrain;
495 y1 = d_yt[jindex_terrain];
496 y2 = d_yt[jindex_terrain+1];
505 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];
509 for (
int it = 0; it < ny && jindex_terrain == -1; it++) {
511 jindex_terrain = it-1;
514 if (jindex_terrain == -1) {
515 jindex_terrain = ny-1;
518 for (
int it = 0; it < nx && iindex_terrain == -1; it++) {
520 iindex_terrain = it-1;
523 if (iindex_terrain == -1) {
524 iindex_terrain = nx-1;
528 x1 = d_xt[iindex_terrain];
529 x2 = d_xt[iindex_terrain+1];
530 y1 = d_yt[jindex_terrain];
531 y2 = d_yt[jindex_terrain+1];
537 ind11 = iindex_terrain * ny + jindex_terrain;
538 ind21 = std::min(iindex_terrain+1,nx-1) * ny + jindex_terrain;
540 ind12 = iindex_terrain * ny + std::min(jindex_terrain+1,ny-1);
541 ind22 = std::min(iindex_terrain+1,nx-1) * ny + std::min(jindex_terrain+1,ny-1);
547 ind11 = jindex_terrain * nx + iindex_terrain;
548 ind12 = std::min(jindex_terrain+1,ny-1) * nx + iindex_terrain;
550 ind21 = jindex_terrain * nx + std::min(iindex_terrain+1,nx-1);
551 ind22 = std::min(jindex_terrain+1,ny-1) * nx + std::min(iindex_terrain+1,nx-1);
554 if (iindex_terrain == nx-1 && jindex_terrain == ny-1)
556 z_arr(i,j,klo) = d_zt[ind11];
558 else if (iindex_terrain != nx-1 && jindex_terrain == ny-1)
563 z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_21*d_zt[ind21])/denom;
565 else if (iindex_terrain == nx-1 && jindex_terrain != ny-1)
570 z_arr(i,j,klo) = (w_11*d_zt[ind11] + w_12*d_zt[ind12])/denom;
579 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 Box zbx
Definition: ERF_SetupDiff.H:9