3d/operators/prolong3.f

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     Author: Ralf Deiterding
c-----------------------------------------------------------------------

      subroutine prolong3(uc,mcx,mcy,mcz,lbc,ubc,
     &     uf,mfx,mfy,mfz,lbf,ubf,
     &     lbr,ubr,shaper,args,argc)
      implicit none
      common /ghosts/ mbc
      integer mbc,mcx,mcy,mcz,mfx,mfy,mfz

      integer  shaper(3)
      integer  args, argc            
      
      real*8   uf(args,mfx,mfy,mfz), uc(args,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, 
     &     refine, stridec, stridef,
     &     ifine, ics, jfine, jcs, kfine, kcs, meqn
      
      real*8   eta1, eta2, eta3

      meqn = args

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)
      refine = stridec/stridef
      
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
         kc = kc + stridec*0.5
         mkc = kc - (kc/stridec)*stridec
         eta3 = (mkc+0.5d0*stridef) / stridec
         
         do 10 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
            jc = jc + stridec*0.5
            mjc = jc - (jc/stridec)*stridec
            eta2 = (mjc+0.5d0*stridef) / stridec
            
            do 10 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
               ic = ic + stridec*0.5
               mic = ic - (ic/stridec)*stridec
               eta1 = (mic+0.5d0*stridef) / stridec
 
!               if(ics+1 .gt. mcx .or. 
!     &              jcs+1 .gt. mcy .or. 
!     &              kcs+1 .gt. mcz .or. 
!     &              ics .lt. 1 .or. jcs .lt. 1 .or. kcs .lt. 1) then
!                  uf(m,ifine,jfine,kfine) = -100000.0
!                  write(0,*)'ERROR in prolongation: ',ics,jcs,kcs
!               endif

               do 10 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)))
               
 10   continue

      return
      end