C     Last change:  LMY  29 Oct 2001   11:42 am
c   Subroutine RDTAPE2.  This routine reads TAPE2 data written by
c   program Parmela.
c
c Written by Lloyd M. Young, 
c Los Alamos National Laboratory,
c Los Alamos, New Mexico 87545, USA.
**********************************************************************
* Copyright, 1965-2002, The Regents of the University of California. *
*                                                                    *
*  This software was produced under a U.S. Government contract       *
*  (W-7405-ENG-36) by Los Alamos National Laboratory, which is       *
*  operated by the University of California for the U.S. Department  *
*  of Energy.  Neither the Government nor the University makes any   *
*  warranty, express or implied, or assumes any liability or         *
*  responsibility for the use of this software.                      *
*                                                                    *
*  This code may not be reproduced or distributed without express    *
*  written consent of the Los Alamos Accelerator Code Group.  It     *
*  may not be rented or leased to others.                            *
**********************************************************************
c
c  This code requires fortran 90 unless it is rewritten so that it does
c  not use the module  readtapeArrays.
c  If you are a programmer writing your own post processors for Parmela
c  you may include this Lahey FORTRAN source file in your own code.
c  You are free to modify this file as needed for your own use.
c  You must use the -nhed option when compiling this code to read
c  the files written by Parmela.
c
c  This file is provided as is without additional documentation except
c  that provided herein.
c
c  Los Alamos National Laboratory cannot provide consulting for source-
c  code developers.
c
      subroutine rdtape2(num,ipt,ieof)
      USE readtapeArrays
      implicit none
      integer*4 num,i,ieof,ipt !,ifile
c
c
c  This routine reads two binary output files written by program
c  Parmela.  The files, whose default names in Parmela are TAPE2
c  or TAPE3 contain the coordinates of particles at the end of
c  beam-line elements (TAPE2) or the coordinates of particles at
c  the end of specified time steps (TAPE3).  The actual file names
c  are controlled by settings in the [Parmela] section of file
c  LANL.INI.
c
c  The file to read is specified by calling parameter ELEMENTOUTNAME
c  and the value of NUM indicates whether the file is a TAPE2 or TAPE3
c  file.
c
c  If NUM is not equal to 5, then the file is of type TAPE2 and
c  coordinates for element number NE will appear in array CORD.
c  The coordinates and other data at each element are stored in array
c  CORD as follows:
c
c  CORD(1,I)=x
c  CORD(2,I)=bgx
c  CORD(3,I)=y
c  CORD(4,I)=bgy
c  CORD(5,I)=Phase (degree)
c  CORD(6,I)=Energy (MeV)
c  CORD(7,I)=original index of the particle in the COR array
c
c  Note that CORD(7,I) is always an integer number even though it is
c  stored in a real array as a real number.
c
c  If NUM is equal to 5, then the file is of type TAPE3 and
c  coordinates for time step number NE will appear in array CORD.
c  The coordinates at each time step are stored in array CORD as
c  follows:
c
c  CORD(1,I)=x
c  CORD(2,I)=y
c  CORD(3,I)=z
c  CORD(4,I)=original index of the particle in the COR array
c  CORD(5,I)=Kinetic Energy (keV)
c
c  Note that CORD(4,I) is always an integer number even though it is
c  stored in a real array as a real number.
c
c for out5 time steping
c ipt=-1 means process all particles
c ipt=0  means other particles i.e. cord(5*n-1)=0
c ipt= 1 means process only beam particles i.e. cord(5*n-1).ne.0
c cord (1-3) is x,y,z : cord(4) is orginal particle number :cord (5) is KE.
c use orginial particle number to identify particles type and charge.
c
c
      integer*4 irun,ngood,npoints,ltype,istat,j!,iyear,month,iday
      integer*4 nsize,nsteps,nsc,nout,left,ibigin,irec,numfile,numf
      REAL*4 ac(16),cc(13) ,rfsize
       common/ncords/npoints,ngood
      real*4 pi,twopi,radian,clight,erest(3),brhof,freq,
     *             wavel,z0,w0,beami,dcon,zlimit
      common/const/pi,twopi,radian,clight,erest,brhof,freq,
     *             wavel,z0,w0,beami,dcon,zlimit,ltype
      character title*80
      common/char/title
      common/misc/irun
      real*8 dwt,wt0,dum,swt,cwt
      common/tstep/dwt,wt0,dum,swt,cwt,nsteps,nsc,nout
      equivalence (ac(1),pi), (cc(1),dwt)
      CHARACTER (LEN=4) :: line
      CHARACTER (LEN=80) :: tempfilename
      LOGICAL timeopen,fileexist
c
      if (.not. ALLOCATED(map) ) then
      inquire (FILE=ElementOutName,FLEN=ifsize(1),EXIST=fileexist)
c      call iosdirinfo(' ',TRIM(ElementOutName),filename,numfil,ifdate,
c     &               ifsize)
c      IF (NUMFIL.NE.0) THEN
c          DO IFILE = 1,NUMFIL
c            CALL IOsFileDate(IFDATE(IFILE),IYEAR,MONTH,IDAY)
c          END DO
c      END IF
      IF(.not.fileexist)then
      WRITE(*,*)' file ',elementoutname,' does not exist.'
      STOP ' file does not exist'
      endif
      open (2,FILE=TRIM(ElementOutName),ERR=6,ACCESS='transparent',
     & STATUS='old', ACTION='read,denynone')
      do i=1,4
      READ(2,rec=i) line(i:i)
      end do
      ireclen=TRANSFER(line,ireclen)
      do i=1,4
      READ(2,rec=i+4) line(i:i)
      end do
      numrec=TRANSFER(line,numrec)
      do i=1,4
      READ(2,rec=i+8) line(i:i)
      end do
      npoints=TRANSFER(line,npoints)
      nplots=(ifsize(1)/ireclen)/numrec -1
      do i=1,4
      READ(2,rec=i+12) line(i:i)
      end do
      irun=TRANSFER(line,irun)
      close (2)
      open(2,file=TRIM(ElementOutName),ERR=7,ACCESS='direct',
     & RECL=ireclen , ACTION='read,denynone')
      lastifsize(1)=ifsize(1)
c
      ALLOCATE(cor(6,npoints),cord(7,npoints),iicol(0:npoints),
     # ip(npoints),icol(npoints), STAT=istat)
      if (istat.ne.0) then
      STOP 'Not enough memory for number of particles'
      end if
      nels=MAX(1000,irun)
      ALLOCATE  (map(0:nels),zloc(0:nels),
     * el(11,0:nels), STAT=istat)
      if (istat.ne.0) then
      STOP'Not enough memory for number of elements specfied in Parmela
     *input file'
      end if
c
          READ(2,rec=1,ERR=14)ireclen,numrec,npoints,irun,title,
     &     ac,nel,cc,ngood
        left=(ireclen/4-57-3*(nel+1))/7   ! nxt write has 52 *4 bytes other than zloc, el and cor(
        ibigin=1
       do irec=1,numrec
         IF(numrec.eq.1)then
          READ(2,rec=irec,ERR=14)ireclen,numrec,npoints,irun,title,
     &     ac,nel,cc,ngood,(zloc(i),el(2,i),map(i),i=0,nel),
     * ((cor(i,j),i=1,6),j=1,npoints),(ip(j),j=1,npoints),
     * (iicol(j),j=1,npoints)
         elseif (irec.eq.1)then
          read(2,rec=irec,ERR=14)ireclen,numrec,npoints,irun,title,
     &     ac,nel,cc,ngood,(zloc(i),el(2,i),map(i),i=0,nel),
     * ((cor(i,j),i=1,6),j=ibigin,left),(ip(j),j=ibigin,left),
     * (iicol(j),j=ibigin,left)
          ibigin=left+1
          left=MIN(left+(ireclen/4)/7,npoints)
         else
         IF(ibigin.gt.left)exit
         read(2,rec=irec,ERR=14)
     * ((cor(i,j),i=1,6),j=ibigin,left),(ip(j),j=ibigin,left),
     * (iicol(j),j=ibigin,left)
          ibigin=left+1
          left=MIN(left+(ireclen/4)/7,npoints)
         endif
       end do
c check for multiple tape2 files and open if they exist.
      maxj=1
      do i=1,nel
      if (map(i).ne.0) then
      maxj=maxj+1
      end if
      end do
       rfsize=REAL(maxj)*REAL(numrec)*REAL(ireclen)
      numfile=rfsize/2.09e9 +1
      maxj=rfsize/(REAL(numrec)*REAL(ireclen))/numfile+1
c
      numfile=MIN(numfile,10)
      do numf=2, numfile
        tempfilename=TRIM(ElementOutName)//CHAR(numf+63)
      open(itape35(numf),file=TRIM(tempfileName),ERR=8,ACCESS='direct',
     & RECL=ireclen , ACTION='read,denynone')
      inquire (FILE=tempfileName,FLEN=ifsize(1),EXIST=fileexist)
c      call iosdirinfo(' ',TRIM(tempfileName),filename,numfil,ifdate,
c     &               ifsize)
      lastifsize(numf)=ifsize(1)
      end do
c
      go to 15
   6  write(*,*)' error opening file ',ElementOutName
      stop 'error opening file'
   8  write(*,*)' error opening file ',tempfileName
      stop 'error opening file'
   14 write(*,*)' error reading cor !!!'
      STOP 'error reading cor'
   7  WRITE(*,"(' error opening file ',a,
     #         ' for direct access,',/,' The file exists but is not a '
     #        ,'valid Parmela V3 tape2 type file.')") ElementOutName
      STOP 'bad "tape2" file'
      end if
   15 continue
c
      if(num.eq.5)then
        inquire (4,OPENED=timeopen)
        if (.not.timeopen) then
         open(4,file=TimeStepOutName,status='OLD',err=151,
     *   form='unformatted',action='read,denynone')
        end if
      endif
c
c  The IP array contains particle numbers (1, 2, or 3) defined on the
c  CHARGE line. The format of the CHARGE line is:
c  CHARGE Charge, ParticleType
c  where Charge is an interger from -127 to +127 and ParticeleType is 1, 2, or 3
c  The charge spcifies the particle type for the following INPUT line in Parmela
c  input file. a different mass maybe associated with each ParticleType
c
c  Program PARGRAF uses a different color for each of the six particle
c  types.  Multiply values in the IP array by 4 subtract 1 and subtract the
c  corresponding value from array IC, to get the following color
c  numbers for the particles.
c
c   ColorNumber       IP    IC
c     2               1    +1
c     4               1    -1
c     6               2    +1
c     8               2    -1
c     10              3    +1
c     12              3    -1
c 
c
      nsize=npoints
      call readbuf(ieof,nsize,num,ipt)
       if(ieof.ne.0 )then
       if (num.eq.5) then
       WRITE(*,*)'time step',ne,'not in file' 
       else
       WRITE(*,*) 'element',ne,' requested not in file'
c if ieof < zero element requested not in map of elements
c if ieof > zero element requested in map but not in file.
       end if
       endif
      return
  151 WRITE(*,*) 'The time-step file does not exist.'
      return
      end
c
      subroutine readbuf(ieof,nsize,num,ipt)
      USE readtapeArrays
      common/ncords/npoints,ngood
      INTEGER*1 ic(npoints)
      integer*4 nsteps,nsc,nout
      real*8 dwt,wt0,dum,swt,cwt
      common/tstep/dwt,wt0,dum,swt,cwt,nsteps,nsc,nout
      real temp
      CHARACTER (LEN=80) :: tempfilename
      if(num.eq.5)then
c for out5 time steping
c ipt=-1 means process all particles
c ipt=0  means other particles i.e. cord(4.n)=0
c ipt= 1 means process only beam particles i.e. cord(4,n).ne.0
c cord (1-3) is x,y,z : cord(4) is particle number :cord (5,n) is KE.
c orginial beam particles have particle number.
      ieof=0
      read (4,err=30,iostat=ieof,end=30)nstep,nbuf,wt,
     * ((cord(i,j),i=1,5),j=1,nbuf),(ic(i),i=1,nbuf)
     *  ,(ip(i),i=1,nbuf),
     &   ax,bx,xbx,ybx,exrms,xmaxrms,
     &   ay,by,xby,yby,eyrms,ymaxrms,
     &   az,bz,xbz,ybz,ezrms,zmaxrms
      nn=0
      do n=1,nbuf
      if(cord(4,n).eq.0. )then    !particle number 0
        IF(ipt.le.0)then
         nn=nn+1
         icol(nn)=4*ip(n)-ic(n)-1
         do j=1,5
         cord(j,nn)=cord(j,n)
         end do
        endif
      ELSEIF (cord(4,n)-1. .EQ. 0.) THEN !reference particle cord(4,n)=1.
        IF(ipt.ne.0)then
         nn=nn+1
c      put reference particle in location 1 and move particle 1 to nn
c      if ipt<1 all particles are being processed and nn=n so save it first.
         do j=1,5
         temp=cord(j,nn)
         cord(j,nn)=cord(j,1)
c         now put reference particle in location 1
         cord(j,1)=cord(j,n)
         cord(j,n)=temp
         end do
         icol(nn)=icol(1)
         icol(1)=4*ip(n)-ic(n)
        endif
      else
        IF(ipt.ne.0)then
         nn=nn+1
         icol(nn)=4*ip(n)-ic(n) !normal beam particle
         do j=1,5
         cord(j,nn)=cord(j,n)
         end do
        endif
      endif
      enddo
      nbuf=nn
      return
      endif
      do i=1,nel
      if (map(i).eq.ne) then
      nnx=i
      GO TO 1
      end if
      end do
      ieof=-1
      return
  1    numf=nnx/maxj +1     !determin which file
       if (numf.gt.10) then
       ieof=1
       return
       end if
       irecs=(nnx-(numf-1)*maxj)*numrec+1
       irece=irecs+numrec-1
       itry=0
       IF(numf.gt.1)then
       tempfilename=TRIM(elementOutName)//CHAR(numf+63)
       else
       tempfilename=elementOutName
       endif
 19    inquire (FILE=TRIM(tempfileName),FLEN=ifsize(1))
c      call iosdirinfo(' ',TRIM(ElementOutName),filename,numfil,ifdate,
c     &               ifsize)
      nplots=((ifsize(1))/ireclen)/numrec
      if (ifsize(1).ne.lastifsize(numf)) then
      CLOSE(itape35(numf))
      open(itape35(numf),file=TRIM(tempfileName),ERR=20,ACCESS='direct',
     & RECL=ireclen , ACTION='read,denynone')
      end if
      lastifsize(numf)=ifsize(1)
        if (nnx-(numf-1)*maxj.gt.nplots) then
         ieof=1
         return
        else
         ieof=0
        end if
c
        READ(itape35(numf),rec=irecs,ERR=30) nbuf
        left=MIN((ireclen/4-4)/8,nbuf)
        ibigin=1
c
        itry=0
        do irec=irecs,irece
  299     IF(numrec.eq.1)then
          READ(itape35(numf),rec=irec,ERR=300)  nbuf,neo,pr,wr,
     *                      ((cord(i,j),i=1,7),j=1,nbuf),
     *                      (ic(i),i=1,nbuf)
         elseif (irec.eq.irecs)then
          READ(itape35(numf),rec=irec,ERR=300) nbuf,neo,pr,wr,
     *                   ((cord(i,j),i=1,7),j=ibigin,left),
     *                     (ic(i),i=ibigin,left)
          ibigin=left+1
          left=MIN(left+(ireclen/4)/8,nbuf)
         else
         if(ibigin.gt.left)exit
         READ(itape35(numf),rec=irec,ERR=300)
     *                  ((cord(i,j),i=1,7),j=ibigin,left),
     *                     (ic(i),i=ibigin,left)
          ibigin=left+1
          left=MIN(left+(ireclen/4)/8,nbuf)
         endif
         itry=0
         cycle
 300    itry=itry+1
        if (itry.lt.2) then
         GO TO 299
         else
         GO TO 30
        end if
        end do
      if (neo.eq.ne) then
      do n=1,nbuf
      icol(n)=iicol(int(cord(7,n)))-ic(n)
      enddo
       return
      else
       write (*,*)" neo and ne disagree, neo,ne ",neo,ne
       return
      end if
c
   30 ieof=1
       return
   20 if (itry.lt.2) then
      itry=itry+1
         READ(2,rec=numrec+1,ERR=30) nbuf
      GO TO 19
      else
      ieof=1
      return
      end if
      end



