Program Project02

      INTEGER   Ndx, Nfx, Ngraph, Nbound, Nmovie, Nrat
      PARAMETER ( Nfx = 1600 )
      REAL*8    DT, DX, Delta, Time, Timpt, Q(-1:1,-Nfx:Nfx)

      OPEN(7, FILE = 'Project02.dat' , FORM = 'formatted'  )
         READ(7, 200, ERR = 1000 ) DT, DX, Delta, Time, Timpt
 200     FORMAT(4(11X,E10.4 / ),11X,E10.4)
         READ(7, 210, ERR = 1000 ) Ngraph, Nmovie, Nbound
 210     FORMAT(2(11X,I5 / ),11X,I5)
      CLOSE(7)

      Ndx = INT(1.0D0/DX)
      IF ( DT .GE. DX ) THEN
         PRINT *, 'DT is too large.'
         STOP
      ENDIF

      CALL BODY(DT,DX,Delta,Time,Timpt,Ngraph,Nmovie,Nbound,Q,Ndx)
      STOP

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

C***********************************************************************
      
      SUBROUTINE BODY(DT,DX,Delta,Time,Timpt,Ngraph,Nmovie,Nbound,Q,Ndx)

      INTEGER   i, j, Ngraph, Nbound, Nmovie, Ndx, Nnon0, Ntime, Nslice
      REAL*8    DT, DX, Delta, Q(-1:1,-Ndx:Ndx), Qdum, Time, Timpt, pi
      PARAMETER ( pi = 3.141593D0 )
      CHARACTER*64 fname

C Set your initial conditions

      Nnon0 = INT(Ndx*0.2D0)
      DO j = -Ndx, -Nnon0-1
         Q(1,j) = 0.0D0
      ENDDO 
      
      DO j = -Nnon0, Nnon0
         IF ( Ngraph .EQ. 1 ) THEN
            Q(1,j) = COS(pi/Delta*(DX*j))/2.0D0 + 0.5D0
         ELSE
            Q(1,j) = 0.04D0*EXP( -(DX*j)**2 / Delta**2 )
         ENDIF
      ENDDO 

      DO j = Nnon0+1, Ndx
         Q(1,j) = 0.0D0
      ENDDO 

      IF ( Nmovie .NE. 0 .OR. Timpt .EQ. 0.0D0 ) THEN
         open(unit=98,file='slice.0')
            DO j = -Ndx, Ndx 
               Qdum = Q(1,j)
               IF ( ABS(Qdum) .LT. 1.0D-89 ) Qdum = 0.0D0
               WRITE(98,220) j*DX, Qdum
 220           FORMAT(1X,F9.6,2X,E17.10)
            ENDDO
         close(unit=98)
      ENDIF

C Now determine your backwards step using the analytical solution.

 400  DO j = -Ndx, -Nnon0-1
         Q(0,j) = 0.0D0
      ENDDO

      DO j = -Nnon0, Nnon0-1
         IF ( Ngraph .EQ. 1 ) THEN
            Q(0,j) = COS(pi/Delta*(DX*j + DT))/2.0D0 + 0.5D0
         ELSE
            Q(0,j) = 0.04D0*EXP( -(DX*j + DT)**2 / Delta**2 )
         ENDIF
      ENDDO 

      DO j = Nnon0, Ndx
         Q(0,j) = 0.0D0
      ENDDO

C Now evolve the current grid.

      Ntime = INT(Time/DT)
      IF ( Nmovie .EQ. 1 ) PRINT *, 'Ntime = ', Ntime
      DO i = 1, Ntime
         DO j = -Ndx, Ndx
            Q(-1,j) = Q(0,j)
            Q(0,j) = Q(1,j)
         ENDDO
         DO j = -Ndx+1, Ndx-1
            Q(1,j) = 2.0D0*Q(0,j) - Q(-1,j) + ( Q(0,j+1) 
     /               - 2.0D0*Q(0,j) + Q(0,j-1) )*( DT/DX )**2
         ENDDO
         IF ( Nbound .EQ. -1 ) THEN
            Q(1,Ndx) = 2.0D0*Q(0,Ndx) - Q(-1,Ndx) + ( Q(0,-Ndx+1) 
     /                 - 2.0D0*Q(0,Ndx) + Q(0,Ndx-1) )*( DT/DX )**2
            Q(1,-Ndx) = Q(1,Ndx)
         ELSEIF ( Nbound .EQ. 0 ) THEN
            Q(1,Ndx) = 0.0D0
            Q(1,-Ndx) = 0.0D0
         ELSEIF ( Nbound .EQ. 1 ) THEN
            Q(1,Ndx) = ( 4.0D0*Q(0,Ndx) - Q(-1,Ndx) + ( DT/DX )
     /                 *( 4.0D0*Q(1,Ndx-1) - Q(1,Ndx-2) ) ) /
     /                 ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
            Q(1,-Ndx) = ( 4.0D0*Q(0,-Ndx) - Q(-1,-Ndx) + ( DT/DX )
     /                  *( 4.0D0*Q(1,-Ndx+1) - Q(1,-Ndx+2) ) ) /
     /                  ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
         ELSE
            PRINT *, 'Invalid value of Nbound' 
            STOP
         ENDIF
       
         Nslice = INT(Timpt/DT)
         IF ( Nmovie .EQ. 0 ) THEN
            IF ( i .NE. Nslice ) GOTO 600
         ENDIF

   	 if (i .lt. 10) then
	   write(fname,901) i
         end if
         if (i .gt. 9 .and. i .lt. 100) then
            write(fname,902) i
         end if
         if (i .gt. 99 .and. i .lt. 1000) then
            write(fname,903) i
         end if
         if (i .gt. 999 .and. i .lt. 10000) then
            write(fname,904) i
         end if

         open(unit=99,file=fname)
            DO j = -Ndx, Ndx 
            Qdum = Q(1,j)
            IF ( ABS(Qdum) .LT. 1.0D-89 ) Qdum = 0.0D0
               WRITE(99,220) j*DX, Qdum
            ENDDO
         close(unit=99)

 600  ENDDO

 901  FORMAT('slice','.',I1) 
 902  FORMAT('slice','.',I2) 
 903  FORMAT('slice','.',I3) 
 904  FORMAT('slice','.',I4) 
      RETURN
      END