Program Converge09

      IMPLICIT NONE

      INTEGER   i, j, k, l, Nf1, Np1, Np2, Np4, Nmax1, Nmax2, Nmax4, 
     /          Nb1, Nb2, Nb4, Ntime, Nmovie, Npick, Ndir
      PARAMETER ( Np1 = 400, Np2 = 800, Np4 = 1600, Nmax1 = 401, 
     /            Nmax2 = 801, Nmax4 = 1601 )
      REAL*8    DT1, DX1, DT2, DX2, DT4, DX4, M, thetarr, thetatt, 
     /          thetar, thetat, Delta, Time, factor, Work(Nmax4),
     /          Q1(-1:1,0:Np4), Mat1(Nmax4,Nmax4), Vec1(Nmax4),
     /          Q2(-1:1,0:Np2), Mat2(Nmax2,Nmax2), Vec2(Nmax2),
     /          Q4(-1:1,0:Np1), Mat4(Nmax1,Nmax1), Vec4(Nmax1),
     /          test
      INTEGER   Ord1(Nmax4), Ord2(Nmax2), Ord4(Nmax1)

      OPEN(7, FILE = 'Converge09.dat' , FORM = 'formatted'  )
         READ(7, 200, ERR = 1000 ) DT4, DX4, M, thetarr, thetatt, 
     /                             thetar, thetat, Delta, Time
 200     FORMAT(8(11X,E11.4 / ),11X,E11.4)
         READ(7, 210, ERR = 1000 ) Nmovie, Npick, Ndir
 210     FORMAT(2(11X,I5 / ),11X,I5)
      CLOSE(7)

      Nf1 = NINT(4.0D0/DX4)
      DT2 = DT4 / 2.0D0
      DX2 = DX4 / 2.0D0
      DT1 = DT4 / 4.0D0
      DX1 = DX4 / 4.0D0

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

      IF ( M .GT. 0.5 .AND. M .LT. 2.5 ) THEN
         Nb1 = INT( ( 2.0D0*M - 1.0D0 )/DX4 ) + 1
         test = ABS( NINT( ( 2.0D0*M - 1.0D0 )/DX4 ) 
     /               - ( 2.0D0*M - 1.0D0 )/DX4 )
         IF ( test .LE. 1.0D-7 ) Nb1 = NINT( ( 2.0D0*M - 1.0D0 )/DX4 )

         Nb2 = INT( ( 2.0D0*M - 1.0D0 )/DX2 ) + 1
         test = ABS( NINT( ( 2.0D0*M - 1.0D0 )/DX2 ) 
     /               - ( 2.0D0*M - 1.0D0 )/DX2 )
         IF ( test .LE. 1.0D-7 ) Nb2 = NINT( ( 2.0D0*M - 1.0D0 )/DX2 )

         Nb4 = INT( ( 2.0D0*M - 1.0D0 )/DX1 ) + 1
         test = ABS( NINT( ( 2.0D0*M - 1.0D0 )/DX1 ) 
     /               - ( 2.0D0*M - 1.0D0 )/DX1 )
         IF ( test .LE. 1.0D-7 ) Nb4 = NINT( ( 2.0D0*M - 1.0D0 )/DX1 )
      ELSE
         Nb1 = 0
         Nb2 = 0
         Nb4 = 0
      ENDIF

C Set up your initial conditions.

      CALL INITIAL(DT1,DX1,M,Nb4,4*Nf1,Q1,Delta,Ndir)
      CALL INITIAL(DT2,DX2,M,Nb2,2*Nf1,Q2,Delta,Ndir)
      CALL INITIAL(DT4,DX4,M,Nb1,Nf1,Q4,Delta,Ndir)

C Calculate and invert your matrices.

      CALL MATRIX(DT1,DX1,M,thetarr,thetatt,thetar,thetat,Nb4,4*Nf1+1,
     /            Mat1,Work,Ord1)
      CALL MATRIX(DT2,DX2,M,thetarr,thetatt,thetar,thetat,Nb2,2*Nf1+1,
     /            Mat2,Work,Ord2)
      CALL MATRIX(DT4,DX4,M,thetarr,thetatt,thetar,thetat,Nb1,Nf1+1,
     /            Mat4,Work,Ord4)

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 EVOL(DT1,DX1,M,thetarr,thetatt,thetar,thetat,Nb4,
     /                   4*Nf1,Q1,4*Nf1+1,Mat1,Vec1,Ord1)
            ENDDO
            CALL EVOL(DT2,DX2,M,thetarr,thetatt,thetar,thetat,Nb2,2*Nf1,
     /                Q2,2*Nf1+1,Mat2,Vec2,Ord2)
         ENDDO
         CALL EVOL(DT4,DX4,M,thetarr,thetatt,thetar,thetat,Nb1,Nf1,Q4,
     /             Nf1+1,Mat4,Vec4,Ord4)

         CALL FACTORING(i,DX1,4*Nf1,Q1,2*Nf1,Q2,Nf1,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,Nf4,Q1,Nf2,Q2,Nf1,Q4,factor,Nmovie,
     /                     Npick)

      IMPLICIT NONE

      INTEGER i, j, Nf4, Nf2, Nf1, Nmovie, Npick
      REAL*8  Q1(-1:1,0:Nf4), Q2(-1:1,0:Nf2), Q4(-1:1,0:Nf1), 
     /        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 = 0, Nf4
            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 + 1.0D0, trunc
            ENDIF
         ENDDO
      ELSEIF ( Npick .EQ. 2 ) THEN
         DO j = 0, Nf2
            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 + 1.0D0, 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 = 0, Nf4
         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 / ( Nf1 + 1.0D0 ) )
      denom = SQRT( denom / ( Nf2 + 1.0D0 ) )
      factor = numer / denom

      RETURN
      END