From e9a2f565396e102aa7abbe467e116e9ebe2c6d8c Mon Sep 17 00:00:00 2001 From: Max Katz Date: Sun, 15 Oct 2023 17:05:29 -0400 Subject: [PATCH] Convert RadInterpBndryData routines to C++ --- Source/radiation/_interpbndry/Make.package | 4 - .../_interpbndry/RADINTERPBNDRYDATA_F.H | 186 ----- .../_interpbndry/RadInterpBndryData.cpp | 233 +++++- .../_interpbndry/RadInterpBndryData_1d.F90 | 106 --- .../_interpbndry/RadInterpBndryData_2d.F90 | 351 -------- .../_interpbndry/RadInterpBndryData_3d.F90 | 753 ------------------ 6 files changed, 197 insertions(+), 1436 deletions(-) delete mode 100644 Source/radiation/_interpbndry/RADINTERPBNDRYDATA_F.H delete mode 100644 Source/radiation/_interpbndry/RadInterpBndryData_1d.F90 delete mode 100644 Source/radiation/_interpbndry/RadInterpBndryData_2d.F90 delete mode 100644 Source/radiation/_interpbndry/RadInterpBndryData_3d.F90 diff --git a/Source/radiation/_interpbndry/Make.package b/Source/radiation/_interpbndry/Make.package index 28b5c8ac99..5120ca1ec5 100644 --- a/Source/radiation/_interpbndry/Make.package +++ b/Source/radiation/_interpbndry/Make.package @@ -4,7 +4,3 @@ CA_BASE=EXE C$(CA_BASE)_sources += RadInterpBndryData.cpp RadBndryData.cpp C$(CA_BASE)_headers += RadInterpBndryData.H RadBndryData.H RadBoundCond.H - -ca_F90$(CA_BASE)_sources += RadInterpBndryData_$(DIM)d.F90 - -ca_F$(CA_BASE)_headers += RADINTERPBNDRYDATA_F.H diff --git a/Source/radiation/_interpbndry/RADINTERPBNDRYDATA_F.H b/Source/radiation/_interpbndry/RADINTERPBNDRYDATA_F.H deleted file mode 100644 index bed644330a..0000000000 --- a/Source/radiation/_interpbndry/RADINTERPBNDRYDATA_F.H +++ /dev/null @@ -1,186 +0,0 @@ -#ifndef CASTRO_RADINTERPBNDRY_F_H -#define CASTRO_RADINTERPBNDRY_F_H - -/*************************************************************/ - -#if defined(RAD_INTERP) - -#if defined(BL_LANG_FORT) - -# define FORT_BDINTERPXLO bdintrpxlo2 -# define FORT_BDINTERPXHI bdintrpxhi2 -# if (AMREX_SPACEDIM > 1) -# define FORT_BDINTERPYLO bdintrpylo2 -# define FORT_BDINTERPYHI bdintrpyhi2 -# endif -# if (AMREX_SPACEDIM > 2) -# define FORT_BDINTERPZLO bdintrpzlo2 -# define FORT_BDINTERPZHI bdintrpzhi2 -# endif - -#else - -#if defined(BL_FORT_USE_LOWERCASE) -# define FORT_BDINTERPXLO bdintrpxlo2 -# define FORT_BDINTERPXHI bdintrpxhi2 -# if (AMREX_SPACEDIM > 1) -# define FORT_BDINTERPYLO bdintrpylo2 -# define FORT_BDINTERPYHI bdintrpyhi2 -# endif -# if (AMREX_SPACEDIM > 2) -# define FORT_BDINTERPZLO bdintrpzlo2 -# define FORT_BDINTERPZHI bdintrpzhi2 -# endif -#elif defined(BL_FORT_USE_UPPERCASE) -# define FORT_BDINTERPXLO BDINTRPXLO2 -# define FORT_BDINTERPXHI BDINTRPXHI2 -# if (AMREX_SPACEDIM > 1) -# define FORT_BDINTERPYLO BDINTRPYLO2 -# define FORT_BDINTERPYHI BDINTRPYHI2 -# endif -# if (AMREX_SPACEDIM > 2) -# define FORT_BDINTERPZLO BDINTRPZLO2 -# define FORT_BDINTERPZHI BDINTRPZHI2 -# endif -#elif defined(BL_FORT_USE_UNDERSCORE) -# define FORT_BDINTERPXLO bdintrpxlo2_ -# define FORT_BDINTERPXHI bdintrpxhi2_ -# if (AMREX_SPACEDIM > 1) -# define FORT_BDINTERPYLO bdintrpylo2_ -# define FORT_BDINTERPYHI bdintrpyhi2_ -# endif -# if (AMREX_SPACEDIM > 2) -# define FORT_BDINTERPZLO bdintrpzlo2_ -# define FORT_BDINTERPZHI bdintrpzhi2_ -# endif -#endif - -#include - -typedef void BDInterpFunc(amrex::Real* bndr, ARLIM_P(blo), ARLIM_P(bhi), - const int* lo, const int* hi, - ARLIM_P(cblo), ARLIM_P(cbhi), - const int* nvar, const int* ratio, - const int* not_covered, - const int* mask, ARLIM_P(mlo), ARLIM_P(mhi), - const amrex::Real* crse, ARLIM_P(clo), ARLIM_P(chi), - amrex::Real* derives); -#ifdef __cplusplus -extern "C" -{ -#endif - - BDInterpFunc FORT_BDINTERPXLO; - BDInterpFunc FORT_BDINTERPXHI; - -#if (AMREX_SPACEDIM > 1) - BDInterpFunc FORT_BDINTERPYLO; - BDInterpFunc FORT_BDINTERPYHI; -#endif - -#if (AMREX_SPACEDIM > 2) - BDInterpFunc FORT_BDINTERPZLO; - BDInterpFunc FORT_BDINTERPZHI; -#endif -#ifdef __cplusplus -} -#endif - -#endif - -/*************************************************************/ - -#else /*RAD_INTERP*/ - -/*************************************************************/ - -#if defined(BL_LANG_FORT) - -# define FORT_BDINTERPXLO bdintrpxlo -# define FORT_BDINTERPXHI bdintrpxhi -# if (AMREX_SPACEDIM > 1) -# define FORT_BDINTERPYLO bdintrpylo -# define FORT_BDINTERPYHI bdintrpyhi -# endif -# if (AMREX_SPACEDIM > 2) -# define FORT_BDINTERPZLO bdintrpzlo -# define FORT_BDINTERPZHI bdintrpzhi -# endif - -#else - -#if defined(BL_FORT_USE_LOWERCASE) -# define FORT_BDINTERPXLO bdintrpxlo -# define FORT_BDINTERPXHI bdintrpxhi -# if (AMREX_SPACEDIM > 1) -# define FORT_BDINTERPYLO bdintrpylo -# define FORT_BDINTERPYHI bdintrpyhi -# endif -# if (AMREX_SPACEDIM > 2) -# define FORT_BDINTERPZLO bdintrpzlo -# define FORT_BDINTERPZHI bdintrpzhi -# endif -#elif defined(BL_FORT_USE_UPPERCASE) -# define FORT_BDINTERPXLO BDINTRPXLO -# define FORT_BDINTERPXHI BDINTRPXHI -# if (AMREX_SPACEDIM > 1) -# define FORT_BDINTERPYLO BDINTRPYLO -# define FORT_BDINTERPYHI BDINTRPYHI -# endif -# if (AMREX_SPACEDIM > 2) -# define FORT_BDINTERPZLO BDINTRPZLO -# define FORT_BDINTERPZHI BDINTRPZHI -# endif -#elif defined(BL_FORT_USE_UNDERSCORE) -# define FORT_BDINTERPXLO bdintrpxlo_ -# define FORT_BDINTERPXHI bdintrpxhi_ -# if (AMREX_SPACEDIM > 1) -# define FORT_BDINTERPYLO bdintrpylo_ -# define FORT_BDINTERPYHI bdintrpyhi_ -# endif -# if (AMREX_SPACEDIM > 2) -# define FORT_BDINTERPZLO bdintrpzlo_ -# define FORT_BDINTERPZHI bdintrpzhi_ -# endif -#endif - -#include - -typedef void BDInterpFunc(amrex::Real* bndr, ARLIM_P(blo), ARLIM_P(bhi), - const int* lo, const int* hi, - ARLIM_P(cblo), ARLIM_P(cbhi), - const int* nvar, const int* ratio, - const int* not_covered, - const int* mask, ARLIM_P(mlo), ARLIM_P(mhi), - const amrex::Real* crse, ARLIM_P(clo), ARLIM_P(chi), - amrex::Real* derives); -#ifdef __cplusplus -extern "C" -{ -#endif - - BDInterpFunc FORT_BDINTERPXLO; - BDInterpFunc FORT_BDINTERPXHI; - -#if (AMREX_SPACEDIM > 1) - BDInterpFunc FORT_BDINTERPYLO; - BDInterpFunc FORT_BDINTERPYHI; -#endif - -#if (AMREX_SPACEDIM > 2) - BDInterpFunc FORT_BDINTERPZLO; - BDInterpFunc FORT_BDINTERPZHI; -#endif -#ifdef __cplusplus -} -#endif - -#endif - -/*************************************************************/ - -#endif /*RAD_INTERP*/ - -/*************************************************************/ - -#endif diff --git a/Source/radiation/_interpbndry/RadInterpBndryData.cpp b/Source/radiation/_interpbndry/RadInterpBndryData.cpp index 2e24aa7061..dd921c01a5 100644 --- a/Source/radiation/_interpbndry/RadInterpBndryData.cpp +++ b/Source/radiation/_interpbndry/RadInterpBndryData.cpp @@ -2,35 +2,9 @@ #include #include #include -#include using namespace amrex; -static BDInterpFunc* bdfunc[2*AMREX_SPACEDIM]; -static int bdfunc_set = 0; - -static void bdfunc_init() -{ - Orientation xloface(0,Orientation::low); - Orientation xhiface(0,Orientation::high); - - bdfunc[xloface] = FORT_BDINTERPXLO; - bdfunc[xhiface] = FORT_BDINTERPXHI; -#if (AMREX_SPACEDIM > 1) - Orientation yloface(1,Orientation::low); - Orientation yhiface(1,Orientation::high); - bdfunc[yloface] = FORT_BDINTERPYLO; - bdfunc[yhiface] = FORT_BDINTERPYHI; -#endif -#if (AMREX_SPACEDIM > 2) - Orientation zloface(2,Orientation::low); - Orientation zhiface(2,Orientation::high); - bdfunc[zloface] = FORT_BDINTERPZLO; - bdfunc[zhiface] = FORT_BDINTERPZHI; -#endif - -} - #if (AMREX_SPACEDIM == 2) #define NUMDERIV 2 #endif @@ -111,8 +85,6 @@ RadInterpBndryData::setBndryValues(BndryRegister& crse, int c_start, int bnd_start, int num_comp, IntVect& ratio, const BCRec& bc) { - if (! bdfunc_set) bdfunc_init(); - // check that boxarrays are identical BL_ASSERT( grids.size() ); BL_ASSERT( grids == fine.boxArray() ); @@ -159,24 +131,213 @@ RadInterpBndryData::setBndryValues(BndryRegister& crse, int c_start, const Mask& mask = *masks[face][grd]; const int* mlo = mask.loVect(); const int* mhi = mask.hiVect(); - const int* mdat = mask.dataPtr(); + Array4 const mask_arr = mask.array(); const FArrayBox& crse_fab = crse[face][mfi]; const int* clo = crse_fab.loVect(); const int* chi = crse_fab.hiVect(); - const Real* cdat = crse_fab.dataPtr(c_start); + Array4 const crse = crse_fab.array(c_start); + + int iclo = 0; + int ichi = 0; + int jclo = 0; + int jchi = 0; + int kclo = 0; + int kchi = 0; + + if (face.coordDir() != 0) { + iclo = cblo[0]; + ichi = cbhi[1]; + } + +#if AMREX_SPACEDIM >= 2 + if (face.coordDir() != 1) { + jclo = cblo[1]; + jchi = cbhi[1]; + } +#endif + +#if AMREX_SPACEDIM == 3 + if (face.coordDir() != 2) { + kclo = cblo[2]; + kchi = cbhi[2]; + } +#endif + + int ratiox = 1; + int ratioy = 1; + int ratioz = 1; + + if (face.coordDir() != 0) { + ratiox = ratio[0]; + } + +#if AMREX_SPACEDIM >= 2 + if (face.coordDir() != 1) { + ratioy = ratio[1]; + } +#endif + +#if AMREX_SPACEDIM == 3 + if (face.coordDir() != 2) { + ratioz = ratio[2]; + } +#endif FArrayBox& bnd_fab = bndry[face][mfi]; const int* blo = bnd_fab.loVect(); const int* bhi = bnd_fab.hiVect(); - Real* bdat = bnd_fab.dataPtr(bnd_start); + Array4 const bdry = bnd_fab.array(bnd_start); int is_not_covered = RadBndryData::not_covered; - bdfunc[face](bdat,ARLIM(blo),ARLIM(bhi), - lo,hi,ARLIM(cblo),ARLIM(cbhi), - &num_comp,ratio.getVect(),&is_not_covered, - mdat,ARLIM(mlo),ARLIM(mhi), - cdat,ARLIM(clo),ARLIM(chi),derives); + + for (int n = 0; n < num_comp; ++n) { + + // Note that only two of these three loops will do something + // nontrivial, depending on which face we are working on. + + for (int koff = 0; koff < ratioz; ++koff) { + Real zz = (koff - 0.5_rt * ratioz + 0.5_rt) / ratioz; + + for (int kc = kclo; kc <= kchi; ++kc) { + int k = ratioz * kc + koff; + + for (int joff = 0; joff < ratioy; ++joff) { + Real yy = (joff - 0.5_rt * ratioy + 0.5_rt) / ratioy; + + for (int jc = jclo; jc <= jchi; ++jc) { + int j = ratioy * jc + joff; + + for (int ioff = 0; ioff < ratiox; ++ioff) { + Real xx = (ioff - 0.5_rt * ratiox + 0.5_rt) / ratiox; + + for (int ic = iclo; ic <= ichi; ++ic) { + int i = ratiox * ic + ioff; + + Real xderiv = 0.0_rt; + Real yderiv = 0.0_rt; + Real zderiv = 0.0_rt; + + Real xxderiv = 0.0_rt; + Real yyderiv = 0.0_rt; + Real zzderiv = 0.0_rt; + + Real xyderiv = 0.0_rt; + Real xzderiv = 0.0_rt; + Real yzderiv = 0.0_rt; + + if (face.coordDir() != 0) { + xderiv = 0.5_rt * (crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n)); + xxderiv = 0.5_rt * (crse(ic+1,jc,kc,n) - 2.0_rt * crse(ic,jc,kc,n) + crse(ic-1,jc,kc,n)); + } + +#if AMREX_SPACEDIM >= 2 + if (face.coordDir() != 1) { + yderiv = 0.5_rt * (crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n)); + yyderiv = 0.5_rt * (crse(ic,jc+1,kc,n) - 2.0_rt * crse(ic,jc,kc,n) + crse(ic,jc-1,kc,n)); + } +#endif + +#if AMREX_SPACEDIM == 3 + if (face.coordDir() != 2) { + zderiv = 0.5_rt * (crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n)); + zzderiv = 0.5_rt * (crse(ic,jc,kc+1,n) - 2.0_rt * crse(ic,jc,kc,n) + crse(ic,jc,kc-1,n)); + } +#endif + +#if AMREX_SPACEDIM == 3 + if (face.coordDir() == 0) { + yzderiv = 0.25_rt * (crse(ic,jc+1,kc+1,n) - crse(ic,jc-1,kc+1,n) + + crse(ic,jc-1,kc-1,n) - crse(ic,jc+1,kc-1,n)); + } + + if (face.coordDir() == 1) { + xzderiv = 0.25_rt * (crse(ic+1,jc,kc+1,n) - crse(ic-1,jc,kc+1,n) + + crse(ic-1,jc,kc-1,n) - crse(ic+1,jc,kc-1,n)); + } + + if (face.coordDir() == 2) { + xyderiv = 0.25_rt * (crse(ic+1,jc+1,kc,n) - crse(ic-1,jc+1,kc,n) + + crse(ic-1,jc-1,kc,n) - crse(ic+1,jc-1,kc,n)); + } +#endif + + if (mask_arr(i-1,j,k) != not_covered) { + xderiv = crse(ic+1,jc,kc,n) - crse(ic,jc,kc,n); + xxderiv = 0.0_rt; + } + if (mask_arr(i+ratiox,j,k) != not_covered) { + xderiv = crse(ic,jc,kc,n) - crse(ic-1,jc,kc,n); + xxderiv = 0.0_rt; + } + if (mask_arr(i-1,j,k) != not_covered && mask_arr(i+ratiox,j,k) != not_covered) { + xderiv = 0.0_rt; + } + +#if AMREX_SPACEDIM >= 2 + if (mask_arr(i,j-1,k) != not_covered) { + yderiv = crse(ic,jc+1,kc,n) - crse(ic,jc,kc,n); + yyderiv = 0.0_rt; + } + if (mask_arr(i,j+ratioy,k) != not_covered) { + yderiv = crse(ic,jc,kc,n) - crse(ic,jc-1,kc,n); + yyderiv = 0.0_rt; + } + if (mask_arr(i,j-1,k) != not_covered && mask_arr(i,j+ratioy,k) != not_covered) { + yderiv = 0.0_rt; + } +#endif + +#if AMREX_SPACEDIM == 3 + if (mask_arr(i,j,k-1) != not_covered) { + yderiv = crse(ic,jc,kc+1,n) - crse(ic,jc,kc,n); + yyderiv = 0.0_rt; + } + if (mask_arr(i,j,k+ratioz) != not_covered) { + yderiv = crse(ic,jc,kc,n) - crse(ic,jc,kc-1,n); + yyderiv = 0.0_rt; + } + if (mask_arr(i,j,k-1) != not_covered && mask_arr(i,j,k+ratioz) != not_covered) { + yderiv = 0.0_rt; + } +#endif + +#if AMREX_SPACEDIM == 3 + if ((mask_arr(i,j+ratioy,k+ratioz) != not_covered) || + (mask_arr(i,j-1 ,k+ratioz) != not_covered) || + (mask_arr(i,j+ratioy,k-1 ) != not_covered) || + (mask_arr(i,j-1 ,k-1 ) != not_covered)) { + yzderiv = 0.0_rt; + } + + if ((mask_arr(i+ratiox,j,k+ratioz) != not_covered) || + (mask_arr(i-1 ,j,k+ratioz) != not_covered) || + (mask_arr(i+ratiox,j,k-1 ) != not_covered) || + (mask_arr(i-1 ,j,k-1 ) != not_covered)) { + xzderiv = 0.0_rt; + } + + if ((mask_arr(i+ratiox,j+ratioy,k) != not_covered) || + (mask_arr(i-1 ,j+ratioy,k) != not_covered) || + (mask_arr(i+ratiox,j-1 ,k) != not_covered) || + (mask_arr(i-1 ,j-1 ,k) != not_covered)) { + xyderiv = 0.0_rt; + } +#endif + + bdry(i,j,k,n) = crse(ic,jc,kc,n) + + xx * xderiv + xx * xx * xxderiv + + yy * yderiv + yy * yy * yyderiv + + zz * zderiv + zz * zz * zzderiv + + xx * yy * xyderiv + xx * zz * xzderiv + yy * zz * yzderiv; + } + } + } + } + } + } + + } // num_comp } else { // physical bndry, copy from ghost region of // corresponding grid diff --git a/Source/radiation/_interpbndry/RadInterpBndryData_1d.F90 b/Source/radiation/_interpbndry/RadInterpBndryData_1d.F90 deleted file mode 100644 index 293d146038..0000000000 --- a/Source/radiation/_interpbndry/RadInterpBndryData_1d.F90 +++ /dev/null @@ -1,106 +0,0 @@ -#include -#include -#include -#include -#include - -#define SDIM 1 - -! --------------------------------------------------------------- -! :: FORT_BDINTERPXLO : Interpolation on Xlo Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(2) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPXLO(bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(2), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(mask) - integer DIMDEC(crse) - integer DIMDEC(cb) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(1) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - integer i, ic, n - - ic = ARG_L1(cb)-1 - i = lo(1)-1 - - do n = 1, nvar - ! interpolate to fine grid - bdry(i,n) = crse(ic,n) - enddo - - return -end subroutine FORT_BDINTERPXLO - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPXHI : Interpolation on Xhi Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(2) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPXHI (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(2), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(mask) - integer DIMDEC(cb) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(1) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - integer i, ic, n - - ic = ARG_H1(cb)+1 - i = hi(1)+1 - - do n = 1, nvar - ! interpolate to fine grid - bdry(i,n) = crse(ic,n) - enddo - - return -end subroutine FORT_BDINTERPXHI diff --git a/Source/radiation/_interpbndry/RadInterpBndryData_2d.F90 b/Source/radiation/_interpbndry/RadInterpBndryData_2d.F90 deleted file mode 100644 index e2198b178a..0000000000 --- a/Source/radiation/_interpbndry/RadInterpBndryData_2d.F90 +++ /dev/null @@ -1,351 +0,0 @@ -#include -#include -#include -#include -#include - -#define SDIM 2 -#define NUMDERIV 2 -#define XDER 1 -#define X2DER 2 - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPXLO : Interpolation on Xlo Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(2) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPXLO (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(2), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(mask) - integer DIMDEC(crse) - integer DIMDEC(cb) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM2(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - REAL_T xx, ainterp - integer i, j, ic, jc, off, n - integer jclo, jchi, ratioy - - ratioy = ratios(2) - - jclo = ARG_L2(cb) - jchi = ARG_H2(cb) - ic = ARG_L1(cb)-1 - i = lo(1)-1 - - do n = 1, nvar - ! define interp coefs - do jc = jclo, jchi - j = ratioy*jc - derives(jc,XDER) = half * (crse(ic,jc+1,n) - crse(ic,jc-1,n)) - derives(jc,X2DER) = half * (crse(ic,jc+1,n) & - - crse(ic,jc ,n) * two & - + crse(ic,jc-1,n)) - - if (mask(i,j-1) .ne. not_covered) then - derives(jc,XDER) = crse(ic,jc+1,n) - crse(ic,jc ,n) - derives(jc,X2DER) = zero - endif - if (mask(i,j+ratioy) .ne. not_covered) then - derives(jc,XDER) = crse(ic,jc ,n) - crse(ic,jc-1,n) - derives(jc,X2DER) = zero - endif - if (mask(i,j-1) .ne. not_covered .and. & - mask(i,j+ratioy) .ne. not_covered) then - derives(jc,XDER) = zero - endif - enddo - - ! interpolate to fine grid - do off = 0, ratioy - 1 - xx = (off - half * ratioy + half)/ratioy - do jc = jclo, jchi - j = ratioy*jc + off - bdry(i,j,n) = crse(ic,jc,n) & - + derives(jc,XDER) *xx & - + derives(jc,X2DER)*xx**2 - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPXLO - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPXHI : Interpolation on Xhi Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(2) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPXHI (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(2), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(mask) - integer DIMDEC(cb) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM2(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - REAL_T xx, ainterp - integer i, j, ic, jc, off, n - integer jclo, jchi, ratioy - - ratioy = ratios(2) - - jclo = ARG_L2(cb) - jchi = ARG_H2(cb) - ic = ARG_H1(cb)+1 - i = hi(1)+1 - - do n = 1, nvar - ! define interp coefs - do jc = jclo, jchi - j = ratioy*jc - derives(jc,XDER) = half * (crse(ic,jc+1,n) - crse(ic,jc-1,n)) - derives(jc,X2DER) = half * (crse(ic,jc+1,n) & - - crse(ic,jc ,n) * two & - + crse(ic,jc-1,n)) - - if (mask(i,j-1) .ne. not_covered) then - derives(jc,XDER) = crse(ic,jc+1,n) - crse(ic,jc ,n) - derives(jc,X2DER) = zero - endif - if (mask(i,j+ratioy) .ne. not_covered) then - derives(jc,XDER) = crse(ic,jc ,n) - crse(ic,jc-1,n) - derives(jc,X2DER) = zero - endif - if (mask(i,j-1) .ne. not_covered .and. & - mask(i,j+ratioy) .ne. not_covered) then - derives(jc,XDER) = zero - endif - enddo - - ! interpolate to fine grid - do off = 0, ratioy - 1 - xx = (off - half * ratioy + half)/ratioy - do jc = jclo, jchi - j = ratioy*jc + off - bdry(i,j,n) = crse(ic,jc,n) & - + derives(jc,XDER) *xx & - + derives(jc,X2DER)*xx**2 - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPXHI - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPYLO : Interpolation on Ylo Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(2) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPYLO (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(2), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(mask) - integer DIMDEC(cb) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM1(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - REAL_T xx, ainterp - integer i, j, ic, jc, off, n - integer iclo, ichi, ratiox - - ratiox = ratios(1) - - iclo = ARG_L1(cb) - ichi = ARG_H1(cb) - jc = ARG_L2(cb)-1 - j = lo(2)-1 - - do n = 1, nvar - ! define interp coefs - do ic = iclo, ichi - i = ratiox*ic - derives(ic,XDER) = half * (crse(ic+1,jc,n) - crse(ic-1,jc,n)) - derives(ic,X2DER) = half * (crse(ic+1,jc,n) & - - crse(ic ,jc,n) * two & - + crse(ic-1,jc,n)) - if (mask(i-1,j) .ne. not_covered) then - derives(ic,XDER) = crse(ic+1,jc,n) - crse(ic ,jc,n) - derives(ic,X2DER) = zero - endif - if (mask(i+ratiox,j) .ne. not_covered) then - derives(ic,XDER) = crse(ic ,jc,n) - crse(ic-1,jc,n) - derives(ic,X2DER) = zero - endif - if (mask(i-1,j) .ne. not_covered .and. & - mask(i+ratiox,j) .ne. not_covered) then - derives(ic,XDER) = zero - endif - enddo - - ! interpolate to fine grid - do off = 0, ratiox - 1 - xx = (off - half * ratiox + half)/ratiox - do ic = iclo, ichi - i = ratiox*ic + off - bdry(i,j,n) = crse(ic,jc,n) & - + derives(ic,XDER) *xx & - + derives(ic,X2DER)*xx**2 - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPYLO - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPYHI : Interpolation on Yhi Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(2) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPYHI (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(2), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(mask) - integer DIMDEC(cb) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM1(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - REAL_T xx, ainterp - integer i, j, ic, jc, off, n - integer iclo, ichi, ratiox - - ratiox = ratios(1) - - iclo = ARG_L1(cb) - ichi = ARG_H1(cb) - jc = ARG_H2(cb)+1 - j = hi(2)+1 - - do n = 1, nvar - ! define interp coefs - do ic = iclo, ichi - i = ratiox*ic - derives(ic,XDER) = half * (crse(ic+1,jc,n) - crse(ic-1,jc,n)) - derives(ic,X2DER) = half * (crse(ic+1,jc,n) & - - crse(ic ,jc,n) * two & - + crse(ic-1,jc,n)) - if (mask(i-1,j) .ne. not_covered) then - derives(ic,XDER) = crse(ic+1,jc,n) - crse(ic ,jc,n) - derives(ic,X2DER) = zero - endif - if (mask(i+ratiox,j) .ne. not_covered) then - derives(ic,XDER) = crse(ic ,jc,n) - crse(ic-1,jc,n) - derives(ic,X2DER) = zero - endif - if (mask(i-1,j) .ne. not_covered .and. & - mask(i+ratiox,j) .ne. not_covered) then - derives(ic,XDER) = zero - endif - enddo - - ! interpolate to fine grid - do off = 0, ratiox - 1 - xx = (off - half * ratiox + half)/ratiox - do ic = iclo, ichi - i = ratiox*ic + off - bdry(i,j,n) = crse(ic,jc,n) & - + derives(ic,XDER) *xx & - + derives(ic,X2DER)*xx**2 - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPYHI - diff --git a/Source/radiation/_interpbndry/RadInterpBndryData_3d.F90 b/Source/radiation/_interpbndry/RadInterpBndryData_3d.F90 deleted file mode 100644 index 95b2250a76..0000000000 --- a/Source/radiation/_interpbndry/RadInterpBndryData_3d.F90 +++ /dev/null @@ -1,753 +0,0 @@ -#include -#include -#include -#include -#include - -#define SDIM 3 -#define NUMDERIV 5 -#define XDER 1 -#define YDER 2 -#define X2DER 3 -#define Y2DER 4 -#define XYDER 5 - -! --------------------------------------------------------------- -! :: FORT_BDINTERPXLO : Interpolation on Xlo Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(3) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array for derivatives -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPXLO (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(3), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(cb) - integer DIMDEC(mask) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM23(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - REAL_T xx, yy, ainterp - integer i, j, k, ic, jc, kc, joff, koff, n - integer jclo, jchi, kclo, kchi, ratioy, ratioz - - ratioy = ratios(2) - ratioz = ratios(3) - - kclo = ARG_L3(cb) - kchi = ARG_H3(cb) - jclo = ARG_L2(cb) - jchi = ARG_H2(cb) - ic = ARG_L1(cb)-1 - i = lo(1)-1 - - do n = 1, nvar - ! define interp coefs - do kc = kclo, kchi - k = ratioz*kc - do jc = jclo, jchi - j = ratioy*jc - derives(jc,kc,XDER) = half*(crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n)) - derives(jc,kc,X2DER) = half*(crse(ic,jc+1,kc,n) - two*crse(ic,jc,kc,n) & - + crse(ic,jc-1,kc,n)) - derives(jc,kc,YDER) = half*(crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n)) - derives(jc,kc,Y2DER) = half*(crse(ic,jc,kc+1,n) - two*crse(ic,jc,kc,n) & - + crse(ic,jc,kc-1,n)) - derives(jc,kc,XYDER) = fourth*(crse(ic,jc+1,kc+1,n) - crse(ic,jc-1,kc+1,n) & - + crse(ic,jc-1,kc-1,n) - crse(ic,jc+1,kc-1,n)) - - if (mask(i,j-1,k) .ne. not_covered) then - derives(jc,kc,XDER) = crse(ic,jc+1,kc,n) - crse(ic,jc,kc,n) - derives(jc,kc,X2DER) = zero - endif - if (mask(i,j+ratioy,k) .ne. not_covered) then - derives(jc,kc,XDER) = crse(ic,jc,kc,n) - crse(ic,jc-1,kc,n) - derives(jc,kc,X2DER) = zero - endif - if (mask(i,j-1,k) .ne. not_covered .and. & - mask(i,j+ratioy,k) .ne. not_covered) then - derives(jc,kc,XDER) = zero - endif - - if (mask(i,j,k-1) .ne. not_covered) then - derives(jc,kc,YDER) = crse(ic,jc,kc+1,n) - crse(ic,jc,kc,n) - derives(jc,kc,Y2DER) = zero - endif - if (mask(i,j,k+ratioz) .ne. not_covered) then - derives(jc,kc,YDER) = crse(ic,jc,kc,n) - crse(ic,jc,kc-1,n) - derives(jc,kc,Y2DER) = zero - endif - if (mask(i,j,k-1) .ne. not_covered .and. & - mask(i,j,k+ratioz) .ne. not_covered) then - derives(jc,kc,YDER) = zero - endif - if (( mask(i,j+ratioy,k+ratioz) .ne. not_covered ) .or. & - ( mask(i,j-1,k+ratioz) .ne. not_covered ) .or. & - ( mask(i,j+ratioy,k-1) .ne. not_covered ) .or. & - ( mask(i,j-1,k-1) .ne. not_covered ) ) then - - derives(jc,kc,XYDER) = zero - endif - enddo - enddo - - ! interpolate to fine grid - do koff = 0, ratioz - 1 - yy = (koff - half * ratioz + half)/ratioz - do kc = kclo,kchi - k = ratioz*kc + koff - do joff = 0, ratioy - 1 - xx = (joff - half * ratioy + half)/ratioy - do jc = jclo, jchi - j = ratioy*jc + joff - bdry(i,j,k,n) = crse(ic,jc,kc,n) + xx*derives(jc,kc,XDER) & - + derives(jc,kc,X2DER)*xx**2 + yy*derives(jc,kc,YDER) & - + derives(jc,kc,Y2DER)*yy**2 + xx*yy*derives(jc,kc,XYDER) - enddo - enddo - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPXLO - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPXHI : Interpolation on Xhi Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(3) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array for derivatives -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPXHI (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(3), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(cb) - integer DIMDEC(mask) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM23(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - REAL_T xx, yy, ainterp - integer i, j, k, ic, jc, kc, joff, koff, n - integer jclo, jchi, kclo, kchi, ratioy, ratioz - - ratioy = ratios(2) - ratioz = ratios(3) - - kclo = ARG_L3(cb) - kchi = ARG_H3(cb) - jclo = ARG_L2(cb) - jchi = ARG_H2(cb) - ic = ARG_H1(cb)+1 - i = hi(1)+1 - - do n = 1, nvar - ! define interp coefs - do kc = kclo, kchi - k = ratioz*kc - do jc = jclo, jchi - j = ratioy*jc - derives(jc,kc,XDER) = half*(crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n)) - derives(jc,kc,X2DER) = half*(crse(ic,jc+1,kc,n) - two*crse(ic,jc,kc,n) & - + crse(ic,jc-1,kc,n)) - derives(jc,kc,YDER) = half*(crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n)) - derives(jc,kc,Y2DER) = half*(crse(ic,jc,kc+1,n) - two*crse(ic,jc,kc,n) & - + crse(ic,jc,kc-1,n)) - derives(jc,kc,XYDER) = fourth*(crse(ic,jc+1,kc+1,n) - crse(ic,jc-1,kc+1,n) & - + crse(ic,jc-1,kc-1,n) - crse(ic,jc+1,kc-1,n)) - - if (mask(i,j-1,k) .ne. not_covered) then - derives(jc,kc,XDER) = crse(ic,jc+1,kc,n) - crse(ic,jc,kc,n) - derives(jc,kc,X2DER) = zero - endif - if (mask(i,j+ratioy,k) .ne. not_covered) then - derives(jc,kc,XDER) = crse(ic,jc,kc,n) - crse(ic,jc-1,kc,n) - derives(jc,kc,X2DER) = zero - endif - if (mask(i,j-1,k) .ne. not_covered .and. & - mask(i,j+ratioy,k) .ne. not_covered) then - derives(jc,kc,XDER) = zero - endif - - if (mask(i,j,k-1) .ne. not_covered) then - derives(jc,kc,YDER) = crse(ic,jc,kc+1,n) - crse(ic,jc,kc,n) - derives(jc,kc,Y2DER) = zero - endif - if (mask(i,j,k+ratioz) .ne. not_covered) then - derives(jc,kc,YDER) = crse(ic,jc,kc,n) - crse(ic,jc,kc-1,n) - derives(jc,kc,Y2DER) = zero - endif - if (mask(i,j,k-1) .ne. not_covered .and. & - mask(i,j,k+ratioz) .ne. not_covered) then - derives(jc,kc,YDER) = zero - endif - if (( mask(i,j+ratioy,k+ratioz) .ne. not_covered ) .or. & - ( mask(i,j-1,k+ratioz) .ne. not_covered ) .or. & - ( mask(i,j+ratioy,k-1) .ne. not_covered ) .or. & - ( mask(i,j-1,k-1) .ne. not_covered ) ) then - - derives(jc,kc,XYDER) = zero - endif - enddo - enddo - - ! interpolate to fine grid - do koff = 0, ratioz - 1 - yy = (koff - half * ratioz + half)/ratioz - do kc = kclo,kchi - k = ratioz*kc + koff - do joff = 0, ratioy - 1 - xx = (joff - half * ratioy + half)/ratioy - do jc = jclo, jchi - j = ratioy*jc + joff - bdry(i,j,k,n) = crse(ic,jc,kc,n) + xx*derives(jc,kc,XDER) & - + derives(jc,kc,X2DER)*xx**2 + yy*derives(jc,kc,YDER) & - + derives(jc,kc,Y2DER)*yy**2 + xx*yy*derives(jc,kc,XYDER) - enddo - enddo - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPXHI - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPYLO : Interpolation on Ylo Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(3) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array for derivatives -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPYLO (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(3), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(cb) - integer DIMDEC(mask) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM13(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - REAL_T xx, yy, ainterp - integer i, j, k, ic, jc, kc, ioff, koff, n - integer iclo, ichi, kclo, kchi, ratiox, ratioz - - ratiox = ratios(1) - ratioz = ratios(3) - - kclo = ARG_L3(cb) - kchi = ARG_H3(cb) - iclo = ARG_L1(cb) - ichi = ARG_H1(cb) - jc = ARG_L2(cb)-1 - j = lo(2)-1 - - do n = 1, nvar - ! define interp coefs - do kc = kclo, kchi - k = ratioz*kc - do ic = iclo, ichi - i = ratiox*ic - derives(ic,kc,XDER) = half*(crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n)) - derives(ic,kc,X2DER) = half*(crse(ic+1,jc,kc,n) - two*crse(ic,jc,kc,n) & - + crse(ic-1,jc,kc,n)) - derives(ic,kc,YDER) = half*(crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n)) - derives(ic,kc,Y2DER) = half*(crse(ic,jc,kc+1,n) - two*crse(ic,jc,kc,n) & - + crse(ic,jc,kc-1,n)) - derives(ic,kc,XYDER) = fourth*(crse(ic+1,jc,kc+1,n) - crse(ic-1,jc,kc+1,n) & - + crse(ic-1,jc,kc-1,n) - crse(ic+1,jc,kc-1,n)) - - if (mask(i-1,j,k) .ne. not_covered) then - derives(ic,kc,XDER) = crse(ic+1,jc,kc,n) - crse(ic,jc,kc,n) - derives(ic,kc,X2DER) = zero - endif - if (mask(i+ratiox,j,k) .ne. not_covered) then - derives(ic,kc,XDER) = crse(ic,jc,kc,n) - crse(ic-1,jc,kc,n) - derives(ic,kc,X2DER) = zero - endif - if (mask(i-1,j,k) .ne. not_covered .and. & - mask(i+ratiox,j,k) .ne. not_covered) then - derives(ic,kc,XDER) = zero - endif - - if (mask(i,j,k-1) .ne. not_covered) then - derives(ic,kc,YDER) = crse(ic,jc,kc+1,n) - crse(ic,jc,kc,n) - derives(ic,kc,Y2DER) = zero - endif - if (mask(i,j,k+ratioz) .ne. not_covered) then - derives(ic,kc,YDER) = crse(ic,jc,kc,n) - crse(ic,jc,kc-1,n) - derives(ic,kc,Y2DER) = zero - endif - if (mask(i,j,k-1) .ne. not_covered .and. & - mask(i,j,k+ratioz) .ne. not_covered) then - derives(ic,kc,YDER) = zero - endif - if (( mask(i+ratiox,j,k+ratioz) .ne. not_covered ) .or. & - ( mask(i-1,j,k+ratioz) .ne. not_covered ) .or. & - ( mask(i+ratiox,j,k-1) .ne. not_covered ) .or. & - ( mask(i-1,j,k-1) .ne. not_covered ) ) then - - derives(ic,kc,XYDER) = zero - endif - enddo - enddo - - ! interpolate to fine grid - do koff = 0, ratioz - 1 - yy = (koff - half * ratioz + half)/ratioz - do kc = kclo,kchi - k = ratioz*kc + koff - do ioff = 0, ratiox - 1 - xx = (ioff - half * ratiox + half)/ratiox - do ic = iclo, ichi - i = ratiox*ic + ioff - bdry(i,j,k,n) = crse(ic,jc,kc,n) + xx*derives(ic,kc,XDER) & - + derives(ic,kc,X2DER)*xx**2 + yy*derives(ic,kc,YDER) & - + derives(ic,kc,Y2DER)*yy**2 + xx*yy*derives(ic,kc,XYDER) - enddo - enddo - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPYLO - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPYHI : Interpolation on Yhi Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(3) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array for derivatives -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPYHI (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(3), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(cb) - integer DIMDEC(mask) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM13(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - REAL_T xx, yy, ainterp - integer i, j, k, ic, jc, kc, ioff, koff, n - integer iclo, ichi, kclo, kchi, ratiox, ratioz - - ratiox = ratios(1) - ratioz = ratios(3) - - kclo = ARG_L3(cb) - kchi = ARG_H3(cb) - iclo = ARG_L1(cb) - ichi = ARG_H1(cb) - jc = ARG_H2(cb)+1 - j = hi(2)+1 - - do n = 1, nvar - ! define interp coefs - do kc = kclo, kchi - k = ratioz*kc - do ic = iclo, ichi - i = ratiox*ic - derives(ic,kc,XDER) = half*(crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n)) - derives(ic,kc,X2DER) = half*(crse(ic+1,jc,kc,n) - two*crse(ic,jc,kc,n) & - + crse(ic-1,jc,kc,n)) - derives(ic,kc,YDER) = half*(crse(ic,jc,kc+1,n) - crse(ic,jc,kc-1,n)) - derives(ic,kc,Y2DER) = half*(crse(ic,jc,kc+1,n) - two*crse(ic,jc,kc,n) & - + crse(ic,jc,kc-1,n)) - derives(ic,kc,XYDER) = fourth*(crse(ic+1,jc,kc+1,n) - crse(ic-1,jc,kc+1,n) & - + crse(ic-1,jc,kc-1,n) - crse(ic+1,jc,kc-1,n)) - - if (mask(i-1,j,k) .ne. not_covered) then - derives(ic,kc,XDER) = crse(ic+1,jc,kc,n) - crse(ic,jc,kc,n) - derives(ic,kc,X2DER) = zero - endif - if (mask(i+ratiox,j,k) .ne. not_covered) then - derives(ic,kc,XDER) = crse(ic,jc,kc,n) - crse(ic-1,jc,kc,n) - derives(ic,kc,X2DER) = zero - endif - if (mask(i-1,j,k) .ne. not_covered .and. & - mask(i+ratiox,j,k) .ne. not_covered) then - derives(ic,kc,XDER) = zero - endif - - if (mask(i,j,k-1) .ne. not_covered) then - derives(ic,kc,YDER) = crse(ic,jc,kc+1,n) - crse(ic,jc,kc,n) - derives(ic,kc,Y2DER) = zero - endif - if (mask(i,j,k+ratioz) .ne. not_covered) then - derives(ic,kc,YDER) = crse(ic,jc,kc,n) - crse(ic,jc,kc-1,n) - derives(ic,kc,Y2DER) = zero - endif - if (mask(i,j,k-1) .ne. not_covered .and. & - mask(i,j,k+ratioz) .ne. not_covered) then - derives(ic,kc,YDER) = zero - endif - - if (( mask(i+ratiox,j,k+ratioz) .ne. not_covered ) .or. & - ( mask(i-1,j,k+ratioz) .ne. not_covered ) .or. & - ( mask(i+ratiox,j,k-1) .ne. not_covered ) .or. & - ( mask(i-1,j,k-1) .ne. not_covered ) ) then - - derives(ic,kc,XYDER) = zero - endif - enddo - enddo - - ! interpolate to fine grid - do koff = 0, ratioz - 1 - yy = (koff - half * ratioz + half)/ratioz - do kc = kclo,kchi - k = ratioz*kc + koff - do ioff = 0, ratiox - 1 - xx = (ioff - half * ratiox + half)/ratiox - do ic = iclo, ichi - i = ratiox*ic + ioff - bdry(i,j,k,n) = crse(ic,jc,kc,n) + xx*derives(ic,kc,XDER) & - + derives(ic,kc,X2DER)*xx**2 + yy*derives(ic,kc,YDER) & - + derives(ic,kc,Y2DER)*yy**2 + xx*yy*derives(ic,kc,XYDER) - enddo - enddo - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPYHI - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPZLO : Interpolation on Zlo Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(3) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array for derivatives -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPZLO (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(3), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(cb) - integer DIMDEC(mask) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM12(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - REAL_T xx, yy, ainterp - integer i, j, k, ic, jc, kc, ioff, joff, n - integer iclo, ichi, jclo, jchi, ratiox, ratioy - - ratiox = ratios(1) - ratioy = ratios(2) - - jclo = ARG_L2(cb) - jchi = ARG_H2(cb) - iclo = ARG_L1(cb) - ichi = ARG_H1(cb) - kc = ARG_L3(cb)-1 - k = lo(3)-1 - - do n = 1, nvar - ! define interp coefs - do jc = jclo, jchi - j = ratioy*jc - do ic = iclo, ichi - i = ratiox*ic - - derives(ic,jc,XDER) = half*(crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n) ) - derives(ic,jc,X2DER) = half*(crse(ic+1,jc,kc,n) - two*crse(ic,jc,kc,n) + crse(ic-1,jc,kc,n) ) - derives(ic,jc,YDER) = half*(crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n) ) - derives(ic,jc,Y2DER) = half*(crse(ic,jc+1,kc,n) - two*crse(ic,jc,kc,n) + crse(ic,jc-1,kc,n) ) - derives(ic,jc,XYDER) = fourth*(crse(ic+1,jc+1,kc,n) - crse(ic-1,jc+1,kc,n) & - + crse(ic-1,jc-1,kc,n) - crse(ic+1,jc-1,kc,n)) - - if (mask(i-1,j,k) .ne. not_covered) then - derives(ic,jc,XDER) = crse(ic+1,jc,kc,n) - crse(ic,jc,kc,n) - derives(ic,jc,X2DER) = zero - endif - if (mask(i+ratiox,j,k) .ne. not_covered) then - derives(ic,jc,XDER) = crse(ic,jc,kc,n) - crse(ic-1,jc,kc,n) - derives(ic,jc,X2DER) = zero - endif - if (mask(i-1,j,k) .ne. not_covered .and. & - mask(i+ratiox,j,k) .ne. not_covered) then - derives(ic,jc,XDER) = zero - endif - - if (mask(i,j-1,k) .ne. not_covered) then - derives(ic,jc,YDER) = crse(ic,jc+1,kc,n) - crse(ic,jc,kc,n) - derives(ic,jc,Y2DER) = zero - endif - if (mask(i,j+ratioy,k) .ne. not_covered) then - derives(ic,jc,YDER) = crse(ic,jc,kc,n) - crse(ic,jc-1,kc,n) - derives(ic,jc,Y2DER) = zero - endif - if (mask(i,j-1,k) .ne. not_covered .and. & - mask(i,j+ratioy,k) .ne. not_covered) then - derives(ic,jc,YDER) = zero - endif - - if (( mask(i+ratiox,j+ratioy,k) .ne. not_covered ) .or. & - ( mask(i-1,j+ratioy,k) .ne. not_covered ) .or. & - ( mask(i+ratiox,j-1,k) .ne. not_covered ) .or. & - ( mask(i-1,j-1,k) .ne. not_covered ) ) then - - derives(ic,jc,XYDER) = zero - endif - enddo - enddo - - ! interpolate to fine grid - do joff = 0, ratioy - 1 - yy = (joff - half * ratioy + half)/ratioy - do jc = jclo,jchi - j = ratioy*jc + joff - do ioff = 0, ratiox - 1 - xx = (ioff - half * ratiox + half)/ratiox - do ic = iclo, ichi - i = ratiox*ic + ioff - bdry(i,j,k,n) = crse(ic,jc,kc,n) + xx*derives(ic,jc,XDER) & - + derives(ic,jc,X2DER)*xx**2 + yy*derives(ic,jc,YDER) & - + derives(ic,jc,Y2DER)*yy**2 + xx*yy*derives(ic,jc,XYDER) - enddo - enddo - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPZLO - - -! --------------------------------------------------------------- -! :: FORT_BDINTERPZHI : Interpolation on Zhi Face -! :: Quadratic Interpolation from crse data -! :: in directions transverse to face of grid -! :: -! :: Inputs/Outputs: -! :: bdry <= fine grid bndry data strip -! :: DIMS(bdry) => index limits of bdry -! :: lo,hi => index limits of grd interior -! :: DIMS(cb) => index limits of coarsened grid interior -! :: nvar => number of variables to interpolate -! :: ratios(3) => refinement ratios -! :: not_covered => mask is set to this value if cell is not -! :: covered by another fine grid and not outside the domain. -! :: mask => fine grid mask bndry strip -! :: DIMS(mask) => index limits of mask array -! :: crse => crse grid bndry data strip -! :: DIMS(crse) => index limits of crse array -! :: derives => crse grid tmp array for derivatives -! --------------------------------------------------------------- - -subroutine FORT_BDINTERPZHI (bdry,DIMS(bdry), & - lo,hi,DIMS(cb),nvar,ratios,not_covered, & - mask,DIMS(mask),crse,DIMS(crse),derives) - - integer nvar, ratios(3), not_covered - integer lo(SDIM), hi(SDIM) - integer DIMDEC(bdry) - integer DIMDEC(cb) - integer DIMDEC(mask) - integer DIMDEC(crse) - REAL_T bdry(DIMV(bdry),nvar) - REAL_T derives(DIM12(cb),NUMDERIV) - integer mask(DIMV(mask)) - REAL_T crse(DIMV(crse),nvar) - - REAL_T xx, yy, ainterp - integer i, j, k, ic, jc, kc, ioff, joff, n - integer iclo, ichi, jclo, jchi, ratiox, ratioy - - ratiox = ratios(1) - ratioy = ratios(2) - - jclo = ARG_L2(cb) - jchi = ARG_H2(cb) - iclo = ARG_L1(cb) - ichi = ARG_H1(cb) - kc = ARG_H3(cb)+1 - k = hi(3)+1 - - do n = 1, nvar - ! define interp coefs - do jc = jclo, jchi - j = ratioy*jc - do ic = iclo, ichi - i = ratiox*ic - derives(ic,jc,XDER) = half*(crse(ic+1,jc,kc,n) - crse(ic-1,jc,kc,n)) - derives(ic,jc,X2DER) = half*(crse(ic+1,jc,kc,n) - two*crse(ic,jc,kc,n) & - + crse(ic-1,jc,kc,n)) - derives(ic,jc,YDER) = half*(crse(ic,jc+1,kc,n) - crse(ic,jc-1,kc,n)) - derives(ic,jc,Y2DER) = half*(crse(ic,jc+1,kc,n) - two*crse(ic,jc,kc,n) & - + crse(ic,jc-1,kc,n)) - derives(ic,jc,XYDER) = fourth*(crse(ic+1,jc+1,kc,n) - crse(ic-1,jc+1,kc,n) & - + crse(ic-1,jc-1,kc,n) - crse(ic+1,jc-1,kc,n)) - - if (mask(i-1,j,k) .ne. not_covered) then - derives(ic,jc,XDER) = crse(ic+1,jc,kc,n) - crse(ic,jc,kc,n) - derives(ic,jc,X2DER) = zero - endif - if (mask(i+ratiox,j,k) .ne. not_covered) then - derives(ic,jc,XDER) = crse(ic,jc,kc,n) - crse(ic-1,jc,kc,n) - derives(ic,jc,X2DER) = zero - endif - if (mask(i-1,j,k) .ne. not_covered .and. & - mask(i+ratiox,j,k) .ne. not_covered) then - derives(ic,jc,XDER) = zero - endif - - if (mask(i,j-1,k) .ne. not_covered) then - derives(ic,jc,YDER) = crse(ic,jc+1,kc,n) - crse(ic,jc,kc,n) - derives(ic,jc,Y2DER) = zero - endif - if (mask(i,j+ratioy,k) .ne. not_covered) then - derives(ic,jc,YDER) = crse(ic,jc,kc,n) - crse(ic,jc-1,kc,n) - derives(ic,jc,Y2DER) = zero - endif - if (mask(i,j-1,k) .ne. not_covered .and. & - mask(i,j+ratioy,k) .ne. not_covered) then - derives(ic,jc,YDER) = zero - endif - if (( mask(i+ratiox,j+ratioy,k) .ne. not_covered ) .or. & - ( mask(i-1,j+ratioy,k) .ne. not_covered ) .or. & - ( mask(i+ratiox,j-1,k) .ne. not_covered ) .or. & - ( mask(i-1,j-1,k) .ne. not_covered ) ) then - - derives(ic,jc,XYDER) = zero - endif - enddo - enddo - - ! interpolate to fine grid - do joff = 0, ratioy - 1 - yy = (joff - half * ratioy + half)/ratioy - do jc = jclo,jchi - j = ratioy*jc + joff - do ioff = 0, ratiox - 1 - xx = (ioff - half * ratiox + half)/ratiox - do ic = iclo, ichi - i = ratiox*ic + ioff - bdry(i,j,k,n) = crse(ic,jc,kc,n) + xx*derives(ic,jc,XDER) & - + derives(ic,jc,X2DER)*xx**2 + yy*derives(ic,jc,YDER) & - + derives(ic,jc,Y2DER)*yy**2 + xx*yy*derives(ic,jc,XYDER) - enddo - enddo - enddo - enddo - enddo - - return -end subroutine FORT_BDINTERPZHI - -#undef NUMDERIV -#undef XDER -#undef YDER -#undef X2DER -#undef Y2DER -#undef XYDER -