!> bound.F
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> ibound(Bound,Var,ioff,joff,koff)
!!
      subroutine ibound(Bound,Var,ioff,joff,koff)
!
! ... Local variables ..................................................
!
#include "param.inc"
c
c     needs to be done in the order i, j, k to fill in all the spaces
c
c     ioff used for things that are dimensioned ni+1 instead of ni
#include "global_dims.inc"
!
! ... Parameter variables ..............................................
!
      dimension Bound(ilow:ihigh,jlow:jhigh,klow:khigh)
!> Description of "Bound"
      dimension Var(ilow:ihigh,jlow:jhigh,klow:khigh)
!> Description of "Var"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in bound.F::ibound(Bound,Var,ioff,joff,koff)"
#endif

      if ( .not. MHD_PE ) return
c
      if (i_global_off .eq. 0 ) then
*                                           low i boundary
*             mirrors about i=1
         do k=1,nl_k+koff
            do j=1,nl_j+joff
               do i=1,no2
                  Var(-i+1,j,k) = Bound(i,j,k)
               enddo
            enddo
         enddo
       endif
       if ( i_global_off+nl_i .eq. ni_global ) then
*                                              high i boundary
*               straight copy from (ni-no2) to ni+
         do k=1,nl_k+koff
            do j=1,nl_j+joff
               do i=1,no2
                  Var(nl_i+i+ioff,j,k) = Bound(nl_i-no2+i,j,k)
               enddo
            enddo
         enddo
       endif
!
! ... End subroutine ibound(Bound,Var,ioff,joff,koff) ..................
!
      return
      end
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!> kbound(Bound,Var,ioff,joff,koff)
!!
!! @param koff - used for things that are dimensioned nk+1 instead of nk
      subroutine kbound(Bound,Var,ioff,joff,koff)
!
! ... Local variables ..................................................
!
#include "param.inc"
#include "global_dims.inc"
!
! ... Parameter variables ..............................................
!
      dimension Bound(ilow:ihigh,jlow:jhigh,klow:khigh)
!> Description of "Bound"
      dimension Var(ilow:ihigh,jlow:jhigh,klow:khigh)
!> Description of "Var"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in bound.F::kbound(Bound,Var,joff,koff)"
#endif

      if ( .not. MHD_PE ) return
c
      if (k_global_off .eq. 0 ) then
*                                           low k boundary
         do k=1,no2
            do j=1,nl_j+joff
               do i=1,nl_i+ioff
                  Var(i,j,-no2+k) = Bound(i,j,k)
               enddo
            enddo
         enddo
       endif
       if ( k_global_off+nl_k .eq. nk_global ) then
*                                              high k boundary
*             high boundary has nl_k => nk 
         do k=1,no2
            do j=1,nl_j+joff
               do i=1,nl_i+ioff
                  Var(i,j,k+nl_k+koff) = Bound(i,j,nl_k-no2+k)
               enddo
            enddo
         enddo
       endif
!
! ... End subroutine kbound(Bound,Var,ioff,joff,koff)...................
!
      return
      end
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!> jbound(Bound,Var,ioff,joff,koff)
!!
!! @param joff -  used for things that are dimensioned nj+1 instead of nj

      subroutine jbound(Bound,Var,ioff,joff,koff)
!
! ... Local variables ..................................................
!
#include "param.inc"
#include "global_dims.inc"
!
! ... Parameter variables ..............................................
!
      dimension Bound(ilow:ihigh,jlow:jhigh,klow:khigh)
!> Description of "Bound"
      dimension Var(ilow:ihigh,jlow:jhigh,klow:khigh)
!> Description of "Var"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in bound.F::jbound(Bound,Var,ioff,joff,koff)"
#endif

      if ( .not. MHD_PE ) return
      if (j_global_off .eq. 0 ) then
*                                           low j boundary
         do k=1,nl_k+koff
            do j=1,no2
               do i=1,nl_i+ioff
                  Var(i,1-j,k) = Bound(i,j,k)
               enddo
            enddo
         enddo
       endif
       if ( j_global_off+nl_j .eq. nj_global ) then
*                                              high j boundary
*             high boundary has nl_j => nj
         do k=1,nl_k+koff
            do j=1,no2
               do i=1,nl_i+ioff
                  Var(i,nl_j+j+joff,k) = Bound(i,nl_j+1-j,k)
               enddo
            enddo
         enddo
       endif
!
! ... End subroutine jbound(Bound,Var,ioff,joff,koff)...................
!
      return
      end
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!> edge_fix
!!   edges needed for the magnetic field calculation
!!
!!   kluge for single grid code
!!
      subroutine edge_fix
!
! ... Local variables ..................................................
!
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in bound.F::edge_fix()"
#endif

      nl_k2 = nl_k/2
c
	if ( .not. mhd_pe ) return
c
c        ij edge
c
      if ( ilow_bound .and. jlow_bound ) then
         do k=1,nl_k/2
            do j=1,no2
               do i=1,no2
                  rho_g(-no2+i,1-j,k) = rho_g(-no2+i,j,k+nl_k2)
                  vx_g(-no2+i,1-j,k) = vx_g(-no2+i,j,k+nl_k2)
                  vy_g(-no2+i,1-j,k) = vy_g(-no2+i,j,k+nl_k2)
                  vz_g(-no2+i,1-j,k) = vz_g(-no2+i,j,k+nl_k2)
                  rho_g(-no2+i,1-j,k+nl_k2) = rho_g(-no2+i,j,k)
                  vx_g(-no2+i,1-j,k+nl_k2) = vx_g(-no2+i,j,k)
                  vy_g(-no2+i,1-j,k+nl_k2) = vy_g(-no2+i,j,k)
                  vz_g(-no2+i,1-j,k+nl_k2) = vz_g(-no2+i,j,k)
               enddo
            enddo
         enddo
       endif
      if ( ihigh_bound .and. jlow_bound ) then
         do k=1,nl_k/2
            do j=1,no2
               do i=1,no2
                  rho_g(nl_i+i,1-j,k) = rho_g(nl_i+i,j,k+nl_k2)
                  vx_g(nl_i+i,1-j,k) = vx_g(nl_i+i,j,k+nl_k2)
                  vy_g(nl_i+i,1-j,k) = vy_g(nl_i+i,j,k+nl_k2)
                  vz_g(nl_i+i,1-j,k) = vz_g(nl_i+i,j,k+nl_k2)
                  rho_g(nl_i+i,1-j,k+nl_k2) = rho_g(nl_i+i,j,k)
                  vx_g(nl_i+i,1-j,k+nl_k2) = vx_g(nl_i+i,j,k)
                  vy_g(nl_i+i,1-j,k+nl_k2) = vy_g(nl_i+i,j,k)
                  vz_g(nl_i+i,1-j,k+nl_k2) = vz_g(nl_i+i,j,k)
               enddo
            enddo
         enddo
       endif
      if ( ilow_bound .and. jhigh_bound ) then
         do k=1,nl_k/2
            do j=1,no2
               do i=1,no2
                  rho_g(-no2+i,nl_j+j,k) =
     $                 rho_g(-no2+i,nl_j+1-j,k+nl_k2)
                  vx_g(-no2+i,nl_j+j,k) = vx_g(-no2+i,nl_j+1-j,k+nl_k2)
                  vy_g(-no2+i,nl_j+j,k) = vy_g(-no2+i,nl_j+1-j,k+nl_k2)
                  vz_g(-no2+i,nl_j+j,k) = vz_g(-no2+i,nl_j+1-j,k+nl_k2)
                  rho_g(-no2+i,nl_j+j,k+nl_k2) =
     $                 rho_g(-no2+i,nl_j+1-j,k)
                  vx_g(-no2+i,nl_j+j,k+nl_k2) = vx_g(-no2+i,nl_j+1-j,k)
                  vy_g(-no2+i,nl_j+j,k+nl_k2) = vy_g(-no2+i,nl_j+1-j,k)
                  vz_g(-no2+i,nl_j+j,k+nl_k2) = vz_g(-no2+i,nl_j+1-j,k)
               enddo
            enddo
         enddo
       endif
      if ( ihigh_bound .and. jlow_bound ) then
         do k=1,nl_k/2
            do j=1,no2
               do i=1,no2
                  rho_g(nl_i+i,nl_j+j,k) =
     $                 rho_g(nl_i+i,nl_j+1-j,k+nl_k2)
                  vx_g(nl_i+i,nl_j+j,k) = vx_g(nl_i+i,nl_j+1-j,k+nl_k2)
                  vy_g(nl_i+i,nl_j+j,k) = vy_g(nl_i+i,nl_j+1-j,k+nl_k2)
                  vz_g(nl_i+i,nl_j+j,k) = vz_g(nl_i+i,nl_j+1-j,k+nl_k2)
                  rho_g(nl_i+i,nl_j+j,k+nl_k2) =
     $                 rho_g(nl_i+i,nl_j+1-j,k)
                  vx_g(nl_i+i,nl_j+j,k+nl_k2) = vx_g(nl_i+i,nl_j+1-j,k)
                  vy_g(nl_i+i,nl_j+j,k+nl_k2) = vy_g(nl_i+i,nl_j+1-j,k)
                  vz_g(nl_i+i,nl_j+j,k+nl_k2) = vz_g(nl_i+i,nl_j+1-j,k)
               enddo
            enddo
         enddo
       endif
c
c       jk edges
c
       if ( jlow_bound .and. klow_bound ) then
       do k=1,no2
          do j=1,no2
             do i =1,nl_i
                rho_g(i,-no2+j,-no2+k) = rho_g(i,-no2+j,nl_k-no2+k)
                vx_g(i,-no2+j,-no2+k) = vx_g(i,-no2+j,nl_k-no2+k)
                vy_g(i,-no2+j,-no2+k) = vy_g(i,-no2+j,nl_k-no2+k)
                vz_g(i,-no2+j,-no2+k) = vz_g(i,-no2+j,nl_k-no2+k)
              enddo
          enddo
       enddo
       endif
       if ( jhigh_bound .and. klow_bound ) then
       do k=1,no2
          do j=1,no2
             do i =1,nl_i
                rho_g(i,nl_j+j,-no2+k) = rho_g(i,nl_j+j,nl_k-no2+k)
                vx_g(i,nl_j+j,-no2+k) = vx_g(i,nl_j+j,nl_k-no2+k)
                vy_g(i,nl_j+j,-no2+k) = vy_g(i,nl_j+j,nl_k-no2+k)
                vz_g(i,nl_j+j,-no2+k) = vz_g(i,nl_j+j,nl_k-no2+k)
              enddo
          enddo
       enddo
       endif
c
       if ( jlow_bound .and. khigh_bound ) then
       do k=1,no2
          do j=1,no2
             do i =1,nl_i
                rho_g(i,-no2+j,nl_k+k) = rho_g(i,-no2+j,k)
                vx_g(i,-no2+j,nl_k+k) = vx_g(i,-no2+j,k)
                vy_g(i,-no2+j,nl_k+k) = vy_g(i,-no2+j,k)
                vz_g(i,-no2+j,nl_k+k) = vz_g(i,-no2+j,k)
              enddo
          enddo
       enddo
       endif
       if ( jhigh_bound .and. khigh_bound ) then
       do k=1,no2
          do j=1,no2
             do i =1,nl_i
                rho_g(i,nl_j+j,nl_k+k) = rho_g(i,nl_j+j,k)
                vx_g(i,nl_j+j,nl_k+k) = vx_g(i,nl_j+j,k)
                vy_g(i,nl_j+j,nl_k+k) = vy_g(i,nl_j+j,k)
                vz_g(i,nl_j+j,nl_k+k) = vz_g(i,nl_j+j,k)
              enddo
          enddo
       enddo
       endif
c
c           ki edges
c
       if (klow_bound .and. ilow_bound ) then
       do k=1,no2
          do j=1,nl_j
             do i=1,no2
                rho_g(-no2+i,j,-no2+k) = rho_g(-no2+i,j,nl_k-no2+k)
                vx_g(-no2+i,j,-no2+k) = vx_g(-no2+i,j,nl_k-no2+k)
                vy_g(-no2+i,j,-no2+k) = vy_g(-no2+i,j,nl_k-no2+k)
                vz_g(-no2+i,j,-no2+k) = vz_g(-no2+i,j,nl_k-no2+k)
             enddo
          enddo
       enddo
       endif
       if (klow_bound .and. ihigh_bound ) then
       do k=1,no2
          do j=1,nl_j
             do i=1,no2
                rho_g(nl_i+i,j,-no2+k) = rho_g(nl_i+i,j,nl_k-no2+k)
                vx_g(nl_i+i,j,-no2+k) = vx_g(nl_i+i,j,nl_k-no2+k)
                vy_g(nl_i+i,j,-no2+k) = vy_g(nl_i+i,j,nl_k-no2+k)
                vz_g(nl_i+i,j,-no2+k) = vz_g(nl_i+i,j,nl_k-no2+k)
             enddo
          enddo
       enddo
       endif
c
       if (khigh_bound .and. ilow_bound ) then
       do k=1,no2
          do j=1,nl_j
             do i=1,no2
                rho_g(-no2+i,j,nl_k+k) = rho_g(-no2+i,j,k)
                vx_g(-no2+i,j,nl_k+k) = vx_g(-no2+i,j,k)
                vy_g(-no2+i,j,nl_k+k) = vy_g(-no2+i,j,k)
                vz_g(-no2+i,j,nl_k+k) = vz_g(-no2+i,j,k)
             enddo
          enddo
       enddo
       endif
       if (klow_bound .and. ihigh_bound ) then
       do k=1,no2
          do j=1,nl_j
             do i=1,no2
                rho_g(nl_i+i,j,nl_k+k) = rho_g(nl_i+i,j,k)
                vx_g(nl_i+i,j,nl_k+k) = vx_g(nl_i+i,j,k)
                vy_g(nl_i+i,j,nl_k+k) = vy_g(nl_i+i,j,k)
                vz_g(nl_i+i,j,nl_k+k) = vz_g(nl_i+i,j,k)
             enddo
          enddo
       enddo
       endif
!
! ... End subroutine edge_fix() ........................................
!
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!