Program ConvVar

      IMPLICIT NONE

      INTEGER   i, j, k, l, Nx1, Np1, Np2, Np4, Nmax1, Nmax2, Nmax4, 
     /          Ntime, Nbound, Nmovie, Nuse1, Nuse2, Nuse4, Npick, 
     /          Ndir
      PARAMETER ( Np1 = 200, Np2 = 400, Np4 = 800, Nmax1 = 401, 
     /            Nmax2 = 801, Nmax4 = 1601 )
      REAL*8    DT1, DX1, DT2, DX2, DT4, DX4, thetax, thetat, thetaf,
     /          acon, Delta, Time, factor, Work(Nmax4),
     /          Q1(-1:1,-Np4:Np4), Mat1(Nmax4,Nmax4), Vec1(Nmax4),
     /          Q2(-1:1,-Np2:Np2), Mat2(Nmax2,Nmax2), Vec2(Nmax2),
     /          Q4(-1:1,-Np1:Np1), Mat4(Nmax1,Nmax1), Vec4(Nmax1)
      INTEGER   Ord1(Nmax4), Ord2(Nmax2), Ord4(Nmax1)

      OPEN(7, FILE = 'ConvVar.dat' , FORM = 'formatted'  )
         READ(7, 200, ERR = 1000 ) DT4, DX4, thetax, thetat, thetaf,
     /                             acon, Delta, Time
 200     FORMAT(7(11X,E11.4 / ),11X,E11.4)
         READ(7, 210, ERR = 1000 ) Nbound, Nmovie, Npick, Ndir
 210     FORMAT(3(11X,I5 / ),11X,I5)
      CLOSE(7)

      Nx1 = NINT(1.0D0/DX4)
      Nuse4 = 8*Nx1 + Nbound
      Nuse2 = 4*Nx1 + Nbound
      Nuse1 = 2*Nx1 + Nbound
      DT2 = DT4 / 2.0D0
      DX2 = DX4 / 2.0D0
      DT1 = DT4 / 4.0D0
      DX1 = DX4 / 4.0D0

      IF ( Nx1 .GT. Np1 ) THEN
         PRINT *, 'The coarsest grid is too fine.'
         STOP
      ENDIF
      Ntime = NINT(Time/DT4)

C Set up your initial conditions

      CALL INITIAL(DT1,DX1,4*Nx1,Q1,Delta,acon,Ndir)
      CALL INITIAL(DT2,DX2,2*Nx1,Q2,Delta,acon,Ndir)
      CALL INITIAL(DT4,DX4,Nx1,Q4,Delta,acon,Ndir)


C Now evolve the solutions to calculate the norm

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

      DO i = 1, Ntime
         DO k = 2*i-1, 2*i
            DO l = 2*k-1, 2*k
               CALL MATRIX(DT1,DX1,l-1,thetax,thetat,thetaf,acon,
     /                     Nbound,Nuse4,Mat1,Work,Ord1)
               CALL EVOL(DT1,DX1,l-1,thetax,thetat,thetaf,acon,
     /                   Nbound,4*Nx1,Q1,Nuse4,Mat1,Vec1,Ord1)
            ENDDO
            CALL MATRIX(DT2,DX2,k-1,thetax,thetat,thetaf,acon,Nbound,
     /                  Nuse2,Mat2,Work,Ord2)
            CALL EVOL(DT2,DX2,k-1,thetax,thetat,thetaf,acon,Nbound,
     /                2*Nx1,Q2,Nuse2,Mat2,Vec2,Ord2)
         ENDDO
         CALL MATRIX(DT4,DX4,i-1,thetax,thetat,thetaf,acon,Nbound,Nuse1,
     /               Mat4,Work,Ord4)
         CALL EVOL(DT4,DX4,i-1,thetax,thetat,thetaf,acon,Nbound,Nx1,Q4,
     /             Nuse1,Mat4,Vec4,Ord4)

         CALL FACTORING(i,DX1,4*Nx1,Q1,2*Nx1,Q2,Nx1,Q4,factor,Nmovie,
     /                  Npick)
         PRINT 500, DT4*i, factor
         WRITE(13,500) DT4*i, factor
 500     FORMAT(1X,F9.6,1P,2X,E13.6)
 510     FORMAT(1X,F9.6,1P,2(2X,E13.6))
      ENDDO

      CLOSE(13)
      STOP

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

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

      SUBROUTINE FACTORING(i,DX1,Nx4,Q1,Nx2,Q2,Nx1,Q4,factor,Nmovie,
     /                     Npick)

      IMPLICIT NONE

      INTEGER i, j, Nx4, Nx2, Nx1, Nmovie, Npick
      REAL*8  Q1(-1:1,-Nx4:Nx4), Q2(-1:1,-Nx2:Nx2), Q4(-1:1,-Nx1:Nx1), 
     /        numer, denom, numpt, denpt, factor, trunc, DX1
      CHARACTER*64 fname

      IF ( Nmovie .EQ. 0 ) GOTO 500

      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 )
      IF ( Npick .EQ. 1 ) THEN
         DO j = -Nx4, Nx4
            IF ( MOD(j,2) .EQ. 0 ) THEN
               trunc = ( Q2(1,j/2) - Q1(1,j) ) / 4.0D0
               IF ( ABS(trunc) .LT. 1.0D-89 ) trunc = 0.0D0
               WRITE(99,320) j*DX1, trunc
            ENDIF
         ENDDO
      ELSEIF ( Npick .EQ. 2 ) THEN
         DO j = -Nx2, Nx2
            IF ( MOD(j,2) .EQ. 0 ) THEN
               trunc = ( Q4(1,j/2) - Q2(1,j) ) / 4.0D0
               IF ( ABS(trunc) .LT. 1.0D-89 ) trunc = 0.0D0
               WRITE(99,320) j*DX1, trunc
            ENDIF
         ENDDO
      ELSE
         PRINT *, 'Invalid value for Npick.'
         STOP
      ENDIF

 320  FORMAT(1X,F9.6,2X,E17.10)
      CLOSE( UNIT = 99 )

 500  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.0D0*Nx1 + 1.0D0 ) )
      denom = SQRT( denom / ( 2.0D0*Nx2 + 1.0D0 ) )
      factor = numer / denom

      RETURN
      END