c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
c WARNING                                                              c
c Warning for change of System and Operating system we have to do      c
c one  change of the variable kvers in main program ligne 282          c
c utilised in ddate                                                    c
c----------------------------------------------------------------------c
c                                                                      c
c                            VERSION V5.03                             c
c                                                                      c
c     October 90  Orsay    V4.00                                       c
c     Release Mai 92       V4.10                                       c
c     Release June 92      V4.20                                       c
c     Release October 92   V4.21                                       c
c     Release December 92  V4.22                                       c
c     Release July 93      V4.30                                       c
c     Release March 94     V4.31                                       c
c     Release April 94     V4.32                                       c
c     Release June  95     V5.00                                       c
c     Release June  96     V5.01                                       c
c     Release July  96     V5.02    VMS / UNIX                         c
c     Release January 97   V5.03    VMS / UNIX                         c
c     Release September 2000 V5.03  Multi plateforme sans cernlib      c
c     Release April 2010   V5.03    devient version officielle LAL     c
c                                                                      c
c                                                                      c
c     B. Mouton  LAL Bat 200 91405 Orsay France                        c
c                                                                      c
c     MOUTON@LAL.IN2P3.FR                                              c
c                                                                      c
c     Version ALPHA                                                    c
c     This version of PARMELA is a pachwork of 3 versions              c
c     one of from Lloyd Young (version of July 88) of Los Alamos,      c
c     another of KT Mc Donald of Brookhaven and one of SLAC with       c
c     some local modifications, with modifications and news            c
c     routines from H Liu of IHEP Beijing,with news routines from      c
c     CEBAF.                                                           c
c                                                                      c
c---PhAse and Radial Motion in Electron Linear Accelerators            c
c---(and associated transport systems)                                 c
c                                                                      c
c   on VAX VMS you have to increase your Paging file quota to 20000    c
c                                                                      c
c   $for/noop parmelav503                                              c
c   $link/exe=parmelav5 parmvelav503,'lib$'                            c
c                                                                      c
c   # Makefile for PARMELA cas un seul fichier                         c
c   # on Alpha  OSF1                                                   c
c   #                                                                  c
c   parmela : parmelav5.f                                              c
c           f77 -g -O0 parmelav5.f -o parmelav5  ${CERNLIBS}           c
c                                                                      c
c   # on HP                                                            c
c   parmela : parmelav5.f                                              c
c     f77 +ppu -K -g -w +E1 -O parmelav503.f -o parmelav5hp ${CERNLIBS}c
c                                                                      c
c----------------------------------------------------------------------c
c                                                                      c
c   Pour le cas de la version eclatee 07/2000                          c
c   #!//bin/ksh                                                        c
c   # gmake sans cernlib                                               c
c   programe="parmela"                                                 c
c                                                                      c
c   for files in *.f                                                   c
c   do                                                                 c
c   objname="${files:%.f}"                                             c
c   echo $objname                                                      c
c   objname=${objname}.o                                               c
c   liste=${liste}" ${files}"                                          c
c   #echo $liste                                                       c
c   echo ${objname}:${files} >>temp1                                   c
c   echo "\t \$(FF) -c \$(FFLAGS) "${files} >>temp1                    c
c   echo " ">>temp1                                                    c
c   done                                                               c
c                                                                      c
c   echo "${programe}: \$(OBJ)" >>temp1                                c
c   #echo "\t \$(FF)  \$(FFLAGS) \$(OBJ) -o "${programe}"\$(VER)"      c
c          \${CERNLIBS} >>temp1                                        c
c   echo "\t \$(FF)  \$(FFLAGS) \$(OBJ) -o "${programe}"\$(VER)"       c
c          >>temp1                                                     c
c   echo " " >>temp1                                                   c
c   echo "clean:" >>temp1                                              c
c   echo "\t rm *.o">>temp1                                            c
c                                                                      c
c   # ecriture dans le fichier de makefile definitif                   c
c                                                                      c
c   ## Pour HP                                                         c
c   echo FF=f77 >>makefile                                             c
c   #echo FFLAGS=+ppu +U77 -K -g -w +E1 -O >>makefile                  c
c   echo FFLAGS= -K -w +E1 +Oall >>makefile                            c
c   echo " " >>makefile                                                c
c   echo "VER=alonehpopt" >>makefile                                   c
c                                                                      c
c   ## Pour OSF                                                        c
c   #echo FF=f77 >>makefile                                            c
c   #echo FFLAGS= -g -O2 >>makefile                                    c
c                                                                      c
c   echo " " >>makefile                                                c
c   echo SRC=${liste} >>makefile                                       c
c   echo " " >>makefile                                                c
c   echo "OBJ=\$(SRC:.f=.o)" >>makefile                                c
c   echo " " >>makefile                                                c
c   cat temp1 >> makefile                                              c
c   rm -f temp1                                                        c
c   #c'est fini                                                        c
c----------------------------------------------------------------------c
c                                                                      c
c          modification des call scheff pour input0 02/92              c
c          modification dans emit90 pour kk=0 05/92 --> version 4.10   c
c          modification dans elips pour negative square root 05/92     c
c          outlal, emit90 and elips are in double precision   05/92    c
c          out6 is in double precision 06/92                           c
c          version 4.20  02/06/92                                      c
c          modification dans out6dp boucle 9996 remplace cor par cord  c
c          version 4.21                                                c
c          version 4.22 ip.eq.999 for imp.dat cell and trwave          c
c          modification dans cellfeld rcir en rcur  04/93              c
c          version 4.30 suppression of calls hbook                     c
c          version 4.31 24/03/94                                       c
c          pvol is in double precision pvoldp 03/94                    c
c          pvoldp supressed icall 03/94                                c
c          intro de cord(6,1)=cor(6,1) in main after 125 continue 03/94c
c          scheff2 loop 20 i=2,nr ----> i=2,nr1 03/94                  c
c          in out6dp change in format 220 zmax = f8.4 in               c
c          zmax = f8.2 03/94                                           c
c          version 4.32 18/04/94                                       c
c          introduction scaling factor on foclal card 04/94            c
c          modif /pcord/ intro nupar 05/94                             c
c          nspl=10000  ----> 10001          11/94                      c
c          nsplc=1497 ----> 30005=10001*3+2 11/94                      c
c          in zoutlal replaced in data bcd : design by trwave 11/94    c
c          renumber subroutine input 27/03/95                          c
c          input 11 27/03/95                                           c
c          histogram1 and histogram2 28/03/95                          c
c          renumber main program 29/03/95                              c
c          subroutine poisson 30/03/95                                 c
c          input 12 31/03/95                                           c
c          subroutine sextupole 11/04/95                               c
c          modification of space charge card 11/04/95                  c
c               new parameter iprog ;                                  c
c               scheff beami iprog point and subroutine dependent      c
c          modification extension infld in sffld 01/06/95              c
c          zout introduction de angle pour signe arctangente 04/09/95  c
c          modification input11 write ngood et tmax=2.5 sigma 12/09/95 c
c          modif dans input11 du calcul de degmax                      c
c                         remplace tmax par tfac   28/09/95            c
c          intro dans parmela routine gpindp plus dans cernlib 28/09/95c
c          deplacement de l'open de parmtrch cas ip=3 19/02/96         c
c          modification input 4 26/06/96                               c
c                               - 8 parameters one first card input 4  c
c          modification input11 05/07/96 two news parameters angle of  c
c            injection of laser light and distane focal point cathode  c
c          modification input10 23/07/96 two news parameters angle of  c
c            injection of laser light and distane focal point cathode  c
c          13/12/96 modif zout pour trwave (phase)                     c
c          13/12/96 modif fieldlal pour nchp.lt.nspl                   c
c          13/01/97 modif fieldlal pour free format                    c
c          13/01/97 modif common/image/                                c
c          25/02/97 intro dans subroutine input pour input15           c
c            dppi,dpepsout double precision et modif de l'appel gpindp c
c          03/09/97 mise sous Unix des statistiques elapsed cpu time   c
c          04/09/97 version pour HP (c_hp)                             c
c          14/11/97 nettoyage des variables non utilisees              c
c                   intro du call date_and_time sous VMS pour an 2000  c
c                   tous les open avec status='scratch' sans filename  c
c          03/02/98 modifation de zout et zoutlal pour cas solenoid    c
c          18/09/98 introduction de iadp a la place ngood dans out6dp  c
c                   write(ndes2), write(ndes1)                         c
c          29/03/00 reprise version eclate a partir du source non      c
c                   non eclate du 31/05/99  pour correction d'un bug   c
c                   dans une version ephemere dite v503b               c
c          18/07/00 abandon de la cernlib                              c
c                   suppression de call histo dans main et input       c
c          20/07/00 Modification de ddate.f pour essayer une methode   c
c                   universelle de date, modification du main          c
c                   introduction de kvers, introduction de la          c
c                   subroutine datime.f machine dependante             c
c                   suppression de l'utilisation de la cernlib         c
c          13/09/00 Introduction des include 'filename.h'              c
c                   Les fichiers .h  contiennent les parameters et les c
c                   common.                                            c
c                   Modification des noms de variable dans scwall      c
c                   Modification npoinz en npoints dans scheff2        c
c                   Declaration real zloc, gam, beami, brhof, radian , c
c                   twopi, z0 dans out6dp                              c
c                   Dans les appels a scheff ajouter 2 variables mt    c
c                        et wtsauv                                     c
c                   Modification de swap pour /pcord/                  c
c                   Suppression des instructions                       c
c                      dimension ac(14) cc(9) bc(1+(1+lmx)*13          c
c                      equivalence (ac(1),pi) (bc(1),nel) (cc(1),wt0)  c
c                   Ajouter mt aux arguments de la subroutine output   c
c                   Debug version modif main pardyn ouvrir trwave      c
c                   new common debug , new possibility for ip          c
c                   negative value means keep file parmdebug           c
c                   common/debug/idebug et tous les idebug             c
c    octobre 2002   augmentation du nombre de cartes cell numc 20-->60 c
c                                                                      c
c    29/10/2002     nom de l'exe sous OSF et HP-UX : parmela           c
c                                                                      c
c            nom de l'exe sous OSF et HP-UX : parmela                  c
c                                                                      c
c    19/04/2007    introduction de la polarisation ksi1,ksi2,ksi3      c
c                  modification lecture input4                         c
c                  modification out6dp  et sortie                      c
c                  18 variables dans ndes1 et ndes2 dans out6          c
c                                                                      c
c Last revision 14/11/2007 mise sous Linux                             c
c             utilisation du debugger                                  c
c             gdb /exp/sera/mouton/Parmela/V503/Linux/parmelalinux     c
c             pour cela faire cp toto.in parmin avant de lancer gdb    c
c                                                                      c
c    14/11/2007     modification des open de file.f, pour utilisation  c
c                   du disque local de la machine /tmp                 c
c                   dans le parmela.sh il faut :                       c
c                                      mv /tmp/parmdesz toto.desz      c
c Last revision 18/06/2008                                             c
c             un seul source de parmela v503   compatible OSF LINUX    c
c             un seul nom d'exe parmela                                c
c                                                                      c
c     Release octobre 2009  V5.03                                      c
c         modification du calcul de dzz dans fieldlal                  c
c              dzz=z(2)-z(1) ----> dzz=(zmax-zmin)/(nchp-1)            c
c         augmentation du nombre de points pour le  champ magnetique   c 
c         nspl 1001 ---> 10001                                         c
c         nouveau parametre nptcb = 10000 dans param_sz.h              c
c         creation d'une nouvelle carte nfieldb avec comme parametre   c
c          nptcbu nombre de points que l'on veut utiliser pour le      c 
c          champ magnetique < ou = 10000                               c
c    04/02/2010 retrait de la carte nfieldb                            c
c               modification de la 1er carte coil                      c
c               des tests dans field, solnoid si nb pt champ b > nptcb c
c                                                                      c 
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
c-----------------------------------------------------------------------
c---->program parmela
      subroutine parmela()
c-----------------------------------------------------------------------
c label n --> n*100
c go to  7000 --> 7008
c do     8000 --> 8017
c format 9000 --> 9025
c-----------------------------------------------------------------------
c
      include 'param_sz.h'
      include 'var_char.h'
      include 'coordcom.h'
      include 'constcom.h'
      include 'debug.h'
      include 'bfieldcom.h'
      include 'flagcom.h'
      include 'misccom.h'
      include 'ncordscom.h'
      include 'outcom.h'
      include 'phizoutcom.h'
      include 'pipecom.h'
      include 'sccom.h'
      include 'syscom.h'
      include 'tstepcom.h'
      include 'ucom.h'
      include 'upawcom.h'
      character*25 kvers
      common/aswap/acor(7,imaa)
      common/back/bmax,byz0,byzs,pzowmin,pzowmax,iback
      common/cathcur/rcur
      common/chars/title
      common/cout7/corout7(6,imaa),bgs,nelal,nbufl,iout7
      common/ct/var,rp,dzmin,dzmax,wrm,dwbm,xrmsbg,yrmsbg,
     1xsmax,xpsmax,ysmax,ypsmax,istat,nebm
      common/ddat/ctime,ddat
      common/errors/eraln(5,0:lmx)
      common/image/jj,ij
      common/intype/intype ! mcd
      common/jones/ajones,zjones,zcath(imaa) ! mcd
      common/out6/nbphase,nbremesh
      common/outbuf/outcor(8,imaa*imb)
      common/outpar/noutp,outpar(0:lmx,10),moutp(10)
      common/part/nparticle(imaa)
      common/pnum/ipnum(imaa)
      common/save/zrsav
      common/sbload1/tsbparm(5,50)
      common/tapes/ifl(2)
      common/version/kvers
      common /wtstep/wtstep,nstep ! mcd
c--------------------------------------------------------------------------
      character Istart_date*9,Iend_date*9,Istart_time*8,Iend_time*8
      real*4    time_begin

c_with cernlib
c      call timest(1.0e9)
c      call timed(t0)
c

      write(*,*) ' Go Parmela'
c_for ibm
c      time_begin = secnds(0.0)

      call datime(Istart_date,Istart_time)
c*
c---set constants
c
       kvers='V503 on LINUX'
c       kvers='V503 on OSF'
c
      pi=4.*atan(1.)
      twopi=2.*pi
      radian=pi/180.
      clight=29979.29   ! 10.e6 cm/s
      erest=.511        ! MeV
      brhof=1704.51     ! gauss.cm moc/e
c     brhof = (9.109389E-31*299792458/1.602177E-19)*(100)*(10000) !1704.5
      zlimit=1000000.   ! cm
c---set flags
      idebug=0          ! for debuging in subroutine ouvrir
      ifld=0            ! flag for external electric field
      poiflag=.false.   ! flag for external magnetic field
      nel=0             ! number of element
      ip=0              ! flag of output type
      isc=0             ! flag for space charge calculation
      istat=0           ! flag for statistics
      ifl(1)=0
      ifl(2)=0
      iback=0
      ltype=1           ! by default disk and washer cells linac
      ifoclal=0         ! bm bzlal
      npipe=0           ! mcd
      nstep=0           ! mcd
      wtstep=0.         ! mcd
      ajones=0.         ! mcd
      nout=0            ! mcd
      do 8001 n=1,10    ! mcd
      moutp(n)=0        ! mcd
      do 8001 i=1,lmx   ! mcd
      outpar(i,n)=0.    ! mcd
 8001 continue          ! mcd
      icfield = 0 ! cbm
      ictrw= 0    ! cbm
      nbphase=0   ! cbm nombre de phase en sortie
      nbremesh=0  ! cbm nombre de remesh
c----set polarization to zero
      do 7999 n=1,imaa
         ksi1(n)=0.
         ksi2(n)=0.
         ksi3(n)=0.
 7999 continue
c---
c---
      open(unit=nin,file='parmin',status='old',
     *     access='sequential',form='formatted')
      open(unit=nnout,file='parmdebug',status='unknown',
     *        access='sequential',form='formatted')
      call ouvrir
c---
      call setwall
c---set misalignments to zero.
      do 8002 n=1,100
      do 8002 i=1,5
      eraln(i,n)=0.
 8002 continue
c---read input data
c---the read(jj)  reads one card with a label and all
c---subsequent cards without labels.  the stop subroutine compares the label
c---with a list of words and assigns a value to jj corresponding to the
c---position of the word in the dictionary.  the word list and corresponding
c---values of jj are:
c---(drift     ,1), (solenoid  ,2), (quad      ,3), (bend      ,4),
c---(buncher   ,5), (chopper   ,6), (cell      ,7), (tank      ,8),
c---(trwave    ,9), (coil      ,10),(run       ,11),(input     ,12),
c---(output    ,13),(title     ,14),(scheff    ,15),(zout      ,16),
c---(adjust    ,17),(start     ,18),(restart   ,19),(continue  ,20),
c---(save      ,21),(end       ,22),(zlimit    ,23),(errors    ,24),
c---(change    ,25),(rotate    ,26),(sbload    ,27),(cfield    ,28),
c---(dpout     ,29),(cathode   ,30),(design    ,31),(pipe      ,32),
c---(foclal    ,33),(backb     ,34),(wiggler   ,35),(alpham    ,36),
c---(stat      ,37),(poisson   ,38),(sextupole ,39)
c---the subroutine reads numbers from the card in a free format.
c---the values are stored in vv, and there are nn of them.
 99   continue
cbm      call cread(jj)
      call cread
      go to (100,100,100,100,100,100,100,100,100,1000,
     1       1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,
     2       2100,2200,2300,2400,2500,100,100,2800,100,100,
     3       3100,3200,3300,3400,100,100,3700,3800,100),jj
c---after transport element or linac cell
c---after drift:     vv(n)=l,aper,iout
c---after solenoid:  vv(n)=l,aper,iout,h
c---after quad:      vv(n)=l,aper,iout,hp
c---after chopper:   vv(n)=0,aper,iout,freq,phi,dphi
c---after cathode:   vv(n)=0,aper,iout,radius of curvate of cathode
c---after wiggler:   vv(n)=l,aper,iout,periods,steps,b0,ky/kw
c---after sextupole: vv(n)=l,aper,iout,Bpole/ape**2
 100  continue
      if(nn.gt.11.and.jj.ne.7.and.jj.ne.9.and.jj.ne.27)go to 7000
      if(nel.ge.lmx)go to 7001
      if(nn.lt.2)go to 7002
      nel=nel+1
      ntype(nel)=jj
      do 8003 n=1,nn
      el(n,nel)=0.      ! mcd
      if(n.gt.nn)go to 8003
      el(n,nel)=vv(n)
 8003 continue
      zloc(nel)=el(1,nel)
      if(nel.gt.1)zloc(nel)=zloc(nel)+zloc(nel-1)
      if (ip.eq.1 .or. ip.eq.2) then
        iout = vv(3)
        write(nnout,9001) nel,jj,vv(1),zloc(nel),vv(2),iout
 9001   format (
     . ' element no. [ne] =',t40,i8/
     . ' element type [ntype(ne)] =',t40,i8/
     . ' element length [el(1,ne)] =',t40,1pg12.4,' cm'/
     . ' z coord at downstream end [zloc(ne)] =',t40,g12.4,' cm'/
     . ' radial aperture at exit [el(2,ne)] =',t40,g12.4,' cm'/
     . ' output at end of element [el(3,ne)] =',t40,i8)
        if((jj.eq.3) .or. (jj.eq.39)) write(nnout,9002) vv(4)
 9002   format(
     . ' magnetic field gradient = ',t40,g12.3,' gauss/cm')
      endif
      if(jj.eq.2) ifld=1
      if(jj.eq.26) go to 2600
      if(jj.eq.29) go to 99
      if(jj.eq.27) go to 2700
      if(jj.eq.30) rcur=vv(4)
      if(jj.eq.35) go to 3500
      if(jj.eq.36) go to 3600
      if(jj.eq.39) go to 99
      go to (99,99,99,400,500,99,700,800,900), jj
c---after bend:  vv(n)=l,aper,iout,wr,alpha,beta1,beta2,psi1,psi2,r1,r2
  400 continue
      do 8004 n=5,9
      el(n,nel)=el(n,nel)*radian
 8004 continue
      go to 99
c---after buncher:  vv(n)=0.,aper,iout,wmax,freq,phi,w
  500 continue
      w=w0
      if(nn.ge.7)w=el(7,nel)
      g=w/erest
      bg=sqrt(g*(2.+g))
      cay=twopi*el(5,nel)/(bg*clight)
      el(8,nel)=cay
      go to 99
c---after trwave;
  900 continue
      ictrw = ictrw+1
      call trwfield
      go to 99
c---after cell:  vv(n)=cl,aper,iout,phi,e0,nc,(dwtmax)
  700 continue
      el(5,nel)=.01*el(5,nel)
      if (nn.lt.8) el(8,nel)=0.
      phideg(nel)=el(4,nel)
c      phi=el(4,nel)*radian this is now defined in cellfld
c      el(9,nel)=sin(phi)    and this
c      el(10,nel)=cos(phi)  and this also
c---if unspecified, set maximum step size in cell to 10 degrees.
      if(nn.lt.7)el(7,nel)=10.
c---get field information for this cell
      cl=el(1,nel)
      if (el(8,nel).ne.0.) cl=2.*cl
      nc=el(6,nel)
c      iel(6,nel)=nc
      if (ip.eq.1 .or. ip.eq.2) then
        write(nnout,9003)vv(4),vv(5),nc,el(7,nel),el(11,nel),el(8,nel)
 9003   format (
     . ' phase of cell [phi0] =',t40,g12.3,' degrees'/
     . ' average axial electric field [el(5,ne)] =',t40,g12.3,' MV/m'/
     . ' cell no. [nc] =',t40,i8/
     . ' max step size [dwtmax] =',t40,g12.3,' degrees'/
     . ' axial static magnetic field =',t40,g12.3,' gauss'/
     . ' cell symmetry [el(8,ne)] =',t40,g12.3)
        if (el(8,nel).eq.0.) write(nnout,9004)
 9004   format (' cell is symmetric about midplane in z')
        if (el(8,nel).gt.0.) write(nnout,9005)
 9005   format (' half cell, with zero field at entrance')
        if (el(8,nel).lt.0.) write(nnout,9006)
 9006   format (' half cell, with zero field at exit')
      endif
      call cellfld(cl,nc)
      go to 99
c---after tank:  tl,aper,iout,e0t,nc,p,phi
  800 continue
      el(8,nel)=el(1,nel)/el(5,nel)
c---convert e0t from mv/m to mv/cm
      el(4,nel)=.01*el(4,nel)
c---calculate beta and k**2 for the tank
      beta=2.*el(8,nel)/(el(6,nel)*wavel)
      cl=el(8,nel)
      caysq=((wavel/(2.*cl))**2 - 1.)*(twopi/wavel)**2
      el(9,nel)=beta
      el(10,nel)=caysq
      go to 99
c---after field:  vv parameters depend on field subroutine
 1000 continue
      icfield=icfield+1
      call field
c---define field structure after coil
c--- vv(n)=z,r,ni,(zbeg,zend,nptcbu,order)
c---in ( ) required on first coil card ignored on rest
      go to 99
c---after errors.  generate errors in geners
 2400 continue
      call geners
      go to 99
c---after run:  vv(n)=irun,ip,freq,z0,w0,ltype
 1100 continue
c     process run card         ! mcd
      if (nn.ge.1) irun=vv(1)  ! mcd
c--debug
c-- Si ip negatif sur carte run on a la version de debug
      if (nn.ge.2) then
         if(vv(2).lt.0.) then
            idebug=1
            write(*,*) '-----------------------------------'
            write(*,*) 'You have ask the version with debug'
            write(*,*) '-----------------------------------'
            write(nnout,*) '-----------------------------------'
            write(nnout,*) 'You have ask the version with debug'
            write(nnout,*) '-----------------------------------'
         else
         endif
         ip=abs(vv(2))
      endif
c--debug
      if (nn.ge.3) freq=vv(3)  ! mcd
      if (nn.lt.4) go to 99 
      if (nn.ge.4) z0=vv(4)    ! mcd
      if (nn.lt.5) go to 99
      if (nn.ge.5) w0=vv(5)    ! mcd
      if (nn.ge.6) ltype=vv(6) ! mcd 
      if(freq.gt.0.) wavel=clight/freq
      write(nnout,9007)irun
 9007 format (/' '/
     . ' run no. =',t40,i8/)
      if((ip.gt.4).and.(ip.ne.999)) then
      ip=1
      write(nnout,9008)ip
 9008 format(
     .' print option modified =',t40,i8)
      else
      write(nnout,9009)ip
 9009 format(
     . ' print option =',t40,i8)
      endif
c---
      if (ip.eq.3) then
      open(unit=nuch,file='parmtrch',access='sequential',
     *     status='unknown')
      write(nuch,*)'  nc  sym  r  z  ez  er  bt  '
      endif
      if (ip.eq.999) open(unit=nimp,file='parmimp',status='unknown')
c---
      write(nnout,9010)freq,wavel,z0,w0,ltype
 9010 format(
     . ' rf frequency =',t40,f9.2,'    MHz'/
     . ' wavelength =',t40,1pg12.3,' cm'/
     , ' z0 =',t40,g12.3,' cm'/
     , ' reference particle energy =',t40,g12.3,' MeV'/
     . ' linac type =',t40,i8)
      do 8005 i=1,4
      corout7(i,1)=0.
      cor(i,1)=0.
 8005 continue
      cor(5,1)=z0
      corout7(5,1)=z0
      corout7(6,1)=w0
      g=w0/erest
      gam(1)=g+1
      cor(6,1)=sqrt(g*(2.+g))
      cord(6,1)=cor(6,1)
      if (npoints.lt.1)npoints=1
      ngood=npoints
      nparticle(1)=1
      phizero(1)=0.
      if(nn.lt.6)go to 99
      go to 99
c---after input:  vv parameters depend on input subroutine
 1200 continue
      call input
      ngood=npoints
      do 8006 i=1,ngood
      nparticle(i)=i
      ipnum(i)=i
 8006 continue
      go to 99
c---after output:  vv parameters stored in optcon array
c---after output:  vv parameters stored in outpar array
c
 1300 continue
      noutp=vv(1)+1
      if(noutp.gt.10 .or. noutp.lt.1) then
      write(ndiag,9011)
 9011 format(' invalid OUTPUT type')
      go to 99
      endif
      moutp(noutp)=nn
      do 8007 n=1,nn
      optcon(n)=vv(n)
      outpar(n,noutp)=vv(n)
 8007  continue
      if(vv(1).gt.1..and.vv(1).lt.5.)then
      imbf=imb-1
      else
      imbf=imb
      endif
      go to 99
c---after title:  read title information from next card
 1400 continue
      read(nin,9012)title
 9012 format(a80)
      write(nnout,*)title
      go to 99
c---after scheff:  vv(1)=beam current; vv(2)= program number
c--- rest depend on space-charge subroutines
 1500 continue
      do 8008 n=1,12
      sce(n)=0.
      if (n.le.nn) sce(n)=vv(n)
 8008 continue
      beami=vv(1)
      iprog=vv(2)
      if(isc.eq.0) then
        if(iprog.eq.1) call scheff1(0.,-99,0.)
        if(iprog.eq.2) call scheff2(0.,-99,0.)
        if(iprog.eq.3) call scheff3(0.,-99,0.)
         isc=1
      else          ! isc
      endif         ! isc
      go to 99
c---after zout.  print locations of elements
 1600 continue
      if(ifoclal.eq.1) then
        call zoutlal
      else
        call zout
      endif
      go to 99
c---after adjust:  call arbitrary subroutine named adjust
 1700 continue
      call adjust
      go to 99
c---after change.  call change subroutine
 2500 continue
      call change
      go to 99
c---after rotate
 2600 continue
      el(4,nel)=el(4,nel)*radian
      go to 99
c---  after sbload paramaters are l,aperature,iout,zlen,ncells,a1,a2,a3
c---  ,a4,a5,ta1,ta2,ta3,ta4,ta5. if ncells is negative then the paramaters
c---  ta1 to ta4 are required. ta5 is optional. if ncells is negative the
c---  transvers wake will be applied also.
 2700 continue
      if(nn.gt.11)then
      numsb=numsb+1
      tsbparm(1,numsb)=vv(11)
      el(11,nel)=numsb
      tsbparm(2,numsb)=vv(12)
      tsbparm(3,numsb)=vv(13)
      tsbparm(4,numsb)=vv(14)
      if(el(5,nel).lt.0..and.nn.lt.14)go to 7002
      if(nn.gt.14)tsbparm(5,numsb)=vv(15)
      endif
      go to 99
c---after start: start dynamics calculations
 1800 continue
      if(nn.lt.5)go to 7002
      wt0=vv(1)
c---  For input type 10 and 11, WTSTEP is the half width of the laser pulse.
c     The clock must be set back by this amount to generate the time
c     distribution of photoelectrons.
      wt0=wt0-wtstep
      dwt=vv(2)
      nsteps=vv(3)
      nsc=vv(4)
      nout=vv(5)
      go to 7003
c---after restart.  set the starting time (phase), initialize cord
c---from cor, and find largest element number (len) that any particle is in.
 1900 if(nn.lt.4)go to 7002
      dwt=vv(1)
      nsteps=vv(2)
      nsc=vv(3)
      nout=vv(4)
c	-----------------
      if(jj.eq.19) then				
  		npoints=-vv(6)
      endif
c	-------------------------
      if(nn.gt.5.and.(vv(6).gt.0..and.vv(6).le.26))then
        i=vv(6)-1
      call orest1(i)
      else
      call orest2
      endif
      read(nsav) cor,wt0,ngs
      close(nsav)
      ngood=ngs
C	----------------------------
      IF(NGS.NE.NPOINTS) then
        call appendparm
        print *,'ngs=',ngs,'npoints=',npoints
        stop ' Abnormal stop ngs.ne.npoints main '
      endif
C	-------------------------------
 7003 continue
      if (ip.eq.1 .or. ip.eq.2)
     .   write (nnout,9013) wt0,dwt,nsteps,vv(5),nsc,nout
 9013 format (
     . ' starting phase =',t40,1pg12.3,' degrees'/
     . ' integration step size =',t40,g12.3,' degrees'/
     . ' no. of steps in integration =',t40,i8/
     . ' shift of reference particle =',t40,g12.3,' cm'/
     . ' call SCHEFF every',t40,i8,'     step(s)'/
     . ' call OUTPUT every',t40,i8,'     step(s)'/)
      wt=wt0
      do 8009 i=1,maxbuf
      neb(i)=0
      nbuf(i)=0
 8009 continue
c---if nn.gt.4, the 5th parameter is dist to shift ref particle
c   only for restart. not for start
      if (jj.eq.19.and.nn.gt.4) then
        cor(5,1)=cor(5,1)+vv(5)
      endif
c	--------------------------------------------------
c--------- vv(5)= dz
      if(jj.eq.19.and.vv(5).ne.0.) then
            cor(1,1)=0.  ! for energy change of ref part
            cor(2,1)=0.
            cor(3,1)=0.
            cor(4,1)=0.
            print *,'vv(7) =',vv(7)
            cor(6,1)=sqrt((1.+vv(7)/0.511)**2-1.)
      endif
c	---------------------------------------------------
      do 8010 n=1,ngood
      do 8011 i=1,6
      cord(i,n)=cor(i,n)
 8011  continue
      gam(n)=sqrt(1.+cord(2,n)**2+cord(4,n)**2+cord(6,n)**2)
c---find which element each particle is in
      ne=0
      z=cord(5,n)
      if(z.lt.0.)go to 7004
      do 8012 i=1,nel
      ne=i
      if(z.lt.zloc(i))go to 7004
 8012 continue
      ne=nel+1
 7004 cord(7,n)=ne
c     if generating from a photocathode, put ne=0
      if(jj.eq.19) go to 8010		
      if (intype.eq.10 .or. intype.eq.11 .or. intype.eq.15) 
     .   cord(7,n)=0.
 8010 continue
c---if nn.gt.4, the 5th parameter is dist to shift ref particle
c   only for restart. not for start
c---start dynamics calculation
      call pardyn
      go to 99
c---after continue:  used to change step size and continue computations
 2000 continue
      if(nn.lt.4)go to 7002
      dwt=vv(1)
      nsteps=vv(2)
      nsc=vv(3)
      nout=vv(4)
c---if nn.gt.4, the 5th parameter is dist to shift ref particle
      if (nn.gt.4)cord(5,1)=cord(5,1)+vv(5)
      call pardyn
      go to 99
c---after save
c---save cord in cor
 2100 continue
      ngs=ngood
      do 8013 n=1,ngood
      do 8013 i=1,6
      cor(i,n)=cord(i,n)
 8013 continue
      write(nnout,*) 'zr of save ',zrsav
      wt0=wt
      if(nn.ge.1.and.(vv(1).gt.0..and.vv(1).le.26))then
      i=vv(1)-1
      call osave1(i)
      if(ifl(2).ne.0) then
      close(nsav)
      ifl(2)=0
      endif
      if(idebug.eq.0) then
        close(unit=nnout,status='delete')
      else
        close(unit=nnout,status='keep')
      endif
      close(nnout)
      open(unit=nnout,file='parmout',access='sequential',
     *     status='unknown',form='formatted')
 7005 continue
      read(nnout,'(10a8)',err=7006)tem
      go to 7005
 7006 continue
      backspace (nnout)
      else
      open (unit=nsav,file='savecor',access='sequential',
     *      form='unformatted',status='new')
      endif
      write(nsav) cor,wt0,ngs,bgs
      close(nsav)
      go to 99
c---after zlimit:  lower limit on lagging particles
 2300 continue
      zlimit=vv(1)
      go to 99
c---after cfield:  read fields from sfo7 as written on tape7
c--- first parameter on cfield is starting cell number
c--- next parameter is ending cell number and third parameter is
c--- increment of cell number as in a do loop
c---line after cfield is filename of sfo7 data
 2800 continue
C modif glm 29/10/2012      read(nin,'(a12)')filesuper
      read(nin,*) filesuper
      write(nnout,'(1x,a12)')filesuper
      if(nn.lt.1)then
      nc=1
      ne=1
      ns=1
      elseif(nn.eq.1)then
      nc=vv(1)
      ne=nc
      ns=1
      elseif(nn.eq.2)then
      nc=vv(1)
      ne=vv(2)
      ns=1
      elseif(nn.ge.3)then
      nc=vv(1)
      ne=vv(2)
      ns=vv(3)
      endif
      call cellfeld(filesuper,nc,ne,ns)
      go to 99
c---after design
 3100 continue
      call design
      go to 99
c---after pipe, vv(n) = radius,zlow,zhigh
 3200 continue
      if (npipe.eq.100) then
        write(nnout,9014)
 9014   format (' Too many PIPE cards')
        go to 99
      endif
      if (vv(1).le.0.) then
        write(nnout,9015)
 9015   format (' Bad radius on PIPE card')
        go to 99
      endif
      npipe = npipe + 1
      rpipe(npipe) = vv(1)
      zlpipe(npipe) = vv(2)
      zhpipe(npipe) = vv(3)
      if (ip.eq.1 .or. ip.eq.2)
     .write(nnout,9016) npipe,vv(1),vv(2),vv(3)
 9016 format (
     . ' Beam pipe for point-by-point image calculation'/
     . ' pipe no. = ',t40,i8/
     . ' pipe radius = ',t40,1pg12.3,' cm'/
     . ' zlow = ',t40,g12.3,' cm'/
     . ' zhigh = ',t40,g12.3,' cm')
      go to 99
c---after foclal vv(n)=zmin,zmax,dzz,npchb,facbz,opt
 3300 continue
      iopt=vv(6)
      if(iopt.eq.0) then
         ifld=0
         ifoclal=0
         read(nin,'(a12)') filebz
         write(nnout,*) 'foclal not used'
      else
         ifld=1
         ifoclal=1
         read(nin,'(a12)') filebz
         call fieldlal(filebz)
      endif
      go to 99
c---after backb back bombardment 
 3400 continue
      iback=1
      bmax=vv(1)
      byz0=vv(2)
      byzs=vv(3)
      iope=0
      if(nn.gt.3) then
           pzowmin=vv(4)
           pzowmax=vv(5)
      else
           pzowmin=0.
           pzowmax=0.
      endif
      if(pzowmax.ne.0.) iope=1
      call opeback(iope)
      go to 99
c---after wiggler l,aper,iout,periods,steps,b0,kx,ky,kw,iprint 
 3500 continue
      if(nn.lt.7) go to 7002
      xlam=el(1,nel)/el(4,nel)
      cayw=twopi/xlam
      cayy=el(7,nel)*cayw
      if(el(7,nel).le.0. .or. el(7,nel).gt.1.)then
         write(nnout,9018)
 9018    format(' error: ky/kw on wiggler must be between 0 and 1')
         jj=100 ! for cread
         go to 99
      endif
      cayx=0.
      cayxsq=cayw**2-cayy**2
      if(cayxsq.gt.0.) cayx=sqrt(cayxsq)
      el(7,nel)=cayx
      el(8,nel)=cayy
      el(9,nel)=cayw
      if(ip.eq.1 .or. ip.eq.2)
     .write(nnout,9019)el(4,nel),el(5,nel),el(6,nel),el(7,nel)
 9019 format(
     .' number of wiggler periods =',t40,g12.3/
     .' step of wiggler =',t40,g12.3/
     .' wiggler field =',t40,g12.3,' gauss'/
     .' ratio of Ky to Kw =',t40,g12.3)
      go to 99
c---after alpha magnet
 3600 continue
      if(ip.eq.1 .or. ip.eq.2)
     .write(nnout,9020) el(1,nel),el(4,nel),el(5,nel),el(7,nel)
 9020 format(
     .' distance travelled by the r.p. =',t40,g12.3,' cm'/
     .' reference particle energy =',t40,g12.3,' MeV'/
     .' alpha magnet angle =',t40,g12.3,' Degrees'/
     .' absolute field gradient =',t40,g12.3,' Gauss/cm')
      go to 99
c---after stat
 3700 continue
      var=vv(1)
      istat=1
      open(unit=nstat,file='parmstat',form='unformatted',
     *     status='unknown')
      write(nstat)ij
      go to 99
c---after Poisson
 3800 continue
C modif glm 29/10/2012      read (nin,'(a12)')filepois
      read (nin,*)filepois
      write(nnout,'(1x,a12)')filepois
      if(nn.ge.1)then
      zoffset=vv(1)
      else
      zoffset=0.
      endif
      if(nn.ge.2)then
      rmult=vv(2)
      else
      rmult=1.
      endif
      call poisson(filepois,zoffset,rmult)
      go to 99
c
cbm retrait de nfieldb 04/02/2010 apres modif de coil
c--- after nfieldb number of points used for the magnetic field
c    < 10000
c 4000 continue
c      nptcbu=vv(1)
c      write(nnout,*)
c     *' Nombre de points maximal pour le champ magnetique ',nptcb
c      write(nnout,*)
c     *' Number of used values for magnetic field',nptcbu
c      if(nptcbu.gt.nptcb) then
c        write(nnout,*)
c     *  ' Number of points for magnetic field gt ',nptcb
c      endif
c      go to 99 
c---after-end output any used buffer
 2200 continue
c---check output buffers.
c---element neo?
      inb=imaa*imbf/npoints
      if(inb.gt.maxbuf)inb=maxbuf
      do 8014 i=1,inb
      if(nbuf(i).eq.0)go to 8014
      neo=neb(i)
c---yes.  empty buffer.
      call output(neo,nbuf(i),outcor(1,(i-1)*npoints+1),npoints)
 8014 continue
      do 8015 j=0,nel
      sum=0.
      sumi=0.
      do 8016 i=ngood+1,npoints
      if(cord(7,i).eq.j)then
         sum=sum+1
         sumi=sumi+(gam(i)-1.)*erest
      endif
 8016 continue
      if(sum.ne.0.)then
         sumi=sumi/sum
         write(nnout,9021) j,sum/npoints,sumi
 9021 format(' fraction lost in element',i3,' = ',f7.5,
     *' average energy =',f9.5,' MeV')
      endif
 8015 continue
      if(ngood.gt.2)then
         sumi=0.
         do 8017 i=1,ngood
         sumi=sumi+(gam(i)-1.)*erest
 8017    continue
         sumi=sumi/ngood
         write(nnout,9022) float(ngood)/npoints,sumi
 9022    format(' fraction good at end of calculation = ',f7.5,
     *   ' average energy = ',f9.5,' MeV ')
      endif
c-----
cbm 18/07/2k      if(lolun1.ne.0) call closehisto
c-----
      write(ndiag,*) '   '
      write(ndiag,*) '  Statistique sur ce run'
      write(ndiag,*) '   '

      call datime(Iend_date,Iend_time)
      
c-with cernlib
c      call timed(t1)
c      CPU_time=t1
c_vms/unix
c_vms
c      write(ndiag,698) 
c     * Istart_date(1:4),Istart_date(5:6),Istart_date(7:8),Istart_time,
c     * Iend_date(1:4),Iend_date(5:6),Iend_date(7:8),Iend_time,
c     * secnds(time_begin),CPU_time
c 698  format(//,' --  on ALPHA VMS -- LAL --',28('-')/
c     * '   Start Execution :  ',a4,'/',a2,'/',a2,' ',a8,/
c     * '   End   Execution :  ',a4,'/',a2,'/',a2,' ',a8,/
c     * '   Time used for calculation = ',f10.2,' [sec]',/
c     * '   CPU Time                  = ',f10.2,' [sec]',/
c     * ' -------------------------------',20('-'))
c_unix
      if (idebug.eq.0) then
      else
        write(ndiag,*) '          ---- Version de DEBUG ----- '
      endif
      write(ndiag,698) kvers,
     *   Istart_date,Istart_time,
     *   Iend_date, Iend_time
c     *   ,secnds(time_begin),CPU_time
 698  format(//,' -- LAL -- ',a,16('-')/
     *          '   Start Execution : ',a9,' ',a8,/
     *          '   End   Execution : ',a9,' ',a8,/
     *          ' -------------------------------',20('-'))
c les 2 lignes supprimees apres suppression de la cernlib
c retrait des appel a la fonction secnds et a la subroutine timed
c     *          '   Time used for calculation = ',f10.2,' [sec]',/
c     *          '   CPU Time                  = ',f10.2,' [sec]',/

c_vms/unix
c------------------------------------------------------------------------
      call appendparm
      if(idebug.eq.0) then
        close(unit=nnout,status='delete')
      else
        close(unit=nnout,status='keep')
        write(*,*) '  '
        write(*,*) ' CLOSE statement keep file parmdebug '
        write(*,*) '  '
      endif
      print *, ' parmela normal end'
      stop ' normal end'
c---error messages
 7000 continue
      write(nnout,9023)
 9023 format(' too many parameters on element card')
      jj=100
      go to 99
 7001 continue
      write(nnout,9024) lmx
 9024 format(' too many elements max is ',i4)
      go to 2200
 7002 write(nnout,9025)
 9025 format(' not enough parameters')
      jj=100
      go to 99
      end

c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*
      subroutine cellfld(cl,nc)
c     cell initialization
c----------------------------------------------------------------------
      save
c
      include 'param_sz.h'
      include 'cfldscom.h'
      include 'constcom.h'
      include 'misccom.h'
      include 'syscom.h'
      include 'ucom.h'
c
      common/jones/ajones,zjones,zcath(imaa)
      common/fcoef/ncoef(100),noff(100),ncer(100)
c--------------------------------------------------------------------------
c*
      hcll(nc)=.5*cl
      beta=2.*cl/wavel
      cellfreq(nc)=freq
      lctype=ltype
      if(nn.ge.9)then
      if(vv(9).ne.0.)cellfreq(nc)=vv(9)
      beta=beta*cellfreq(nc)/freq
      endif
      if(nn.ge.10.and.vv(10).ne.0.)lctype=vv(10)
      phi=el(4,nel)*radian
      el(9,nel)=sin(phi)
      el(10,nel)=cos(phi)
      el(4,nel)=cellfreq(nc)
      do 30 i = 1,numcf
      fc(i,nc) = 0.
   30 continue
      if (ip.eq.1 .or. ip.eq.2) then
        write(nnout,20) lctype,cellfreq(nc),hcll(nc)
   20   format (
     . ' cell type =',t40,i8/
     . ' cell frequency =',t40,1pg12.3,' MHz'/
     . ' "half-cell" length =',t40,g12.3,' cm')
        if (nel.eq.1.and.ajones.gt.0.) then
          write(nnout,21)
   21     format (' rf fields from analytic form of M. Jones')
          return
        endif
      endif
      if(nn.gt.11) then
      np=nn
      ncoef(nel)=vv(12)
      noff(nel)=vv(13)
      ncer(nel)=vv(14)
      nmax=14+ncoef(nel)
      if(np.gt.nmax)np=nmax
      do 10 i=15,nmax
      fc(i-14,nc)=vv(i)
   10 continue
        if (ip.eq.1 .or. ip.eq.2) then
          write(nnout,39) ncoef(nel),noff(nel)
   39     format (' No. of Fourier coefs =',t40,i8/
     .    ' odd/even cosines =',t40,i8)
          write(nnout,40) nmax
   40     format (' Using vv(15 to ',i2,') as Fourier coefficents of'
     .    ,' this cell')
          write(nnout,41) ncer(nel)
   41     format (' No. of Fourier coefs used in er calculation =',
     .    t50,i8)
        endif
        go to 100
      endif
c        get fourier coefficients for special cells from tables
      ncoef(nel) = 14
      noff(nel)  =  1
      ncer(nel)  = 14
      iflag(nc)=1
      if(iflag(nc).ne.0)return
      if(lctype.eq.8.and.abs(hcll(nc)-2.21).lt..01)beta=1.02
      if(lctype.eq.6)call rtmfc(beta,fc(1,nc))
      if(lctype.eq.2)call sccfc(beta,fc(1,nc))
      if(lctype.ne.1.and.lctype.ne.2)call genfc(beta,fc(1,nc))
      return
c         make lookup table for rf field
  100 continue
      call genfld(nc)
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*