Program ProjImp

      IMPLICIT NONE

      INTEGER   i, j, Nx, Nfx, Nmovie, Nmax, Ntime, Nbound, Ncloak,
     /          Nslice, Nuse, Nseg
      PARAMETER ( Nfx = 1600, Nmax = 3201 )
      REAL*8    DT, DX, thetax, thetat, beta, Delta, Time, Timpt, 
     /          Timseg, Q(-1:1,-Nfx:Nfx), Mat(Nmax,Nmax), Vec(Nmax),
     /          Work(Nmax), Cloak(-Nfx:Nfx)
      INTEGER   Ord(Nmax)

      OPEN(7, FILE = 'Project.dat' , FORM = 'formatted'  )
         READ(7, 200, ERR = 1000 ) DT, DX, thetax, thetat, beta, 
     /                             Delta, Time, Timpt, timseg
 200     FORMAT(8(11X,E11.4 / ),11X,E11.4)
         READ(7, 210, ERR = 1000 ) Nmovie, Nbound, Ncloak
 210     FORMAT(2(11X,I5 / ),11X,I5)
      CLOSE(7)

      Nx = NINT(1.0D0/DX)
      Ntime = NINT(Time/DT)
      Nseg = NINT(Timseg/DT)
      Nuse = 2*Nx + Nbound

      IF ( Nmovie .EQ. 1 ) PRINT *, '# of steps = ', Ntime / Nseg
      Nslice = NINT(Timpt/DT)

C Set up your initial conditions

      CALL INITIAL(DT,DX,Nx,Q,Delta,beta,Cloak)
      IF ( Nmovie .EQ. 1 .OR. Timpt .EQ. 0.0D0 ) 
     /   CALL OUTPUT(0,DX,Nx,Cloak)
      CALL MATRIX(DT,DX,thetax,thetat,beta,Nbound,Nuse,Mat,Work,Ord)

C Now evolve the solution over time.

      DO I = 1, Ntime
         CALL EVOL(DT,DX,thetax,thetat,beta,Nbound,Nx,Q,Nuse,Mat,
     /             Vec,Ord,Cloak,Ncloak)

         IF ( Nmovie .EQ. 1 .AND. MOD(I,Nseg) .EQ. 0 ) THEN
            CALL OUTPUT(I/Nseg,DX,Nx,Cloak)
         ELSEIF ( I .EQ. Nslice ) THEN
            CALL OUTPUT(I,Dx,Nx,Cloak)
         ENDIF
      ENDDO
      STOP

 1000 PRINT *, 'Error in reading file'
      STOP
      END

C***********************************************************************

      SUBROUTINE OUTPUT(I,DX,Nx,Cloak)

      IMPLICIT NONE

      INTEGER      I, j, Nx
      REAL*8       DX, Cloak(-Nx:Nx), Qdum
      CHARACTER*64 fname

      IF ( I .LT. 10 ) THEN
         WRITE(fname,901) I
 901     FORMAT('slice','.',I1) 
      ELSEIF ( I .GT. 9 .AND. I .LT. 100 ) THEN
         WRITE(fname,902) I
 902     FORMAT('slice','.',I2) 
      ELSEIF ( I .GT. 99 .AND. I .LT. 1000 ) THEN
         WRITE(fname,903) I
 903     FORMAT('slice','.',I3) 
      ELSEIF ( I .GT. 999 .AND. I .LT. 10000 ) THEN
         WRITE(fname,904) I
 904     FORMAT('slice','.',I4) 
      ELSE
         PRINT *, 'I is too high.'
         STOP
      ENDIF

      OPEN( UNIT = 99, FILE = fname )
      DO j = -Nx, Nx 
         Qdum = Cloak(j)
         IF ( ABS(Qdum) .LT. 1.0D-89 ) Qdum = 0.0D0
         WRITE(99,320) j*DX, Qdum
 320     FORMAT(1X,F9.6,2X,E17.10)
      ENDDO
      CLOSE( UNIT = 99 )

      RETURN
      END