PROGRAM ptmass
c
c--Takes a BINARY output 'P file' and outputs a text file for each
c     sink particle that contains its mass versus time
c     This code also re-orders sinks for MPI dumps according to iunique
c
      IMPLICIT DOUBLE PRECISION (D)
      IMPLICIT INTEGER (I-N)
      IMPLICIT REAL (A-C, E-H, O-Z)

      PARAMETER (idim = 200)
      PARAMETER (nlines = 1000000)

      DIMENSION listmerge(idim)
      DIMENSION pmasses(idim,nlines),fftimes(nlines)

      INTEGER*8 iunique, listunique(idim)
      LOGICAL*1 ifound(idim)

      CHARACTER*21 infile, outfile
      CHARACTER*19 outbase
      CHARACTER*3 num

      WRITE (*,*) 'Enter input filename:' 
      READ (*,99001) infile
99001 FORMAT(A21)

      WRITE (*,*) 'Enter base name for output files:' 
      READ (*,99002) outbase
99002 FORMAT(A19)
c
c--Compute
c
      ilines = 0
      nptmass = 0
      nptmasscurrent = 0
      nmerge = 0
      oldtime = 0.

      DO i = 1, idim
         listmerge(i) = i
      END DO

      OPEN (15, FILE=infile, STATUS='unknown', FORM='unformatted')
 100  CONTINUE
c
c--Read time and number of point masses
c
 110  nptmassold = nptmass
 111  READ (15, END=200, ERR=200) 
     &     fftime, time, nptmass

      IF (time.LT.oldtime) THEN
         WRITE (*,*) 'ERROR - time.LT.oldtime ',time,oldtime
         STOP
      ELSE
         oldtime = time
      ENDIF

c      write (*,*) fftime, time, nptmass

      IF (nptmass.LT.0) THEN
         READ (15) imerged1, imerged2
         WRITE (*,*) 'Merged ',imerged1, imerged2
         nmerge = nmerge + 1
         nptmassold = nptmassold - 1
         DO i = imerged2, idim
            listmerge(i) = listmerge(i+1)
         END DO
         GOTO 111
 112     nptmass = nptmassold
         GOTO 200
      ENDIF

      ifindmerge = 0
      IF (nptmassold .GT. nptmass) THEN
         WRITE (*,*) 'MERGER not done via -ve nptmass ',
     &        nptmass, nptmassold, fftime, time
         STOP
         nptmassold = nptmassold - 1
         ifindmerge = 1
         DO i = 1, nptmasscurrent
            ifound(i) = .FALSE.
         END DO
      ENDIF

      IF (nptmass+nmerge.GT.idim) THEN
      WRITE (*,*) 'GREATER THAN ',idim,' POINT MASSES ',nptmass,nmerge
         STOP
      ENDIF
c
c--Read each time dump
c
      ilines = ilines + 1
      IF (ilines.GT.nlines) THEN
         WRITE (*,*) 'ERROR - number of lines too big'
         STOP
      ENDIF
      fftimes(ilines) = fftime
      DO k = 1, nptmass
         READ (15, END=200, ERR=200) iunique, 
     &        xi, yi, zi, vxi, vyi, vzi, pmassi, 
     &        rho4, nactotal2, ptmass4, 
     &        spinxi, spinyi, spinzi

         IF (iunique.GT.2000000000) THEN
            WRITE (*,*) 'ERROR - iunique > 2^31 ',iunique
            STOP
         ENDIF
c
c--Search for iunique in existing sink list
c
         DO iii = 1, nptmasscurrent
            IF (iunique.EQ.listunique(iii)) THEN
               i = listmerge(iii)
               ifound(iii) = .TRUE.
               GOTO 555
            ENDIF
         END DO
         nptmasscurrent = nptmasscurrent + 1
         WRITE (*,*) 'New sink found: ',iunique,' assigned ',
     &        nptmasscurrent,' at time ',fftime
         listunique(nptmasscurrent) = iunique
         i = listmerge(nptmasscurrent)

 555     CONTINUE

         pmasses(i,ilines) = pmassi
 
      END DO
c
c--If merger identified through iunique
c
      IF (ifindmerge.EQ.1) THEN
         DO imerged = 1, nptmasscurrent
            IF (.NOT.ifound(imerged)) THEN
               WRITE (*,*) 'Merged ',imerged
               nmerge = nmerge + 1
               DO i = imerged, idim-1
                  listmerge(i) = listmerge(i+1)
               END DO
               
               GOTO 876
            END IF
         END DO
      ENDIF
 876  CONTINUE
c
c--Read next time
c
      GOTO 100
c
c--End of P-file
c
 200  CLOSE(15)

c
c--For each point mass, output evolution file
c
      DO i = 1, nptmass
         IF (i.GT.99) THEN
            WRITE (num, 88000) i
88000       FORMAT (I3)
         ELSEIF (i.GT.9) THEN
            WRITE (num, 88001) i
88001       FORMAT ('0',I2)
         ELSE 
            WRITE (num, 88002) i
88002       FORMAT ('00',I1)
         ENDIF
         DO junk = LEN(outbase), 1, -1
            IF (outbase(junk:junk).NE.' ') THEN
               iplace = junk
               GOTO 150
            ENDIF
         END DO
 150     outfile = outbase(1:junk) // num

         OPEN (23, FILE=outfile, POSITION='APPEND')
         DO j = 1, ilines
            IF (pmasses(i,j).GT.0.0) THEN
               WRITE(23, 99005) fftimes(j), pmasses(i,j)
99005          FORMAT (2(1PE15.8,1X),I10)
            ENDIF
         END DO
         CLOSE (23)
      END DO

      STOP
      END