!> test-arrays.F: FORTRAN helper routines to use for debugging purposes.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> test_arrays
!! Call this subroutine anywhere in the code to gain access to global
!! variables.  This is extremely useful within the TotalView debugger
!! for inspecting how variables change
!!
      subroutine test_arrays
! ... Local variables ..................................................
c      basic dimensioning for MHD
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
***#include "var.inc"
***#include "var1.inc"
#include "var2.inc"
***#include "meter.inc"
#include "twod.inc"
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::test-arrays(...)"
#endif

#ifdef DEBUGTESTARRAY
c    >jgl test
        write (9,1910) mype
 1910   format ('In test_arrays  ',i3)
#endif
      if ( .not. MHD_PE ) return
      aa = 2.0
      if ( aa .eq. 1. ) aa = 4.0
! ... End subroutine test_arrays .......................................
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> write_cg_to_file
!! Write the sound speed (cg) to file for debugging purposes.
!!
      subroutine write_cg_to_file(offset)
! ... Local variables ..................................................
c      basic dimensioning for MHD
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
***#include "var.inc"
***#include "var1.inc"
#include "var2.inc"
***#include "meter.inc"
#include "twod.inc"
!
      integer file_unit
      integer offset
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::write_cg_to_file(...)"
#endif
! Write sound speed c_g to file
      file_unit = 20 + mype + offset
      write(file_unit) rho_g
! write dimension info for c_g to file
      file_unit = 30 + mype + offset
      write(file_unit, *) ilow, ihigh
      write(file_unit, *) jlow, jhigh
      write(file_unit, *) klow, khigh

! ... End subroutine write_cg_to_file...................................
      return
      end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> fort_divcheck(...)
!!
!!
       subroutine fort_divcheck(dv,av,ibad)
! ... Local variables ..................................................
c      basic dimensioning for MHD
#include "global_dims.inc"
#include "param.inc"
***#include "var.inc"
***#include "var1.inc"
***#include "var2.inc"
***#include "meter.inc"
***#include "twod.inc"
! ... Parameter variables ..............................................
       dimension dv(ilow:ihigh,jlow:jhigh,klow:khigh)
       dimension av(ilow:ihigh,jlow:jhigh,klow:khigh)
       dimension ibad(ilow:ihigh,jlow:jhigh,klow:khigh)
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::fort_divcheck(...)"
#endif
#ifdef DEBUGTESTARRAY
c    >jgl test
        write (9,1910) mype
 1910   format ('In test_arrays  ',i3)
#endif
      if ( .not. MHD_PE ) return
      aa = 2.0
      if ( aa .eq. 1. ) aa = 4.0
! ... End subroutine fort_divcheck .....................................
      return
      end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> fdivcheck(...)
!!
!!
      subroutine fdivcheck(zero,flux,i2)
! ... Local variables ..................................................
#include "global_dims.inc"
#include "param.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "var2.inc"
#include "bzero.inc"
#include "run-time.inc"
      real div_0(ni,nj,nk)
      real, dimension(:,:,:), allocatable :: bii, bjj,bkk,div
      real, dimension(:,:,:), allocatable :: div_abs, div_ratio
      real, dimension(:,:,:), allocatable :: f_abs, flux_tot
      data first/.true./
      logical first
! ... Parameter variables ..............................................
      logical flux, zero, i2
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::fdivcheck(...)"
#endif
      if (mhd_pe) then

      allocate(bii(nl_i+1,nl_j+1,nl_k+1))
      allocate(bjj(nl_i+1,nl_j+1,nl_k+1))
      allocate(bkk(nl_i+1,nl_j+1,nl_k+1))
      allocate(div(nl_i,nl_j,nl_k))
      allocate(div_abs(nl_i,nl_j,nl_k))
      allocate(div_ratio(nl_i,nl_j,nl_k))
      allocate(f_abs(nl_i,nl_j,nl_k))
      allocate(flux_tot(nl_i,nl_j,nl_k))
c
      if (.not. flux) call ave2flux
c
      if (zero) then
         do k=1,nl_k
            do j=1,nl_j
               do i=1,nl_i+1
                  if (i2) then
                  bii(i,j,k) = bi2(i,j,k)+face_i(i,j,k)*bnqfi(i,j,k)
                  else
                  bii(i,j,k) = bi_g(i,j,k)+face_i(i,j,k)*bnqfi(i,j,k)
                  endif
               enddo
            enddo
         enddo
         do k=1,nl_k
            do j=1,nl_j+1
               do i=1,nl_i
                  if (i2) then
                  bjj(i,j,k) = bj2(i,j,k)+face_j(i,j,k)*bnqfj(i,j,k)
                  else
                  bjj(i,j,k) = bj_g(i,j,k)+face_j(i,j,k)*bnqfj(i,j,k)
                  endif
               enddo
            enddo
         enddo
         do k=1,nl_k+1
            do j=1,nl_j
               do i=1,nl_i
                  if (i2) then
                  bkk(i,j,k) = bk2(i,j,k)+face_k(i,j,k)*bnqfk(i,j,k)
                  else
                  bkk(i,j,k) = bk_g(i,j,k)+face_k(i,j,k)*bnqfk(i,j,k)
                  endif
               enddo
            enddo
         enddo
      else
         do k=1,nl_k
            do j=1,nl_j
               do i=1,nl_i+1
                  if ( i2 ) then
                  bii(i,j,k) = bi2(i,j,k)
                  else
                  bii(i,j,k) = bi_g(i,j,k)
                  endif
               enddo
            enddo
         enddo
         do k=1,nl_k
            do j=1,nl_j+1
               do i=1,nl_i
                  if ( i2 ) then
                  bjj(i,j,k) = bj2(i,j,k)
                  else
                  bjj(i,j,k) = bj_g(i,j,k)
                  endif
               enddo
            enddo
         enddo
         do k=1,nl_k+1
            do j=1,nl_j
               do i=1,nl_i
                  if (i2) then
                  bkk(i,j,k) = bk2(i,j,k)
                  else
                  bkk(i,j,k) = bk_g(i,j,k)
                  endif
               enddo
            enddo
         enddo
      endif
c
      divmax = 0.
      imax = 1
      jmax = 1
      kmax = 1
      do k=1,nl_k
         do j=1,nl_j
            do i=1,nl_i
               div_abs(i,j,k) = abs(
     $              bii(i+1,j,k)-bii(i,j,k)+bjj(i,j+1,k)-bjj(i,j,k)+
     $              bkk(i,j,k+1)-bkk(i,j,k))
               f_abs(i,j,k) = face_i(i,j,k)+face_i(i+1,j,k)+
     $              face_j(i,j,k)+
     $              face_j(i,j+1,k)+face_k(i,j,k)+face_k(i,j,k+1)
               if (first) then
                  div_0(i,j,k) = (div_abs(i,j,k)+1.e-10)
                  div_ratio(i,j,k) = 1.0
               else
                  div_ratio(i,j,k) = div_abs(i,j,k)/div_0(i,j,k)
               endif
               flux_tot(i,j,k) =
     $              abs(bii(i+1,j,k))+abs(bii(i,j,k))+
     $              abs(bjj(i,j+1,k))+abs(bjj(i,j,k))+
     $              abs(bkk(i,j,k+1))+abs(bkk(i,j,k))
               div(i,j,k) = div_abs(i,j,k)/(flux_tot(i,j,k)+
     $              1.e-6*face_k(i,j,k))
               if (div(i,j,k) .gt. divmax) then
                  divmax = div(i,j,k)
                  imax =i
                  jmax = j
                  kmax =k
*                  if ( mype .eq. 1 ) then
*                     write (6,*) '>>>>>>>>>>>>>>>>>>>>>>>>>>>  ',
*     $                    i,j,k
*                     write (6,911) div(i,j,k),bii(i:i+1,j,k),
*     $                    bjj(i,j:j+1,k),bkk(i,j,k:k+1)
 911                 format(1pe12.4/6(1pe13.4))
*                  endif
               endif
            enddo
         enddo
      enddo
c
      do lush=1,num_mhd
         iunit = 70+lush
         if (lush .eq. mype) then
      write (6,901) mype, imax,jmax,kmax, divmax,
     $     div_ratio(imax,jmax,kmax),div_abs(imax,jmax,kmax),
     $     div_abs(imax,jmax,kmax)/f_abs(imax,jmax,kmax),
     $     flux_tot(imax,jmax,kmax)/f_abs(imax,jmax,kmax),
     $     imax+i_global_off,jmax+j_global_off,kmax+k_global_off
      write (iunit,902) lstep, mype, imax,jmax,kmax, divmax,
     $     div_ratio(imax,jmax,kmax),div_abs(imax,jmax,kmax),
     $     div_abs(imax,jmax,kmax)/f_abs(imax,jmax,kmax),
     $     flux_tot(imax,jmax,kmax)/f_abs(imax,jmax,kmax),
     $     imax+i_global_off,jmax+j_global_off,kmax+k_global_off
 901  format(4i6,3x,5(1pe14.6),2x,3i6)
 902  format(5i6,3x,5(1pe14.6),2x,3i6)
      call flush(iunit)
          endif
      enddo
      if (.not. flux) call flux2ave
      deallocate(bii)
      deallocate(bjj)
      deallocate(bkk)
      deallocate(div)
      deallocate(div_abs)
      deallocate(div_ratio)
      deallocate(f_abs)
      deallocate(flux_tot)
      first = .false.
c
      endif ! end of if around mhd pe code
! ... End subroutine fdivcheck(zero,flux,i2) ............................
      return
      end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> add_bzero
!!
!!
      subroutine add_bzero
! ... Local variables ..................................................
#include "global_dims.inc"
#include "param.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "var2.inc"
#include "bzero.inc"
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::add_bzero()"
#endif
      if (mhd_pe) then

         do k=1,nl_k
            do j=1,nl_j
               do i=1,nl_i+1
                  bi_g(i,j,k) = bi_g(i,j,k)+face_i(i,j,k)*bnqfi(i,j,k)
               enddo
            enddo
         enddo
         do k=1,nl_k
            do j=1,nl_j+1
               do i=1,nl_i
                  bj_g(i,j,k) = bj_g(i,j,k)+face_j(i,j,k)*bnqfj(i,j,k)
               enddo
            enddo
         enddo
         do k=1,nl_k+1
            do j=1,nl_j
               do i=1,nl_i
                  bk_g(i,j,k) = bk_g(i,j,k)+face_k(i,j,k)*bnqfk(i,j,k)
               enddo
            enddo
         enddo
         endif
! ... End subroutine add_bzero ..........................................
         return
         end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> subtract_bzero
!!
!!

      subroutine subtract_bzero
! ... Local variables ..................................................
#include "global_dims.inc"
#include "param.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "var2.inc"
#include "bzero.inc"
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::subtract_bzero()"
#endif
      if (mhd_pe) then

         do k=1,nl_k
            do j=1,nl_j
               do i=1,nl_i+1
                  bi_g(i,j,k) = bi_g(i,j,k)-face_i(i,j,k)*bnqfi(i,j,k)
               enddo
            enddo
         enddo
         do k=1,nl_k
            do j=1,nl_j+1
               do i=1,nl_i
                  bj_g(i,j,k) = bj_g(i,j,k)-face_j(i,j,k)*bnqfj(i,j,k)
               enddo
            enddo
         enddo
         do k=1,nl_k+1
            do j=1,nl_j
               do i=1,nl_i
                  bk_g(i,j,k) = bk_g(i,j,k)-face_k(i,j,k)*bnqfk(i,j,k)
               enddo
            enddo
         enddo
         endif
! ... End subroutine subtract_bzero .....................................
         return
         end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> ion_test_arrays
!!
!!

#if defined(ION_ON) && !defined(USE_MIX)      
      subroutine ion_test_arrays
! ... Local variables .................................................. 
c      basic dimensioning for MHD
#include "global_dims.inc"
#include "ion-param.inc"
#include "ion_P++_pointers.inc"
#include "ion_fortran_P++_arrays.inc"
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::ion_test_arrays()"
#endif
#ifdef DEBUGTESTARRAY
c    >jgl test
        write (9,1910) mype
 1910   format ('In test_arrays  ',i3)
#endif
      if ( .not. ION_PE ) return
      aa = 2.0
      if ( aa .eq. 1. ) aa = 4.0

! ... End subroutine ion_test_arrays ....................................

      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#endif

!> test_ppp(XOX)
!!
!!
      subroutine test_ppp(XOX)
! ... Local variables ..................................................
c      basic dimensioning for MHD
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
***#include "var.inc"
***#include "var1.inc"
#include "var2.inc"
***#include "meter.inc"
#include "twod.inc"
      dimension XOX(ilow:ihigh,jlow:jhigh,klow:khigh)
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::test_ppp(...)"
#endif
      if ( .not. MHD_PE ) return
      aa = 2.0
      if ( XOX(1,1,1) .eq. 1. ) aa = 4.0
c
      write (6,*) '>>>>>>>>> rho_g', rho_g(1,1,1),rho_g(nl_i,nl_j,nl_k)
! ... End subroutine test_ppp(XOX) .....................................
      return
      end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> load_b
!!
!!

      subroutine load_b
! ... Local variables ..................................................
c      basic dimensioning for MHD
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
***#include "var.inc"
***#include "var1.inc"
#include "var2.inc"
***#include "meter.inc"
#include "twod.inc"
      common /dbzz/ dbi(nip1,nj,nk),dbj(ni,njp1,nk),
     $     dbk(ni,nj,nkp1)
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::load_b()"
#endif
#ifdef DEBUG
c    >jgl test
        write (9,1910) mype
 1910   format ('In test_arrays: load_b  ',i3)
#endif
      if ( .not. MHD_PE ) return
      do k=1,nl_k
         do j=1,nl_j
            do i=1,nl_i+1
               dbi(i,j,k) = bi_g(i,j,k)
            enddo
         enddo
      enddo
      do k=1,nl_k
         do j=1,nl_j+1
            do i=1,nl_i
               dbj(i,j,k) = bj_g(i,j,k)
            enddo
         enddo
      enddo
      do k=1,nl_k+1
         do j=1,nl_j
            do i=1,nl_i
               dbk(i,j,k) = bk_g(i,j,k)
            enddo
         enddo
      enddo
! ... End subroutine load_b .............................................
      return
      end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> exam_b
!!
!!
      subroutine exam_b
! ... Local variables ..................................................
c      basic dimensioning for MHD
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
***#include "var.inc"
***#include "var1.inc"
#include "var2.inc"
***#include "meter.inc"
#include "twod.inc"
      common /dbzz/ dbi(nip1,nj,nk),dbj(ni,njp1,nk),
     $     dbk(ni,nj,nkp1)
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in test-arrays.F::exam_b()"
#endif
#ifdef DEBUGTESTARRAY
c    >jgl test
        write (9,1910) mype
 1910   format ('In test_arrays  ',i3)
#endif
      if ( .not. MHD_PE ) return
      do k=1,nl_k
         do j=1,nl_j
            do i=1,nl_i+1
               dbi(i,j,k) = bi_g(i,j,k) - dbi(i,j,k)
            enddo
         enddo
      enddo
      do k=1,nl_k
         do j=1,nl_j+1
            do i=1,nl_i
               dbj(i,j,k) = bj_g(i,j,k) -dbj(i,j,k)
            enddo
         enddo
      enddo
      do k=1,nl_k+1
         do j=1,nl_j
            do i=1,nl_i
               dbk(i,j,k) = bk_g(i,j,k) - dbk(i,j,k)
            enddo
         enddo
      enddo
! ... End subroutine exam_b .............................................
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!