Program Converge

      INTEGER   Ntime, Np1, Np4, Np2, Ngraph, Nbound, Nx1
      PARAMETER ( Np1 = 400, Np2 = 800, Np4 = 1600 )
      REAL*8    DT, DX, Delta, Q1(-1:1,-Np4:Np4), Q2(-1:1,-Np2:Np2),
     /          Q4(-1:1,-Np1:Np1), Time

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

      Nx1 = INT(1.0D0 / DX )
      IF ( Nx1 .GT. Np1 ) THEN
         PRINT *, 'The coarsest grid is too fine.'
         STOP
      ENDIF

      CALL BODY(DT,DX,Delta,Time,Ngraph,Nbound,Q1,Q2,Q4,Nx1,2*Nx1,4*Nx1)
      STOP

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

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

      SUBROUTINE BODY(DT,DX,Delta,Time,Ngraph,Nbound,Q1,Q2,Q4,Nx1,Nx2, 
     /                Nx4)

      INTEGER   i, j, k, l, Ntime, Nx1, Nx4, Nx2, Nnon0, Ngraph, Nbound
      REAL*8    DT, DX, Delta, Q1(-1:1,-Nx4:Nx4), Q2(-1:1,-Nx2:Nx2),
     /          Q4(-1:1,-Nx1:Nx1), Time, numer, denom, numpt, denpt,
     /          factor, Dfx

C Initialize all Data. 
      
      Nnon0 = INT(Nx4*0.2D0)
      DO j = -Nx4, -Nnon0-1
         Q1(1,j) = 0.0D0
      ENDDO 
      
      Dfx = DX / 4.0D0
      DO j = -Nnon0, Nnon0
         IF ( Ngraph .EQ. 1 ) THEN
            Q1(1,j) = COS(pi/Delta*(Dfx*j))/2.0D0 + 0.5D0
         ELSE
            Q1(1,j) = 0.04D0*EXP( -(Dfx*j)**2 / Delta**2 )
         ENDIF
      ENDDO 

      DO j = Nnon0+1, Nx4
         Q1(1,j) = 0.0D0
      ENDDO 

      DO j = -Nx2, Nx2
         Q2(1,j) = Q1(1,2*j)
      ENDDO
      DO j = -Nx1, Nx1
         Q4(1,j) = Q1(1,4*j)
      ENDDO

C Take a backwards step.  Then evolve the finest grid.

      DO j = -Nx4, Nx4
         Q1(0,j) = Q1(1,j)
      ENDDO

      OPEN( UNIT = 11, file = 'Factor.dat' )

      Ntime = INT(Time/DT)
      DO i = 1, Ntime

         DO k = 4*i-3, 4*i
            DO j = -Nx4, Nx4
               Q1(-1,j) = Q1(0,j)
               Q1(0,j) = Q1(1,j)
            ENDDO
            DO j = -Nx4+1, Nx4-1
               Q1(1,j) = 2.0D0*Q1(0,j) - Q1(-1,j) + ( Q1(0,j+1) 
     /                   - 2.0D0*Q1(0,j) + Q1(0,j-1) )*( DT/DX )**2
            ENDDO

            IF ( Nbound .EQ. -1 ) THEN
               Q1(1,Nx4) = 2.0D0*Q1(0,Nx4) - Q1(-1,Nx4) + ( Q1(0,-Nx4+1) 
     /            - 2.0D0*Q1(0,Nx4) + Q1(0,Nx4-1) )*( DT/DX )**2
               Q1(1,-Nx4) = Q1(1,Nx4)
            ELSEIF ( Nbound .EQ. 0 ) THEN
               Q1(1,Nx4) = 0.0D0
               Q1(1,-Nx4) = 0.0D0
            ELSEIF ( Nbound .EQ. 1 ) THEN
               Q1(1,Nx4) = ( 4.0D0*Q1(0,Nx4) - Q1(-1,Nx4) + ( DT/DX )
     /                       *( 4.0D0*Q1(1,Nx4-1) - Q1(1,Nx4-2) ) ) /
     /                     ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
               Q1(1,-Nx4) = ( 4.0D0*Q1(0,-Nx4) - Q1(-1,-Nx4) + ( DT/DX )
     /                       *( 4.0D0*Q1(1,-Nx4+1) - Q1(1,-Nx4+2) ) ) /
     /                      ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
            ELSE
               PRINT *, 'Invalid value of Nbound' 
               STOP
            ENDIF

            IF ( k .EQ. 1 ) THEN
               DO j = -Nx2, Nx2
                  Q2(0,j) = Q1(1,2*j)
               ENDDO
            ENDIF
            IF ( k .EQ. 3 ) THEN
               DO j = -Nx1, Nx1
                  Q4(0,j) = Q1(1,4*j)
               ENDDO
            ENDIF
         ENDDO

C Now evolve the other grids.

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

         DO j = -Nx1, Nx1
            Q4(-1,j) = Q4(0,j)
            Q4(0,j) = Q4(1,j)
         ENDDO

         DO j = -Nx1+1, Nx1-1
            Q4(1,j) = 2.0D0*Q4(0,j) - Q4(-1,j) + ( Q4(0,j+1) 
     /               - 2.0D0*Q4(0,j) + Q4(0,j-1) )*( DT/DX )**2
         ENDDO

         IF ( Nbound .EQ. -1 ) THEN
            Q4(1,Nx1) = 2.0D0*Q4(0,Nx1) - Q4(-1,Nx1) + ( Q4(0,-Nx1+1) 
     /                  - 2.0D0*Q4(0,Nx1) + Q4(0,Nx1-1) )*( DT/DX )**2
            Q4(1,-Nx1) = Q4(1,Nx1)
         ELSEIF ( Nbound .EQ. 0 ) THEN
            Q4(1,Nx1) = 0.0D0
            Q4(1,-Nx1) = 0.0D0
         ELSEIF ( Nbound .EQ. 1 ) THEN
            Q4(1,Nx1) = ( 4.0D0*Q4(0,Nx1) - Q4(-1,Nx1) + ( DT/DX )
     /                    *( 4.0D0*Q4(1,Nx1-1) - Q4(1,Nx1-2) ) ) /
     /                  ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
            Q4(1,-Nx1) = ( 4.0D0*Q4(0,-Nx1) - Q4(-1,-Nx1) + ( DT/DX )
     /                     *( 4.0D0*Q4(1,-Nx1+1) - Q4(1,-Nx1+2) ) ) /
     /                   ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
         ELSE
            PRINT *, 'Invalid value of Nbound' 
            STOP
         ENDIF

C Now calculate the convergence factor

         numer = 0.0D0
         denom = 0.0D0
         DO j = -Nx4, Nx4
            IF ( MOD(j,2) .EQ. 0 ) THEN
               denpt = ( Q2(1,j/2) - Q1(1,j) )**2
               denom = denom + denpt
            ENDIF
               IF ( MOD(j,4) .EQ. 0 ) THEN
               numpt = ( Q4(1,j/4) - Q2(1,j/2) )**2
               numer = numer + numpt
            ENDIF
         ENDDO

         numer = SQRT( numer / ( 2*Nx1 + 1 ) )
         denom = SQRT( denom / ( 2*Nx2 + 1 ) )
         factor = numer / denom

         WRITE(11,500) DT*i, factor
 500     FORMAT(1X,F9.6,2X,1P,E13.6)

      ENDDO

      CLOSE(11)
      RETURN
      END