subroutine pardyn
c---calculate dynamics for nsteps in time (phase), each step of length dwt.
c---calculate space-charge impulse at beginning of each nsc steps.  call
c---output at end of every nout steps.
c---insert common blocks here
c----------------------------------------------------------------------
      save
      include 'param_sz.h'
      include 'var_char.h'
      include 'cfldscom.h'
      include 'constcom.h'
      include 'coordcom.h'
      include 'debug.h'
      include 'flagcom.h'
      include 'misccom.h'
      include 'ncordscom.h'
      include 'outcom.h'
      include 'pcordcom.h'
      include 'sccom.h'
      include 'syscom.h'
      include 'tstepcom.h'
      include 'ucom.h'
c
      common/aswap/acor(7,imaa)
      common/back/bmax,byz0,byzs,pzowmin,pzowmax,iback
      common/chars/title
      common/ct/var,rp,dzmin,dzmax,wrm,dwbm,xrmsbg,yrmsbg,
     1          xsmax,xpsmax,ysmax,ypsmax,istat,nebm
      common/errors/eraln(5,0:lmx)
      common/intype/intype ! mcd
      common/jones/ajones,zjones,zcath(imaa) ! mcd
      common/out6/nbphase,nbremesh
      common/outbuf/outcor(8,imaa*imb)
      common/part/nparticle(imaa)
      common/pnum/ipnum(imaa)
      common/save/zrsav
      common/wtstep/wtestep,nstep ! mcd
      real xyz(7)
      logical outflag
      equivalence (xyz(1),x)
      data outflag/.true./
c--------------------------------------------------------------------------
c*
      if(ngood.le.0)return
      inb=imaa*imbf/npoints
      if(inb.gt.maxbuf)inb=maxbuf
      dcon=wavel/360.         ! conversion degrees to cm
      dwtsc=nsc*dwt*dcon
      mt=0
c     For input type 10 and 11 nstep = no. of extra steps to be added
c     at beginning of loop to allow for generation at the photocathode
c     nstep is define in input 10 and 11
      mt=mt-nstep ! mcd
      msbl=0
      mdpoutl=0
c---start loop on time steps
   10 continue
      mt=mt+1
      swt=sin(wt*radian)
      cwt=cos(wt*radian)
      wt1=wt+dwt
c---set space-charge and output flags
      if(nsc.le.0)msc=1
c      if(nsc.gt.0)msc=mod(mt-1,nsc)
      if(nsc.gt.0)msc=mod(mt,nsc)
      if(beami.eq.0.)msc=1
      if(nout.le.0)mout=1
      if(nout.gt.0.and.mt.ne.0)mout=mod(mt,nout)
c---calculate space-charge fields from present particle distribution.
c---and apply space-charge impulse to particles.
      if(msc.eq.0) then
        if(iprog.eq.1) call scheff1(dwtsc,mt,wt1)
        if(iprog.eq.2) call scheff2(dwtsc,mt,wt1)
        if(iprog.eq.3) call scheff3(dwtsc,mt,wt1)
      else          ! isc
      endif         ! isc
c---check single-bunch beam loadin flag
      if(msbl.ne.0)call sbload(nesb)
      if(mdpoutl.ne.0)call dpout(nedpout)
      mdpoutl=0
      msbl=0
      np=0
      nle=0
c---start loop on particles
   20 continue
      np=np+1
cbm      itoto=np
      itoto=0
      nupar=np
c---get particle coordinates
   25 do 26 i=1,7
        xyz(i)=cord(i,np) !xyz(1,2,...,7) equiv to x,bgx,..,ne
   26 continue
      ne=rne
      gamma=gam(np)
      wtp=wt
      dwtp=dwt
      phinought=phizero(np)
c Upstream of first element (ne.le.0)
      if((intype.eq.10 .or. intype.eq.11 .or. intype.eq.15)
     .  .and. ne.le.0) then
c        5/10/87 Revised treatment of particles upstream of first element.
        bz=bgz/gamma
      if(bz.le.0.0) then   ! back-bombardment
            if(iback.eq.0) go to 155
            if(ne.gt.1 .or. (np.eq.1 .or. z.lt.z0)) go to 155
      endif
c---calculate remaining distance this particle would travel in this
c---step at this velocity
        dz=bz*dcon*dwtp
        zn = z + dz
c---set "end of element" flag to off
        iend=0
c        zcat = z coord of cathode for emission of this particle
        zcat = zcath(np)
        if(zn.lt.zcat) then
c         particle remains upstream of first element
c         advance the particle by dz
          z = zn
          x = x + dz*bgx/bgz
          y = y + dz*bgy/bgz
          go to 150
        endif
c         particle crosses into the first element
c---calculate time spent before reaching the first element.
        dz = zcat - z
        dwtp = dz/(bz*dcon)
c        move particle to beginning of first element
        z = zcat
        x = x + dz*bgx/bgz
        y = y + dz*bgy/bgz
        ne=1
        rne=1.
        nt=ntype(ne)
c        go to 150 -> 180 -> 30 -> 35 -> 70 to apply cell impulse and drift
        go to 150
      endif
c End upstream of first element
cmcd
      if(ne.gt.nel)go to 40
   30 nt=ntype(ne)
      if(nt.ne.27.and.nt.ne.29) go to 35
c---sbload element. set flag,zlen and and cells if np=1
c-- it is now in sbload
      if(np.ne.1) go to 34
      if(nt.eq.27)then
         msbl=1
         nesb=ne
      else
         mdpoutl=1
         nedpout=ne
      endif
   34 if(ibeg.eq.0)then
         ne=ne+1
         rne=rne+1
      else
         ne=ne-1
         rne=rne-1
      endif
      go to 30
c---apply cell impulse
   35 continue
      if(nt.eq.7) call cellimp(wtp,dwtp)
c---apply trwave impulse
      if(nt.eq.9) call trwimp(wtp,dwtp)
   40 continue
      bz=bgz/gamma
c------allow particles to go backward
c      if(bz.le.0.)go to 155
      if(bz.le.0.0) then   
            if(iback.eq.0) go to 155
            if(ne.gt.1 .or. (np.eq.1 .or. z.lt.z0)) go to 155
      endif
c---calculate remaining distance this particle would travel in this
c---step at this velocity
      dz=bz*dcon*dwtp
      zn=z+dz
c---set "end of element" flag to off
      iend=0
c-----allow particle to go backwards, set "begining of element" flag off
      ibeg=0
      if(ne.gt.0) go to 50
c---particle is upstream of first element.  assume drift.
      ne=0.
      rne=0.
      nt=ntype(1)
      go to 70
   50 continue
      if(ne.gt.nel) go to 70
      if(nt.eq.8)go to 130
      if(nt.eq.2)go to 80
      if(nt.eq.4)go to 100
      if(nt.eq.35)go to 146 
      if(nt.eq.36)go to 105
c-----allow particle to go backwards
c-----check if particle will cross into previous element during this step
c-----ntype 1,2,5,7, and 9 perform this check so skip it.
      if(nt.eq.1 .or. nt.eq.2 .or. nt.eq.7 .or. nt.eq.9
     *  .or. nt.eq.5)go to 65
      if(dz.lt.0 .and. ne.ge.1 .and. zn.lt.zloc(ne-1))then
c---particle will cross into previous element during this step
      zn=zloc(ne-1)
      dwtp=(zn-z)/(bz*dcon)
      ibeg=1
      go to 64
      endif
      if(zn.lt.zloc(ne)) go to 65
c---particle will cross into next element during this step
      zn=zloc(ne)
      dwtp=(zn-z)/(bz*dcon)
      iend=1
   64 continue
      if(nt.eq.26)go to 140
      if(nt.eq.30)go to 145
      if(nt.eq.39)go to 147
   65 continue
      go to (70,80,90,100,110,120,75,130,70), nt
c---after drift, linac cell impulse, or trwave impulse
   70 continue
      call drift(dwtp,iend,ibeg)
      go to 150
   75 continue
      if(el(11,ne).eq.0)go to 70
      call solnoid(dwtp,iend,ibeg,el(11,ne))
      go to 150
c---after solenoid
   80 continue
      call solnoid(dwtp,iend,ibeg,0.)
      go to 150
c---after quad
   90 continue
      call quad(zn)
      go to 150
c---after bend
  100 continue
      call bend(dwtp,iend)
      go to 150
c---after alpha magnet
  105 continue
      call alpham(dwtp,iend,np)
      go to 150
c---after buncher
  110 continue
      phase=(wtp*el(5,ne)/freq+el(6,ne))*radian
      call buncher(phase,iend,ibeg)
      dwtp=0.
      go to 150
c---after chopper
  120 continue
      phase= wtp*el(4,ne)/freq - el(5,ne)
      if(phase.gt.0.)phase=amod(phase,360.)
      if(phase.lt.0.)phase=360.-amod(-phase,360.)
      if(phase.gt.180.)phase=phase-360.
      if(abs(phase).gt.el(6,ne))go to 155
      go to 150
c---after tank
  130 continue
      call tank(wtp,dwtp,iend,ibeg)
      bz=bgz/gamma
      go to 150
c---after rotate
  140 continue
      call rotate
      go to 150
c---after cathode
  145 continue
      call cathode(dwtp,iend,ibeg)
      go to 150
c---after wiggler 
  146 continue
      call wiggler(dwtp,iend)
      go to 150      
c---after sextupole
  147 continue
      call sextupole(zn)
      go to 150
c---check for loss.
  150 continue
c-------allow particle to go backwards
      if(bz.le.0.0) then   
            if(iback.eq.0) go to 155
            if(ne.gt.1 .or. (np.eq.1 .or. z.lt.0)) go to 155
      endif
      if(ne.le.0 .or. ne.gt.nel)go to 152
      apsq=el(2,ne)**2    ! aperture
      rsq=x**2+y**2
      if(rsq.gt.apsq)go to 155
  152 continue
      wtp=wtp+dwtp
      dwtp=wt1-wtp
c---problem with number of significant digits so set dwtp=0 if very small.
      if(abs(dwtp).lt.1.e-4)dwtp=0.
c---check for end of element
      if(iend.eq.0.and.ibeg.eq.0) go to 180
c---particle is at end or begining of element.
c---if particle is at begining forget it
      if(ibeg.ne.0)go to 170
c---this is special code for pepperpot.
c      if(ne.eq.13.and.np.ne.1)then
c--- at end of element 13 so apply pepper pot.
c     x0=-1.1
c     do 151 i=1,10
c     x0=x0+.2
c     y0=-1.1
c     do 151 j=1,10
c     y0=y0+.2
c     if((x-x0)**2+(y-y0)**2.lt..00064516)go to 153
c151  continue
c--- particle stopped at pepperpot
c     go to 155
c     endif
  153 continue
      if(np.ne.1)go to 160
c---this is the reference particle.  save phase and energy.
      pr(ne)=wtp
      wr(ne)=(gamma-1.)*erest
      if(ne.eq.nel)then
c--- this is end of last element so output the cordinates of all the
c--- particles for FELEX  in the format: X,Bgx,Y,Bgy,Z,Bgz,Q(coul)
      zrsav=cord(5,1)
      ifelex=0
      if(ifelex.eq.1) then
       if(outflag)then
       open(unit=ndpout,file='dpout',status='new',form='formatted')
       outflag=.false.
       endif
      write(ndpout,'(6(1x,e15.9),1x,e15.9)')((cord(i,j),i=1,6),weight(j)
     *   ,j=1,ngood)
      endif
      endif
      go to 160
c---particle has been lost
  155 continue
      if(np.eq.1)go to 220
      call swap(np,wtp)
      if(np.lt.ngood)np=np-1
      ngood=ngood-1
      go to 187
c---check for output at end of element.
  160 continue
c      if(iel(3,ne).ne.0) then
      nuel=el(3,ne)
      if(nuel.ne.0) then
      call sortie(np,ngood,wtp)
      endif
c      if(iel(3,ne).eq.0)go to 170
      if(nuel.eq.0)go to 170
c---check if referance particle particle is less than zlimit ahead of
c---this particle if it is ignore it.
      if(cord(5,1)-cord(5,np).gt.abs(zlimit)) go to 170   
c---save coordinates in output buffer.
c---does a buffer already exist for this element
cliu      if(ne.eq.0) go to 170
      do 162 i=1,inb
      k=i
      if (ne.eq.neb(i)) go to 168
  162 continue
c---no buffer. are there any unused buffers
      do 164 i=1,inb
      k=i
      if (neb(i).eq.0) go to 166
  164 continue
c---no unused buffers. forget it.
      write(ndiag,*)' no buffer avaliable for element',ne,' particle',np
      go to 170
c---new buffer.  set neb
  166 continue
      neb(k)=ne
  168 continue
      if(nbuf(k).eq.npoints)then
c---buffer is full so output it. Then store particle in it.
      call output(neb(k),nbuf(k),outcor(1,(k-1)*npoints+1),npoints,mt)
      nbuf(k)=0
      neb(k)=0
      endif
      nbuf(k)=nbuf(k)+1
      nb=(k-1)*npoints+nbuf(k)
      outcor(1,nb)=x
      outcor(2,nb)=bgx/bgz
      outcor(3,nb)=y
      outcor(4,nb)=bgy/bgz
      outcor(5,nb)=wtp
      outcor(6,nb)=(gamma-1.)*erest
      outcor(7,nb)=ipnum(np)
      outcor(8,nb)=weight(np)
c---if there is a background bfield the beam is rotating.
c   therefor transform to the rotation frame of referance.
      if(ifld.ne.0 .or. poiflag)then
         if(ifoclal.eq.1) then
           cay=bfldlal(z)/(brhof*sqrt(gamma**2-1.))
         else
           cay=bfld(z,sqrt(x**2+y**2),brfld)/(brhof*sqrt(gamma**2-1.))
         endif
         outcor(2,nb)=outcor(2,nb)-.5*cay*y
         outcor(4,nb)=outcor(4,nb)+.5*cay*x
      endif
  170 continue
      if(iback.eq.1) then   
           if(iend.eq.1)ne=ne+1
           if(ibeg.eq.1)ne=ne-1     
      else
           ne=ne+1
      endif
      rne=ne
      if(ne.le.0.and.bgz.lt.0.)go to 155
      nt=ntype(ne)
c---check for misalignment errors
      if(iend.ne.0) then
          if(eraln(1,ne).eq.0.)go to 180
c---adjust transverse coordinates to axis of element
          con=bgx**2+bgy**2+bgz**2
          x=x-eraln(2,ne)
          bgx=bgx-bgz*eraln(3,ne)
          y=y-eraln(4,ne)
          bgy=bgy-bgz*eraln(5,ne)
          bgz=sign(sqrt(con-bgx**2-bgy**2),bgz)
      elseif(ibeg.ne.0)then
          if(eraln(1,ne+1).eq.0.)go to 180
          con=bgx**2+bgy**2+bgz**2
          x=x+eraln(2,ne+1)
          bgx=bgx+bgz*eraln(3,ne+1)
          y=y+eraln(4,ne+1)
          bgy=bgy+bgz*eraln(5,ne+1)
          bgz=sign(sqrt(con-bgx**2-bgy**2),bgz)      
      endif
c---does this particle have farther to go ?
  180 continue
      if(dwtp.gt.0.)go to 30
c---check distance limit if zlimit is negative particle is not discarded
      if((cord(5,1)-z).gt.zlimit.and.zlimit.gt.0.)go to 155
c---save particle back in cord array
      do 185 i=1,7
      cord(i,np)=xyz(i)
 185  continue
      gam(np)=gamma
      if(ne.gt.nel)nle=nle+1
      if(iback.eq.0) go to 187
      if(ne.ne.1 .or. pzowmax.eq.0.0)go to 187 
      if(phinought.gt.pzowmin.and.phinought.lt.pzowmax) then 
         if(amod(wt,5.).eq.0.0) write(nback) wtp,bz,x,y,z
      endif 
c---get next particle, if any
  187 continue
      if(np.lt.ngood)go to 20
c---end of one time step.  check for output
  190 continue
      if (ngood.le.0)go to 200
      wt=wt1
      if(mout.eq.0) then
         call output(mout,mout,outcor,npoints,mt)
      endif
c---check output buffer.  have all good particles already passed
c---element neo?
      do 197 i=1,inb
c---check output buffer. if full output it. if reference particle
c---is greater then abs(zlimit) beyond it output it
      if(nbuf(i).eq.npoints)then
          neo=neb(i)
          go to 196
      endif
      if(neb(i).ne.0.and.cord(5,1).gt.zloc(neb(i))+abs(zlimit)) then
          neo=neb(i)
          go to 196
      endif
      if(nbuf(i).eq.0)go to 197
      neo=neb(i)
      do 195 np=1,ngood
      nne=cord(7,np)
      if(nne.le.neo)go to 197
  195 continue
c---yes.  empty buffer.
  196 continue
      call output(neo,nbuf(i),outcor(1,(i-1)*npoints+1),npoints,mt)
      nbuf(i)=0
      neb(i)=0
  197 continue
  199 continue
      if(mt.ge.nsteps)return
      if(nle.ge.ngood)return
      go to 10
  200 continue
      write(ndiag, 210)mt
  210 format(' no more good particles after step',i5/)
      go to 221
  220 continue
      write(nnout, 230)wtp,xyz
      write(*,*) ' Reference particle lost '
  221 continue
      do 225 j=1,inb
      if(neb(j).eq.0) go to 225
      if(neb(j).lt.ne) then
        call output
     1      (neb(j),nbuf(j),outcor(1,(j-1)*npoints+1),npoints,mt)
      endif
      neb(j)=0.
      nbuf(j)=0.
  225 continue
  230 format(' reference particle lost, wt=',f8.1,/
     1'  x, bgx, y, bgy, z, bgz, rne',/    
     2 4f8.4,2f8.2,f5.0,/)
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*