Program Project09

      IMPLICIT NONE

      INTEGER   i, Nb, Nf, Nfx, Nmovie, Nmax, Ntime, Nslice, Nseg, Ndir
      PARAMETER ( Nfx = 3200, Nmax = 3201 )
      REAL*8    DT, DX, M, thetarr, thetatt, thetar, thetat, Delta, 
     /          Time, Timpt, Timseg, Q(-1:1,0:Nfx), Mat(Nmax,Nmax),
     /          Vec(Nmax), Work(Nmax), test
      INTEGER   Ord(Nmax), Ord2(Nmax)

      OPEN(7, FILE = 'Project09.dat' , FORM = 'formatted'  )
         READ(7, 200, ERR = 1000 ) DT, DX, M, thetarr, thetatt, thetar,
     /                             thetat, Delta, Time, Timpt, Timseg
 200     FORMAT(10(11X,E11.4 / ),11X,E11.4)
         READ(7, 210, ERR = 1000 ) Nmovie, Ndir
 210     FORMAT(1(11X,I5 / ),11X,I5)
      CLOSE(7)

      Nf = NINT(4.0D0/DX)
      Ntime = NINT(Time/DT)
      Nseg = NINT(Timseg/DT)
      IF ( M .GT. 0.5 .AND. M .LT. 2.5 ) THEN
         Nb = INT( ( 2.0D0*M - 1.0D0 )/DX ) + 1
         test = ABS( NINT( ( 2.0D0*M - 1.0D0 )/DX ) 
     /               - ( 2.0D0*M - 1.0D0 )/DX )
         IF ( test .LE. 1.0D-7 ) Nb = NINT( ( 2.0D0*M - 1.0D0 )/DX )
      ELSE
         Nb = 0
      ENDIF

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

C Set up your initial conditions

      CALL INITIAL(DT,DX,M,Nb,Nf,Q,Delta,Ndir)
      IF ( Nmovie .EQ. 1 .OR. Timpt .EQ. 0.0D0 ) 
     /   CALL OUTPUT(0,DX,Nf,Q)

      CALL MATRIX(DT,DX,M,thetarr,thetatt,thetar,thetat,Nb,Nf+1,Mat,
     /            Work,Ord)
      DO I = 1, Ntime
         CALL EVOL(DT,DX,M,thetarr,thetatt,thetar,thetat,Nb,Nf,Q,
     /             Nf+1,Mat,Vec,Ord)
         IF ( Nmovie .EQ. 1 .AND. MOD(I,Nseg) .EQ. 0 ) THEN
            CALL OUTPUT(I/Nseg,DX,Nf,Q)
         ELSEIF ( I .EQ. Nslice ) THEN
            CALL OUTPUT(I,Dx,Nf,Q)
         ENDIF
      ENDDO
      STOP

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

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

      SUBROUTINE OUTPUT(I,DX,Nx,Q)

      IMPLICIT NONE

      INTEGER      I, j, Nx
      REAL*8       DX, Q(-1:1,0: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 = 0, Nx 
         Qdum = Q(1,j)
         IF ( ABS(Qdum) .LT. 1.0D-89 ) Qdum = 0.0D0
         WRITE(99,320) j*DX + 1.0D0, Qdum
 320     FORMAT(1X,F9.6,2X,E17.10)
      ENDDO
      CLOSE( UNIT = 99 )

      RETURN
      END