Program Converge02

      INTEGER   Ntime, Np1, Np2, Np4, 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 = 'Converge02.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 = NINT( 1.0D0 / DX )
      IF ( Nx1 .GT. Np1 ) THEN
         PRINT *, 'The coarsest grid is too fine.'
         STOP
      ENDIF
      IF ( DT .GE. DX ) THEN
         PRINT *, 'DT is too large.'
         STOP
      ENDIF

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

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

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

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

      INTEGER   i, j, k, l, Ntime, Nx1, Nx2, Nx4, Nnon0, Ngraph, Nbound, 
     /          Npa
      REAL*8    DT, DX, Delta, Q1(-1:1,-Nx4:Nx4), Q2(-1:1,-Nx2:Nx2),
     /          Q4(-1:1,-Nx1:Nx1), Time, norm1, norm2, norm4, factor42,
     /          factor21, dx2, dx4, dt2, dt4, pi
      PARAMETER ( Npa = 1600, pi = 3.141593D0 )
      REAL*8    QA(-Npa:Npa)

C Initialize all Data. 
      
      Nnon0 = INT(Nx4*0.2D0)
      Dx4 = DX / 4.0D0
      DO j = -Nx4, -Nnon0-1
         Q1(1,j) = 0.0D0
      ENDDO 
      DO j = -Nnon0, Nnon0
         IF ( Ngraph .EQ. 1 ) THEN
            Q1(1,j) = COS(pi/Delta*(Dx4*j))/2.0D0 + 0.5D0
         ELSE
            Q1(1,j) = 0.04D0*EXP( -(Dx4*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 using the analytic solution.

      Dt4 = DT / 4.0D0
      DO j = -Nx4, -Nnon0-1
         Q1(0,j) = 0.0D0
      ENDDO 
      DO j = -Nnon0, Nnon0-1
         IF ( Ngraph .EQ. 1 ) THEN
            Q1(0,j) = COS(pi/Delta*(Dx4*j+Dt4))/2.0D0 + 0.5D0
         ELSE
            Q1(0,j) = 0.04D0*EXP( -(Dx4*j+Dt4)**2 / Delta**2 )
         ENDIF
      ENDDO 
      DO j = Nnon0, Nx4
         Q1(0,j) = 0.0D0
      ENDDO 

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

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

C Now get the norm for each timestep.  Be careful to evolve all three grids to
C the SAME moment in time.

      OPEN( UNIT = 11, file = 'Factor02a.dat' )
      OPEN( UNIT = 12, file = 'Factor02b.dat' )

      Ntime = INT(Time/DT)
      DO i = 1, Ntime
         DO k = 2*i-1, 2*i
            DO l = 2*k-1, 2*k
               CALL GETNORM(l,DT4,DX4,Delta,Ngraph,Nbound,Q1,QA,Nx4,
     /                      norm1)
            ENDDO
            CALL GETNORM(k,DT2,DX2,Delta,Ngraph,Nbound,Q2,QA,Nx2,norm2)
            factor21 = norm2 / norm1
            WRITE(11,500) DT*i, factor21
         ENDDO
         CALL GETNORM(i,DT,DX,Delta,Ngraph,Nbound,Q4,QA,Nx1,norm4)
         factor42 = norm4 / norm2
         WRITE(12,500) DT*i, factor42
 500     FORMAT(1X,F9.6,1P,2X,E13.6)
      ENDDO

      CLOSE(11)
      RETURN
      END

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

      SUBROUTINE GETNORM(i,DT,DX,Delta,Ngraph,Nbound,Q,QA,Nx,norm)

      INTEGER   i, j, Nx, Nnon0f, Nnon0b, Nnon0l, Nnon0t, Ngraph, Nbound
      REAL*8    DT, DX, Delta, Q(-1:1,-Nx:Nx), QA(-Nx:Nx), norm, front,
     /          back, lead, tail, test

      DO j = -Nx, Nx
         Q(-1,j) = Q(0,j)
         Q(0,j) = Q(1,j)
      ENDDO

C Evolve to the next timestep.

      DO j = -Nx+1, Nx-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,Nx) = 2.0D0*Q(0,Nx) - Q(-1,Nx) + ( Q(0,-Nx+1) 
     /             - 2.0D0*Q(0,Nx) + Q(0,Nx-1) )*( DT/DX )**2
         Q(1,-Nx) = Q(1,Nx)
      ELSEIF ( Nbound .EQ. 0 ) THEN
         Q(1,Nx) = 0.0D0
         Q(1,-Nx) = 0.0D0
      ELSEIF ( Nbound .EQ. 1 ) THEN
         Q(1,Nx) = ( 4.0D0*Q(0,Nx) - Q(-1,Nx) + ( DT/DX )
     /             *( 4.0D0*Q(1,Nx-1) - Q(1,Nx-2) ) ) /
     /             ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
         Q(1,-Nx) = ( 4.0D0*Q(0,-Nx) - Q(-1,-Nx) + ( DT/DX )
     /              *( 4.0D0*Q(1,-Nx+1) - Q(1,-Nx+2) ) ) /
     /              ( 3.0D0*( 1.0D0 +( DT/DX ) ) )
      ELSE
         PRINT *, 'Invalid value of Nbound' 
         STOP
      ENDIF

C Calculate the Analytic Solution.:

C The Periodic Boundary Condition

      DO j = -Nx, Nx
         QA(j) = 0.0D0
      ENDDO

      IF ( Nbound .EQ. -1 ) THEN
         front = DMOD((1.2D0 + i*DT),2.0D0) - 1.0D0
         IF ( front .GE. 0.0D0 ) THEN
            Nnon0f = INT(front/DX)
         ELSE
            Nnon0f = INT(front/DX) - 1
         ENDIF
         back = DMOD((0.8D0 + i*DT),2.0D0) - 1.0D0
         IF ( back .LE. 0.0D0 ) THEN
            Nnon0b = INT(back/DX)
         ELSE
            Nnon0b = INT(back/DX) + 1
         ENDIF
         test = ABS( NINT(front/DX) - front/DX )
         IF ( test .LE. 1.0D-7 ) Nnon0f = NINT(front/DX)
         test = ABS( NINT(back/DX) - back/DX )
         IF ( test .LE. 1.0D-7 ) Nnon0b = NINT(back/DX)
         IF ( front .LT. back ) THEN
            DO j = -Nx, Nnon0f
               IF ( Ngraph .EQ. 1 ) THEN
                  QA(j) = COS(pi/Delta*( DX*j - front + 0.2D0 ))/2.0D0
     /                    + 0.5D0
               ELSE
                  QA(j) = 0.04D0*EXP( -( DX*j - front + 0.2D0 )**2
     /                     / Delta**2 )
               ENDIF
            ENDDO
            DO j = Nnon0b, Nx
               IF ( Ngraph .EQ. 1 ) THEN
                  QA(j) = COS(pi/Delta*( DX*j - back - 0.2D0 ))/2.0D0
     /                    + 0.5D0
               ELSE
                  QA(j) = 0.04D0*EXP( -( DX*j - back - 0.2D0 )**2 
     /                      / Delta**2 )
               ENDIF
            ENDDO
         ELSE
            DO j = Nnon0b, Nnon0f
               IF ( Ngraph .EQ. 1 ) THEN
                  QA(j) = COS(pi/Delta*( DX*j - front + 0.2D0 ))/2.0D0
     /                    + 0.5D0
               ELSE
                  QA(j) = 0.04D0*EXP( -( DX*j - front + 0.2D0 )**2
     /                    / Delta**2 )
               ENDIF
            ENDDO
         ENDIF
      ENDIF

C The Fixed Boundary Condition

      IF ( Nbound .EQ. 0 ) THEN
         front = DMOD((2.2D0 + i*DT),4.0D0) - 2.0D0
         IF ( front .GE. 0.0D0 ) THEN
            Nnon0f = INT(front/DX)
         ELSE
            Nnon0f = INT(front/DX) - 1
         ENDIF
         back = DMOD((1.8D0 + i*DT),4.0D0) - 2.0D0
         IF ( back .LE. 0.0D0 ) THEN
            Nnon0b = INT(back/DX)
         ELSE
            Nnon0b = INT(back/DX) + 1
         ENDIF
         test = ABS( NINT(front/DX) - front/DX )
         IF ( test .LE. 1.0D-7 ) Nnon0f = NINT(front/DX)
         test = ABS( NINT(back/DX) - back/DX )
         IF ( test .LE. 1.0D-7 ) Nnon0b = NINT(back/DX)
         IF ( front .LT. -1.0D0 .OR. back .GT. 1.0D0 ) GOTO 400
         IF ( front .LT. -0.6D0 ) THEN
            DO j = -Nx, Nnon0f
               IF ( Ngraph .EQ. 1 ) THEN
                  QA(j) = COS(pi/Delta*( DX*j - front + 0.2D0 ))/2.0D0
     /                    + 0.5D0
               ELSE
                  QA(j) = 0.04D0*EXP( -( DX*j - front + 0.2D0 )**2
     /                      / Delta**2 )
               ENDIF
            ENDDO
         ELSEIF ( back .GT. 0.6D0 ) THEN
            DO j = Nnon0b, Nx
               IF ( Ngraph .EQ. 1 ) THEN
                  QA(j) = COS(pi/Delta*( DX*j - back - 0.2D0 ))
     /                      /2.0D0 + 0.5D0
               ELSE
                  QA(j) = 0.04D0*EXP( -( DX*j - back - 0.2D0 )**2 
     /                      / Delta**2 )
               ENDIF
            ENDDO
         ELSE
            DO j = Nnon0b, Nnon0f
               IF ( Ngraph .EQ. 1 ) THEN
                  QA(j) = COS(pi/Delta*( DX*j - front + 0.2D0 ))/2.0D0
     /                    + 0.5D0
               ELSE
                  QA(j) = 0.04D0*EXP( -( DX*j - front + 0.2D0 )**2
     /                      / Delta**2 )
               ENDIF
            ENDDO
         ENDIF

 400     IF ( (3.8D0 - i*DT) .LT. 0.0D0 ) THEN
            lead = 2.0D0 + DMOD((3.8D0 - i*DT),4.0D0) 
         ELSE
            lead = DMOD((3.8D0 - i*DT),4.0D0) - 2.0D0
         ENDIF
         IF ( lead .LE. 0.0D0 ) THEN
            Nnon0l = INT(lead/DX)
         ELSE
            Nnon0l = INT(lead/DX) + 1
         ENDIF
         IF ( (0.2D0 - i*DT) .LT. 0.0D0 ) THEN
            tail = 2.0D0 + DMOD((0.2D0 - i*DT),4.0D0) 
         ELSE
            tail = DMOD((0.2D0 - i*DT),4.0D0) - 2.0D0
         ENDIF
         IF ( tail .GE. 0.0D0 ) THEN
            Nnon0t = INT(tail/DX)
         ELSE
            Nnon0t = INT(tail/DX) - 1
         ENDIF
         test = ABS( NINT(lead/DX) - lead/DX )
         IF ( test .LE. 1.0D-7 ) Nnon0l = NINT(lead/DX)
         test = ABS( NINT(tail/DX) - tail/DX )
         IF ( test .LE. 1.0D-7 ) Nnon0t = NINT(tail/DX)
         IF ( lead .GT. 1.0D0 .OR. tail .LT. -1.0D0 ) GOTO 500
         IF ( lead .GT. 0.6D0 ) THEN
            DO j = Nnon0l, Nx
               IF ( Ngraph .EQ. 1 ) THEN
                  QA(j) = QA(j)
     /                    - COS(pi/Delta*( DX*j - lead - 0.2D0 ))/2.0D0
     /                    - 0.5D0
               ELSE
                  QA(j) = QA(j)
     /                    - 0.04D0*EXP( -( DX*j - lead - 0.2D0 )**2 
     /                      / Delta**2 )
               ENDIF
            ENDDO
         ELSEIF ( tail .LT. -0.6D0 ) THEN
            DO j = -Nx, Nnon0t
               IF ( Ngraph .EQ. 1 ) THEN
                  QA(j) = QA(j)
     /                    - COS(pi/Delta*( DX*j - tail + 0.2D0 ))/2.0D0
     /                    - 0.5D0
               ELSE
                  QA(j) = QA(j) 
     /                    - 0.04D0*EXP( -( DX*j - tail + 0.2D0 )**2
     /                      / Delta**2 )
               ENDIF
            ENDDO
         ELSE
            DO j = Nnon0l, Nnon0t
               IF ( Ngraph .EQ. 1 ) THEN
                  QA(j) = QA(j)
     /                    - COS(pi/Delta*( DX*j - lead - 0.2D0 ))/2.0D0
     /                    - 0.5D0
               ELSE
                  QA(j) = QA(j)
     /                    - 0.04D0*EXP( -( DX*j - lead - 0.2D0 )**2
     /                      / Delta**2 )
               ENDIF
            ENDDO
         ENDIF
      ENDIF

C The Radiation Boundary Condition

      IF ( Nbound .EQ. 1 ) THEN
         front = 0.2D0 + i*DT
         back = -0.2D0 + i*DT
         IF ( back .GT. 1.0D0 ) GOTO 500
         IF ( back .LE. 0.0D0 ) THEN
            Nnon0b =  INT(back/DX)
         ELSE
            Nnon0b =  INT(back/DX) + 1
         ENDIF
         test = ABS( NINT(back/DX) - back/DX )
         IF ( test .LE. 1.0D-7 ) Nnon0b = NINT(back/DX)
         IF ( front .LE. 1.0D0 ) THEN
            Nnon0f = INT(front/DX)
            test = ABS( NINT(front/DX) - front/DX )
            IF ( test .LE. 1.0D-7 ) Nnon0f = NINT(front/DX)
         ELSE
            Nnon0f = Nx
         ENDIF
         DO j = Nnon0b, Nnon0f
            IF ( Ngraph .EQ. 1 ) THEN
               QA(j) = COS(pi/Delta*(DX*j-DT*i))/2.0D0 + 0.5D0
            ELSE
               QA(j) = 0.04D0*EXP( -(DX*j-DT*i)**2 / Delta**2 )
            ENDIF
         ENDDO
      ENDIF

C Now calculate the norm.

 500  norm = 0.0D0
      DO j = -Nx, Nx
         norm = norm + ( Q(1,j) - QA(j) )**2
      ENDDO
      norm = SQRT( norm / ( 2.0D0*Nx + 1.0D0 ) )

      RETURN
      END