c-----------------------------------------------------------------------
c Three-dimensional prolongation operator for AMROC.
c A fine grid value is replaced by the value of a linear function
c through the neighbouring coarse grid values at the center
c of the particular fine grid cell.
c
c Interface:
c mfx,mfy,mfz := shape of fine grid function
c mcx,mcy,mcz := shape of coarse grid function
c
c uf(,,) := fine grid function
c uc(,,) := coarse grid function
c
c lbc(3) := lower bound for coarse grid
c ubc(3) := upper bound for coarse grid
c lbf(3) := lower bound for fine grid
c ubf(3) := upper bound for fine grid
c lbr(3) := lower bound for region prolongation desired
c ufr(3) := upper bound for region prolongation desired
c shaper(3) := shape of region prolongation desired
c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@cacr.caltech.edu
c-----------------------------------------------------------------------
subroutine prolong3(uc,mcx,mcy,mcz,lbc,ubc,
& uf,mfx,mfy,mfz,lbf,ubf,
& lbr,ubr,shaper,meqn,mbc)
implicit none
integer meqn,mbc,mcx,mcy,mcz,mfx,mfy,mfz
integer shaper(3)
double precision uf(meqn,mfx,mfy,mfz), uc(meqn,mcx,mcy,mcz)
integer lbf(3), ubf(3),
& lbc(3), ubc(3),
& lbr(3), ubr(3),
& getindx
c Local variables
integer i, j, k, m, ic, jc, kc, mic, mjc, mkc,
& stridec, stridef,
& ifine, ics, jfine, jcs, kfine, kcs
double precision eta1, eta2, eta3
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c See definition of member-function extents() in BBox.h
c for calculation of stride
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
stridec = (ubc(1) - lbc(1))/(mcx-1)
stridef = (ubf(1) - lbf(1))/(mfx-1)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Prolongation region is defined on the domain of the fine grid.
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
do 10 k=lbr(3), ubr(3), stridef
kfine = getindx(k, lbf(3), stridef)
kcs = getindx(k, lbc(3), stridec)
kc = k - lbc(3)
mkc = kc - (kc/stridec)*stridec
if(mkc .lt. stridec*0.5) then
kcs = kcs - 1
end if
if(kcs+1.gt.mcz .or. kcs.lt.1 .or.
& kfine.gt.mfz .or. kfine.lt.1) goto 10
kc = kc + stridec*0.5
mkc = kc - (kc/stridec)*stridec
eta3 = (mkc+0.5d0*stridef) / stridec
do 14 j=lbr(2), ubr(2), stridef
jfine = getindx(j, lbf(2), stridef)
jcs = getindx(j, lbc(2), stridec)
jc = j - lbc(2)
mjc = jc - (jc/stridec)*stridec
if(mjc .lt. stridec*0.5) then
jcs = jcs - 1
end if
if(jcs+1.gt.mcy .or. jcs.lt.1 .or.
& jfine.gt.mfy .or. jfine.lt.1) goto 14
jc = jc + stridec*0.5
mjc = jc - (jc/stridec)*stridec
eta2 = (mjc+0.5d0*stridef) / stridec
do 16 i=lbr(1), ubr(1), stridef
ifine = getindx(i, lbf(1), stridef)
ics = getindx(i, lbc(1), stridec)
ic = i - lbc(1)
mic = ic - (ic/stridec)*stridec
if(mic .lt. stridec*0.5) then
ics = ics - 1
end if
if(ics+1.gt.mcx .or. ics.lt.1 .or.
& ifine.gt.mfx .or. ifine.lt.1) goto 16
ic = ic + stridec*0.5
mic = ic - (ic/stridec)*stridec
eta1 = (mic+0.5d0*stridef) / stridec
do 20 m=1, meqn
uf(m,ifine,jfine,kfine) = (1.d0-eta3)*
& ((1.d0-eta2)*((1.d0-eta1)*uc(m,ics, jcs ,kcs ) +
& eta1 *uc(m,ics+1,jcs, kcs ))+
& eta2 *((1.d0-eta1)*uc(m,ics, jcs+1,kcs ) +
& eta1 *uc(m,ics+1,jcs+1,kcs )))+
& eta3 *
& ((1.d0-eta2)*((1.d0-eta1)*uc(m,ics, jcs ,kcs+1) +
& eta1 *uc(m,ics+1,jcs, kcs+1))+
& eta2 *((1.d0-eta1)*uc(m,ics, jcs+1,kcs+1) +
& eta1 *uc(m,ics+1,jcs+1,kcs+1)))
20 continue
16 continue
14 continue
10 continue
return
end