!  Prepared for the GEM Magnetic Reconnection Challenge on November 22, 2000
!  at Lindau (MPAe) by Ken Nishikawa
!
!  Attention: 
!    This version does not include background plasmas. Since the simulation is
!    intedned for 2-dimenstional system, the third dimension is very small.
!    And this code can be extended for a full 3-D simulation by increasing the
!    third demendtion (myy).
!
!    This version runs on ORIGIN2000 with frotran77. For other machines 
!    it may require some modifications.
!
!    The diagnostic tools are not well adjusted for the reconnection study
!    contact Ken Nishikawa at kenichi@physics.rutgers.edu for assistance.
!
!  Any comments are welcome and please send them to Ken Nishikawa
!
!     July 3, 2001
!
c  TRIdimensional STANford code, TRISTAN, fully electromagnetic,
C  with full relativistic particle dynamics. Written during spring
C  1990 by OSCAR BUNEMAN, with help from TORSTEN NEUBERT and
C  KEN NISHIKAWA.
c  Modified for Cray-2 by Ken Nishikawa
c
c  Whistler waves with puls electron beam
c  Modified for a study of magnetotail simulation
c       on Jan. 27, 1995 at Max-Plank-Institut fuer Extraterrische Physik, 
c                              Berlin
c   This code was modified for Harris model
c
c   1. Separate "input" 
c   2. Diagnosis for particle trajectories (dpartw.f < partcl00)
c   3. Diagnosis for local electric and magnetic fields
c                                          (dlocebw.f < locleb00)
c                       with spectra       (dlocebwl.f < locleb00)
c                a hodogram                (dbpol.f    < locleb00)
c   4. Diagnosis of electric and magnetic fields 
c        on the x - y plane at z = z1       (debcw.f < elfdxy00 emfdxy00)
c        on the x - z plane at y = y1       (debcw.f < elfdxz00 emfdxz00)
c        charge density (x-z) at y= y1      (debcw.f < vecptz00)
c        color for charge with magnetic f   (debcl.f < all)
c   5. Diagnosis for 2-d velocity space for electrons and ions
c        including localized beam and background electrons
c                                          (dveldw.f < veldel00 veldio00)
c   6. Diagnosis for electrons and ions in the phase spaces
c                 (x - vx, vy, vz)         (dpshock1d.f < prtl.d)
c 
c   7. Radiating boundaries      (October 30,1990)
c           x,z: radiating waves and reflecting particles
c             y: periodic for waves and particles
c
c            3 =< x < mx - 2
c            3 =< y < my - 2
c   This was changed after consideration
c            3 =< z < mz - 2
c
c   8. Density contours and flux in the x - z plane
c               in the sheet 42 < y < 44          (March 12, 1991)  
c                                           (densfvw.f < densfl00)
c
c   9. Diagnosis for energy (subroutine energy(t))
c              kinetic energy and field energy  (November 6, 1990)
c                                           (diagew.f < denerg00)
c
c  10. Charge density perturbation    (January, 1992)
c                rhox(mxx) bxcx(mxx)  at y=43, z=80
c                rhoz(mzz) bxcz(mzz)  at x=43, y=43
c                          bxsx(mxx)  at y=20  z=80
c                          bxsz(mzz)  at y=20  x=43
c                    profile of waves    (dshprof.f < drhoxz00)
c                    spectra (linear)    (drhospcb.f < drhoxz00)
c                    spectra (local)     (drhofun.f < drhoxz00)
c
c  11.   Conducting walls using 'conduct'  (January, 1992)
c
      program whist
      parameter(nptl=5000000,mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz,
     1      mx4=mxx+4,my4=myy+4)
      parameter(mbeam=9300,mback=9300)
      real me,mi
      logical in
      common /partls/x(nptl),y(nptl),z(nptl),u(nptl),v(nptl),w(nptl)
      common /fields/ex1(ngtl),ey1(ngtl),ez1(ngtl),bx1(ngtl),
     &by1(ngtl),bz1(ngtl)
     &,sm(27),q,qe,qi,qme,qmi,c
     &,ms(27),mx,my,mz,ix,iy,iz,lot,maxptl,maxhlf,ions,lecs,nstep
      common /veldis/ xvne,xvni,nvcell,nvcell1,nvcel2
      common /phas/ mcont(mbeam),mmax,kcont(mback),kmax
      common /para1/ b0x,b0y,b0z,deltm,vithml,vethml
      common /energ/ ntimee,me,mi,elosi,elose,anpo
      common /avec/ az(mx4,my4)                    
      common /jet/ xst,xlt,yst,ylt,zst,zlt,vdre,vdri
     1            ,xdel,ydel,zdel,avdn,xa,rand1
      common /tear/ alhc,zz0,xx1,xx2,xx3,ep0,ep1,ep2,ep3,and0
     1         ,ail
c
C  The sizes of the particle and field arrays must be chosen in
C  accordance with the requirements of the problem and the limitation
C  of the computer memory. The choice '8192' (with the bz1 array
C  curtailed to accommodate the remainder of the fields common) was
C  made to suit the segmented memory of a PC. Any changes in these
C  array sizes must be copied exactly into the COMMON statements of
C  the "surface", "mover" and "depsit" subroutines.
C Fields can be treated as single-indexed or triple-indexed:
C For CRAY-s, the first two field dimensions should be ODD, as here:
c
      dimension ex(mxx,myy,mzz),ey(mxx,myy,mzz),ez(mxx,myy,mzz),
     &bx(mxx,myy,mzz),by(mxx,myy,mzz),bz(mxx,myy,mzz)
c
c Definition of dimension for data on positions of electrons and ions
c for trajectories  dpart(3,40)
c
c Local electrostatic and magnetic fields   dleb(6,9) 
c
      dimension dpart(3,40), dleb(6,9)
c
c Electric and magnetic fields on x - y plane at z = z1
c                                 and on x - z plane at y = y1
c Vector potential (z - com) on x - y plane at z = z1
c
      dimension delzx(mxx,myy), delzy(mxx,myy), 
     &          demzx(mxx,myy), demzy(mxx,myy),
     &          delyx(mxx,mzz), delyz(mxx,mzz),
     &          demyx(mxx,mzz), demyz(mxx,mzz),
     &          dvecz(mxx,mzz)
c   according veclib 
c     integer*8 iseed
c     integer*4 iseed
c     real*8 dran
c
c 2-d velocity distributions for electrons and ions
c
c   Electrons:  suexy(200,200),  suexz(200,200)
c   Ions:       suixy(200,200),  suixz(200,200)
c
c
c density contours and flow velocity
c
c   Electrons:   dense(mxx,myy), fvlex(mxx,myy), fvley(mxx,myy)
c   Ions:        densi(mxx,myy), fvlix(mxx,myy), fvliy(mxx,myy)
c
c
c  charge density
c
      dimension rhox(mxx),rhoz(mxx)
      dimension bxcx(mxx),bxcz(mxx),bxsx(mxx),bxsz(mxx)
      dimension excx(mxx),excz(mxx),exsx(mxx),exsz(mxx)
      dimension ffx(mxx),ffz(mxx),fxx(mxx),fzz(mzz)
c
c--------------------------------------------
c
c     dimension aux(36),misc(38)
c    &,rho(65,65,64),tt(33,33,33)
c
      equivalence(ex1(1),ex(1,1,1)),(ey1(1),ey(1,1,1)),(ez1(1),ez(1,
     &1,1)),(bx1(1),bx(1,1,1)),(by1(1),by(1,1,1)),(bz1(1),bz(1,1,1))
c
c
c  ------- for the repeatable run
c
c     write(6,*) 'begining'
      last=1
c     last=2
c     last=10
c     write(6,*) 'begining 1'
      if(last.gt.1) go to 7521
c     last=100
c     last=1000
c     last=2000
c     last=3000
c     last=4000
c     last=5000
c     last=6000
c     last=7000
c     last=8000
c     last=9000
c     last=10000
c     last=11000
c     last=12000
c     last=13000
c     last=14000
c     last=15000
c     last=16000
c     last=17000
c     last=18000
c     last=19000
c     last=20000
c     if(last.gt.64) go to 7521
      open(unit=6,file='kenout5101r',status='new')
c     *  access='sequential')
c     write(6,*) 'begining 2'
      open(unit=7,file='fldsw.d',status='new',form='unformatted')
      open(unit=8,file='prtlw.d',status='new',form='unformatted')
      open(unit=14,file='partcl5101r',status='new',form='unformatted')
      open(unit=15,file='locleb5101r',status='new',form='unformatted')
      open(unit=16,file='elfdxy5101r',status='new',form='unformatted')
      open(unit=17,file='emfdxy5101r',status='new',form='unformatted')
      open(unit=18,file='elfdxz5101r',status='new',form='unformatted')
      open(unit=19,file='emfdxz5101r',status='new',form='unformatted')
      open(unit=20,file='vecptz5101r',status='new',form='unformatted')
      open(unit=21,file='veldel5101r',status='new',form='unformatted')
      open(unit=22,file='veldio5101r',status='new',form='unformatted')
      open(unit=23,file='phaseb5101r',status='new',form='unformatted')
      open(unit=24,file='phaseg5101r',status='new',form='unformatted')
      open(unit=25,file='densfv5101r',status='new',form='unformatted')
      open(unit=26,file='denerg5101r',status='new',form='unformatted')
      open(unit=27,file='drhoxz5101r',status='new',form='unformatted')
      write(6,*) 'before input0'
c   
c--------------------------------------------------------------------
c For initialization
c
      elosi=0.0
      elose=0.0
      anpo=0.0
C     Begin time stepping
      nstep=0
      ntimee=0
c
c     write(6,*) 'before input'
      call input
c     write(6,*) 'after input'

c
C Initialise the fields, typically to uniform components, such as
C just a uniform magnetic field parallel to the z-axis:
c
c  Initializing the magnetotail magnetic field  
c     

      do k=1,mz
      do j=1,my
      do i=1,mx
      ark=float(k)
      ffz(i)=(ark-zz0)/alhc
      ex(i,j,k)=0.0
      ey(i,j,k)=0.0
      ez(i,j,k)=0.0
      bx(i,j,k)=b0z*tanh(ffz(i))
      by(i,j,k)=0.0 
      bz(i,j,k)=0.0
      enddo
      enddo
      enddo
C (Remember that bx,by,bz are really c*Bx, c*By, c*Bz.)
C Initial fields, both electric and magnetic, must be divergence-free.
c
      go to 7816
 7521 continue
      open(unit=6,file='kenout5101r',status='unknown')
      open(unit=7,file='fldsw.d',status='old',form='unformatted')
      open(unit=8,file='prtlw.d',status='old',form='unformatted')
      open(unit=14,file='partcl5101r',status='unknown',
     1   form='unformatted')
      open(unit=15,file='locleb5101r',status='unknown',
     1   form='unformatted')
      open(unit=16,file='elfdxy5101r',status='unknown',
     1   form='unformatted')
      open(unit=17,file='emfdxy5101r',status='unknown',
     1   form='unformatted')
      open(unit=18,file='elfdxz5101r',status='unknown',
     1   form='unformatted')
      open(unit=19,file='emfdxz5101r',status='unknown',
     1   form='unformatted')
      open(unit=20,file='vecptz5101r',status='unknown',
     1   form='unformatted')
      open(unit=21,file='veldel5101r',status='unknown',
     1   form='unformatted')
      open(unit=22,file='veldio5101r',status='unknown',
     1   form='unformatted')
      open(unit=23,file='phaseb5101r',status='unknown',
     1   form='unformatted')
      open(unit=24,file='phaseg5101r',status='unknown',
     1   form='unformatted')
      open(unit=25,file='densfv5101r',status='unknown',
     1   form='unformatted')
      open(unit=26,file='denerg5101r',status='unknown',
     1   form='unformatted')
      open(unit=27,file='drhoxz5101r',status='unknown',
     1   form='unformatted')
c   
c--------------------------------------------------------------------
c
      read(7)ex,ey,ez,bx,by,bz,sm,ms,c,mx,my,mz,ix,iy,iz,lot,
     &     qe,qi,qme,qmi,nstep,b0x,b0y,b0z,az
     1   ,xst,xlt,yst,ylt,zst,zlt,vdre,vdri
     1            ,xdel,ydel,zdel,avdn,xa,rand1
      read(8)ions,lecs,maxptl,maxhlf,x,y,z,u,v,w,
     1       mcont,mmax,kcont,kmax,deltm,vithml,vethml,
     1       ntimee,me,mi,elosi,elose,anpo,
     1       xvne,xvni,nvcell,nvcell1,nvcel2
     1    , alhc,zz0,xx1,xx2,xx3,ep0,ep1,ep2,ep3,and0
      rewind(7)
      rewind(8)
 7816 continue
c
c perturbations

      alxx=mx-5.0
      alzz=mz-5.0
      psi0=0.1*b0z*ail
c     psi0=0.
      pi=3.141592654
      cxx=2.0*pi/alxx
      czz=pi/alzz
      b0xx=-pi*psi0/alzz
      b0zz=2.0*pi*psi0/alxx
c
c
c for dpart
      ncallp=1000 
c for dloceb
      ncalll=1000 
c for deba
      ncallf=1000
c     ncallf=5
c for deba (vector potential)
      ncallv=ncallf
c for dveld, dphas, densfv
      ncalld=1000
c     ncalld=5
c for energy
      ncalle=1000
c
c Initailization for dvecz(mxx,myy)
c
c     do 100 i = 1, mx 
c     do 100 j = 1, my 
c 100 dvecz(i,j)=0.0
c     write(6,*) 'after intit for dvecz'
c        
   6  nstep=nstep+1
c Attention
c the modification of current due to this perturbation is not included
c   July 5, 2001
      do k=3,mz-2
      do j=3,my-2
      do i=3,mx-2
      ari=float(i) -3.0
      ark=float(k) -3.0
      ffz(i)=(ark-zz0)/alhc
! attention
      fzz(k)=czz*ark -0.5*pi
      fxx(i)=cxx*ari -pi
      bx(i,j,k)=bx(i,j,k) +b0xx*cos(fxx(i))*sin(fzz(k))
      bz(i,j,k)=bz(i,j,k) +b0zz*sin(fxx(i))*cos(fzz(k))
      enddo
      enddo
      enddo
C Before moving particles, the magnetic field is Maxwell-advanced
C by half a timestep:
c     call preledge(by,bz,bx,ey,ez,ex,iy,iz,ix,my,mz,mx,1,c)
c     call preledge(bz,bx,by,ez,ex,ey,iz,ix,iy,mz,mx,my,1,c)
c     call preledge(bx,by,bz,ex,ey,ez,ix,iy,iz,mx,my,mz,1,c)
c 
      do 7 i=3,mx-3
      do 7 j=3,my-3
      do 7 k=1,mz-1
      bx(i,j,k)=bx(i,j,k) + (.5*c) *
     &          (ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k))
      by(i,j,k)=by(i,j,k) + (.5*c) *
     &          (ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k))
  7   bz(i,j,k)=bz(i,j,k) + (.5*c) *
     &          (ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k))
c
      t=deltm*nstep
c This was the (mx-3)*(my-3)*(mz-3) core where particles are initialised.
c They also want to know the fields just outside this core, at i=1,
c j=1, k=1 and i=mx-1, j=my-1, k=mz-1. Use the periodicity to get those
c fields: (Only z-direction is periodic)
      call copylayr(bx,by,bz,mx,my,mz,mx-2,3,my-2,3)
      call copylayr(bx,by,bz,mx,my,mz,2,mx-3,2,my-3)
c     call copylayr(bx,by,bz,mx,my,mz,my-2,3)       
c     call copylayr(bx,by,bz,mx,my,mz,2,my-3)
c
C Now move ions:
      call mover(1,ions,qmi)
C and electrons:
      call mover(1+maxhlf,lecs+maxhlf,qme)
c
c     if(mod(nstep,ncalle).eq.0) call energy(t)
c
C The Maxwell-advance of the fields begins with another half-step
C advance of the magnetic field since for Maxwell's equations the
C B - information and the E - information must be staggered in time.
C In space, their information is also staggered. Here we show the
C locations where field components are recorded:
C
C
C
C                z
C                ^
C                |
C                |
C                |
C                *---------------Ey---------------*
C               /|                                |
C              / |                                |
C             /  |                                |
C            Ex  |                                |
C           /    |                                |
C          /     |                                |
C         /      |                                |
C        *       |                                |
C        |       |                                |
C        |       |                                |
C        |       |                                |
C        |       Ez              Bx               Ez
C        |       |                                |
C        |       |                                |
C        |       |                                |
C        |   By  |                                |
C        |       |                                |
C        |       |                                |
C        |       |                                |
C        Ez      |                                |
C        |       |                                |
C        |       |                                |
C        |       |                                |
C        |       *---------------Ey---------------*-----> y
C        |      /                                /
C        |     /                                /
C        |    /                                /
C        |   Ex              Bz               Ex
C        |  /                                /
C        | /                                /
C        |/                                /
C        *---------------Ey---------------*
C       /
C      /
C     /
C    x
C
C
C     Ex(i,j,k) = value of Ex at x=i+.5, y=j,    z=k
C     Ey(i,j,k) = value of Ey at x=i,    y=j+.5, z=k
C     Ez(i,j,k) = value of Ez at x=i,    y=j,    z=k+.5
C
C     Bx(i,j,k) = value of Bx at x=i,    y=j+.5, z=k+.5
C     By(i,j,k) = value of By at x=i+.5, y=j,    z=k+.5
C     Bz(i,j,k) = value of Bz at x=i+.5, y=j+.5, z=k
C
C     Maxwell's laws are implemented (to central difference
C     accuracy) in the form:
C     
C     Change of flux of B through a cell face
C         = - circulation of E around that face
C
C     Change of flux of E through a cell face of the offset grid
C         = circulation of B around that face - current through it
C
C Second half-advance of magnetic field:
c 
      do 8 i=3,mx-3
      do 8 j=3,my-3
      do 8 k=1,mz-1
      bx(i,j,k)=bx(i,j,k) + (.5*c) *
     &          (ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k))
      by(i,j,k)=by(i,j,k) + (.5*c) *
     &          (ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k))
  8   bz(i,j,k)=bz(i,j,k) + (.5*c) *
     &          (ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k))
C Front,right and top layers of B must be obtained from a special
C boundary routine based on Lindman's method:
c Find out which is unnecessary ?????
c  Update B-fields on three surface:
c     call surface(by,bz,bx,ey,ez,ex,iy,iz,ix,my,mz,mx,1,c)
c     call surface(bz,bx,by,ez,ex,ey,iz,ix,iy,mz,mx,my,1,c)
      call surface(bx,by,bz,ex,ey,ez,ix,iy,iz,mx,my,mz,1,c)
c  And post-process the edge:
c     call postedge(by,bz,bx,ey,ez,ex,iy,iz,ix,my,mz,mx,1,c)
c     call postedge(bz,bx,by,ez,ex,ey,iz,ix,iy,mz,mx,my,1,c)
c     call postedge(bx,by,bz,ex,ey,ez,ix,iy,iz,mx,my,mz,1,c)
c Outer layers by periodicity: (only z)
      call copylayr(bx,by,bz,mx,my,mz,mx-2,3,my-2,3)
      call copylayr(bx,by,bz,mx,my,mz,2,mx-3,2,my-3)
c     call copylayr(bx,by,bz,mx,my,mz,my-2,3)       
c     call copylayr(bx,by,bz,mx,my,mz,2,my-3)        
c
c  Update E-field similarly, beginning with preliminary edge updates:
c     call preledge(ey,ez,ex,by,bz,bx,-iy,-iz,-ix,my,mz,mx,lot,c)
c     call preledge(ez,ex,ey,bz,bx,by,-iz,-ix,-iy,mz,mx,my,lot,c)
c     call preledge(ex,ey,ez,bx,by,bz,-ix,-iy,-iz,mx,my,mz,lot,c)
C Full advance of the electric field:
c 
      do 9 i=3,mx-3
      do 9 j=3,my-3
      do 9 k=2,mz
      ex(i,j,k)=ex(i,j,k) + c *
     &          (by(i,j,k-1)-by(i,j,k)-bz(i,j-1,k)+bz(i,j,k))
      ey(i,j,k)=ey(i,j,k) + c *
     &          (bz(i-1,j,k)-bz(i,j,k)-bx(i,j,k-1)+bx(i,j,k))
  9   ez(i,j,k)=ez(i,j,k) + c *
     &          (bx(i,j-1,k)-bx(i,j,k)-by(i-1,j,k)+by(i,j,k))
c????????????
c????????????
C  ----------
C Boundary values of the E - field must be provided at rear, left
C and bottom faces of the field domain:
c ?????? which is unnecessary?
c  and bottom faces of the field domain:
c     call surface(ey,ez,ex,by,bz,bx,-iy,-iz,-ix,my,mz,mx,lot,c)
c     call surface(ez,ex,ey,bz,bx,by,-iz,-ix,-iy,mz,mx,my,lot,c)
      call surface(ex,ey,ez,bx,by,bz,-ix,-iy,-iz,mx,my,mz,lot,c)
c  And post-processing of edge:
c     call postedge(ey,ez,ex,by,bz,bx,-iy,-iz,-ix,my,mz,mx,lot,c)
c     call postedge(ez,ex,ey,bz,bx,by,-iz,-ix,-iy,mz,mx,my,lot,c)
c     call postedge(ex,ey,ez,bx,by,bz,-ix,-iy,-iz,mx,my,mz,lot,c)
c attention
      call clearlay(ex,ey,ez,mx,my,mz,mx,my)   
      call clearlay(ex,ey,ez,mx,my,mz,mx-1,my-1)     
      call clearlay(ex,ey,ez,mx,my,mz,mx-2,my-2)      
      call clearlay(ex,ey,ez,mx,my,mz,2,2)   
      call clearlay(ex,ey,ez,mx,my,mz,1,1)   
c
C The currents due to the movement of each charge q are applied to the
C E-M fields as decrements of E-flux through cell faces. The movement
C of particles which themselves cross cell boundaries has to be split
C into several separate moves, each only within one cell. Each of
C these moves contributes to flux across twelve faces.
C Ions and electrons are processed in two loops, changing the sign
C of the charge in-between. These loops cannot be vectorised:
C particles get processed on by one. Here is a good place to
C insert the "if" clauses for applying boundary conditions to
C the particles, such as reflection, periodicity, replacement
C by inward moving thermal or streaming particles, etc.
C First clear peripheral layers to receive overflow currents:
c
C The currents due to the movement of each charge q are applied to the
C E-M fields as decrements of E-flux through cell faces. The movement
C of particles which themselves cross cell boundaries has to be split
C into several separate moves, each only within one cell. Each of
C these moves contributes to flux across twelve faces.
C Ions and electrons are processed in two loops, changing the sign
C of the charge in-between. These loops cannot be vectorised:
C particles get processed one by one. Here is a good place to
C insert the statements for applying boundary conditions to
C the particles, such as reflection, periodicity, replacement
C by inward moving thermal or streaming particles, etc.
C Split and deposit ions currents:
C -----------
C  The split routines call the deposit routine. They also return the verdict
C  whether the particles have left the computational domain, whether they are
C  " in " or " .not.in ". In the solar wind application we eliminate particles
C  which have left the domain (which are  .not.in). Such particles leave their
C  charges behind on the boundary surfaces.
c
      q=qi
      vdr=0.0
      vth2=vithml
      n=1
c  
  52  continue
C  Previous position:  for reflected particles
      x0=x(n)-u(n)
      y0=y(n)-v(n)
      z0=z(n)-w(n)
c  
c     if (x(n).ge.3..and.x(n).lt.mx-2..and.
c    &    z(n).ge.3..and.z(n).lt.mz-2.)go to 171
      if (z(n).ge.3..and.z(n).lt.mz-2.)go to 171
C  Find a border crossing:
c     if(x(n).ge.mx-2.) then
c     x(n)=amin1(amax1(x(n),3.),mx-2.000001)
      z(n)=amin1(amax1(z(n),3.),mz-2.000001)
c     t=amin1((x(n)-x0)/u(n),(y(n)-y0)/v(n))
c     tx=(x(n)-x0)/u(n)
c     ty=(y(n)-y0)/v(n)
      tz=(z(n)-z0)/w(n)
c     rabs=dran(iseed)
      rabs=2.0
c     if(tx.gt.tz) go to 705
C  Stop particles at first border encountered:
      y(n)=y0+tz*v(n)
      x(n)=x0+tz*u(n)
C  Randomise their velocity to a typical inward velocity:
c     u(n)=-sign(vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1                +dran(iseed)-2.0),u(n))
c atteintion
      u(n)=-u(n)   
c     v(n)=-v(n)
c     w(n)=-w(n)
c     v(n)=vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1     +dran(iseed)-2.0)
c     w(n)=vdr+vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1     +dran(iseed)-2.0)
c     if(rabs.gt.0.8) then     
      x(n)=x(n)+(1.0-tz)*u(n)
      y(n)=y(n)+(1.0-tz)*v(n)
      z(n)=2.0*zz0-z(n)+(1.0-tz)*w(n)
      z0=z(n)-(1.0-tz)*w(n)
c     endif
c     else
c     rabs=0.1
c     go to 172
c     endif
c
      go to 175
C  Stop particles at first border encountered:
c705  x(n)=x0+tz*u(n)
c     y(n)=y0+tz*v(n)
C  Randomise their velocity to a typical inward velocity:
c     u(n)=-u(n)
c     v(n)=v(n)
c     w(n)=w(n)
c     u(n)=vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1          +dran(iseed)-2.0)
c     v(n)=-sign(vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1                +dran(iseed)-2.0),v(n))
c attention  keep the same velocity after reflection
c     w(n)=vdr+vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1              +dran(iseed)-2.0)
c     if(rabs.gt.0.8) then     
c     x(n)=x(n)+(1.0-tz)*u(n)
c     y(n)=y(n)+(1.0-tz)*v(n)
c     z(n)=2.0*zz0-z(n)+(1.0-tz)*w(n)
c     z0=z(n)-w(n)
c     endif
  175 continue
      go to 172
c
C Make particles periodic in the axial dimension, with period mz-3:
c Attention
 171  continue
      rabs=2.0
 172  continue
      per=sign(.5*(my-5.),y(n)-3.)+sign(.5*(my-5.),y(n)-my+2.)
c attention
      y(n)=y(n)-per
      y0=y0-per   
c
      per=sign(.5*(mx-5.),x(n)-3.)+sign(.5*(mx-5.),x(n)-mx+2.)
c attention
      x(n)=x(n)-per
      x0=x0-per   
C Split particles which cross cell boundaries and deposit currents:
 153  call xsplit(x(n),y(n),z(n),x0,y0,z0)
c     if(rabs.gt.0.8) go to 58
c     if(rabs.gt.1.5) go to 58
c  checking the lost kinetic energy
c     elosi=elosi+0.5*mi*(u(n)*u(n)+v(n)*v(n)+w(n)*w(n))
C  Replace by the ion from the top of the stack:
c     x(n)=x(ions)
c     y(n)=y(ions)
c     z(n)=z(ions)
c     u(n)=u(ions)
c     v(n)=v(ions)
c     w(n)=w(ions)
c     ions=ions-1
c     n=n-1
  58  n=n+1
      if(n.le.ions)go to 52
      q=qe
      n=maxhlf+1
      vdr=0.0
      vth2=vethml
c attention
c  
  53  continue
C  Previous position:  for reflected particles
      x0=x(n)-u(n)
      y0=y(n)-v(n)
c attention
      z0=z(n)-w(n)
c  
c     if (x(n).ge.3..and.x(n).lt.mx-2..and.
c    &    z(n).ge.3..and.z(n).lt.mz-2.)go to 271
      if (z(n).ge.3..and.z(n).lt.mz-2.)go to 271
C  Find a border crossing:
c     if(x(n).ge.mx-2.) then
c     x(n)=amin1(amax1(x(n),3.),mx-2.000001)
      z(n)=amin1(amax1(z(n),3.),mz-2.000001)
c     t=amin1((x(n)-x0)/u(n),(y(n)-y0)/v(n))
c     tx=(x(n)-x0)/u(n)
      tz=(z(n)-z0)/w(n)
c     rabs=dran(iseed)
      rabs=2.0
c     if(tx.gt.tz) go to 805
C  Stop particles at first border encountered:
      y(n)=y0+tz*v(n)
      x(n)=x0+tz*u(n)
C  Randomise their velocity to a typical inward velocity:
c     u(n)=-sign(vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1                +dran(iseed)-2.0),u(n))
      u(n)=-u(n)  
c     v(n)=-v(n)
c     w(n)=-w(n)
c     v(n)=vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1          +dran(iseed)-2.0)
c     w(n)=vdr+vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1              +dran(iseed)-2.0)
c     if(rabs.gt.0.5) then     
      x(n)=x(n)+(1.0-tz)*u(n)
      y(n)=y(n)+(1.0-tz)*v(n)
      z(n)=2.0*zz0-z(n)+(1.0-tz)*w(n)
      z0=z(n)-(1.0-tz)*w(n)
c     endif
c     else
c     rabs=0.1       
c     go to 272 
c     endif
      go to 375
C  Stop particles at first border encountered:
c805  x(n)=x0+tz*u(n)
c     y(n)=y0+tz*v(n)
c  Randomise their velocity to a typical inward velocity:
c     u(n)=-u(n)
c     v(n)=v(n)
c     w(n)=w(n)
c     u(n)=vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1          +dran(iseed)-2.0)
c     v(n)=-sign(vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1                +dran(iseed)-2.0),v(n))
c attention  keep the same velocity after reflection
c     w(n)=vdr+vth2*(dran(iseed)+dran(iseed)+dran(iseed)
c    1              +dran(iseed)-2.0)
c     if(rabs.gt.0.5) then     
c     x(n)=x(n)+(1.0-tz)*u(n)
c     y(n)=y(n)+(1.0-tz)*v(n)
c     z(n)=2.0*zz0-z(n)+(1.0-tz)*w(n)
c     z0=z(n)-w(n)
c     endif
  375 continue
      go to 272
c
C Make particles periodic in the axial dimension, with period mz-3:
c Attention
 271  continue
      rabs=2.0
 272  continue
      per=sign(.5*(my-5.),y(n)-3.)+sign(.5*(my-5.),y(n)-my+2.)
c attention
      y(n)=y(n)-per
      y0=y0-per   
c
      per=sign(.5*(mx-5.),x(n)-3.)+sign(.5*(mx-5.),x(n)-mx+2.)
c attention
      x(n)=x(n)-per
      x0=x0-per   
C Split particles which cross cell boundaries and deposit currents:
 253  call xsplit(x(n),y(n),z(n),x0,y0,z0)
c     if(rabs.gt.0.5) go to 59
c     if(rabs.gt.1.5) go to 59
c  checking the lost kinetic energy
c     elose=elose+0.5*me*(u(n)*u(n)+v(n)*v(n)+w(n)*w(n))
C  Replace by the electron from the top of the stack:
c     x(n)=x(maxhlf+lecs)
c     y(n)=y(maxhlf+lecs)
c     z(n)=z(maxhlf+lecs)
c     u(n)=u(maxhlf+lecs)
c     v(n)=v(maxhlf+lecs)
c     w(n)=w(maxhlf+lecs)
c     lecs=lecs-1
c     n=n-1
  59  n=n+1
      if(n.le.maxhlf+lecs)go to 53
C  We spread out any surface charges, and let surface electrons recombine
C  with surface ions by making the surfaces slightly conducting:
c     call conduct(ix,iy,iz,mx,my,mz,ey,ez,.25)
c     call conduct(iy,iz,ix,my,mz,mx,ez,ex,.25)
c for periodic with z-direction the following line is not necessary
c     call conduct(iz,ix,iy,mz,mx,my,ex,ey,.25)
c
 173  call addlayer(ex,ey,ez,mx,my,mz,5,mx,5,my)     
      call addlayer(ex,ey,ez,mx,my,mz,4,mx-1,4,my-1)       
      call addlayer(ex,ey,ez,mx,my,mz,3,mx-2,3,my-2)       
      call addlayer(ex,ey,ez,mx,my,mz,mx-3,2,my-3,2)        
      call addlayer(ex,ey,ez,mx,my,mz,mx-4,1,my-4,1)       
c
c
      call copylayr(ex,ey,ez,mx,my,mz,mx-2,3,my-2,3)       
      call copylayr(ex,ey,ez,mx,my,mz,mx-1,4,my-1,4)        
      call copylayr(ex,ey,ez,mx,my,mz,mx,5,my,5)    
      call copylayr(ex,ey,ez,mx,my,mz,2,mx-3,2,my-3)       
      call copylayr(ex,ey,ez,mx,my,mz,1,mx-4,1,my-4)       
c
c
C  ---------------
C  Inject solar wind, according to the density, velocity, and time step
C  
      if(mod(nstep,ncalld).eq.0) then
c
c charge density record:
c
      nxpt=25 
      nypt=my/2 +1
      nzpt=zz0
         do 275 i=4,mx-3
         bxcx(i-3)=0.5*(bx(i,nypt,nzpt)+bx(i,nypt-1,nzpt-1))
         excx(i-3)=0.5*(by(i,nypt,nzpt)+by(i-1,nypt,nzpt-1))
         bxcz(i-3)=0.5*(bz(i,nypt,nzpt)+bz(i-1,nypt-1,nzpt))
         rhoz(i-3)=0.5*(ex(i-1,nypt,nzpt)+ex(i,nypt,nzpt))
         bxsx(i-3)=0.5*(ey(i,nypt-1,nzpt)+ey(i,nypt,nzpt))
         excz(i-3)=0.5*(ez(i,nypt,nzpt-1)+ez(i,nypt,nzpt))
         bxsz(i-3)=0.5*((by(i,nypt,nzpt)+by(i-1,nypt,nzpt))
     1             -(by(i,nypt,nzpt-1)+by(i-1,nypt,nzpt-1))
     1             -(bz(i,nypt,nzpt)+bz(i-1,nypt,nzpt))
     1             +(bz(i,nypt-1,nzpt)+bz(i-1,nypt-1,nzpt)))
         exsx(i-3)=0.5*((bz(i,nypt,nzpt)+bz(i,nypt-1,nzpt))
     1             -(bz(i-1,nypt,nzpt)+bz(i-1,nypt-1,nzpt))
     1             -(bx(i,nypt,nzpt)+bx(i,nypt-1,nzpt))
     1             +(bx(i,nypt,nzpt-1)+bx(i,nypt-1,nzpt-1)))
         exsz(i-3)=0.5*((bx(i,nypt,nzpt)+bx(i,nypt,nzpt-1))
     1             -(bx(i,nypt-1,nzpt)+bx(i,nypt-1,nzpt-1))
     1             -(by(i,nypt,nzpt)+by(i,nypt,nzpt-1))
     1             +(by(i-1,nypt,nzpt)+by(i-1,nypt,nzpt-1)))
 275     rhox(i-3)=ex(i,nypt,nzpt)-ex(i-1,nypt,nzpt)
     &            +ey(i,nypt,nzpt)-ey(i,nypt-1,nzpt)
     &            +ez(i,nypt,nzpt)-ez(i,nypt,nzpt-1)
c        do 276 k=4,mz-3
c        bxcz(k-3)=0.5*(bx(nxpt,nypt,k)+bx(nxpt+1,nypt,k))
c276     rhoz(k-3)=ex(nxpt,nypt,k)-ex(nxpt-1,nypt,k)
c    &            +ey(nxpt,nypt,k)-ey(nxpt,nypt-1,k)
c    &            +ez(nxpt,nypt,k)-ez(nxpt,nypt,k-1)
c  for the outside of the beam
c      nypt=5
c        do 277 i=4,mx-3
c277     bxsx(i-3)=0.5*(bx(i,nypt,nzpt)+bx(i+1,nypt,nzpt))
c        do 278 k=4,mz-3
c278     bxsz(k-3)=0.5*(bx(nxpt,nypt,k)+bx(nxpt+1,nypt,k))
c
      write(27) rhox,rhoz,bxcx,bxcz,bxsx,bxsz
     1                   ,excx,excz,exsx,exsz
C charge density record:
c        do 75 k=2,mz-2
c        do 75 j=2,my-2
c        do 75 i=2,mx-2
c 75     rho(i-1,j-1,k-1)=ex(i,j,k)-ex(i-1,j,k)+ey(i,j,k)-ey(i,j-1,k)
c    &            +ez(i,j,k)-ez(i,j,k-1)
c     do 76 k=1,64
c     call vfht(6,rho(1,1,k),1,65,64)
c 76  call vfht(6,rho(1,1,k),65,1,64)
c     do 77 j=1,64
c 77  call vfht(6,rho(1,j,1),4225,1,64)
c     do 78 k=-16,16
c     do 78 j=-16,16
c     do 78 i=-16,16
c 78  tt(i+17,j+17,k+17)=rho(33+i-isign(32,i),33+j-isign(32,j),33+k
c    &-isign(32,k))
c        write(8)tt
   54 continue

      endif
c
c  Velocity distributions
c 
c     if(mod(nstep,ncalld).eq.0.or.nstep.eq.1)
      if(mod(nstep,ncalld).eq.0)
     .call veldist(t)
c
c Phase space for beam and background electrons
c 
c     if(mod(nstep,ncalld).eq.0.or.nstep.eq.1)
c    .call diffus(t)
c
c Density and flow velocity
c
      if(mod(nstep,ncalld).eq.0.or.nstep.eq.1)
     .call densfv(t)
c
c  Attention
c
      kz1 = 0.5*mz+1
c  attention
c     jy1 = 43-1
      jy1 = 0.5*my+1
c
c     if(mod(nstep,ncallf).ne.0) go to 93
      if(mod(nstep,ncallf).eq.0.or.nstep.eq.1) then      
c
c Buffer out the data for local electric and magnetic fields
c
      do 91 i = 2, mxx-1
      do 91 j = 2, myy-1
      delzx(i,j)=0.5*(ex(i,j,kz1) +ex(i-1,j,kz1))    
      delzy(i,j)=0.5*(ey(i,j,kz1) +ey(i,j-1,kz1))    
      demzx(i,j)=0.5*(bx(i,j,kz1) +bx(i,j-1,kz1-1))    
      demzy(i,j)=0.5*(by(i,j,kz1) +by(i-1,j,kz1-1))    
   91 continue
c
      do 92 i = 2, mxx-1 
      do 92 k = 2, mzz-1
      delyx(i,k)=0.5*(ex(i,jy1,k) +ex(i-1,jy1,k))    
      delyz(i,k)=0.5*(ez(i,jy1,k) +ez(i,jy1,k-1))    
      demyx(i,k)=0.5*(bx(i,jy1,k) +bx(i,jy1-1,k-1))
      demyz(i,k)=0.5*(bz(i,jy1,k) +bz(i-1,jy1-1,k))
   92 continue
      do 921 k=1,mzz
      demyz(1,k)=0.0 
  921 continue
      delzx(1,1)=t
c
      write(16) delzx,delzy
c
      write(17) demzx,demzy
c
      write(18) delyx,delyz
c
      write(19) demyx,demyz
c
      endif
c  93 continue
c
c     do 95 i = 1, mx 
c     do 95 j = 1, my 
c     dvecz(i,j)=dvecz(i,j) +0.5*(ez(i,j,kz1)+ez(i,j,kz1-1))
c  95 continue
c
      if(mod(nstep,ncallv).eq.0.or.nstep.eq.1) then
      do 95 k=2,mz-2
      do 95 i=2,mx-2
 95   dvecz(i,k)=ex(i,jy1,k)-ex(i-1,jy1,k)
     &          +ey(i,jy1,k)-ey(i,jy1-1,k)
     &          +ez(i,jy1,k)-ez(i,jy1,k-1)
c    
      write(20) dvecz
      endif
c
c
c     if(mod(nstep,ncallp).ne.0) go to 68
      if(mod(nstep,ncallp).eq.0.or.nstep.eq.1) then     
c
c Buffer out the data for trajectories dpart(m,n)
c   m=1: x,    m=2: y,   m=3: z
c   n=1, 20: ions,     n=21, 40 electrons
c
c     write(6,*) 'x(162)=',x(162), 'z(5186)=',z(5186)  
      nint=0
c     ndelp=0.018*nptl
      ndelp=0.018*ions
      do 70 i = 1, 40
      nn = nint + i * ndelp
      if(i.ge.21)
     . nn = maxhlf+(i-20)*ndelp
      dpart(1,i) = x(nn)
      dpart(2,i) = y(nn)
  70  dpart(3,i) = z(nn)
c
      write(14) dpart
      endif
c 68  continue
c
c     if(mod(nstep,ncalle).eq.0) call energy(t)
c
c     if(mod(nstep,ncalll).ne.0) go to 158
      if(mod(nstep,ncalll).eq.0.or.nstep.eq.1) then      
c
c Buffer out the local electrostatic and elecromagnetic fields
c
      jz=0.5*mz+1
      jy=0.5*my+1
      ii=0
      do 72 i = 1, 9
      jx = 1 + 10*(i-1)
      do 71 k = 1, 1
      ii=ii+1
      dleb(1,ii) = 0.5*(ex(jx,jy,jz)+ex(jx-1,jy,jz))
      dleb(2,ii) = 0.5*(ey(jx,jy,jz)+ey(jx,jy-1,jz))
      dleb(3,ii) = 0.5*(ez(jx,jy,jz)+ez(jx,jy,jz-1))
      dleb(4,ii) = 0.5*(bx(jx,jy,jz)+bx(jx,jy-1,jz-1))
      dleb(5,ii) = 0.5*(by(jx,jy,jz)+by(jx-1,jy,jz-1))
      dleb(6,ii) = 0.5*(bz(jx,jy,jz)+bz(jx-1,jy-1,jz))
   71 continue
   72 continue
c
      write(15) dleb
c
      endif
C Countdown:
 158  continue
c     write(6,*) 'nstep= ',nstep
      if (nstep.le.last) go to 6
      write(6,*) 'last=',last
      write(7)ex,ey,ez,bx,by,bz,sm,ms,c,mx,my,mz,ix,iy,iz,lot,
     &     qe,qi,qme,qmi,nstep,b0x,b0y,b0z,az
     1   ,xst,xlt,yst,ylt,zst,zlt,vdre,vdri
     1            ,xdel,ydel,zdel,avdn,xa,rand1
      write(8)ions,lecs,maxptl,maxhlf,x,y,z,u,v,w,
     1       mcont,mmax,kcont,kmax,deltm,vithml,vethml,
     1       ntimee,me,mi,elosi,elose,anpo,
     1       xvne,xvni,nvcell,nvcell1,nvcel2
     1    , alhc,zz0,xx1,xx2,xx3,ep0,ep1,ep2,ep3,and0
      close(6)
      close(7)
      close(8)
      close(14)
      close(15)
      close(16)
      close(17)
      close(18)
      close(19)
      close(20)
      close(21)
      close(22)
      close(23)
      close(24)
      close(25)
      close(26)
      close(27)
c
c
c for steps with energy histries
c
      write (6,*) 'ntimee= ',ntimee
c
C  The user must decide what information is to be written out at
C  each timestep, what only occasionally, and what only at the end.
C  It may be wise to write out both COMMONS before stopping: they
C  can then be read in again for a continuation of the run.
  101 format(1h ,'nstep= ',i7)
 1000 continue
      stop
      end
c
c------------------------------------------------------------------------
      subroutine input
c
      parameter(nptl=5000000,mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz,
     1      mx4=mxx+4,my4=myy+4)
      parameter (mbeam=9300,mback=9300)
      real me,mi
      common /partls/x(nptl),y(nptl),z(nptl),u(nptl),v(nptl),w(nptl)
      common /fields/ex1(ngtl),ey1(ngtl),ez1(ngtl),bx1(ngtl),
     &by1(ngtl),bz1(ngtl)
     &,sm(27),q,qe,qi,qme,qmi,c
     &,ms(27),mx,my,mz,ix,iy,iz,lot,maxptl,maxhlf,ions,lecs,nstep
      common /veldis/ xvne,xvni,nvcell,nvcell1,nvcel2
c
c  Diagnosis for localized beam electrons and background electrons
c
      common /phas/ mcont(mbeam),mmax,kcont(mback),kmax
      common /para1/ b0x,b0y,b0z,deltm,vithml,vethml
      common /energ/ ntimee,me,mi,elosi,elose,anpo
      common /avec/ az(mx4,my4)                    
      common /jet/ xst,xlt,yst,ylt,zst,zlt,vdre,vdri
     1            ,xdel,ydel,zdel,avdn,xa,rand1
      common /tear/ alhc,zz0,xx1,xx2,xx3,ep0,ep1,ep2,ep3,and0
     1         ,ail
c     integer*8 iseed
c     real*8 dran
c
c
c  This may cause the problem
c     equivalence(ex1(1),ex(1,1,1)),(ey1(1),ey(1,1,1)),(ez1(1),ez(1,
c    &1,1)),(bx1(1),bx(1,1,1)),(by1(1),by(1,1,1)),(bz1(1),bz(1,1,1))
c
c maxptl=nptl   mxx=mx  myy=my  mzz=mz
c
      maxptl=nptl     
      maxhlf=maxptl/2
      mx=mxx 
      my=myy
      mz=mzz
C Strides for single-indexed field arrays:
      ix=1
      iy=mx
      iz=iy*my
      lot=iz*mz
c   Parameters for collisionless reconnection
c   L
      c=.5
C Miscellaneuos constants:
      pi=3.141592654
      pi2=2.0*pi
c
c     qe=-.0625
      qe=-.007
c     qe=-.01
c     qe=-.025
c     qi=.0625
      qi=-qe
c     me=.0625
      me=.15   
!     me=1.0  
c     me=0.015625
c  mi/me=64
c     mi=4.
c     mi=2.4      
c   mi/me=16
      mi=25.0*me
c  mi/me=128.
c     mi=8.
c amie=mi/me
      amie=mi/me
      qme=qe/me
      qmi=qi/mi
c attention the number should be based on simulations
      an0=50.0
c     an0=1.0
      wpe2=qe*qe*an0/me
      wpe=sqrt(wpe2)
c
      deltm=wpe
      wpi2=qi*qi*an0/mi
      wpi=sqrt(wpi2)
      ail=c/wpi
c for the same temperature for electron and ions
c  vithml/root square(mi)
c
      vethml4=0.05        
c     vethml4=0.025   
c     vithml4=0.046875/2./4.      
      vithml4=vethml4/sqrt(5.0)
c     vithml4=0.025*0.25          
      vethml=vethml4
      vithml=vithml4
      te=me*vethml*vethml
      ti=mi*vithml*vithml
c tie=ti/te
      tie=ti/te
      debye=vethml/wpe
      skind=c/wpe
      alhc = 0.5*ail
      nzz0 = 0.5*mz +1.0
      zz0 = nzz0
      xx1 = 6.25*alhc 
      xx2 = 15.6*alhc
      xx3 = mx        
      xx3 = xx3 -2.0
      ep0 = 0.00
      ep1 = 0.15
      ep2 = 0.05
      ep3 = (ep0 +ep2)/(xx3 -xx2)
c
      thet0=3.141*55.0/180.0
      b0x=0.00
c     b0x=0.20
      b0y=0.00
c attention
      b0z=0.3475
c  attention
c
c     b0z=sqrt(8.0*3.1417*an0*(ti+te))*c
      wce=qe*sqrt(b0x**2+b0y**2+b0z**2)/(me*c)
      wci=qi*sqrt(b0x**2+b0y**2+b0z**2)/(mi*c)
      ratiom=mi/me
      wecp=wce/wpe
      wicp=wci/wpi
c
      valfc=wci/wpi
      cveth=c/vethml
      beta=2.0/(cveth**2*wecp**2)
      rhoe=vethml/wce
      rhoi=vithml/wci
c
c  for velocity distribution
c
      amage=100./10.
      amagi=100./8.0
      xvne=amage/vethml
      xvni=amagi/vithml
      kpara=100
      kpara2=2*kpara
      nvcell=kpara2
      nvcell1=nvcell+1
      nvcel2=kpara
c
c List of Parameters 
c
      write(6,100) mx, my, mz, qe, qi, me, mi, qme, qmi, ratiom, c
      write(6,200) wpe, wpi, wce, wci, wecp, wicp, b0x, b0y, 
     1    b0z, beta
      write(6,201) vethml, vithml, valfc, cveth
      write(6,202) deltm, debye, skind, rhoe, rhoi
      write(6,203) te,ti,tie,amie
      write(6,204) ail,zz0,alhc
c
c    
C Our finite difference equations imply delta_t = delta_x =
C delta_y = delta_z = 1. So c must satisfy the Courant condition.
C The bx,by and bz arrays are really records of c*Bx, c*By,
C c*Bz: this makes for e <-----> b symmetry in Maxwell's equations.
C Otherwise, units are such that epsilon_0  is 1.0 and hence mu_0
C is 1/c**2. This means that for one electron per cell (see example
C of particle initialisation below) omega_p-squared is qe**2/me.
c
c------------------------------------------------------------------------
C For use in the boundary field calculation:
c     rs=(1.-c)/(1.+c)
c     tsq=.1666667
c     ps=(1.-tsq)*c/(1.+c)
c     os=.5*(1.+tsq)*c/(1.+c)
c
c-----------------------------------------------------------------------
C Data for smoothing: the currents fed into Maxwell's equations
C are smoothed by convolving with the sequence  .25, .5, .25 in
C each dimension. Generate the 27 weights ("sm") and index
C displacements ("ms"):
      n=1
      do 1 nz=-1,1
      do 1 ny=-1,1
      do 1 nx=-1,1
      sm(n)=.015625*(2-nx*nx)*(2-ny*ny)*(2-nz*nz)
      ms(n)=ix*nx+iy*ny+iz*nz
  1   n=n+1
      write(6,*) (sm(i),i=1,27)
      write(6,*) (ms(i),i=1,27)
c
c-----------------------------------------------------------------------
c DIstribution of electrons and ions
c
C In the particle arrays the ions are at the bottom, the electrons
C are stacked against the top, the total number not exceeding maxptl:
C The number of ions, "ions", need not be the same as the number of
C electrons, "lecs".
C The code treats unpaired electrons as having been initially
C dissociated from infinitely heavy ions which remain in situ.
C Initialise the particles: Place electrons in same locations as ions
C for zero initial net charge density. Keep particles 2 units away from
C the lower boundaries of the field domain, 3 units away from the upper
C boundaries. For instance, fill the interior uniformily:
c
c----- Ions ------------------------------------------------------------------
c
c           ions=0
c           do 80 k=1,mz-5
c           do 80 j=1,my-5
c           do 80 i=1,mx-5
c           ions=ions+1
c           x(ions)= 2.5+i+(dran(iseed)-0.5)
c           y(ions)= 2.5+j+(dran(iseed)-0.5)
c 80        z(ions)= 2.5+k+(dran(iseed)-0.5)
c
c---------- Electrons ------------------------------
c
C Put electrons in the same places:
c     lecs=ions
c     do 4 n=1,lecs
c     x(n+maxhlf)=x(n)
c     y(n+maxhlf)=y(n)
c 4   z(n+maxhlf)=z(n)
c
c
c
c-------------------------------------------------------------------
c  Initial velocities
c
C Initialise velocities: these should not exceed  c  in magnitude!
C For thermal distributions, add three or four random numbers for
C each component and scale.
C  Initialise random number generator:
c-------- Ions ---------------------------------------------
c     do 85 n=1,ions
c     u(n)=vithml4*(dran(iseed)+dran(iseed)+dran(iseed)
c    1             +dran(iseed)-2.0) 
c     w(n)=vithml4*(dran(iseed)+dran(iseed)+dran(iseed)
c    1             +dran(iseed)-2.0) 
c     u(maxhlf+n)=vethml4*(dran(iseed)+dran(iseed)+dran(iseed)
c    1                    +dran(iseed)-2.0) 
c     v(maxhlf+n)=vethml4*(dran(iseed)+dran(iseed)+dran(iseed)
c    1                    +dran(iseed)-2.0) 
c 85  w(maxhlf+n)=vethml4*(dran(iseed)+dran(iseed)+dran(iseed)
c    1                    +dran(iseed)-2.0) 
c
c------------------------------------------------------------------------
c
c Initializing the current sheet in  of x - z plane
c
c   therm vel of electron: vethml
c   therm vel of ions    : vithml
c  velocity of electrons: vdre
c  velocity of ions     : vdri
      vdri=2.0*vithml*vithml/(alhc*wci)
      vdre=-te/ti*vdri   
c attention   delt = 1.0 
c
      lx = mxx-5
      ly = myy-5
      lz = mzz-5
c dens: the number density in a cell
      denx=3.8
      deny=3.8
      denz=3.8
!     denx=3.83
!     deny=3.83
!     denz=3.83
!     denx=3.835
!     deny=3.835
!     denz=3.835
!     denx=3.745
!     deny=3.745
!     denz=3.745
      dens=denx*deny*denz
      dlx = 1.0/denx
      dly = 1.0/deny
      dlz = 1.0/denz
      inx = lx*(1.0-dlx)*denx
      iny = ly*(1.0-dly)*deny
      inz = lz*(1.0-dlz)*denz
      write(6,*) 'denx,deny,denz,dens,dlx,dly,dlz',
     *   denx,deny,denz,dens,dlx,dly,dlz
      lx1 = lx + inx
      ly1 = ly + iny
      lz1 = lz + inz
      dlx = float(lx)/float(lx1)
      dly = float(ly)/float(ly1)
      dlz = float(lz)/float(lz1)
      dens= float((lx1*ly1*lz1)/(lx*ly*lz))
      dlxh = 0.5*dlx
      dlyh = 0.5*dly
      dlzh = 0.5*dlz
      xst = 3.0 -dlxh
      yst = 3.0 -dlyh
      zst = 3.0 -dlzh
      densg=dens
      write(6,*)'after adjusted denx,deny,denz,dens,dlx,dly,dlz',
     1  denx,deny,denz,dens,dlx,dly,dlz
c first set up ions in the cloud
      ions=0
      lecs=0
      do k=1,lz1
      do j=1,ly1
      do i=1,lx1
      xx=xst+dlx*i
      yy=yst+dly*j
      zz=zst+dlz*k
      rand1=1.0/cosh((zz-zz0)/alhc)**2
      aran = rand()
      if(aran.le.rand1) then  
      ions=ions+1
C  Initializing the density and velocity
C  
      x(ions)=xx+dlx*(rand()-0.5)
      y(ions)=yy+dly*(rand()-0.5)
      z(ions)=zz+dlz*(rand()-0.5)
      u(ions)=vithml*(rand()+rand()
     1                +rand()+rand()+rand()-2.5)
c     u(ions)=vdri+vithml*sqrt(-2.0*log(rand()))
      v(ions)=vdri+vithml*(rand()+rand()+rand()
     1                +rand()+rand()-2.5)
      w(ions)=vithml*(rand()+rand()+rand()
     1                +rand()+rand()-2.5)
c 
      lecs=lecs+1
      x(lecs+maxhlf)=x(ions)
      y(lecs+maxhlf)=y(ions)
      z(lecs+maxhlf)=z(ions)
      u(lecs+maxhlf)=vethml*(rand()+rand()
     1               +rand()+rand()+rand()-2.5)
c     u(lecs+maxhlf)=vdre+vethml*sqrt(-2.0*log(rand()))
      v(lecs+maxhlf)=vdre+vethml*(rand()+rand()+rand()
     1               +rand()+rand()-2.5)
      w(lecs+maxhlf)=vethml*(rand()+rand()+rand()
     1               +rand()+rand()-2.5)
      endif
      enddo
      enddo
      enddo

c   attention
      mmax=1
      kmax=1
      mcont(1)=1
      kcont(1)=2

      write(6,*) 'ions=',ions,'  elct=',lecs 
c attention 
      xa=xst
      write(6,*)mmax, kmax
C Initialise the fields, in this case we initialize the electric field zero:
c     do 20 k=1,mz
c     do 20 j=1,my
c     do 20 i=1,mx
c     ex(i,j,k)=0.
c     ey(i,j,k)=0.
c  20 ez(i,j,k)=0.
C (Remember that bx,by,bz are really c*Bx, c*By, c*Bz.)
C Initial fields, both electric and magnetic, must be divergence-free.
C Part of the Earth's magnetic field would be ok. If the Earth is included
C in the field domain, its magnetic dipole field is readily established
C by maintaining a steady ring current in the Earth's core.
c      
c----------------
c attention
c   here   a**2=rad2, xc=xcn, yc=ycn, a=rad
C  Using a vector potential component  az  we initialise the
C  magnetic field to that of an axial current with Bennett profile.
C  This profile has density proportional to 1./(a**2+r**2)**2 and
C  creates a vector potential az proportional to log(a**2+r**2).
C  Here r**2 = (x-xc)**2 + (y-yc)**2 and  a  is the radius
C  outside of which there is as much current as inside. We choose:
c    we have already chosen
c     xc=35.
c     yc=35.
c     a=8.
C  The number density of ions and electrons (particles per cell) at
C  the center is as it would be if the current were uniform inside
C  a cylinder of radius  a  and zero outside. Assuming the total
C  number of pairs to be maxhlf and the length of the column mz-3,
C  this leads to a central density:
c     dns=float(maxhlf)/((mz-3.)*3.141593*a**2)
c     dns=float(mmax)/((mz-5.)*3.141593*rad2)
C  We choose the drift velocities of the ions and electrons:
c we have already chosen
c     vdri=.25*c
c     vdre=-.25*c
C  and hence a relative velocity:
      vrel=vdri-vdre
C  This creates a central current density qi*dns*vrel and a
C  vector potential az=az0*log(a**2+r**2) where
c     az0=-qi*dns*vrel*a**2/(4.*c)
c     az0=-qi*dns*vrel*rad2/(4.*c)
C  We calculate az first only in one octant:
c  attention
c attention
c     do 28 j=1,myy
c     do 28 i=1,mxx
c     azt(i,j)=az(i,j)
c  28 continue
C  Now we derive bx and by by curling az. Thus we ensure that our
C  initial magnetic field is divergenceless. We also create an axial
C  component bz (overwriting the az array!) of value:
c-----------------------------------------------------------------------
c
c     vdrirt=-vdre*float(mmax)/float(kmax)
c
c  Return current
c
c     do 211 i=1,kmax
c     k=kcont(i)
c     w(k)=w(k)+vdrirt
c 211 continue
c     write(6,555) vdre, vethmb, kmax, mmax, beamr 
      write(6,555) vdre, vdri, kmax, mmax, beamr, vdrirt
  555 format(/1h ,'vdre= ',f9.6,1x,'vdri= ',f9.6,1x,'kmax  =',
     1       i6,1x,'mmax= ',i6,1x,'be/tl= ',f6.4//
     2       1h ,'vdriftrt= ', f9.6)
c    
c
  100 format(1h , 2x,'mx= ',i5,2x,'my= ',i5,2x,'mz= ',i5//5x,
     1   'qe= ',f9.4,2x,'qi= ',f9.4,2x,'me= ',f9.4,2x,'mi= ',f9.4//5x
     2   'qme= ',f9.4,2x,'qmi= ',f9.4,2x,'mi/me=',f8.3,2x,'c= ',f6.4)
c
  200 format(/1h ,2x,'wpe= ',f7.5,2x,'wpi= ',f8.6,2x,'wce= ',f8.5,2x,
     1        'wci= ',f9.6,2x,/1h ,2x,'wecp= ',f9.6,2x,'wicp= ',f9.6,
     2        2x,/1h ,'b0x= ',f6.3,2x,'b0y= ',f6.3,2x,
     3         'b0z= ',f6.3,2x,'beta= ',e10.4)
c
  201 format(/1h ,2x,'vet= ',f7.5,2x,'vit= ',f8.6,2x,'vac= ',f8.5,2x,
     1         'cvet= 'f8.4)
c
  202 format(/1h ,2x,'deltm= ',f6.4,2x,'debye= ',f6.4,2x,'skind= ',f7.3,
     1  2x,'re= ',f8.5,2x,'ri= ',f8.5)
  203 format(/1h ,2x,'te= ',e12.4,2x,'ti= ',e12.4,2x,'ti/te= ',e12.4,
     1         2x,'mi/me= ',f7.2)
      write(6,204) ail,zz0,alhc
  204 format(/1h ,2x,'ionin= ',e12.4,2x,'zz0= ',e12.4,2x,'alhc= ',e12.4)
 2489 continue
      return
      end
C--------
      subroutine copylayr(bx,by,bz,mx,my,mz,lt,ls,nt,ns)      
      dimension bx(mx,my,mz),by(mx,my,mz),bz(mx,my,mz)
      do 2 j=1,my
      do 2 k=1,mz
      bx(lt,j,k)=bx(ls,j,k)
      by(lt,j,k)=by(ls,j,k)
   2  bz(lt,j,k)=bz(ls,j,k)
      do 1 i=1,mx
      do 1 k=1,mz
      bx(i,nt,k)=bx(i,ns,k)
      by(i,nt,k)=by(i,ns,k)
   1  bz(i,nt,k)=bz(i,ns,k)
      return
      end
C--------
      subroutine clearlay(ex,ey,ez,mx,my,mz,m,l)     
      dimension ex(mx,my,mz),ey(mx,my,mz),ez(mx,my,mz)
      do 1 j=1,my
      do 1 k=1,mz
      ex(m,j,k)=0.
      ey(m,j,k)=0.
   1  ez(m,j,k)=0.
      do 2 i=1,mx
      do 2 k=1,mz
      ex(i,l,k)=0.
      ey(i,l,k)=0.
   2  ez(i,l,k)=0.
      return
      end
C--------
      subroutine addlayer(ex,ey,ez,mx,my,mz,mt,ms,lt,ls)        
      dimension ex(mx,my,mz),ey(mx,my,mz),ez(mx,my,mz)
      do 1 j=1,my
      do 1 k=1,mz
      ex(mt,j,k)=ex(mt,j,k)+ex(ms,j,k)
      ey(mt,j,k)=ey(mt,j,k)+ey(ms,j,k)
   1  ez(mt,j,k)=ez(mt,j,k)+ez(ms,j,k)
      do 2 i=1,mx
      do 2 k=1,mz
      ex(i,lt,k)=ex(i,lt,k)+ex(i,ls,k)
      ey(i,lt,k)=ey(i,lt,k)+ey(i,ls,k)
   2  ez(i,lt,k)=ez(i,lt,k)+ez(i,ls,k)
      return
      end
c
c----------------------------------------------
      subroutine densfv(t)
      parameter(nptl=5000000,mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      common /partls/x(nptl),y(nptl),z(nptl),u(nptl),v(nptl),w(nptl)
      common /fields/ex1(ngtl),ey1(ngtl),ez1(ngtl),bx1(ngtl),
     &by1(ngtl),bz1(ngtl)
     &,sm(27),q,qe,qi,qme,qmi,c
     &,ms(27),mx,my,mz,ix,iy,iz,lot,maxptl,maxhlf,ions,lecs,nstep
      dimension dense(mxx,myy), densi(mxx,myy),
     1   fvlex(mxx,myy), fvley(mxx,myy), fvlix(mxx,myy), fvliy(mxx,myy)
c
      do 1 i=1, mxx
      do 1 j=1, myy
      dense(i,j)=0.0
      densi(i,j)=0.0
      fvlex(i,j)=0.0
      fvley(i,j)=0.0
      fvlix(i,j)=0.0
      fvliy(i,j)=0.0
    1 continue
c
c  The sheet between zmin and zmax
c
      zmin=0.5*mz-1
      zmax=0.5*mz+1
c
c  The sheet between ymin and ymax
c
      ymin=my*0.5-1
      ymax=my*0.5+1
      azmax=myy
c
c   Attention
c
c       Actually Density and flux on the x - z plane
c
      do 2 n=1,ions
c     if(z(n).lt.zmin.or.z(n).gt.zmax) go to 2
      if(y(n).lt.ymin.or.y(n).gt.zmax) go to 2
      if(z(n).gt.azmax) go to 2
      i=ifix(x(n))
c     j=ifix(y(n)) 
      j=ifix(z(n)) 
      densi(i,j)=densi(i,j)+1.0
      fvlix(i,j)=fvlix(i,j)+u(n)
c     fvliy(i,j)=fvliy(i,j)+v(n)
      fvliy(i,j)=fvliy(i,j)+w(n)
    2 continue
      do 3 n=maxhlf+1,lecs+maxhlf
c     if(z(n).lt.zmin.or.z(n).gt.zmax) go to 3
      if(y(n).lt.ymin.or.y(n).gt.zmax) go to 3
      if(z(n).gt.azmax) go to 3
      i=ifix(x(n))
c     j=ifix(y(n)) 
      j=ifix(z(n)) 
      dense(i,j)=dense(i,j)+1.0
      fvlex(i,j)=fvlex(i,j)+u(n)
c     fvley(i,j)=fvley(i,j)+v(n)
      fvley(i,j)=fvley(i,j)+w(n)
    3 continue
c     do 10 i=1,mxx
c     do 10 j=1,myy
c     if(densi(i,j).lt.0.000000001) densi(i,j)=1.0
c     if(dense(i,j).lt.0.000000001) dense(i,j)=1.0
c     fvlix(i,j)=fvlix(i,j)/densi(i,j)
c     fvliy(i,j)=fvliy(i,j)/densi(i,j)
c     fvlex(i,j)=fvlex(i,j)/dense(i,j)
c     fvley(i,j)=fvley(i,j)/dense(i,j)
c  10 continue
      densi(1,1)=t
      write(25) dense,fvlex,fvley,densi,fvlix,fvliy
      return
 2489 write(6,*) 'bad buffer out at unit 25'
      end
c
c---------------------
      subroutine diffus(t)
      parameter(nptl=5000000,mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      parameter (mbeam=9300,mback=9300)
      real me,mi
      common /partls/x(nptl),y(nptl),z(nptl),u(nptl),v(nptl),w(nptl)
      common /fields/ex1(ngtl),ey1(ngtl),ez1(ngtl),bx1(ngtl),
     &by1(ngtl),bz1(ngtl)
     &,sm(27),q,qe,qi,qme,qmi,c
     &,ms(27),mx,my,mz,ix,iy,iz,lot,maxptl,maxhlf,ions,lecs,nstep
      common /veldis/ xvne,xvni,nvcell,nvcell1,nvcel2
c
c  Diagnosis for localized beam electrons and background electrons
c
      common /phas/ mcont(mbeam),mmax,kcont(mback),kmax
      dimension xeb(mbeam), yeb(mbeam), zeb(mbeam), 
     1          vxb(mbeam), vyb(mbeam), vzb(mbeam)
      dimension xeg(mback), yeg(mback), zeg(mback),
     1          vxg(mback), vyg(mback), vzg(mback)
c
      do 100 m=1, mmax
!     i=mcont(m)
      xeb(m)=x(i)
      yeb(m)=y(i)
      zeb(m)=z(i)
      vxb(m)=u(i)
      vyb(m)=v(i)
      vzb(m)=w(i)
  100 continue
c
      do 101 k=1, kmax
c     i=kcont(k)
      xeg(k)=x(i)
      yeg(k)=y(i)
      zeg(k)=z(i)
      vxg(k)=u(i)
      vyg(k)=v(i)
      vzg(k)=w(i)
  101 continue
c
      vzb(mbeam)=t
      vzb(mbeam-1)=mmax
      vzb(mbeam-2)=kmax
c
c
c  Buffer out data for beam and background electrons
c  
      write(23) xeb,yeb,zeb,vxb,vyb,vzb
c     write(24) xeg,yeg,zeg,vxg,vyg,vzg
 2489 continue
      return
      end
c
c
c---------------------
      subroutine energy(t)
      parameter(nptl=5000000,mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      parameter(mbeam=9300,mback=9300)
      real mi,me
      common /partls/x(nptl),y(nptl),z(nptl),u(nptl),v(nptl),w(nptl)
      common /fields/ex1(ngtl),ey1(ngtl),ez1(ngtl),bx1(ngtl),
     &by1(ngtl),bz1(ngtl)
     &,sm(27),q,qe,qi,qme,qmi,c
     &,ms(27),mx,my,mz,ix,iy,iz,lot,maxptl,maxhlf,ions,lecs,nstep
c
c  Diagnosis for localized beam electrons and background electrons
c
      common /phas/ mcont(mbeam),mmax,kcont(mback),kmax
      common /para1/ b0x,b0y,b0z,deltm,vithml,vethml
      common /energ/ ntimee,me,mi,elosi,elose,anpo
      dimension ex(mxx,myy,mzz),ey(mxx,myy,mzz),ez(mxx,myy,mzz),
     &bx(mxx,myy,mzz),by(mxx,myy,mzz),bz(mxx,myy,mzz)
      dimension hist(5,3) 
      equivalence(ex1(1),ex(1,1,1)),(ey1(1),ey(1,1,1)),(ez1(1),ez(1,
     &1,1)),(bx1(1),bx(1,1,1)),(by1(1),by(1,1,1)),(bz1(1),bz(1,1,1))
c
      ntimee=ntimee+1
c
      akine=0.0
      akini=0.0
      akpab=0.0
      akpeb=0.0
      drifb=0.0
      efx=0.0
      efy=0.0
      efz=0.0
      bfx=0.0
      bfy=0.0
      bfz=0.0
      anxp=0.0
      anxn=0.0
      anyp=0.0
      anyn=0.0
c
      do 1 i=1, ions
c     j=maxptl-lecs+i
      akini=akini+u(i)*u(i)+v(i)*v(i)+w(i)*w(i)
    1 continue
      do 107 j=1+maxhlf,lecs+maxhlf
      akine=akine+u(j)*u(j)+v(j)*v(j)+w(j)*w(j)
  107 continue
      akini=0.5*qi*akini/qmi
      akine=0.5*qe*akine/qme
c
c  Diagnosis for localized beam electrons and background electrons
c
      do 100 m=1, mmax
!     i=mcont(m)
      drifb=drifb+w(i)
      akpeb=akpeb+u(i)*u(i)+v(i)*v(i)
      akpab=akpab+w(i)*w(i)
  100 continue
      drifb=drifb/float(mmax)
      akpab=0.5*qe*akpab/qme
      akpeb=0.5*qe*akpeb/qme
c
      do 2 i=3, mx-2
      do 2 j=3, my-2 
      do 2 k=2, mz-2
      efx=efx+ex(i,j,k)*ex(i,j,k)
      efy=efy+ey(i,j,k)*ey(i,j,k)
      efz=efz+ez(i,j,k)*ez(i,j,k)
      bfx=bfx+(bx(i,j,k)-b0x)*(bx(i,j,k)-b0x)
      bfy=bfy+(by(i,j,k)-b0y)*(by(i,j,k)-b0y)
      bfz=bfz+(bz(i,j,k)-b0z)*(bz(i,j,k)-b0z)
    2 continue
c     bfx=bfx/c/c
c     bfr=bfy/c/c
c     bfz=bfz/c/c
c
      do 3 j=3,my-2
      do 3 k=3,mz-3
      anxp=anxp-0.125*(ey(3,j,k)   +ey(4,j,k+1)
     1                +ey(4,j,k)   +ey(3,j,k+1))
     1               *(bz(3,j,k)   +bz(3,j,k+1))
     2         +0.125*(ey(mx-3,j,k)+ey(mx-2,j,k+1)
     2                +ey(mx-2,j,k)+ey(mx-3,j,k+1))
     3               *(bz(mx-3,j,k)+bz(mx-3,j,k+1))
    3 continue
c
      do 4 j=3,my-3
      do 4 k=3,mz-2
      anxn=anxn-0.125*(ez(3,j,k)   +ez(4,j+1,k)
     1                +ez(4,j,k)   +ez(3,j+1,k))
     1               *(by(3,j,k)   +by(3,j+1,k))
     2         +0.125*(ez(mx-3,j,k)+ez(mx-2,j+1,k)
     2                +ez(mx-2,j,k)+ez(mx-3,j+1,k))
     3               *(by(mx-3,j,k)+by(mx-3,j+1,k))
    4 continue
c
      do 5 i=3,mx-3
      do 5 k=3,mz-2
      anyp=anyp-0.125*(ez(i,3,k)   +ez(i+1,4,k)
     1                +ez(i,4,k)   +ez(i+1,3,k))
     1               *(bx(i,3,k)   +bx(i+1,3,k))
     2         +0.125*(ez(i,my-3,k)+ez(i+1,my-2,k)
     2                +ez(i,my-2,k)+ez(i+1,my-3,k))
     3               *(bx(i,my-3,k)+bx(i+1,my-3,k))
    5 continue
c
      do 6 i=3,mx-2
      do 6 k=3,mz-3
      anyn=anyn-0.125*(ex(i,3,k)   +ex(i,4,k+1)
     1                +ex(i,4,k)   +ex(i,3,k+1))
     1               *(bz(i,3,k)   +bz(i,3,k+1))
     2         +0.125*(ex(i,my-3,k)+ex(i,my-2,k+1)
     2                +ex(i,my-2,k)+ex(i,my-3,k+1))
     3               *(bz(i,my-3,k)+bz(i,my-3,k+1))
    6 continue
c
      hist(1,1)=akini
      hist(1,2)=akine
      hist(1,3)=akpab
      hist(2,1)=akpeb
      hist(2,2)=drifb
      hist(2,3)=t
c
      hist(3,1)=efx
      hist(3,2)=efy
      hist(3,3)=efz
      hist(4,1)=bfx
      hist(4,2)=bfy
      hist(4,3)=bfz
      hist(5,1)=elosi
      hist(5,2)=elose
      hist(5,3)=anxp-anxn+anyp-anyn
c     anpo     =anpo +anxp-anxn+anyp-anyn
c     hist(5,3)=anpo
c
      write(26) hist
      write(6,*) 'ntimee= ',ntimee
c     write(6,*) 'hist(1,1)= ',hist(1,1),'hist(2,2)= ',hist(2,2),
c    1    'hist(3,2)= ',hist(3,2),'hist(4,3)= ',hist(4,3)
      write(6,*) ((hist(m,n),m=1,5),n=1,3)
      return
 1000 stop
      end
c
c ----------
c 
      subroutine veldist(t)
      parameter(nptl=5000000,mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      common /partls/x(nptl),y(nptl),z(nptl),u(nptl),v(nptl),w(nptl)
      common /fields/ex1(ngtl),ey1(ngtl),ez1(ngtl),bx1(ngtl),
     &by1(ngtl),bz1(ngtl)
     &,sm(27),q,qe,qi,qme,qmi,c
     &,ms(27),mx,my,mz,ix,iy,iz,lot,maxptl,maxhlf,ions,lecs,nstep
      common /veldis/ xvne,xvni,nvcell,nvcell1,nvcel2
      dimension kkex(65),kkey(65),kkez(65),kkix(65),kkiy(65),kkiz(65)
      dimension suexy(200,200), suexz(200,200), suixy(200,200),
     &          suixz(200,200)
      integer cvmgm
      external cvmgm
c
      nbe=lecs/64
      nbi=ions/64
      laste=lecs - 64*nbe
      if(laste.ne.0) nbe=nbe+1
      if(laste.eq.0) laste=64
      lasti=ions - 64*nbi
      if(lasti.ne.0) nbi=nbi+1
      if(lasti.eq.0) lasti=64
c     write(6,*) 'Begining of subroutine veldist'
      do 5057 j=1,nvcell
      do 5057 i=1,nvcell
      suexy(i,j)=0.0
      suexz(i,j)=0.0
      suixy(i,j)=0.0
      suixz(i,j)=0.0
 5057 continue
c
c  Ions
c
      i2=0
      i3=64
      do 20 kk=1,nbi
      i1=i2+1
      if(kk.eq.nbi) i3=lasti
      i2=i2+i3
      do 30 i=i1,i2
      k=i-i1+1
      kkix(k)=u(i)*xvni+nvcel2
      kkiy(k)=v(i)*xvni+nvcel2
      kkiz(k)=w(i)*xvni+nvcel2
      kkix(k)=cvmgm(kkix(k),nvcell1,kkix(k)-nvcell)
      kkix(k)=cvmgm(nvcell1,kkix(k),kkix(k)-1)
      kkiy(k)=cvmgm(kkiy(k),nvcell1,kkiy(k)-nvcell)
      kkiy(k)=cvmgm(nvcell1,kkiy(k),kkiy(k)-1)      
      kkiz(k)=cvmgm(kkiz(k),nvcell1,kkiz(k)-nvcell)
   30 kkiz(k)=cvmgm(nvcell1,kkiz(k),kkiz(k)-1)
      do 20 k=1,i3
c     if(kkix(k).ge.200) kkix(k)=200
c     if(kkiy(k).ge.200) kkiy(k)=200
c     if(kkiz(k).ge.200) kkiz(k)=200
      suixy(kkix(k),kkiy(k))=suixy(kkix(k),kkiy(k))+1
      suixz(kkix(k),kkiz(k))=suixz(kkix(k),kkiz(k))+1
   20 continue
c
c     write(6,*) 'After ions'
c
c  Electrons
c
c     i2=maxptl - lecs
      i2=maxhlf+1
      i3=64
      do 21 kk=1,nbe
      i1=i2+1
      if(kk.eq.nbe) i3=laste
      i2=i2+i3
      do 32 i=i1,i2
      k=i-i1+1
      kkex(k)=u(i)*xvne+nvcel2
      kkey(k)=v(i)*xvne+nvcel2
      kkez(k)=w(i)*xvne+nvcel2
      kkex(k)=cvmgm(kkex(k),nvcell1,kkex(k)-nvcell)
      kkex(k)=cvmgm(nvcell1,kkex(k),kkex(k)-1)
      kkey(k)=cvmgm(kkey(k),nvcell1,kkey(k)-nvcell)
      kkey(k)=cvmgm(nvcell1,kkey(k),kkey(k)-1)      
      kkez(k)=cvmgm(kkez(k),nvcell1,kkez(k)-nvcell)
   32 kkez(k)=cvmgm(nvcell1,kkez(k),kkez(k)-1)
      do 21 k=1,i3
c     if(kkex(k).ge.200) kkex(k)=200
c     if(kkey(k).ge.200) kkey(k)=200
c     if(kkez(k).ge.200) kkez(k)=200
      suexy(kkex(k),kkey(k))=suexy(kkex(k),kkey(k))+1
      suexz(kkex(k),kkez(k))=suexz(kkex(k),kkez(k))+1
   21 continue
c for time
      suexy(1,1)=t
c
c     write(6,*) 'After electrons'
      write(21) suexy,suexz
c
      write(22) suixy,suixz
 1000 continue
      return
      end
c
c
c
C -----------
      subroutine surface(bx,by,bz,ex,ey,ez,ix,iy,iz,mx,my,mz,m00,c)
      parameter(mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      dimension bx(ngtl),by(ngtl),bz(ngtl),ex(ngtl),
     1             ey(ngtl),ez(ngtl)
      rs=2.*c/(1.+c)
      s=.4142136
      os=.5*(1.-s)*rs
      do 3 m=m00+iz*(mz-1),m00+iz*(mz-1)+iy*(my-2),iy
      do 1 n=m,m+ix*(mx-2),ix
   1  bz(n)=bz(n)+.5*c*(ex(n+iy)-ex(n)-ey(n+ix)+ey(n))
cdir$ ivdep
      do 3 n=m+ix,m+ix*(mx-2),ix
   3  bx(n)=bx(n)+rs*(bx(n-iz)-bx(n)+s*(bz(n)-bz(n-ix)))-os*(
     &ez(n+iy)-ez(n))-(os-c)*(ez(n+iy-iz)-ez(n-iz))-c*(ey(n)-ey(n-iz))
      do 2 m=m00+iz*(mz-1),m00+iz*(mz-1)+ix*(mx-2),ix
cdir$ ivdep
      do 4 n=m+iy,m+iy*(my-2),iy
   4  by(n)=by(n)+rs*(by(n-iz)-by(n)+s*(bz(n)-bz(n-iy)))+os*(
     &ez(n+ix)-ez(n))+(os-c)*(ez(n+ix-iz)-ez(n-iz))+c*(ex(n)-ex(n-iz))
      do 2 n=m,m+iy*(my-2),iy
   2  bz(n)=bz(n)+.5*c*(ex(n+iy)-ex(n)-ey(n+ix)+ey(n))
      return
      end
C-------------
      subroutine preledge(bx,by,bz,ex,ey,ez,ix,iy,iz,mx,my,mz,m,c)
      parameter(mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      dimension bx(ngtl),by(ngtl),bz(ngtl),ex(ngtl),
     1             ey(ngtl),ez(ngtl)
      s=.4142136
      t=c/(2.*c+1.+s)
cdir$ ivdep
      do 1 n=m+iy*(my-1)+iz*(mz-1)+ix,m+iy*(my-1)+iz*(mz-1)+ix*(mx-2),ix
   1  bx(n)=bx(n-iy-iz)+(1.-4.*t)*bx(n)+(1.-2.*t)*(bx(n-iy)+bx(n-iz))
     &+s*t*(by(n)-by(n-ix)+by(n-iz)-by(n-ix-iz)
     &+bz(n)-bz(n-ix)+bz(n-iy)-bz(n-ix-iy))
      r=4./(2.*c+2.+s)
cdir$ ivdep
      do 2 n=m+iy*(my-1),m+iy*(my-1)+ix*(mx-2),ix
   2  bz(n)=(1.-c*r)*(bz(n)+bz(n+iz))+bz(n-iy)+bz(n+iz-iy)-r*(by(n)+c*
     &((1.-s)*(ez(n+ix)-ez(n))+(1.+s)*.25*(ey(n+ix)-ey(n)+ey(n+iz+ix)
     &-ey(n+iz)+ey(n+ix-iy)-ey(n-iy)+ey(n+iz+ix-iy)-ey(n+iz-iy))))
cdir$ ivdep
      do 3 n=m+iz*(mz-1),m+iz*(mz-1)+ix*(mx-2),ix
   3  by(n)=(1.-c*r)*(by(n)+by(n+iy))+by(n-iz)+by(n+iy-iz)-r*(bz(n)-c*
     &((1.-s)*(ey(n+ix)-ey(n))+(1.+s)*.25*(ez(n+ix)-ez(n)+ez(n+iy+ix)
     &-ez(n+iy)+ez(n+ix-iz)-ez(n-iz)+ez(n+iy+ix-iz)-ez(n+iy-iz))))
      p=(1.+c)*2./(1.+2.*c*(1.+c*s))
      q=c*s*2./(1.+2.*c*(1.+c*s))
cdir$ ivdep
      do 4 n=m+iy*(my-1)+iz*(mz-1),m+ix*(mx-2)+iy*(my-1)+iz*(mz-1),ix
      temp=bz(n)-.5*c*(1.-s)*(ey(n+ix)-ey(n)+ey(n+ix-iy)-ey(n-iy))
      bz(n)=bz(n-iy)-bz(n)+p*temp+q*by(n)
   4  by(n)=by(n-iz)-by(n)+p*by(n)+q*temp
      return
      end
C-------------
      subroutine postedge(bx,by,bz,ex,ey,ez,ix,iy,iz,mx,my,mz,m,c)
      parameter(mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      dimension bx(ngtl),by(ngtl),bz(ngtl),ex(ngtl),
     1             ey(ngtl),ez(ngtl)
      s=.4142136
      p=(1.+c)*2./(1.+2.*c*(1.+c*s))
      q=c*s*2./(1.+2.*c*(1.+c*s))
cdir$ ivdep
      do 4 n=m+iy*(my-1)+iz*(mz-1),m+ix*(mx-2)+iy*(my-1)+iz*(mz-1),ix
      temp=by(n-iz)-.5*c*(1.-s)*(ez(n+ix)-ez(n)+ez(n+ix-iz)-ez(n-iz))
      bz(n)=bz(n-iy)+bz(n)-q*temp-p*bz(n-iy)
   4  by(n)=by(n-iz)+by(n)-q*bz(n-iy)-p*temp
      t=c/(2.*c+1.+s)
cdir$ ivdep
      do 1 n=m+iy*(my-1)+iz*(mz-1)+ix,m+iy*(my-1)+iz*(mz-1)+ix*(mx-2),ix
   1  bx(n)=bx(n)-(1.-4.*t)*bx(n-iy-iz)-(1.-2.*t)*(bx(n-iy)+bx(n-iz))
     &+s*t*(by(n)-by(n-ix)+by(n-iz)-by(n-ix-iz)
     &+bz(n)-bz(n-ix)+bz(n-iy)-bz(n-ix-iy))
      r=4./(2.*c+2.+s)
cdir$ ivdep
      do 2 n=m+iy*(my-1),m+iy*(my-1)+ix*(mx-2),ix
   2  bz(n)=bz(n)-bz(n+iz)-(1.-c*r)*(bz(n-iy)+bz(n+iz-iy))+r*by(n)
cdir$ ivdep
      do 3 n=m+iz*(mz-1),m+iz*(mz-1)+ix*(mx-2),ix
   3  by(n)=by(n)-by(n+iy)-(1.-c*r)*(by(n-iz)+by(n+iy-iz))+r*bz(n)
      return
      end
C--------
      subroutine mover(n1,n2,qm)
      parameter(nptl=5000000,mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      common /partls/x(nptl),y(nptl),z(nptl),u(nptl),v(nptl),w(nptl)
C (Field components are treated as single-indexed in this subroutine)
      common /fields/ex(ngtl),ey(ngtl),ez(ngtl),bx(ngtl),
     &by(ngtl),bz(ngtl)
     &,sm(27),q,qe,qi,qme,qmi,c
     &,ms(27),mx,my,mz,ix,iy,iz,lot,maxptl,maxhlf,ions,lecs,nstep
c     write(6,*) 'begining of mover','c= ',c
      do 1 n=n1,n2
C Cell index & displacement in cell:
      i=x(n)
      dx=x(n)-i
      j=y(n)
      dy=y(n)-j
      k=z(n)
      dz=z(n)-k
      l=i+iy*(j-1)+iz*(k-1)
C Field interpolations are tri-linear (linear in x times linear in y
C times linear in z). This amounts to the 3-D generalisation of "area
C weighting". A modification of the simple linear interpolation formula
C         f(i+dx) = f(i) + dx * (f(i+1)-f(i))
C is needed since fields are recorded at half-integer locations in certain
C dimensions: see comments and illustration with the Maxwell part of this
C code. One then has first to interpolate from "midpoints" to "gridpoints"
C by averaging neighbors. Then one proceeds with normal interpolation.
C Combining these two steps leads to:
C   f at location i+dx  = half of f(i)+f(i-1) + dx*(f(i+1)-f(i-1))
C where now f(i) means f at location i+1/2. The halving is absorbed
C in the final scaling.
C   E-component interpolations:
      f=ex(l)+ex(l-ix)+dx*(ex(l+ix)-ex(l-ix))
      f=f+dy*(ex(l+iy)+ex(l-ix+iy)+dx*(ex(l+ix+iy)-ex(l-ix+iy))-f)
      g=ex(l+iz)+ex(l-ix+iz)+dx*(ex(l+ix+iz)-ex(l-ix+iz))
      g=g+dy*
     & (ex(l+iy+iz)+ex(l-ix+iy+iz)+dx*(ex(l+ix+iy+iz)-ex(l-ix+iy+iz))-g)
      ex0=(f+dz*(g-f))*(.25*qm)
C   ------------
      f=ey(l)+ey(l-iy)+dy*(ey(l+iy)-ey(l-iy))
      f=f+dz*(ey(l+iz)+ey(l-iy+iz)+dy*(ey(l+iy+iz)-ey(l-iy+iz))-f)
      g=ey(l+ix)+ey(l-iy+ix)+dy*(ey(l+iy+ix)-ey(l-iy+ix))
      g=g+dz*
     & (ey(l+iz+ix)+ey(l-iy+iz+ix)+dy*(ey(l+iy+iz+ix)-ey(l-iy+iz+ix))-g)
      ey0=(f+dx*(g-f))*(.25*qm)
C   ------------
      f=ez(l)+ez(l-iz)+dz*(ez(l+iz)-ez(l-iz))
      f=f+dx*(ez(l+ix)+ez(l-iz+ix)+dz*(ez(l+iz+ix)-ez(l-iz+ix))-f)
      g=ez(l+iy)+ez(l-iz+iy)+dz*(ez(l+iz+iy)-ez(l-iz+iy))
      g=g+dx*
     & (ez(l+ix+iy)+ez(l-iz+ix+iy)+dz*(ez(l+iz+ix+iy)-ez(l-iz+ix+iy))-g)
      ez0=(f+dy*(g-f))*(.25*qm)
C   ---------
C   B-component interpolations:
      f=bx(l-iy)+bx(l-iy-iz)+dz*(bx(l-iy+iz)-bx(l-iy-iz))
      f=bx(l)+bx(l-iz)+dz*(bx(l+iz)-bx(l-iz))+f+dy*
     & (bx(l+iy)+bx(l+iy-iz)+dz*(bx(l+iy+iz)-bx(l+iy-iz))-f)
      g=bx(l+ix-iy)+bx(l+ix-iy-iz)+dz*(bx(l+ix-iy+iz)-bx(l+ix-iy-iz))
      g=bx(l+ix)+bx(l+ix-iz)+dz*(bx(l+ix+iz)-bx(l+ix-iz))+g+dy*
     & (bx(l+ix+iy)+bx(l+ix+iy-iz)+dz*(bx(l+ix+iy+iz)-bx(l+ix+iy-iz))-g)
      bx0=(f+dx*(g-f))*(.125*qm/c)
C   ------------
      f=by(l-iz)+by(l-iz-ix)+dx*(by(l-iz+ix)-by(l-iz-ix))
      f=by(l)+by(l-ix)+dx*(by(l+ix)-by(l-ix))+f+dz*
     & (by(l+iz)+by(l+iz-ix)+dx*(by(l+iz+ix)-by(l+iz-ix))-f)
      g=by(l+iy-iz)+by(l+iy-iz-ix)+dx*(by(l+iy-iz+ix)-by(l+iy-iz-ix))
      g=by(l+iy)+by(l+iy-ix)+dx*(by(l+iy+ix)-by(l+iy-ix))+g+dz*
     & (by(l+iy+iz)+by(l+iy+iz-ix)+dx*(by(l+iy+iz+ix)-by(l+iy+iz-ix))-g)
      by0=(f+dy*(g-f))*(.125*qm/c)
C   ------------
      f=bz(l-ix)+bz(l-ix-iy)+dy*(bz(l-ix+iy)-bz(l-ix-iy))
      f=bz(l)+bz(l-iy)+dy*(bz(l+iy)-bz(l-iy))+f+dx*
     & (bz(l+ix)+bz(l+ix-iy)+dy*(bz(l+ix+iy)-bz(l+ix-iy))-f)
      g=bz(l+iz-ix)+bz(l+iz-ix-iy)+dy*(bz(l+iz-ix+iy)-bz(l+iz-ix-iy))
      g=bz(l+iz)+bz(l+iz-iy)+dy*(bz(l+iz+iy)-bz(l+iz-iy))+g+dx*
     & (bz(l+iz+ix)+bz(l+iz+ix-iy)+dy*(bz(l+iz+ix+iy)-bz(l+iz+ix-iy))-g)
      bz0=(f+dz*(g-f))*(.125*qm/c)
C   ---------
C   First half electric acceleration, with relativity's gamma:
c     write(6,*) 'Begining of electron acceleration'
c
      g=c/sqrt(c**2-u(n)**2-v(n)**2-w(n)**2)
c     write(6,*) 'g= ',g
      u0=g*u(n)+ex0
      v0=g*v(n)+ey0
      w0=g*w(n)+ez0
C   First half magnetic rotation, with relativity's gamma:
      g=c/sqrt(c**2+u0**2+v0**2+w0**2)
      bx0=g*bx0
      by0=g*by0
      bz0=g*bz0
      f=2./(1.+bx0*bx0+by0*by0+bz0*bz0)
      u1=(u0+v0*bz0-w0*by0)*f
      v1=(v0+w0*bx0-u0*bz0)*f
      w1=(w0+u0*by0-v0*bx0)*f
C   Second half mag. rot'n & el. acc'n:
      u0=u0+v1*bz0-w1*by0+ex0
      v0=v0+w1*bx0-u1*bz0+ey0
      w0=w0+u1*by0-v1*bx0+ez0
C   Relativity's gamma:
      g=c/sqrt(c**2+u0**2+v0**2+w0**2)
      u(n)=g*u0
      v(n)=g*v0
      w(n)=g*w0
C   Position advance:
      x(n)=x(n)+u(n)
      y(n)=y(n)+v(n)
  1   z(n)=z(n)+w(n)
      return
      end
C ------------
      subroutine conduct(ix,iy,iz,mx,my,mz,ey,ez,s)
c     dimension ey(1),ez(1)
      parameter(mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      dimension ey(ngtl),ez(ngtl)
C  field arrays are treated as one-dimensional arrays in this subroutine
C  Current density components jy=s*ey, jz=s*ez are spread laterally in the
C  proportion .25,.5,.25 required by our smoothing convention for all field
C  sources.
      do 1 i=1+2*ix,1+ix*(mx-3),ix*(mx-5)
      do 1 j=i+2*iy,i+iy*(my-4),iy
      do 1 k=j+2*iz,j+iz*(mz-4),iz
      temp=.25*s*ey(k)
      ey(k-ix)=ey(k-ix)-temp
      ey(k+ix)=ey(k+ix)-temp
      ey(k)=ey(k)-2.*temp
      temp=.25*s*ez(k)
      ez(k-ix)=ez(k-ix)-temp
      ez(k+ix)=ez(k+ix)-temp
   1  ez(k)=ez(k)-2.*temp
      return
      end
C  --------
      subroutine xsplit(x,y,z,x0,y0,z0)
      if(ifix(x).ne.ifix(x0).and.((x-x0).ne.0.))go to 1
      call ysplit(x,y,z,x0,y0,z0)
      return
  1   x1=.5*(1+ifix(x)+ifix(x0))
      y1=y0+(y-y0)*((x1-x0)/(x-x0))
      z1=z0+(z-z0)*((x1-x0)/(x-x0))
      call ysplit(x,y,z,x1,y1,z1)
      call ysplit(x1,y1,z1,x0,y0,z0)
      return
      end
C --------
      subroutine ysplit(x,y,z,x0,y0,z0)
      if(ifix(y).ne.ifix(y0).and.((y-y0).ne.0.))go to 1
      call zsplit(x,y,z,x0,y0,z0)
      return
  1   y1=.5*(1+ifix(y)+ifix(y0))
      z1=z0+(z-z0)*((y1-y0)/(y-y0))
      x1=x0+(x-x0)*((y1-y0)/(y-y0))
      call zsplit(x,y,z,x1,y1,z1)
      call zsplit(x1,y1,z1,x0,y0,z0)
      return
      end
C --------
      subroutine zsplit(x,y,z,x0,y0,z0)
      if(ifix(z).ne.ifix(z0).and.((z-z0).ne.0.))go to 1
      call depsit(x,y,z,x0,y0,z0)
      return
  1   z1=.5*(1+ifix(z)+ifix(z0))
      x1=x0+(x-x0)*((z1-z0)/(z-z0))
      y1=y0+(y-y0)*((z1-z0)/(z-z0))
      call depsit(x,y,z,x1,y1,z1)
      call depsit(x1,y1,z1,x0,y0,z0)
      return
      end
C ----------------------
      subroutine depsit(x,y,z,x0,y0,z0)
C (Field components are treated as single-indexed in this subroutine)
      parameter(mxx=505,myy=9,mzz=255,ngtl=mxx*myy*mzz)
      common /fields/ex(ngtl),ey(ngtl),ez(ngtl),bx(ngtl),
     &by(ngtl),bz(ngtl)
     &,sm(27),q,qe,qi,qme,qmi,c
     &,ms(27),mx,my,mz,ix,iy,iz,lot,maxptl,maxhlf,ions,lecs,nstep
C   cell indices of half-way point:
      i=.5*(x+x0)
      j=.5*(y+y0)
      k=.5*(z+z0)
C   displacements in cell of half-way point:
      dx=.5*(x+x0) - i
      dy=.5*(y+y0) - j
      dz=.5*(z+z0) - k
      l=i+iy*(j-1)+iz*(k-1)
c   for 2<= z < mz-1
c     l=i+iy*(j-1)
c    1   +iz*(k-(isign(mz-3,k-2)+isign(mz-3,k-mz+1))/2-1)
C   current elements:
      qu=q*(x-x0)
      qv=q*(y-y0)
      qw=q*(z-z0)
      delt=.08333333*qu*(y-y0)*(z-z0)
C  Directive specifically for the CRAY cft77 compiler:
cdir$ ivdep
C  (This will make the compiler use the "gather-scatter" hardware.)
C  If one desires NO smoothing (risking the presence of alias-prone
C  high harmonics), one can replace the statement "do 1 n=1,27" by
C     n=14   (only n=14)
C  and boost the value of  q  by a factor 8.
      do 1 n=1,27
      ex(ms(n)+l+iy+iz)=ex(ms(n)+l+iy+iz)-sm(n)*(qu*dy*dz+delt)
      ex(ms(n)+l+iz)=ex(ms(n)+l+iz)-sm(n)*(qu*(1.-dy)*dz-delt)
      ex(ms(n)+l+iy)=ex(ms(n)+l+iy)-sm(n)*(qu*dy*(1.-dz)-delt)
      ex(ms(n)+l)=ex(ms(n)+l)-sm(n)*(qu*(1.-dy)*(1.-dz)+delt)
      ey(ms(n)+l+iz+ix)=ey(ms(n)+l+iz+ix)-sm(n)*(qv*dz*dx+delt)
      ey(ms(n)+l+ix)=ey(ms(n)+l+ix)-sm(n)*(qv*(1.-dz)*dx-delt)
      ey(ms(n)+l+iz)=ey(ms(n)+l+iz)-sm(n)*(qv*dz*(1.-dx)-delt)
      ey(ms(n)+l)=ey(ms(n)+l)-sm(n)*(qv*(1.-dz)*(1.-dx)+delt)
      ez(ms(n)+l+ix+iy)=ez(ms(n)+l+ix+iy)-sm(n)*(qw*dx*dy+delt)
      ez(ms(n)+l+iy)=ez(ms(n)+l+iy)-sm(n)*(qw*(1.-dx)*dy-delt)
      ez(ms(n)+l+ix)=ez(ms(n)+l+ix)-sm(n)*(qw*dx*(1.-dy)-delt)
  1   ez(ms(n)+l)=ez(ms(n)+l)-sm(n)*(qw*(1.-dx)*(1.-dy)+delt)
      return
      end
      function cvmgm(i1,i2,i3)
      integer i1, i2, i3, cvmgm
      if(i3.lt.0) then
	 cvmgm=i1
	   else
	     cvmgm=i2
	      endif
      return
      end