RFLO_CalcGradConnBc.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
!******************************************************************************
! 
! Purpose: specify gradients at the current patch subject to connecting 
!          boundary condition (block interface or periodic)
!
! Description: if a region side has a patch with connecting bc, treatment
!              applied on whole side to omit complications at inter-patch 
!              junctions; non-conneting patches will be overruled by other bc
!
! Input: region        = data of current region
!        patch         = current patch
!        iConBc        = connecting Bc treatment flag (1:treated, 0: not yet)
!        iBegV, iEndV  = begin and end var index
!        iBegG, iEndG  = begin and end gradient index
!        var           = variables, the gradient of which to be determined
!  
! Output: gradi, gradj, gradk = gradients at the patch faces subject to 
!                connecting boundary conditions
!
! Notes: 1) gradients at edges and corners can be computed more accurately
!           if the required information of face vectors are available
!        2) three other routines contained in this file:
!           - RFLO_SetIndexRange to set index-range to whole region-side
!           - RFLO_CalcSideGrad to compute gradients on region-side
!           - RFLO_CopyPatchEdgeGrad to copy gradients to patch edges
!
!******************************************************************************
!
! $Id: RFLO_CalcGradConnBc.F90,v 1.3 2008/12/06 08:44:06 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE rflo_calcgradconnbc( region,patch,iConBc,iBegV,iEndV,iBegG,iEndG, &
                                var,gradi,gradj,gradk )

  USE moddatatypes
  USE modbndpatch, ONLY   : t_patch
  USE moddatastruct, ONLY : t_region
  USE modinterfaces, ONLY : rflo_getpatchindices, rflo_getpatchdirection, &
        rflo_getcelloffset, rflo_getnodeoffset, rflo_getdimensphys
  USE moderror
  USE modparameters
  IMPLICIT NONE

#include "Indexing.h"

! ... parameters
  TYPE(t_region) :: region
  TYPE(t_patch)  :: patch
  INTEGER        :: iconbc(6),ibegv,iendv,ibegg,iendg
  REAL(RFREAL), POINTER :: var(:,:), gradi(:,:), gradj(:,:), gradk(:,:)

! ... loop variables
  INTEGER :: i, j, k, l, lx, ly, lz

! ... local variables
  INTEGER :: ilev, lbound, bctype
  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend
  INTEGER :: inbeg, inend, jnbeg, jnend, knbeg, knend
  INTEGER :: idir, jdir, kdir
  INTEGER :: icoff,ijcoff,inoff,ijnoff
  INTEGER :: ijkc, im1jkc, ijm1kc, ijkm1c, im1jm1kc, im1jkm1c, ijm1km1c, &
             im1jm1km1c, ijkn, im1jkn, ijm1kn, ijkm1n, nvar

  REAL(RFREAL), POINTER :: aci(:,:), acj(:,:), ack(:,:)
  REAL(RFREAL), POINTER :: si(:,:), sj(:,:), sk(:,:), vol(:)
  REAL(RFREAL)          :: fnx, fny, fnz, rvol, fface(ibegv:iendv), tvar

!******************************************************************************

  CALL registerfunction( region%global,'RFLO_CalcGradConnBc',&
       'RFLO_CalcGradConnBc.F90' )
 
! get cell indices and geometrical parameters --------------------------------

  ilev   = region%currLevel
  lbound = patch%lbound
  bctype = patch%bcType

  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
                           jpcbeg,jpcend,kpcbeg,kpcend )

! complete the outer layers --------------------------------------------------

  IF ((bctype>=bc_tra_peri   .AND. bctype<=bc_tra_peri+  bc_range) .OR. &
      (bctype>=bc_rot_peri   .AND. bctype<=bc_rot_peri+  bc_range) .OR. &
      (bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range)) THEN

    CALL rflo_getpatchindices( region,patch,ilev,ibeg,iend,jbeg,jend,kbeg,kend )
    CALL rflo_getpatchdirection( patch,idir,jdir,kdir )
    CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
    CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )

    nvar = iendv - ibegv + 1

! - get dimensions and pointers

    aci => region%levels(ilev)%grid%c2eCoI
    acj => region%levels(ilev)%grid%c2eCoJ
    ack => region%levels(ilev)%grid%c2eCoK
    si  => region%levels(ilev)%grid%si
    sj  => region%levels(ilev)%grid%sj
    sk  => region%levels(ilev)%grid%sk
    vol => region%levels(ilev)%grid%vol

    IF (iconbc(lbound)==0) THEN
      IF (lbound==1 .OR. lbound==2) THEN
        jbeg=jpcbeg
        jend=jpcend
        kbeg=kpcbeg
        kend=kpcend
      ELSE IF (lbound==3 .OR. lbound==4) THEN
        ibeg=ipcbeg
        iend=ipcend
        kbeg=kpcbeg
        kend=kpcend
      ELSE IF (lbound==5 .OR. lbound==6) THEN
        ibeg=ipcbeg
        iend=ipcend
        jbeg=jpcbeg
        jend=jpcend
      ENDIF
      CALL rflo_setindexrange
      CALL rflo_calcsidegrad
      iconbc(lbound)=1
    ENDIF

    CALL rflo_setindexrange
    CALL rflo_copypatchedgegrad

  ENDIF

! finalize --------------------------------------------------------------------

  CALL deregisterfunction( region%global )

! ==============================================================================
! gradient completion subroutine for dummy cells at connecting boundaries
! ==============================================================================

CONTAINS

  SUBROUTINE rflo_setindexrange

! set index-range be used in RFLO_CalcSideGrad and RFLO_CopyPatchEdgeGrad ----
      
    IF (lbound==1 .OR. lbound==3 .OR. lbound==5) THEN
      inbeg = ibeg-idir*(region%nDumCells-1)
      jnbeg = jbeg-jdir*(region%nDumCells-1)
      knbeg = kbeg-kdir*(region%nDumCells-1)
      inend = iend-2*idir+1
      jnend = jend-2*jdir+1
      knend = kend-2*kdir+1
    ELSE
      inbeg = ibeg-2*idir
      jnbeg = jbeg-2*jdir
      knbeg = kbeg-2*kdir
      inend = iend-idir*(region%nDumCells-1)+1
      jnend = jend-jdir*(region%nDumCells-1)+1
      knend = kend-kdir*(region%nDumCells-1)+1
    ENDIF

  END SUBROUTINE rflo_setindexrange

  SUBROUTINE rflo_calcsidegrad

! ... local variables

  INTEGER :: inode,jnode,knode

! compute gradients at sides of current region having connecting bc --------

  inode = 0
  jnode = 0
  knode = 0
  IF( lbound==2 .OR. lbound==4 .OR. lbound==6 ) THEN
    inode = -idir
    jnode = -jdir
    knode = -kdir
  ENDIF

  IF (((knend>knbeg).AND.(jnend>jnbeg).AND.(inend>inbeg)).OR. &
      (lbound==1).OR.(lbound==3).OR.(lbound==5)) THEN

    DO k=knbeg+knode,knend
      DO j=jnbeg+jnode,jnend
        DO i=inbeg+inode,inend  

          ijkc    = indijk(i  ,j  ,k  ,icoff,ijcoff)
          im1jkc  = indijk(i-1,j  ,k  ,icoff,ijcoff)
          ijm1kc  = indijk(i  ,j-1,k  ,icoff,ijcoff)
          ijkm1c  = indijk(i  ,j  ,k-1,icoff,ijcoff)
          im1jm1kc= indijk(i-1,j-1,k  ,icoff,ijcoff)
          im1jkm1c= indijk(i-1,j  ,k-1,icoff,ijcoff)
          ijm1km1c= indijk(i  ,j-1,k-1,icoff,ijcoff)
          ijkn    = indijk(i  ,j  ,k  ,inoff,ijnoff)
          im1jkn  = indijk(i-1,j  ,k  ,inoff,ijnoff)
          ijm1kn  = indijk(i  ,j-1,k  ,inoff,ijnoff)
          ijkm1n  = indijk(i  ,j  ,k-1,inoff,ijnoff)

! - gradients at I-Face

! --- i-face of auxiliary control volume (face vectors pointed inwards/negative)

          fnx = .5_rfreal*(si(xcoord,ijkn)+si(xcoord,im1jkn))
          fny = .5_rfreal*(si(ycoord,ijkn)+si(ycoord,im1jkn))
          fnz = .5_rfreal*(si(zcoord,ijkn)+si(zcoord,im1jkn))

          DO l=ibegv,iendv
            lx=l+ibegg-ibegv 
            ly=lx+nvar
            lz=ly+nvar
            tvar=var(l,im1jkc)
            gradi(lx,im1jkn) = gradi(lx,im1jkn) - fnx*tvar
            gradi(ly,im1jkn) = gradi(ly,im1jkn) - fny*tvar
            gradi(lz,im1jkn) = gradi(lz,im1jkn) - fnz*tvar

            gradi(lx,ijkn) = gradi(lx,ijkn) + fnx*tvar
            gradi(ly,ijkn) = gradi(ly,ijkn) + fny*tvar
            gradi(lz,ijkn) = gradi(lz,ijkn) + fnz*tvar
          END DO

! --- j-face of auxiliary control volume (face vectors pointed inwards/negative)

          fnx = .5_rfreal*(sj(xcoord,ijkn)+sj(xcoord,im1jkn))
          fny = .5_rfreal*(sj(ycoord,ijkn)+sj(ycoord,im1jkn))
          fnz = .5_rfreal*(sj(zcoord,ijkn)+sj(zcoord,im1jkn))

          DO l=ibegv,iendv
            lx=l+ibegg-ibegv 
            ly=lx+nvar
            lz=ly+nvar
!            fface(l) = .25_RFREAL*(var(l,im1jkC)+var(l,ijkC)+ &
!                                   var(l,im1jm1kC)+var(l,ijm1kC))

            fface(l) = ack(1,ijkn)*var(l,im1jm1kc)+ack(2,ijkn)*var(l,ijm1kc)+ &
                       ack(3,ijkn)*var(l,ijkc)+ack(4,ijkn)*var(l,im1jkc)
                                 
            gradi(lx,ijm1kn) = gradi(lx,ijm1kn) - fnx*fface(l)
            gradi(ly,ijm1kn) = gradi(ly,ijm1kn) - fny*fface(l)
            gradi(lz,ijm1kn) = gradi(lz,ijm1kn) - fnz*fface(l)

            gradi(lx,ijkn) = gradi(lx,ijkn) + fnx*fface(l)
            gradi(ly,ijkn) = gradi(ly,ijkn) + fny*fface(l)
            gradi(lz,ijkn) = gradi(lz,ijkn) + fnz*fface(l)
          END DO

! --- k-face of auxiliary control volume (face vectors pointed inwards/negative)

          fnx = .5_rfreal*(sk(xcoord,im1jkn)+sk(xcoord,ijkn))
          fny = .5_rfreal*(sk(ycoord,im1jkn)+sk(ycoord,ijkn))
          fnz = .5_rfreal*(sk(zcoord,im1jkn)+sk(zcoord,ijkn))

          DO l=ibegv,iendv
            lx=l+ibegg-ibegv 
            ly=lx+nvar
            lz=ly+nvar
!            fface(l) = .25_RFREAL*(var(l,im1jkC)+var(l,ijkC)+ &
!                                   var(l,im1jkm1C)+var(l,ijkm1C))

            fface(l) = acj(1,ijkn)*var(l,im1jkm1c)+acj(2,ijkn)*var(l,im1jkc)+ &
                       acj(3,ijkn)*var(l,ijkc)+acj(4,ijkn)*var(l,ijkm1c)

            gradi(lx,ijkm1n) = gradi(lx,ijkm1n) - fnx*fface(l)
            gradi(ly,ijkm1n) = gradi(ly,ijkm1n) - fny*fface(l)
            gradi(lz,ijkm1n) = gradi(lz,ijkm1n) - fnz*fface(l)

            gradi(lx,ijkn) = gradi(lx,ijkn) + fnx*fface(l)
            gradi(ly,ijkn) = gradi(ly,ijkn) + fny*fface(l)
            gradi(lz,ijkn) = gradi(lz,ijkn) + fnz*fface(l)
          END DO

! - gradients at J-Face

! --- j-face of auxiliary control volume (face vectors pointed inwards/negative)

          fnx = .5_rfreal*(sj(xcoord,ijkn)+sj(xcoord,ijm1kn))
          fny = .5_rfreal*(sj(ycoord,ijkn)+sj(ycoord,ijm1kn))
          fnz = .5_rfreal*(sj(zcoord,ijkn)+sj(zcoord,ijm1kn))

          DO l=ibegv,iendv
            lx=l+ibegg-ibegv 
            ly=lx+nvar
            lz=ly+nvar
            tvar=var(l,ijm1kc)
            gradj(lx,ijm1kn) = gradj(lx,ijm1kn) - fnx*tvar
            gradj(ly,ijm1kn) = gradj(ly,ijm1kn) - fny*tvar
            gradj(lz,ijm1kn) = gradj(lz,ijm1kn) - fnz*tvar

            gradj(lx,ijkn) = gradj(lx,ijkn) + fnx*tvar
            gradj(ly,ijkn) = gradj(ly,ijkn) + fny*tvar
            gradj(lz,ijkn) = gradj(lz,ijkn) + fnz*tvar
          END DO

! --- i-face of auxiliary control volume (face vectors pointed inwards/negative)

          fnx = .5_rfreal*(si(xcoord,ijkn)+si(xcoord,ijm1kn))
          fny = .5_rfreal*(si(ycoord,ijkn)+si(ycoord,ijm1kn))
          fnz = .5_rfreal*(si(zcoord,ijkn)+si(zcoord,ijm1kn))

          DO l=ibegv,iendv
            lx=l+ibegg-ibegv 
            ly=lx+nvar
            lz=ly+nvar
!            fface(l) = .25_RFREAL*(var(l,im1jkC)+var(l,ijkC)+ &
!                                   var(l,im1jm1kC)+var(l,ijm1kC))

            fface(l) = ack(1,ijkn)*var(l,im1jm1kc)+ack(2,ijkn)*var(l,ijm1kc)+ &
                       ack(3,ijkn)*var(l,ijkc)+ack(4,ijkn)*var(l,im1jkc)
                                 
            gradj(lx,im1jkn) = gradj(lx,im1jkn) - fnx*fface(l)
            gradj(ly,im1jkn) = gradj(ly,im1jkn) - fny*fface(l)
            gradj(lz,im1jkn) = gradj(lz,im1jkn) - fnz*fface(l)

            gradj(lx,ijkn) = gradj(lx,ijkn) + fnx*fface(l)
            gradj(ly,ijkn) = gradj(ly,ijkn) + fny*fface(l)
            gradj(lz,ijkn) = gradj(lz,ijkn) + fnz*fface(l)
          END DO

! --- k-face of auxiliary control volume (face vectors pointed inwards/negative)

          fnx = .5_rfreal*(sk(xcoord,ijm1kn)+sk(xcoord,ijkn))
          fny = .5_rfreal*(sk(ycoord,ijm1kn)+sk(ycoord,ijkn))
          fnz = .5_rfreal*(sk(zcoord,ijm1kn)+sk(zcoord,ijkn))

          DO l=ibegv,iendv
            lx=l+ibegg-ibegv 
            ly=lx+nvar
            lz=ly+nvar
!            fface(l) = .25_RFREAL*(var(l,ijm1kC)+var(l,ijkC)+ &
!                                   var(l,ijm1km1C)+var(l,ijkm1C))

            fface(l) = aci(1,ijkn)*var(l,ijm1km1c)+aci(2,ijkn)*var(l,ijkm1c)+ &
                       aci(3,ijkn)*var(l,ijkc)+aci(4,ijkn)*var(l,ijm1kc)
                                 
            gradj(lx,ijkm1n) = gradj(lx,ijkm1n) - fnx*fface(l)
            gradj(ly,ijkm1n) = gradj(ly,ijkm1n) - fny*fface(l)
            gradj(lz,ijkm1n) = gradj(lz,ijkm1n) - fnz*fface(l)

            gradj(lx,ijkn) = gradj(lx,ijkn) + fnx*fface(l)
            gradj(ly,ijkn) = gradj(ly,ijkn) + fny*fface(l)
            gradj(lz,ijkn) = gradj(lz,ijkn) + fnz*fface(l)
          END DO

! - gradients at K-Face

! --- k-face of auxiliary control volume (face vectors pointed inwards/negative)

          fnx = .5_rfreal*(sk(xcoord,ijkn)+sk(xcoord,ijkm1n))
          fny = .5_rfreal*(sk(ycoord,ijkn)+sk(ycoord,ijkm1n))
          fnz = .5_rfreal*(sk(zcoord,ijkn)+sk(zcoord,ijkm1n))

          DO l=ibegv,iendv
            lx=l+ibegg-ibegv 
            ly=lx+nvar
            lz=ly+nvar
            tvar=var(l,ijkm1c)
            gradk(lx,ijkm1n) = gradk(lx,ijkm1n) - fnx*tvar
            gradk(ly,ijkm1n) = gradk(ly,ijkm1n) - fny*tvar
            gradk(lz,ijkm1n) = gradk(lz,ijkm1n) - fnz*tvar

            gradk(lx,ijkn) = gradk(lx,ijkn) + fnx*tvar
            gradk(ly,ijkn) = gradk(ly,ijkn) + fny*tvar
            gradk(lz,ijkn) = gradk(lz,ijkn) + fnz*tvar
          END DO

! --- j-face of auxiliary control volume (face vectors pointed inwards/negative)

          fnx = .5_rfreal*(sj(xcoord,ijkn)+sj(xcoord,ijkm1n))
          fny = .5_rfreal*(sj(ycoord,ijkn)+sj(ycoord,ijkm1n))
          fnz = .5_rfreal*(sj(zcoord,ijkn)+sj(zcoord,ijkm1n))

          DO l=ibegv,iendv
            lx=l+ibegg-ibegv 
            ly=lx+nvar
            lz=ly+nvar
!            fface(l) = .25_RFREAL*(var(l,ijkm1C)+var(l,ijkC)+ &
!                                   var(l,ijm1km1C)+var(l,ijm1kC))

            fface(l) = aci(1,ijkn)*var(l,ijm1km1c)+aci(2,ijkn)*var(l,ijkm1c)+ &
                       aci(3,ijkn)*var(l,ijkc)+aci(4,ijkn)*var(l,ijm1kc)

            gradk(lx,ijm1kn) = gradk(lx,ijm1kn) - fnx*fface(l)
            gradk(ly,ijm1kn) = gradk(ly,ijm1kn) - fny*fface(l)
            gradk(lz,ijm1kn) = gradk(lz,ijm1kn) - fnz*fface(l)

            gradk(lx,ijkn) = gradk(lx,ijkn) + fnx*fface(l)
            gradk(ly,ijkn) = gradk(ly,ijkn) + fny*fface(l)
            gradk(lz,ijkn) = gradk(lz,ijkn) + fnz*fface(l)
          END DO

! --- i-face of auxiliary control volume (face vectors pointed inwards/negative)

          fnx = .5_rfreal*(si(xcoord,ijkm1n)+si(xcoord,ijkn))
          fny = .5_rfreal*(si(ycoord,ijkm1n)+si(ycoord,ijkn))
          fnz = .5_rfreal*(si(zcoord,ijkm1n)+si(zcoord,ijkn))

          DO l=ibegv,iendv
            lx=l+ibegg-ibegv 
            ly=lx+nvar
            lz=ly+nvar
!            fface(l) = .25_RFREAL*(var(l,im1jkC)+var(l,ijkC)+ &
!                                   var(l,im1jkm1C)+var(l,ijkm1C))

            fface(l) = acj(1,ijkn)*var(l,im1jkm1c)+acj(2,ijkn)*var(l,im1jkc)+ &
                       acj(3,ijkn)*var(l,ijkc)+acj(4,ijkn)*var(l,ijkm1c)

            gradk(lx,im1jkn) = gradk(lx,im1jkn) - fnx*fface(l)
            gradk(ly,im1jkn) = gradk(ly,im1jkn) - fny*fface(l)
            gradk(lz,im1jkn) = gradk(lz,im1jkn) - fnz*fface(l)

            gradk(lx,ijkn) = gradk(lx,ijkn) + fnx*fface(l)
            gradk(ly,ijkn) = gradk(ly,ijkn) + fny*fface(l)
            gradk(lz,ijkn) = gradk(lz,ijkn) + fnz*fface(l)
          END DO

        END DO ! i
      END DO   ! j
    END DO     ! k

  END IF

  IF (lbound==2) THEN

! - finish the last i-face (face vectors pointed inwards/negative)

    DO k=knbeg,knend-1 
      DO j=jnbeg,jnend-1  

        im1jkc  = indijk(inend  ,j  ,k  ,icoff,ijcoff)
        ijkn    = indijk(inend+1,j  ,k  ,inoff,ijnoff)
        im1jkn  = indijk(inend  ,j  ,k  ,inoff,ijnoff)

        fnx = .5_rfreal*(si(xcoord,ijkn)+si(xcoord,im1jkn))
        fny = .5_rfreal*(si(ycoord,ijkn)+si(ycoord,im1jkn))
        fnz = .5_rfreal*(si(zcoord,ijkn)+si(zcoord,im1jkn))
        DO l=ibegv,iendv
          lx=l+ibegg-ibegv 
          ly=lx+nvar
          lz=ly+nvar
          tvar=var(l,im1jkc)
          gradi(lx,im1jkn) = gradi(lx,im1jkn) - fnx*tvar
          gradi(ly,im1jkn) = gradi(ly,im1jkn) - fny*tvar
          gradi(lz,im1jkn) = gradi(lz,im1jkn) - fnz*tvar
        END DO

      END DO ! j
    END DO   ! k

  ELSEIF (lbound==4) THEN

! - finish the last j-face (face vectors pointed inwards/negative)

    DO k=knbeg,knend-1
      DO i=inbeg,inend-1  
 
        ijm1kc  = indijk(i  ,jnend  ,k  ,icoff,ijcoff)
        ijkn    = indijk(i  ,jnend+1,k  ,inoff,ijnoff)
        ijm1kn  = indijk(i  ,jnend  ,k  ,inoff,ijnoff)

        fnx = .5_rfreal*(sj(xcoord,ijkn)+sj(xcoord,ijm1kn))
        fny = .5_rfreal*(sj(ycoord,ijkn)+sj(ycoord,ijm1kn))
        fnz = .5_rfreal*(sj(zcoord,ijkn)+sj(zcoord,ijm1kn))
        DO l=ibegv,iendv
          lx=l+ibegg-ibegv 
          ly=lx+nvar
          lz=ly+nvar
          tvar=var(l,ijm1kc)
          gradj(lx,ijm1kn) = gradj(lx,ijm1kn) - fnx*tvar
          gradj(ly,ijm1kn) = gradj(ly,ijm1kn) - fny*tvar
          gradj(lz,ijm1kn) = gradj(lz,ijm1kn) - fnz*tvar
        END DO

      END DO ! i
    END DO   ! k

  ELSEIF (lbound==6) THEN

! - finish the last k-face (face vectors pointed inwards/negative)

    DO j=jnbeg,jnend-1 
      DO i=inbeg,inend-1  

        ijkm1c  = indijk(i  ,j  ,knend  ,icoff,ijcoff)
        ijkn    = indijk(i  ,j  ,knend+1,inoff,ijnoff)
        ijkm1n  = indijk(i  ,j  ,knend  ,inoff,ijnoff)

        fnx = .5_rfreal*(sk(xcoord,ijkn)+sk(xcoord,ijkm1n))
        fny = .5_rfreal*(sk(ycoord,ijkn)+sk(ycoord,ijkm1n))
        fnz = .5_rfreal*(sk(zcoord,ijkn)+sk(zcoord,ijkm1n))
        DO l=ibegv,iendv
          lx=l+ibegg-ibegv 
          ly=lx+nvar
          lz=ly+nvar
          tvar=var(l,ijkm1c)
          gradk(lx,ijkm1n) = gradk(lx,ijkm1n) - fnx*tvar
          gradk(ly,ijkm1n) = gradk(ly,ijkm1n) - fny*tvar
          gradk(lz,ijkm1n) = gradk(lz,ijkm1n) - fnz*tvar
        END DO

      END DO ! i
    END DO   ! j

  END IF

! finally, division by face averaged volume

  DO k=knbeg,knend
    DO j=jnbeg,jnend
      DO i=inbeg,inend  

        ijkc    = indijk(i  ,j  ,k  ,icoff,ijcoff)
        im1jkc  = indijk(i-1,j  ,k  ,icoff,ijcoff)
        ijm1kc  = indijk(i  ,j-1,k  ,icoff,ijcoff)
        ijkm1c  = indijk(i  ,j  ,k-1,icoff,ijcoff)
        ijkn    = indijk(i  ,j  ,k  ,inoff,ijnoff)
       
        rvol = 2.0_rfreal/(vol(im1jkc)+vol(ijkc))
        DO l=ibegg,iendg
          gradi(l,ijkn) = gradi(l,ijkn)*rvol
        END DO

        rvol = 2.0_rfreal/(vol(ijm1kc)+vol(ijkc))
        DO l=ibegg,iendg
          gradj(l,ijkn) = gradj(l,ijkn)*rvol
        END DO

        rvol = 2.0_rfreal/(vol(ijkm1c)+vol(ijkc))
        DO l=ibegg,iendg
          gradk(l,ijkn) = gradk(l,ijkn)*rvol
        END DO

      END DO ! i
    END DO   ! j
  END DO     ! k

  END SUBROUTINE rflo_calcsidegrad

 
  SUBROUTINE rflo_copypatchedgegrad

! ... local variables

  INTEGER :: ijkn1

! copy gradients to patch edges ---------------------------------------------

  IF ((lbound==1).OR.(lbound==2)) THEN 
    DO i=inbeg,inend
      DO k=knbeg,knend-1
        ijkn    = indijk(i   ,jnbeg  ,k,inoff,ijnoff)
        ijkn1   = indijk(ibeg,jnbeg  ,k,inoff,ijnoff)
        gradj(:,ijkn)=gradj(:,ijkn1)
        ijkn    = indijk(i   ,jnend  ,k,inoff,ijnoff)
        ijkn1   = indijk(ibeg,jnend  ,k,inoff,ijnoff)
        gradj(:,ijkn)=gradj(:,ijkn1)
      END DO  
      DO j=jnbeg,jnend-1
        ijkn    = indijk(i   ,j,knbeg  ,inoff,ijnoff)
        ijkn1   = indijk(ibeg,j,knbeg  ,inoff,ijnoff)
        gradk(:,ijkn)=gradk(:,ijkn1)
        ijkn    = indijk(i   ,j,knend  ,inoff,ijnoff)
        ijkn1   = indijk(ibeg,j,knend  ,inoff,ijnoff)
        gradk(:,ijkn)=gradk(:,ijkn1)
      END DO  
    END DO   ! i
  ELSEIF ((lbound==3).OR.(lbound==4)) THEN
    DO j=jnbeg,jnend
      DO k=knbeg,knend-1
        ijkn    = indijk(inbeg  ,j   ,k,inoff,ijnoff)
        ijkn1   = indijk(inbeg  ,jbeg,k,inoff,ijnoff)
        gradi(:,ijkn)=gradi(:,ijkn1)
        ijkn    = indijk(inend  ,j   ,k,inoff,ijnoff)
        ijkn1   = indijk(inend  ,jbeg,k,inoff,ijnoff)
        gradi(:,ijkn)=gradi(:,ijkn1)
      END DO  
      DO i=inbeg,inend-1
        ijkn    = indijk(i,j   ,knbeg  ,inoff,ijnoff)
        ijkn1   = indijk(i,jbeg,knbeg  ,inoff,ijnoff)
        gradk(:,ijkn)=gradk(:,ijkn1)
        ijkn    = indijk(i,j   ,knend  ,inoff,ijnoff)
        ijkn1   = indijk(i,jbeg,knend  ,inoff,ijnoff)
        gradk(:,ijkn)=gradk(:,ijkn1)
      END DO  
    END DO   ! j
  ELSEIF ((lbound==5).OR.(lbound==6)) THEN
    DO k=knbeg,knend
      DO i=inbeg,inend-1
        ijkn    = indijk(i,jnbeg  ,k   ,inoff,ijnoff)
        ijkn1   = indijk(i,jnbeg  ,kbeg,inoff,ijnoff)
        gradj(:,ijkn)=gradj(:,ijkn1)
        ijkn    = indijk(i,jnend  ,k   ,inoff,ijnoff)
        ijkn1   = indijk(i,jnend  ,kbeg,inoff,ijnoff)
        gradj(:,ijkn)=gradj(:,ijkn1)
      END DO  
      DO j=jnbeg,jnend-1
        ijkn    = indijk(inbeg  ,j,k   ,inoff,ijnoff)
        ijkn1   = indijk(inbeg  ,j,kbeg,inoff,ijnoff)
        gradi(:,ijkn)=gradi(:,ijkn1)
        ijkn    = indijk(inend  ,j,k   ,inoff,ijnoff)
        ijkn1   = indijk(inend  ,j,kbeg,inoff,ijnoff)
        gradi(:,ijkn)=gradi(:,ijkn1)
      END DO  
    END DO   ! k
  END IF     ! lbound 

  END SUBROUTINE rflo_copypatchedgegrad

END SUBROUTINE rflo_calcgradconnbc

!******************************************************************************
!
! RCS Revision history:
! 
! $Log: RFLO_CalcGradConnBc.F90,v $
! Revision 1.3  2008/12/06 08:44:06  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:19  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2004/11/29 21:25:15  wasistho
! lower to upper case
!
! Revision 1.11  2004/08/03 22:48:43  wasistho
! changed cell2edge averaging to grid dependent avg
!
! Revision 1.10  2004/01/22 03:56:27  wasistho
! extrapolated to patch edges from different direction of side faces
!
! Revision 1.9  2003/12/07 17:48:39  jiao
! Chaged in order to compile with Intel compilers on NCSA machines
!
! Revision 1.8  2003/11/20 16:40:34  mdbrandy
! Backing out RocfluidMP changes from 11-17-03
!
! Revision 1.5  2003/05/15 02:57:01  jblazek
! Inlined index function.
!
! Revision 1.4  2003/01/10 17:58:42  jblazek
! Added missing explicit interfaces.
!
! Revision 1.3  2002/09/27 00:57:09  jblazek
! Changed makefiles - no makelinks needed.
!
! Revision 1.2  2002/09/05 17:40:19  jblazek
! Variable global moved into regions().
!
! Revision 1.1  2002/09/02 22:58:54  wasistho
! RFLO grad routines migrated from rocflo to libflo
!
! Revision 1.3  2002/07/22 22:59:11  jblazek
! Some more clean up.
!
! Revision 1.2  2002/07/19 23:43:00  wasistho
! made compliant with CODING RULE
!
!******************************************************************************