!> fixB.F
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Result:  bi_g,bj_g,bk_g stores the fluxes after fields are subtracted out
      subroutine FixBout
! ... 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 "run-time.inc"
!
! ... Begin ............................................................

#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in fixB.F::FixBout()"
#endif
      if ( .not. MHD_PE ) return
c
c     Add an estimate of cartesian fields from subtracted fluxes at the
c     cell centers
      do k=1,nl_k
         do j=1,nl_j
            do i=1,nl_i
	  bx_g(i,j,k) = bx_g(i,j,k)+
     $		0.5*(bxqnk(i,j,k)+bxqnk(i,j,k+1))
	  by_g(i,j,k) = by_g(i,j,k)+
     $		0.5*(byqnk(i,j,k)+byqnk(i,j,k+1))
	  bz_g(i,j,k) = bz_g(i,j,k)+
     $		0.5*(bzqnk(i,j,k)+bzqnk(i,j,k+1))
          enddo
       enddo
      enddo
c     Add the subtracted fields back to the perturbation fields
        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)+bnqfi(i,j,k)
        enddo
        enddo
        enddo
c
        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)+bnqfj(i,j,k)
        enddo
        enddo
        enddo
c
        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)+bnqfk(i,j,k)
        enddo
        enddo
        enddo

        return
        end subroutine FixBout
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> ResetB()
!!
      subroutine ResetB
! ... Local 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 "run-time.inc"
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in fixB.F::ResetB()"
#endif
      if ( .not. MHD_PE ) return
c
      do k=1,nl_k
         do j=1,nl_j
            do i=1,nl_i
	  bx_g(i,j,k) = bx_g(i,j,k)-
     $		0.5*(bxqnk(i,j,k)+bxqnk(i,j,k+1))
	  by_g(i,j,k) = by_g(i,j,k)-
     $		0.5*(byqnk(i,j,k)+byqnk(i,j,k+1))
	  bz_g(i,j,k) = bz_g(i,j,k)-
     $		0.5*(bzqnk(i,j,k)+bzqnk(i,j,k+1))
          enddo
       enddo
      enddo
        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)-bnqfi(i,j,k)
        enddo
        enddo
        enddo
c
        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)-bnqfj(i,j,k)
        enddo
        enddo
        enddo
c
        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)-bnqfk(i,j,k)
        enddo
        enddo
        enddo

        return
        end subroutine ResetB
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!