*********************************************************

      subroutine emfldspar5(t,x,y,z,bx,by,bz,nStts)

      include 'maxdef.h'

      double precision t, x, y, z, bx, by, bz
      integer nStts

c     This next two line lead to wastage of memory 
c     if you really care about memeory usuage
c     just use one of the array that you need.

      double precision bx1d,by1d,bz1d 
      common /fieldslab/ bx1d(NZMAX+2),by1d(NZMAX+2),bz1d

      double precision bx2d,by2d 
      common /field2d/ bx2d(NNMAX),by2d(NNMAX)          !bx2d(NX,NY),by2d(NX,NY)

c
c     Given the time and particle position (x,y,z) evaluate the
c     Cartesian electric and magnetic field components at 
c     ex(t,x,y,z), ey(t,x,y,z), ... etc. This routine uses
c     trilinear interpolation to interpolate from the
c     3D grid to the particle position, and linear interpolation
c     in the time domain.
c
c     INPUT ARGUMENTS
c
c        t -- The current time
c
c        x, y, z -- The (particle) position at which to evaluate 
c                   the fields
c
c     OUTPUT ARGUMENTS
c
c        ex,ey,ez -- The cartesian components of the electric field
c                    evaluated at the position (x,y,z).
c
c        bx,by,bz -- The cartesian components of the magnetic field
c                    evaluated at the position (x,y,z).
c
c
c     COMMON BLOCKS USED:

      integer NX, NY, NZ
      double precision lenx,leny,lenz
      common /box/ lenx,leny,lenz,NX,NY,NZ


      double precision delx,dely,delz
      common /grids/ delx,dely,delz
c
      logical debug
      common /debug/ debug

      double precision  b0, dbob
      common /bfield/  b0, dbob

      integer nkmin1d,nkmax1d
      double precision alpha1,lambda_parx,lambda_pary
      double precision eslab,fracx
      common /slabmodel/ eslab,fracx,lambda_parx,lambda_pary,alpha1,
     &                  nkmin1d,nkmax1d
      
      integer nxhalf, nyhalf
      double precision kmin,kmax,xi
      double precision alpha2,lambda_per
      double precision e2d
      common /twodmodel/ e2d,lambda_per,alpha2,kmin,kmax,
     &                  nxhalf,nyhalf,xi
 
      integer iseed_1,iseed_2,iseed_3
      common /seed/ iseed_1,iseed_2,iseed_3

      integer intinput(100)
      common /intinput/ intinput

      double precision dbleinput(100)
      common /dbleinput/ dbleinput

      integer printb
      common /printb/ printb

      double precision alpha3
      common /particle/ alpha3

      double precision th0, ph0
      common /inisph/ th0, ph0

c     LOCALS
c
      double precision xx,yy,zz
      double precision bxslab,byslab
      double precision t1
      double precision rr0

c     double precision exslab,eyslab,ezslab
      double precision bx2d2d,by2d2d

      double precision PI
      parameter (PI = 3.14159265358979323846d0)
c
c     Apply the periodic boundary conditions to the particle. We map
c     the particle position back into the unit cube. That is we map the
c     position back into the range 0 <= x < lenx, 0 <= y < leny,
c     0 <= z < lenz. This mapped position is used to calculate the 
c     EM fields. The particle position itself is unchanged, i.e. we 
c     replicate the field as necessary outside the primary box.     
c
c     This is the place to add the rest of the parameter
c     that is needed to be loaded into the common block
c     which is read in from int_values and double_values

      printb  = intinput(1)
      iseed_1 = intinput(2)
      iseed_2 = intinput(3)
      iseed_3 = intinput(4)
      nkmin1d = intinput(5)
      nkmax1d = intinput(6)
      nxhalf  = intinput(7)
      nyhalf  = intinput(8)

      eslab       = dbleinput(1) 
      fracx       = dbleinput(2)
      lambda_parx = dbleinput(3)
      lambda_pary = dbleinput(4)
      alpha1      = dbleinput(5)

      e2d         = dbleinput(6)
      lambda_per  = dbleinput(7)
      alpha2      = dbleinput(8)

      xi          = dbleinput(9)
      kmin        = dbleinput(10)
      kmax        = dbleinput(11)

      alpha3      = dbleinput(12)
      th0         = dbleinput(13)
      ph0         = dbleinput(14)
      rr0         = dbleinput(15)

c      write(*,*) iseed_1, iseed_2, iseed_3, nkmin, nkmax
c      write(*,*) eslab, fracx, lambda, alpha
c
c      if ( x .lt. 0.0d0 ) then
c        x = x+(1+int(-x/lenx))*lenx
c      endif
c      x = x-int(x/lenx)*lenx

c      if ( y .lt. 0.0d0 ) then
c        y = y+(1+int(-y/leny))*leny
c      endif
c      y = y-int(y/leny)*leny

c      if ( z .lt. 0.0d0 ) then
c        z = z+(1+int(-z/lenz))*lenz
c      endif
c      z = z-int(z/lenz)*lenz

      t1 = 0.0
c
c     Define xx, yy, zz -- the ptcl position in units
c     of dx, dy, dz in preparation for call to trilnt.
c

      call mhdfldsslab(t1)
      call mhdflds2d(t1)


      zz = z/delz 
      xx = (x-ph0)/delx
      yy = (y-th0)/dely

      if(zz .gt. NZ) then  
          nStts = 0
      else if(zz .lt. 0) then
          nStts = 0
      endif

c      if(zz .gt. NZ) then  
c          zz = zz-int(zz/NZ)*NZ
c      else if(zz .lt. 0) then
c          zz = zz+(1+int(-zz/NZ))*NZ
c      endif

c     Interpolate the fields from their fixed grid positions 
c     (i*dx, j*dy, k*dz) to the (intergrid) particle 
c     position (x,y,z)

      call lint(bx1d,NZ,zz,bxslab)
      call lint(by1d,NZ,zz,byslab)

      call lintarea2(bx2d,NX,NY,xx,yy,bx2d2d)
      call lintarea2(by2d,NX,NY,xx,yy,by2d2d)

      if (eslab .eq. 0.0) then
         bxslab = 0.0
         byslab = 0.0
      endif
      if (e2d .eq. 0.0) then
         bx2d2d = 0.0
         by2d2d = 0.0
      endif

      bx = (bxslab*rr0/z)+bx2d2d
      by = (byslab*rr0/z)+by2d2d
      bz = b0*(rr0*rr0)/(z*z)

      return
      end
