!> hack.F
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> MHD_HACK
!!
!!     This code is meant to become obsolete as rapidly as possible. It
!!     is designed to make the vintage global MHD code run in a parallel
!!     fashion under MPI as quickly as possible -- with few concessions
!!     to memory use or efficiency.
!!
!!     The plan is to use the A++P++ classes to automatically partition
!!     the MHD data and handle communications. The original single
!!     processor code then needs to be fooled into thinking that the
!!     local data arrays are the whole calculation.
!!
!!     There are two major issues that arise. The first is the
!!     synchronization of the P++ local data structure with the
!!     fortran. The second is the handling of local boundary conditions.
!!
!!   1 P++ allows one to define the number of internal guard cells that
!!     local copies have. Suppose the number of guard cells is defined
!!     to be ng, and the number of "active" cells in the 3 directions
!!     are mi, mj, mk. Then the dimension of the local P++ data array is
!!     (mi+2ng,mj+2ng,mk+2ng). The start of the "active" region is at
!!     (ng+1,ng+1,ng+1) [I'm using Fortran index convention]. This is
!!     true for regions without internal boundaries. Even if the global
!!     P++ array is dimensioned (1:NI,1:NJ,1:NK), it acts as though it
!!     is dimensioned (-ng+1:NI+ng,-ng+1:NJ+ng,-ng+1:NK+ng). If a
!!     pointer to the local array is handed to fortran, it suffices to
!!     dimension the array as
!!     (-ng+1:mi+ng,-ng+1:mj+ng,-ng+1:mk+ng). Then (1,1,1) will point to
!!     the origin of the local "active" cells.
!!
!!     Unfortunately, this layout is not consistent with the layout in
!!     the vintage code. The quickest hack around this is to make a copy
!!     of the relevant regions of the P++ array to conform to the memory
!!     layout that the old code expects. (I said this was not going to
!!     be efficient). In addition, it would be efficient to use pointers
!!     ( and dynami!!dimensioning) for the local fortran arrays, but
!!     it's a quicker hack to use stati!!dimensioning for the fortran
!!     via parameter statements. These parameters must be sufficiently
!!     large that the distributed arrays fit inside. Therefore, a
!!     certain amount of information abou the partitioning will need to
!!     be obtained a priori.
!!
!!   2 Boundary conditions are handled in the vintage code by a gaggle of
!!     individual subroutine calls. The P++ arrays automatically contain
!!     the interior boundary information. One way to handle this is make
!!     all the boundary routines in BRD and BRD2 use buffer arrays ( as
!!     they currently do for the solar wind data). Oncce again this is
!!     an inefficient hack as we will need to do more local copies of
!!     the P++ data, but it is probably the simplestway to get a running
!!     code.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       subroutine MHD_HACK
! ... Variables ........................................................
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
#include "help.inc"
#include "bzero.inc"
***#include "var.inc"
***#include "var1.inc"
#include "var2.inc"
***#include "meter.inc"
#include "twod.inc"
#include "run-time.inc"
C
#include "b_dbl.inc"
C
#ifdef DIFFBUG
	dimension drho_2(ni,nj,nk),dvx2(ni,nj,nk),dvy2(ni,nj,nk),dvz2(ni,nj,nk)
	dimension dbvx2(ni,nj,nk),dbvy2(ni,nj,nk),dbvz2(ni,nj,nk)
        dimension dbi2(nip1,nj,nk),dbj2(ni,njp1,nk),dbk2(ni,nj,nkp1)
#endif
        dimension  DUMMY(NIP1,NK,12)
c
        DIMENSION DSUM(NI,NK),RHOSUM(NI,NK),VXSUM(NI,NK),VYSUM(NI,NK),
     $          VZSUM(NI,NK),CSUM(NI,NK),fisum(nip1,nk),bisum(nip1,nk),
     $       fjsum(ni,nk),bjsum(ni,nk),bdvxsum(ni,nk),bdvysum(ni,nk),
     $       bdvzsum(ni,nk)
        dimension bxa(li),bya(li),bza(li),bdot(li),rhob(li)
        dimension vxa(li),vya(li),vza(li)
        equivalence (bxa,dsum),(bya,rhosum),(bza,vxsum),(bdot,vysum),
     $          (rhob,vzsum),(vxa,csum),(vya,fisum),(vza,bisum)
        EQUIVALENCE (DUMMY(1,1,1),DSUM),(DUMMY(1,1,2),RHOSUM),
     $          (DUMMY(1,1,3),VXSUM),(DUMMY(1,1,4),VZSUM),
     $          (DUMMY(1,1,5),CSUM),(dummy(1,1,6),fisum),
     $          (dummy(1,1,7),bisum),(dummy(1,1,8),fjsum),
     $       (dummy(1,1,9),bjsum),(dummy(1,1,10),bdvxsum),
     $       (dummy(1,1,11),bdvysum),(dummy(1,1,12),bdvzsum)
        dimension df(nip1,nk),bjs(ni,nk),bjc(ni,nk),cjsum(ni,nk),
     $       sjsum(ni,nk)
*        dimension diverg(nj,nk2,5)
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in hack.F::MHD_HACK()"
#endif
#ifdef DEBUG
        write (9,*) 'in mhd_hack ',mype, MHD_PE
#endif
        if ( .not. MHD_PE ) return
        nl_ip1 = nl_i+1
        nl_jp1 = nl_j+1
        nl_kp1 = nl_k+1
c
c
c       get extrapolated mid-step values
c
        if ( dt_old .lt. 0. ) then
        tfactor = 0.
        base_factor = 1.0
        else
        tfactor = 0.5*(dt/dt_old)
        base_factor = 1. + tfactor
        endif
c
        call second(time0)
c
        do k =klow,khigh
           do j=jlow,jhigh
              do i=ilow,ihigh
                 rho1_g(i,j,k) = base_factor*rho_g(i,j,k) -
     $                               tfactor*rho1_g(i,j,k)
                 c1_g(i,j,k) = base_factor*c_g(i,j,k) -
     $                               tfactor*c1_g(i,j,k)
                 vx1_g(i,j,k) = base_factor*vx_g(i,j,k) -
     $                               tfactor*vx1_g(i,j,k)
                 vy1_g(i,j,k) = base_factor*vy_g(i,j,k) -
     $                               tfactor*vy1_g(i,j,k)
                 vz1_g(i,j,k) = base_factor*vz_g(i,j,k) -
     $                               tfactor*vz1_g(i,j,k)
                 bx1_g(i,j,k) = base_factor*bx_g(i,j,k) -
     $                               tfactor*bx1_g(i,j,k)
                 by1_g(i,j,k) = base_factor*by_g(i,j,k) -
     $                               tfactor*by1_g(i,j,k)
                 bz1_g(i,j,k) = base_factor*bz_g(i,j,k) -
     $                               tfactor*bz1_g(i,j,k)
                 bi1_g(i,j,k) = base_factor*bi_g(i,j,k) -
     $                               tfactor*bi1_g(i,j,k)
                 bj1_g(i,j,k) = base_factor*bj_g(i,j,k) -
     $                               tfactor*bj1_g(i,j,k)
                 bk1_g(i,j,k) = base_factor*bk_g(i,j,k) -
     $                               tfactor*bk1_g(i,j,k)
              enddo
           enddo
        enddo
c
        call second(time1)
C
C       SET UP THE ACCUMULATOR ARRAYS
C
C
        DO 200 K=1,NL_K
        DO 200 J=1,NL_J
        DO 200 I=1,NL_I
        RHO2(I,J,K) = RHO_g(I,J,K)*VOLUME(I,J,K)
        C2(I,J,K) = RHO2(I,J,K)*(0.5*(VX_g(I,J,K)**2+VY_g(I,J,K)**2
     $          +VZ_g(I,J,K)**2)+ GGM1IN*C_g(I,J,K)**2)
C  these accumulators are set to zero to ease implementation of the
C       alfven correction
C
        VX2(I,J,K) = 0.
        VY2(I,J,K) = 0.
        VZ2(I,J,K) = 0.
        BDVX2(I,J,K) = 0.
        BDVY2(I,J,K) = 0.
        BDVZ2(I,J,K) = 0.
  200   CONTINUE
        do 210 K=1,NL_KP1
        DO 210 J=1,NL_J
        DO 210 I=1,NL_I
  210   BK2(I,J,K) = BK_g(I,J,K)*FACE_k(I,J,K)
        DO 211 K=1,NL_K
        DO 211 J=1,NL_JP1
        DO 211 I=1,NL_I
  211   BJ2(I,J,K) = BJ_g(I,J,K)*FACE_J(I,J,K)
        DO 212 K=1,NL_K
        DO 212 J=1,NL_J
        DO 212 I=1,NL_IP1
  212   BI2(I,J,K) = BI_g(I,J,K)*FACE_I(I,J,K)
C
C
#ifdef DIFFBUG
	do k=1,nk
	do j=1,nj
	do i=1,ni
        drho_2(i,j,k) = rho2(i,j,k)
        dvx2(i,j,k) =vx2(i,j,k)
        dvy2(i,j,k) =vy2(i,j,k)
        dvz2(i,j,k) =vz2(i,j,k)
        dbvx2(i,j,k) =bdvx2(i,j,k)
        dbvy2(i,j,k) =bdvy2(i,j,k)
        dbvz2(i,j,k) =bdvz2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj
	do i=1,ni+1
        dbi2(i,j,k) = bi2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj+1
	do i=1,ni
        dbj2(i,j,k) = bj2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk+1
	do j=1,nj
	do i=1,ni
        dbk2(i,j,k) = bk2(i,j,k)
	enddo
	enddo
	enddo
#endif
#ifdef DEBUG
        write (9,*) 'just above ijmhd  ', mype, 'so what?'
#endif
C
#ifdef PSCTEST
*
*   test code
               myunit = mype+40
               k=25
               level3=25
            if (level3 .eq. 25) then
               write (myunit,*) '>>>>>>>>>>>>>>>>>  level3 =',level3,k
        do j=llow,nl_j+no2
           do i=1,nl_ip1
              write (myunit,978) i,j,k,tran_i(i,j,k,1),
     $             tran_i(i,j,k,2),
     $             tran_i(i,j,k,3)
              call flush(myunit)
 978          format(' t ',3i6,3(1pe13.3))
           enddo
        enddo
            endif
*
*
#endif
C
        CALL IJMHD(DT,RHO_g,VX_g,VY_g,VZ_g,C_g,BX_g,BY_g,BZ_g,
     $       BI_g,BJ_g,BK_g,
     $          RHO1_g,VX1_g,VY1_g,VZ1_g,C1_g,BX1_g,BY1_g,BZ1_g,
     $       BI1_g,BJ1_g,BK1_g,
     $       ilow,ihigh,jlow,jhigh,klow,khigh)
C
*       write(9,*) 'lstep',lstep
*       write(9,*) 'ijmhd'
*       write (9,*) 'vz2',(vz2(1,8,k),k=1,16)
*       write(9,*) 'bdvz2',(bdvz2(1,8,k),k=1,16)
*       write(9,*) 'bi2',(bi2(2,8,k),k=1,16)
*       write(9,*) 'bj2',(bj2(1,8,k),k=1,16)
C
C
#ifdef DIFFBUG
	do k=1,nk
	do j=1,nj
	do i=1,ni
	drho_2(i,j,k) = rho2(i,j,k) -drho_2(i,j,k) 
	dvx2(i,j,k) = vx2(i,j,k) -dvx2(i,j,k) 
	dvy2(i,j,k) = vy2(i,j,k) -dvy2(i,j,k) 
	dvz2(i,j,k) = vz2(i,j,k) -dvz2(i,j,k) 
	dbvx2(i,j,k) = bdvx2(i,j,k) -dbvx2(i,j,k) 
	dbvy2(i,j,k) = bdvy2(i,j,k) -dbvy2(i,j,k) 
	dbvz2(i,j,k) = bdvz2(i,j,k) -dbvz2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj
	do i=1,ni+1
        dbi2(i,j,k) = bi2(i,j,k)-dbi2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj+1
	do i=1,ni
        dbj2(i,j,k) = bj2(i,j,k)-dbj2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk+1
	do j=1,nj
	do i=1,ni
        dbk2(i,j,k) = bk2(i,j,k)-dbk2(i,j,k)
	enddo
	enddo
	enddo
c
c  put them back for the next round
c
	do k=1,nk
	do j=1,nj
	do i=1,ni
	drho_2(i,j,k) = rho2(i,j,k)
	dvx2(i,j,k) = vx2(i,j,k)
	dvy2(i,j,k) = vy2(i,j,k)
	dvz2(i,j,k) = vz2(i,j,k)
	dbvx2(i,j,k) = bdvx2(i,j,k)
	dbvy2(i,j,k) = bdvy2(i,j,k)
	dbvz2(i,j,k) = bdvz2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj
	do i=1,ni+1
        dbi2(i,j,k) = bi2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj+1
	do i=1,ni
        dbj2(i,j,k) = bj2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk+1
	do j=1,nj
	do i=1,ni
        dbk2(i,j,k) = bk2(i,j,k)
	enddo
	enddo
	enddo
#endif	
#ifdef DEBUG
        write (9,*) 'just above jkmhd  ', mype, 'so what?'
#endif
C
        CALL JKMHD(DT,RHO_g,VX_g,VY_g,VZ_g,C_g,BX_g,BY_g,BZ_g,
     $       BI_g,BJ_g,BK_g,
     $          RHO1_g,VX1_g,VY1_g,VZ1_g,C1_g,BX1_g,BY1_g,BZ1_g,
     $       BI1_g,BJ1_g,BK1_g,
     $       ilow,ihigh,jlow,jhigh,klow,khigh)
C
#ifdef DIFFBUG
	do k=1,nk
	do j=1,nj
	do i=1,ni
	drho_2(i,j,k) = rho2(i,j,k) -drho_2(i,j,k) 
	dvx2(i,j,k) = vx2(i,j,k) -dvx2(i,j,k) 
	dvy2(i,j,k) = vy2(i,j,k) -dvy2(i,j,k) 
	dvz2(i,j,k) = vz2(i,j,k) -dvz2(i,j,k) 
	dbvx2(i,j,k) = bdvx2(i,j,k) -dbvx2(i,j,k) 
	dbvy2(i,j,k) = bdvy2(i,j,k) -dbvy2(i,j,k) 
	dbvz2(i,j,k) = bdvz2(i,j,k) -dbvz2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj
	do i=1,ni+1
        dbi2(i,j,k) = bi2(i,j,k)-dbi2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj+1
	do i=1,ni
        dbj2(i,j,k) = bj2(i,j,k)-dbj2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk+1
	do j=1,nj
	do i=1,ni
        dbk2(i,j,k) = bk2(i,j,k)-dbk2(i,j,k)
	enddo
	enddo
	enddo
c
c  put them back for the next round
c
	do k=1,nk
	do j=1,nj
	do i=1,ni
	drho_2(i,j,k) = rho2(i,j,k)
	dvx2(i,j,k) = vx2(i,j,k)
	dvy2(i,j,k) = vy2(i,j,k)
	dvz2(i,j,k) = vz2(i,j,k)
	dbvx2(i,j,k) = bdvx2(i,j,k)
	dbvy2(i,j,k) = bdvy2(i,j,k)
	dbvz2(i,j,k) = bdvz2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj
	do i=1,ni+1
        dbi2(i,j,k) = bi2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj+1
	do i=1,ni
        dbj2(i,j,k) = bj2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk+1
	do j=1,nj
	do i=1,ni
        dbk2(i,j,k) = bk2(i,j,k)
	enddo
	enddo
	enddo
#endif	
C
#ifdef DEBUG
        write (9,*) 'just above kimhd  ', mype, 'so what?'
#endif
C
C
        CALL KIMHD(DT,RHO_g,VX_g,VY_g,VZ_g,C_g,BX_g,BY_g,BZ_g,
     $       BI_g,BJ_g,BK_g,
     $          RHO1_g,VX1_g,VY1_g,VZ1_g,C1_g,BX1_g,BY1_g,BZ1_g,
     $       BI1_g,BJ1_g,BK1_g,
     $       ilow,ihigh,jlow,jhigh,klow,khigh)
C
C
#ifdef DIFFBUG
	do k=1,nk
	do j=1,nj
	do i=1,ni
	drho_2(i,j,k) = rho2(i,j,k) -drho_2(i,j,k) 
	dvx2(i,j,k) = vx2(i,j,k) -dvx2(i,j,k) 
	dvy2(i,j,k) = vy2(i,j,k) -dvy2(i,j,k) 
	dvz2(i,j,k) = vz2(i,j,k) -dvz2(i,j,k) 
	dbvx2(i,j,k) = bdvx2(i,j,k) -dbvx2(i,j,k) 
	dbvy2(i,j,k) = bdvy2(i,j,k) -dbvy2(i,j,k) 
	dbvz2(i,j,k) = bdvz2(i,j,k) -dbvz2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj
	do i=1,ni+1
        dbi2(i,j,k) = bi2(i,j,k)-dbi2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj+1
	do i=1,ni
        dbj2(i,j,k) = bj2(i,j,k)-dbj2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk+1
	do j=1,nj
	do i=1,ni
        dbk2(i,j,k) = bk2(i,j,k)-dbk2(i,j,k)

	enddo
	enddo
	enddo
c
c  put them back for the next round
c
	do k=1,nk
	do j=1,nj
	do i=1,ni
	drho_2(i,j,k) = rho2(i,j,k)
	dvx2(i,j,k) = vx2(i,j,k)
	dvy2(i,j,k) = vy2(i,j,k)
	dvz2(i,j,k) = vz2(i,j,k)
	dbvx2(i,j,k) = bdvx2(i,j,k)
	dbvy2(i,j,k) = bdvy2(i,j,k)
	dbvz2(i,j,k) = bdvz2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj
	do i=1,ni+1
        dbi2(i,j,k) = bi2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk
	do j=1,nj+1
	do i=1,ni
        dbj2(i,j,k) = bj2(i,j,k)
	enddo
	enddo
	enddo
	do k=1,nk+1
	do j=1,nj
	do i=1,ni
        dbk2(i,j,k) = bk2(i,j,k)
	enddo
	enddo
	enddo
#endif	
C
C
C
C     get volume integral of internal energy in c2
C
#ifdef DEBUG
        write (9,*) 'after kimhd  ', mype, 'so what?'
#endif
        DO 700 K=1,NL_K
        DO 700 J=1,NL_J
        DO 700 I=1,NL_I
        DUMMY(i,j,1) = 1./RHO2(I,J,K)
        dummy(i,j,2) = rho_g(i,j,k)*volume(I,j,k)
        VXa(i) = (VX2(I,J,K) + vx_g(i,j,k)*dummy(i,j,2))*DUMMY(i,j,1)
        VYa(i) = (VY2(I,J,K) + vy_g(i,j,k)*dummy(i,j,2))*DUMMY(i,j,1)
        VZa(i) = (VZ2(I,J,K) + vz_g(i,j,k)*dummy(i,j,2))*DUMMY(i,j,1)
        C2(I,J,K) = rho2(i,j,k)*amax1(1.,GGM1*( C2(I,J,K)*DUMMY(i,j,1)
     $          -0.5*(VXa(i)**2+VYa(i)**2+VZa(i)**2) ))
  700   CONTINUE
C
C
C    CONVERT THE CONSERVED QUANTITIES BACK TO PRIMITIVE VARIABLES
C
C
C
        DO 924 K=1,NL_K
        DO 924 J=1,NL_J
        DO 924 I=1,NL_I
        c2(i,j,k) = sqrt( amax1(1.,c2(i,j,k)/rho2(i,j,k)) )
        RHO2(I,J,K) = RHO2(I,J,K)/VOLUME(I,J,K)
  924   CONTINUE
C
        DO 930 K=1,NL_K
        DO 930 J=1,NL_J
        DO 930 I=1,NL_IP1
  930   BI2(I,J,K) = BI2(I,J,K)/FACE_I(I,J,K)
C
        DO 931 K=1,NL_K
        DO 931 J=1,NL_JP1
        DO 931 I=1,NL_I
  931   BJ2(I,J,K) = BJ2(I,J,K)/FACE_J(I,J,K)
C
        DO 932 K=1,NL_KP1
        DO 932 J=1,NL_J
        DO 932 I=1,NL_I
  932   BK2(I,J,K) = BK2(I,J,K)/FACE_K(I,J,K) 
C
C
        do 938 i=1,nl_i
        do 938 k=1,nl_k
        do 938 j=1,nl_j
        rhob(j) = volume(i,j,k)*(rho2(i,j,k)-rho_g(i,j,k))
        vx2(i,j,k) = vx2(i,j,k) - rhob(j)*vx_g(i,j,k)
        vy2(i,j,k) = vy2(i,j,k) - rhob(j)*vy_g(i,j,k)
        vz2(i,j,k) = vz2(i,j,k) - rhob(j)*vz_g(i,j,k)
  938   continue
C
        do 940 i=1,nl_i
        do 940 k=1,nl_k
        do 940 j=1,nl_j
        bxa(j) = bx1_g(i,j,k)+0.5*(bxqfk(i,j,k)+bxqfk(i,j,k+1))
        bya(j) = by1_g(i,j,k)+0.5*(byqfk(i,j,k)+byqfk(i,j,k+1))
        bza(j) = bz1_g(i,j,k)+0.5*(bzqfk(i,j,k)+bzqfk(i,j,k+1))
        bdot(j) = alcon*(bxa(j)*vx2(i,j,k)+bya(j)*vy2(i,j,k)+
     $                bza(j)*vz2(i,j,k))/rho2(i,j,k)
        rhob(j) = volume(i,j,k)*(rho2(i,j,k)+
     $    alcon*(bxa(j)**2+bya(j)**2+bza(j)**2) )
        vx2(i,j,k) = vx_g(i,j,k) +
     $           (vx2(i,j,k)+bdvx2(i,j,k)+bxa(j)*bdot(j))/rhob(j)
        vy2(i,j,k) = vy_g(i,j,k) +
     $           (vy2(i,j,k)+bdvy2(i,j,k)+bya(j)*bdot(j))/rhob(j)
        vz2(i,j,k) = vz_g(i,j,k) +
     $           (vz2(i,j,k)+bdvz2(i,j,k)+bza(j)*bdot(j))/rhob(j)
  940   continue
C 
#ifdef DEBUG
        write (9,*) 'just above kluge  ', mype, 'so what?'
#endif
c     
c     Kluge to fix runaway speeds and temperatures
c 
c 
c     reshuffle the arrays -- note we copy the guard cells from rho to rho1
c     and so on -- cuts down on communication
c
#ifdef DEBUG
        write (9,*) 'just above reshuffle  ', mype, 'so what?'
#endif
         call second(time2)
         do k=klow,khigh
            do j=jlow,jhigh
               do i=ilow,ihigh
                  rho1_g(i,j,k) = rho_g(i,j,k)
                  vx1_g(i,j,k) = vx_g(i,j,k)
                  vy1_g(i,j,k) = vy_g(i,j,k)
                  vz1_g(i,j,k) = vz_g(i,j,k)
                  c1_g(i,j,k)  = c_g(i,j,k)
                  bx1_g(i,j,k) = bx_g(i,j,k)
                  by1_g(i,j,k) = by_g(i,j,k)
                  bz1_g(i,j,k) = bz_g(i,j,k)
                  bi1_g(i,j,k) = bi_g(i,j,k)
                  bj1_g(i,j,k) = bj_g(i,j,k)
                  bk1_g(i,j,k) = bk_g(i,j,k)
               enddo
            enddo
        enddo
c
c
        do k=1,nl_k
           do j=1,nl_j
              do i=1,nl_i
                 rho_g(i,j,k) = rho2(i,j,k)
                 vx_g(i,j,k)  = vx2(i,j,k)
                 vy_g(i,j,k)  = vy2(i,j,k)
                 vz_g(i,j,k)  = vz2(i,j,k)
                 c_g(i,j,k)   = c2(i,j,k)
              enddo
           enddo
         enddo
        do k=1,nl_k
           do j=1,nl_j
              do i=1,nl_i+1
                 bi_g(i,j,k) = bi2(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) = bj2(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) = bk2(i,j,k)
              enddo
           enddo
         enddo
c 
          call second(time4)
#ifdef DEBUG
         if (mype .eq. 1 .or. mype .eq. 4) write (6,*)
     $        "timings in hack for PE: ",mype, time0, time1,
     $        time2,time4
#endif 
      dt_old = dt
#ifdef DEBUG
      write (6,*) ">>>>>>>>>>>>>>>>>>>>> HACK ",mype,lstep,dt_old
#endif
!
! ... End subroutine MHD_HACK() ........................................
!
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> flux_fix(DT)
!!
      subroutine flux_fix(DT)
! ... Local Variables ..................................................
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
#include "help.inc"
#include "var.inc"
#ifdef BZEROS
#include "bzero.inc"
#endif
#include "b_dbl.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in hack.F::flux_fix(DT)"
#endif
      nl_ip1 = nl_i+1
      nl_jp1 = nl_j+1
      nl_kp1 = nl_k+1
c
      do k=1,nl_k
         do j=1,nl_j
            do i=1,nl_ip1
               bi_dbl(i,j,k) = bi_dbl(i,j,k) + DBLE(dt)*(
     $               DBLE(Ek_g(i,j+1,k)) - DBLE(Ek_g(i,j,k))
     $              -DBLE(Ej_g(i,j,k+1))  + DBLE(Ej_g(i,j,k)) )
            enddo
         enddo
      enddo
            
c
      do k=1,nl_k
         do j=1,nl_jp1
            do i=1,nl_i
               bj_dbl(i,j,k) = bj_dbl(i,j,k) + DBLE(dt)*(
     $               DBLE(Ei_g(i,j,k+1)) - DBLE(Ei_g(i,j,k))
     $              -DBLE(Ek_g(i+1,j,k))  + DBLE(Ek_g(i,j,k)) )
            enddo
         enddo
      enddo
            
c
      do k=1,nl_kp1
         do j=1,nl_j
            do i=1,nl_i
               bk_dbl(i,j,k) = bk_dbl(i,j,k) + DBLE(dt)*(
     $               DBLE(Ej_g(i+1,j,k)) - DBLE(Ej_g(i,j,k))
     $              -DBLE(Ei_g(i,j+1,k))  + DBLE(Ei_g(i,j,k)) )
            enddo
         enddo
      enddo
c
c
      call ave2flux
c
      do k=1,nl_kp1
         do j = 1,nl_j
            do i=1,nl_i
               bk_g(i,j,k) = bk_dbl(i,j,k)
            enddo
         enddo
      enddo
      do k=1,nl_k
         do j = 1,nl_jp1
            do i=1,nl_i
               bj_g(i,j,k) = bj_dbl(i,j,k)
            enddo
         enddo
      enddo
      do k=1,nl_k
         do j = 1,nl_j
            do i=1,nl_ip1
               bi_g(i,j,k) = bi_dbl(i,j,k)
            enddo
         enddo
      enddo
c
      call flux2ave
c
#ifdef BZEROS
      do k=1,nl_kp1
         do j = 1,nl_j
            do i=1,nl_i
               bk_g(i,j,k) = bk_g(i,j,k)-bnqfk(i,j,k)
            enddo
         enddo
      enddo
      do k=1,nl_k
         do j = 1,nl_jp1
            do i=1,nl_i
               bj_g(i,j,k) = bj_g(i,j,k)-bnqfj(i,j,k)
            enddo
         enddo
      enddo
      do k=1,nl_k
         do j = 1,nl_j
            do i=1,nl_ip1
               bi_g(i,j,k) = bi_g(i,j,k)-bnqfi(i,j,k)
            enddo
         enddo
      enddo
#endif
!
! ... End subroutine FLUX_FIX(DT) ......................................
!
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!