Program Project1

      INTEGER   Ndx, Nfx, Ngraph, Nbound, Nmovie, Nrat
      PARAMETER ( Nfx = 1600 )
      REAL*8    DT, DX, Delta, Time, Timpt, Q(-1:1,-Nfx:Nfx),
     /          Qf(-1:1,-Nfx:Nfx)
      OPEN(7, FILE = 'Project01.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)

C Initialize all Data. 
      
      Ndx = INT(1.0D0/DX)
      IF ( Ndx .GT. Nfx ) THEN
         PRINT *, 'Your grid is too fine.'
         STOP
      ENDIF

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

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

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

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

      Nrat = Nfx / Ndx
      IF ( MOD(Nfx,Ndx) .NE. 0 ) THEN
         PRINT *, 'You are not using an acceptable grid.'
         STOP
      ENDIF

C Set your initial conditions

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

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

      DO j = -Ndx, Ndx
         Q(1,j) = Qf(1,Nrat*j)
      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 ( Qdum .LT. 1.0E-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 calculate your backwards step by evolving the finest grid.

 400  DO j = -Nfx, Nfx
         Qf(0,j) = Qf(1,j)
      ENDDO

      IF ( Nrat .EQ. 1 ) GOTO 500

      DO i = 1, Nrat-1
 
         DO j = -Nfx, Nfx
            Qf(-1,j) = Qf(0,j)
            Qf(0,j) = Qf(1,j)
         ENDDO
         DO j = -Nfx+1, Nfx-1
            Qf(1,j) = 2.0D0*Qf(0,j) - Qf(-1,j) + ( Qf(0,j+1) 
     /                - 2.0D0*Qf(0,j) + Qf(0,j-1) )*( DT/DX )**2
         ENDDO
         IF ( Nbound .EQ. -1 ) THEN
            Qf(1,Nfx) = 2.0D0*Qf(0,Nfx) - Qf(-1,Nfx) + ( Qf(0,-Nfx+1) 
     /                  - 2.0D0*Qf(0,Nfx) + Qf(0,Nfx-1) )*( DT/DX )**2
            Qf(1,-Nfx) = Qf(1,Nfx)
         ELSEIF ( Nbound .EQ. 0 ) THEN
            Qf(1,Nfx) = 0.0D0
            Qf(1,-Nfx) = 0.0D0
         ELSEIF ( Nbound .EQ. 1 ) THEN
            Qf(1,Nfx) = ( 4.0D0*Qf(0,Nfx) - Qf(-1,Nfx) + ( DT/DX )
     /                  *( 4.0D0*Qf(1,Nfx-1) - Qf(1,Nfx-2) ) ) /
     /                  ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
            Qf(1,-Nfx) = ( 4.0D0*Qf(0,-Nfx) - Qf(-1,-Nfx) + ( DT/DX )
     /                   *( 4.0D0*Qf(1,-Nfx+1) - Qf(1,-Nfx+2) ) ) /
     /                   ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
         ELSE
            PRINT *, 'Invalid value of Nbound' 
            STOP
         ENDIF
      ENDDO

C Step #(Nrat-1) of the finest grid will be #(-1) of the current grid.

 500  DO j = -Ndx, Ndx
         Q(0,j) = Qf(1,Nrat*j)
      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 ( Qdum .LT. 1.0E-89 ) Qdum = 0.0D0
               WRITE(99,220) j*DX, Q(1,j)
            ENDDO
         close(unit=99)

 600  ENDDO

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