!> boundaries.F
!!
!!  boundaries are applied in the order of I, then J, then K in order
!!  to make all the boundaries come out right. I is the only one which
!!  can't be set by symmetry.
!!
!!  one subtle point about the P++ updateGhostBoundaries() is that it
!!  does not update the exterior ghost points. Therefore, we have to set
!!  those here. It also means that update of the interior ghost cells
!!  should be done before these routines are called.
!!
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> KBOUND_FOR
!!
      Subroutine KBOUND_FOR
!
! ... 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 boundaries.F::KBOUND_FOR()"
#endif
c
      if ( .not. mhd_pe ) return
c
c  all processors are k boundaries
c
      do k=1,no2
         do j=jlow,jhigh
            do i=ilow,ihigh
               rho_g(i,j,-no2+k) = rho_g(i,j,nl_k-no2+k)
               vx_g(i,j,-no2+k) = vx_g(i,j,nl_k-no2+k)
               vy_g(i,j,-no2+k) = vy_g(i,j,nl_k-no2+k)
               vz_g(i,j,-no2+k) = vz_g(i,j,nl_k-no2+k)
               c_g(i,j,-no2+k) = c_g(i,j,nl_k-no2+k)
               bx_g(i,j,-no2+k) = bx_g(i,j,nl_k-no2+k)
               by_g(i,j,-no2+k) = by_g(i,j,nl_k-no2+k)
               bz_g(i,j,-no2+k) = bz_g(i,j,nl_k-no2+k)
               bk_g(i,j,-no2+k) = bk_g(i,j,nl_k-no2+k)
c
               rho_g(i,j,nl_k+k) = rho_g(i,j,k)
               vx_g(i,j,nl_k+k) = vx_g(i,j,k)
               vy_g(i,j,nl_k+k) = vy_g(i,j,k)
               vz_g(i,j,nl_k+k) = vz_g(i,j,k)
               c_g(i,j,nl_k+k) = c_g(i,j,k)
               bx_g(i,j,nl_k+k) = bx_g(i,j,k)
               by_g(i,j,nl_k+k) = by_g(i,j,k)
               bz_g(i,j,nl_k+k) = bz_g(i,j,k)
               bk_g(i,j,nl_k+k+1) = bk_g(i,j,k+1)
            enddo
         enddo
      enddo
      do k=1,no2-1
         do j=jlow,jhigh
            do i=ilow,ihigh
               bk_g(i,j,nl_k+k+1) = bk_g(i,j,k+1)
            enddo
         enddo
      enddo
      do k=1,no2
         do j=jlow,jhigh
            do i=ilow,ihigh
               bj_g(i,j,-no2+k) = bj_g(i,j,nl_k-no2+k)
               bj_g(i,j,nl_k+k) = bj_g(i,j,k)
            enddo
         enddo
      enddo
      do k=1,no2
         do j=jlow,jhigh
            do i=ilow,ihigh
               bi_g(i,j,-no2+k) = bi_g(i,j,nl_k-no2+k)
               bi_g(i,j,nl_k+k) = bi_g(i,j,k)
            enddo
         enddo
      enddo
!
! ... End subroutine KBOUND_FOR.........................................
!
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> JBOUND_FOR
!! 
      subroutine JBOUND_FOR
!
! ... 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 boundaries.F::JBOUND_FOR()"
#endif
      if ( .not. mhd_pe ) return
c
      if ( jlow_bound ) then
         do j=jlow,0
         do k=1,nl_k
            koff = mod(k+nl_k/2+1,nl_k)-1
            do i=ilow,ihigh
               rho_g(i,j,k) = rho_g(i,1-j,koff)
               vx_g(i,j,k) = vx_g(i,1-j,koff)
               vy_g(i,j,k) = vy_g(i,1-j,koff)
               vz_g(i,j,k) = vz_g(i,1-j,koff)
               c_g(i,j,k) = c_g(i,1-j,koff)
               bx_g(i,j,k) = bx_g(i,1-j,koff)
               by_g(i,j,k) = by_g(i,1-j,koff)
               bz_g(i,j,k) = bz_g(i,1-j,koff)
               bi_g(i,j,k) = bi_g(i,1-j,koff)
            enddo
         enddo
      enddo
      do j=jlow,0
         do k=1,nl_k
            koff = mod(k+nl_k/2+1,nl_k)-1
            do i=ilow,ihigh
               bj_g(i,j,k) = -bj_g(i,2-j,koff)
            enddo
         enddo
      enddo
      do j=jlow,0
         do k=1,nl_k+1
            koff = mod(k+nl_k/2+1,nl_k)-1
            do i=ilow,ihigh
               bk_g(i,j,k) = -bk_g(i,1-j,koff)
            enddo
         enddo
      enddo
      endif
c
      if ( jhigh_bound ) then
         do j=1,no2
         do k=1,nl_k
            koff = mod(k+nl_k/2+1,nl_k)-1
            do i=ilow,ihigh
              rho_g(i,nl_j+j,k) = rho_g(i,nl_j+1-j,koff)
               vx_g(i,nl_j+j,k) = vx_g(i,nl_j+1-j,koff)
               vy_g(i,nl_j+j,k) = vy_g(i,nl_j+1-j,koff)
               vz_g(i,nl_j+j,k) = vz_g(i,nl_j+1-j,koff)
               c_g(i,nl_j+j,k) = c_g(i,nl_j+1-j,koff)
               bx_g(i,nl_j+j,k) = bx_g(i,nl_j+1-j,koff)
               by_g(i,nl_j+j,k) = by_g(i,nl_j+1-j,koff)
               bz_g(i,nl_j+j,k) = bz_g(i,nl_j+1-j,koff)
               bi_g(i,nl_j+j,k) = bi_g(i,nl_j+1-j,koff)
            enddo
         enddo
      enddo
      do j=1,no2-1
         do k=1,nl_k
            koff = mod(k+nl_k/2+1,nl_k)-1
            do i=ilow,ihigh
               bj_g(i,nl_j+1+j,k) = -bj_g(i,nl_j+1-j,koff)
            enddo
         enddo
      enddo
      do j=1,no2
         do k=1,nl_k+1
            koff = mod(k+nl_k/2+1,nl_k)-1
            do i=ilow,ihigh
               bk_g(i,nl_j+j,k) = -bk_g(i,nl_j+1-j,koff)
            enddo
         enddo
      enddo

      endif
!
! ... End subroutine JBOUND_FOR.........................................
!
      return
      end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> IBOUND_FOR
!!
      subroutine IBOUND_FOR
!
! ... 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 "boundx.inc"
#include "run-time.inc"
#include "var2.inc"
      logical offtable
      data cone_0,cone_0sq,cwidth/0.3,0.09,12./
      data small/1.e-20/
      real*4 time4
!
! ... Begin ............................................................
!
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in boundaries.F::IBOUND_FOR()"
#endif
      time4 = real(time,4)
c
      if ( .not. mhd_pe ) return
c
      if ( ilow_bound ) then
         do k=1,nl_k
            do j=jlow,jhigh
               do i=1,no2
                  rho_g(-no2+i,j,k) = rho_g(no2+1-i,j,k)
                  vx_g(-no2+i,j,k) = vx_g(no2+1-i,j,k)
                  vy_g(-no2+i,j,k) = vy_g(no2+1-i,j,k)
                  vz_g(-no2+i,j,k) = vz_g(no2+1-i,j,k)
                  c_g(-no2+i,j,k) = c_g(no2+1-i,j,k)
                  bx_g(-no2+i,j,k) = bx_g(no2+1-i,j,k)
                  by_g(-no2+i,j,k) = by_g(no2+1-i,j,k)
                  bz_g(-no2+i,j,k) = bz_g(no2+1-i,j,k)
                  bi_g(-no2+i,j,k) = bi_g(no2+1-i,j,k)
*                  bj_g(-no2+i,j,k) = -bj_g(no2+1-i,j,k)
*                  bk_g(-no2+i,j,k) = -bk_g(no2+1-i,j,k)
*  as featured in the old-code: 
                  bj_g(-no2+i,j,k) = bj_g(1,j,k)
                  bk_g(-no2+i,j,k) = bk_g(1,j,k)
               enddo
            enddo
         enddo
c
c  set slip hardwall boundary
#ifndef NO_HARD_WALL
         do k=1,nl_k
            do j=jlow,jhigh
               vn = tran_i(1,j,k,1)*vx_g(1,j,k) +
     $              tran_i(1,j,k,2)*vy_g(1,j,k) +
     $              tran_i(1,j,k,3)*vz_g(1,j,k)
               do i=1,no2
                  vx_g(-no2+i,j,k) = vx_g(no2+1-i,j,k)-
     $                 2.0*vn*tran_i(1,j,k,1)
                  vy_g(-no2+i,j,k) = vy_g(no2+1-i,j,k)-
     $                 2.0*vn*tran_i(1,j,k,2)
                  vz_g(-no2+i,j,k) = vz_g(no2+1-i,j,k)-
     $                 2.0*vn*tran_i(1,j,k,3)
               enddo
            enddo
         enddo
#endif
c     
      endif
c

      if ( ihigh_bound ) then
c
c   set up boundaries
c
#ifdef CONSTTEST
*
*   going to use constant boundaries for test
*
         rho_q = 1.1e-23
         vx_q = -4.e7
         vy_q = 0.
         vz_q = 0.
         c_q = 4.e6
         bx_q = 0.
         by_q = 0.
         bz_q = 5.e-5
*
* end of constants
*
#endif
       if ( by_coef .eq. 0.0 .and. bz_coef .eq. 0.0 ) then
        cos_f = 1.0
        sin_f = 0.0
        cos_b_phi = 1.0
        sin_b_phi = 0.0
        vfrx = vfx
        vfry = vfy
        vfrz = vfz
        bx_0 = by_0*by_coef + bz_0*bz_coef + bx_zero
        bx_1 = by_1*by_coef + bz_1*bz_coef + bx_zero
      else
        bcc = sqrt(by_coef**2+bz_coef**2)
        cos_b_phi = by_coef/bcc
        sin_b_phi = bz_coef/bcc
        btt = sqrt(1.+bcc**2)
        cos_f = 1./btt
        sin_f = bcc/btt
        tan_f = sin_f/cos_f
        bx_0 = by_0*by_coef + bz_0*bz_coef + bx_zero
        bx_1 = by_1*by_coef + bz_1*bz_coef + bx_zero
        vfrx = vfx*cos_f**2
        vfry = -vfx*cos_f*sin_f*cos_b_phi
        vfrz = -vfx*cos_f*sin_f*sin_b_phi
      endif
c
        ind = max(0,min(itable,nint((time4-tzero)/deltaT)))
        if ( ind .eq. 0 ) then
           offtable = .true.
           ind = 1
        else 
           offtable = .false.
        endif
           vffx = vjx(ind)*cos_F**2
           vffy = -vjx(ind)*cos_f*sin_f*cos_b_phi
           vffz = -vjx(ind)*cos_f*sin_f*sin_b_phi
        xfront = (time4-tzero)*vffx
        yfront = (time4-tzero)*vffy
        zfront = (time4-tzero)*vffz
        vfront = sqrt(vffx**2+vffy**2+vffz**2)
        sin_tilt = sin(Zangle(ind))
        cos_tilt = cos(Zangle(ind))
c
         do k=1,nl_k
            do j=jlow,jhigh
               do i=1,no2
                  xdd = 0.125*(x(nl_i+i,j,k)+x(nl_i+1+i,j,k)+
     $                 x(nl_i+i,j+1,k)+x(nl_i+1+i,j+1,k)+
     $                 x(nl_i+i,j,k+1)+x(nl_i+1+i,j,k+1)+
     $                 x(nl_i+i,j+1,k+1)+x(nl_i+1+i,j+1,k+1))
                  ydd = 0.125*(y(nl_i+i,j,k)+y(nl_i+1+i,j,k)+
     $                 y(nl_i+i,j+1,k)+y(nl_i+1+i,j+1,k)+
     $                 y(nl_i+i,j,k+1)+y(nl_i+1+i,j,k+1)+
     $                 y(nl_i+i,j+1,k+1)+y(nl_i+1+i,j+1,k+1))
                  zdd = 0.125*(z(nl_i+i,j,k)+z(nl_i+1+i,j,k)+
     $                 z(nl_i+i,j+1,k)+z(nl_i+1+i,j+1,k)+
     $                 z(nl_i+i,j,k+1)+z(nl_i+1+i,j,k+1)+
     $                 z(nl_i+i,j+1,k+1)+z(nl_i+1+i,j+1,k+1))
                  tm = time4 - tzero - (vffx*xdd +vffy*ydd + vffz*zdd)
     $                 /vfront**2
                  index_now = max(0,min(itable,nint(tm/deltaT)))
                  if ( index_now .eq. 0 ) then
                     offtable = .true.
                     index_now = 1
                  else 
                     offtable = .false.
                  endif
                  xpp = xdd*cos_tilt - zdd*sin_tilt
                  zpp = xdd*sin_tilt + zdd*cos_tilt
                  cone = (ydd**2+zpp**2+small)/
     $                 (xpp**2+ydd**2+zpp**2+small)
             if ( xdd .lt. 0.0 .and. cone .lt. cone_0sq ) then
                  rho_g(nl_i+i,j,k) = rho_g(nl_i,j,k)*drop(i)
                  vx_g(nl_i+i,j,k) = vx_g(nl_i,j,k)
                  vy_g(nl_i+i,j,k) = vy_g(nl_i,j,k)
                  vz_g(nl_i+i,j,k) = vz_g(nl_i,j,k)
                  c_g(nl_i+i,j,k) = C_g(nl_i,j,k)
                  bx_g(nl_i+i,j,k) = bx_g(nl_i,j,k)
                  by_g(nl_i+i,j,k) = By_g(nl_i,j,k)
                  bz_g(nl_i+i,j,k) = Bz_g(nl_i,j,k)

             else
#ifndef CONSTTEST
               if ( offtable ) then
                  bx_g(nl_i+i,j,k) = bx_zero
                  by_g(nl_i+i,j,k) = 0.
                  bz_g(nl_i+i,j,k) = 0.
                  index_now = 1
               else
                  bx_g(nl_i+i,j,k) = BjY(index_now)*by_coef +
     $                 BjZ(index_now)*bz_coef + bx_zero
                  by_g(nl_i+i,j,k) = BjY(index_now)
                  bz_g(nl_i+i,j,k) = BjZ(index_now)
               endif
                  rho_g(nl_i+i,j,k) = Njp(index_now)
                  vx_g(nl_i+i,j,k) = Vjx(index_now)
                  vy_g(nl_i+i,j,k) = Vjy(index_now)
                  vz_g(nl_i+i,j,k) = Vjz(index_now)
                  c_g(nl_i+i,j,k) = Cs(index_now)

#else
                  rho_g(nl_i+i,j,k) = rho_q
                  vx_g(nl_i+i,j,k) = vx_q
                  vy_g(nl_i+i,j,k) = vy_q
                  vz_g(nl_i+i,j,k) = vz_q
                  c_g(nl_i+i,j,k) = C_q
                  bx_g(nl_i+i,j,k) = bx_q
                  by_g(nl_i+i,j,k) = By_q
                  bz_g(nl_i+i,j,k) = Bz_q
#endif
             endif     
               enddo
            enddo
         enddo
c
c   extra k faces
c
         do k=1,nl_k+1
            do j=jlow,jhigh
               do i=1,no2
                  xdd = 0.25*(x(nl_i+i,j,k)+x(nl_i+1+i,j,k)+
     $                 x(nl_i+i,j+1,k)+x(nl_i+1+i,j+1,k))
                  ydd = 0.25*(y(nl_i+i,j,k)+y(nl_i+1+i,j,k)+
     $                 y(nl_i+i,j+1,k)+y(nl_i+1+i,j+1,k))
                  zdd = 0.25*(z(nl_i+i,j,k)+z(nl_i+1+i,j,k)+
     $                 z(nl_i+i,j+1,k)+z(nl_i+1+i,j+1,k))
                  tm = time4 - tzero - (vffx*xdd +vffy*ydd + vffz*zdd)
     $                 /vfront**2
                  index_now = max(0,min(itable,nint(tm/deltaT)))
                  if ( index_now .eq. 0 ) then
                     offtable = .true.
                     index_now = 1
                  else 
                     offtable = .false.
                  endif
               if ( offtable ) then
                  bxdd = bx_zero
                  bydd = 0.
                  bzdd = 0.
                  index_now = 1
               else
                  bxdd = BjY(index_now)*by_coef +
     $                 BjZ(index_now)*bz_coef + bx_zero
                  bydd = BjY(index_now)
                  bzdd = BjZ(index_now)
               endif
#ifndef CONSTTEST
                  bk_g(nl_i+i,j,k) = (
     $                 tran_k(nl_i+i,j,k,1)*bxdd +
     $                 tran_k(nl_i+i,j,k,2)*bydd +
     $                 tran_k(nl_i+i,j,k,3)*bzdd )
#else
                  bk_g(nl_i+i,j,k) = (
     $                 tran_k(nl_i+i,j,k,1)*bx_q +
     $                 tran_k(nl_i+i,j,k,2)*by_q +
     $                 tran_k(nl_i+i,j,k,3)*bz_q )
#endif
               enddo
            enddo
         enddo
c
c   extra j faces
c
         do k=1,nl_k+1
            do j=jlow,jhigh
               do i=1,no2
                  xdd = 0.25*(x(nl_i+i,j,k)+x(nl_i+1+i,j,k)+
     $                 x(nl_i+i,j,k+1)+x(nl_i+1+i,j,k+1))
                  ydd = 0.25*(y(nl_i+i,j,k)+y(nl_i+1+i,j,k)+
     $                 y(nl_i+i,j,k+1)+y(nl_i+1+i,j,k+1))
                  zdd = 0.25*(z(nl_i+i,j,k)+z(nl_i+1+i,j,k)+
     $                 z(nl_i+i,j,k+1)+z(nl_i+1+i,j,k+1))
                  tm = time4 - tzero - (vffx*xdd +vffy*ydd + vffz*zdd)
     $                 /vfront**2
                  index_now = max(0,min(itable,nint(tm/deltaT)))
                  if ( index_now .eq. 0 ) then
                     offtable = .true.
                     index_now = 1
                  else 
                     offtable = .false.
                  endif
               if ( offtable ) then
                  bxdd = bx_zero
                  bydd = 0.
                  bzdd = 0.
                  index_now = 1
               else
                  bxdd = BjY(index_now)*by_coef +
     $                 BjZ(index_now)*bz_coef + bx_zero
                  bydd = BjY(index_now)
                  bzdd = BjZ(index_now)
               endif
#ifdef CONSTTEST
                  bj_g(nl_i+i,j,k) = (
     $                 tran_j(nl_i+i,j,k,1)*bx_q +
     $                 tran_j(nl_i+i,j,k,2)*by_q +
     $                 tran_j(nl_i+i,j,k,3)*bz_q )
#else
                  bj_g(nl_i+i,j,k) = (
     $                 tran_j(nl_i+i,j,k,1)*bxdd +
     $                 tran_j(nl_i+i,j,k,2)*bydd +
     $                 tran_j(nl_i+i,j,k,3)*bzdd )
#endif
               enddo
            enddo
         enddo
c
c   k electric field
c
         do k=1,nl_k
            do j=1,nl_j+1
               xdd = 0.5*(x(nl_i+1,j,k)+x(nl_i+1,j,k+1))
               ydd = 0.5*(y(nl_i+1,j,k)+y(nl_i+1,j,k+1))
               zdd = 0.5*(z(nl_i+1,j,k)+z(nl_i+1,j,k+1))
                  tm = time4 - tzero - (vffx*xdd +vffy*ydd + vffz*zdd)
     $                 /vfront**2
                  index_now = max(0,min(itable,nint(tm/deltaT)))
                  if ( index_now .eq. 0 ) then
                     offtable = .true.
                     index_now = 1
                  else 
                     offtable = .false.
                  endif
                  xpp = xdd*cos_tilt - zdd*sin_tilt
                  zpp = xdd*sin_tilt + zdd*cos_tilt
                  cone = sqrt((ydd**2+zpp**2+small)/
     $                 (xpp**2+ydd**2+zpp**2+small))
#ifndef CONSTTEST
               vxdd = VjX(index_now)
               vydd = VjY(index_now)
               vzdd = VjZ(index_now)
               if ( offtable ) then
                  bxdd = bx_zero
                  bydd = 0.
                  bzdd = 0.
                  index_now = 1
               else
                  bxdd = BjY(index_now)*by_coef +
     $                 BjZ(index_now)*bz_coef + bx_zero
                  bydd = BjY(index_now)
                  bzdd = BjZ(index_now)
               endif
               exd = vydd*bzdd - bydd*vzdd
               eyd = vzdd*bxdd - vxdd*bzdd
               ezd = vxdd*bydd - vydd*bxdd
#else
               vxdd = vx_q
               vydd = vy_q
               vzdd = vz_q
                  bxdd = bx_q
                  bydd = by_q
                  bzdd = bz_q
               exd = vydd*bzdd - bydd*vzdd
               eyd = vzdd*bxdd - vxdd*bzdd
               ezd = vxdd*bydd - vydd*bxdd
#endif
             if ( xdd .lt. 0.0) then
               bmix_ij(j,k) = 0.5*(1.-tanh((cone-cone_0)*cwidth))
             else
                bmix_ij(j,k) = 0.0
             endif
               ek_g(nl_I+1,j,k) = edge_k(nl_i+1,j,k)* (
     $              exd*vec_k(nl_i+1,j,k,7) +
     $              eyd*vec_k(nl_i+1,j,k,8) +
     $              ezd*vec_k(nl_i+1,j,k,9) )
            enddo
         enddo
c
c   j electric field
c
         do k=1,nl_k+1
            do j=1,nl_j
               xdd = 0.5*(x(nl_i+1,j,k)+x(nl_i+1,j+1,k))
               ydd = 0.5*(y(nl_i+1,j,k)+y(nl_i+1,j+1,k))
               zdd = 0.5*(z(nl_i+1,j,k)+z(nl_i+1,j+1,k))
                  tm = time4 - tzero - (vffx*xdd +vffy*ydd + vffz*zdd)
     $                 /vfront**2
                  index_now = max(0,min(itable,nint(tm/deltaT)))
                  if ( index_now .eq. 0 ) then
                     offtable = .true.
                     index_now = 1
                  else 
                     offtable = .false.
                  endif
                  xpp = xdd*cos_tilt - zdd*sin_tilt
                  zpp = xdd*sin_tilt + zdd*cos_tilt
                  cone = sqrt((ydd**2+zpp**2+small)/
     $                 (xpp**2+ydd**2+zpp**2+small))         
#ifndef CONSTTEST
               vxdd = VjX(index_now)
               vydd = VjY(index_now)
               vzdd = VjZ(index_now)
               if ( offtable ) then
                  bxdd = bx_zero
                  bydd = 0.
                  bzdd = 0.
                  index_now = 1
               else
                  bxdd = BjY(index_now)*by_coef +
     $                 BjZ(index_now)*bz_coef + bx_zero
                  bydd = BjY(index_now)
                  bzdd = BjZ(index_now)
               endif
               exd = vydd*bzdd - bydd*vzdd
               eyd = vzdd*bxdd - vxdd*bzdd
               ezd = vxdd*bydd - vydd*bxdd
#else
               vxdd = Vx_q
               vydd = Vy_q
               vzdd = Vz_q
                  bxdd = bx_q
                  bydd = By_q
                  bzdd = Bz_q
               exd = vydd*bzdd - bydd*vzdd
               eyd = vzdd*bxdd - vxdd*bzdd
               ezd = vxdd*bydd - vydd*bxdd
#endif
             if ( xdd .lt. 0.0) then
               bmix_ki(j,k) = 0.5*(1.-tanh((cone-cone_0)*cwidth))
             else
                bmix_ki(j,k) = 0.0
             endif
               ej_g(nl_I+1,j,k) = edge_j(nl_i+1,j,k)* (
     $              exd*vec_j(nl_i+1,j,k,7) +
     $              eyd*vec_j(nl_i+1,j,k,8) +
     $              ezd*vec_j(nl_i+1,j,k,9) )
            enddo
         enddo
c
         endif
!
! ... End subroutine IBOUND_FOR(...)d...................................
!
      return
      end

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