!-- File  verify.f90  ------------------------------------------------------
!
!  Performs several verification calculations given a file of grid spacings 
!  and some observed quantity corresponding to each grid spacing.
!
!  Computes:
!  - order of convergence
!  - Richardson extrapolation to zero grid spacing
!  - grid convergence indices (GCI)
!
!--------------------------------------------------------------------------

   PROGRAM verify

   IMPLICIT NONE

   INTEGER :: n
   INTEGER :: nd

   INTEGER, PARAMETER :: ndmax = 10

   REAL :: ratio
   REAL :: fexact
   REAL :: fsafe
   REAL :: p

   REAL, DIMENSION(ndmax) :: f
   REAL, DIMENSION(ndmax) :: gcif
   REAL, DIMENSION(ndmax) :: r
   REAL, DIMENSION(ndmax) :: x

!--------------------------------------------------------------------------

!..Write Header.

   WRITE(*,*) ' '
   WRITE(*,*) '--- VERIFY: Performs verification calculations ---'

!..Read in the file, determine the number of data points, and output.

   DO n = 1, ndmax
     READ(*,*,end=10) x(n), f(n)
   ENDDO

10 CONTINUE

   nd = n - 1

   WRITE(*,*) ' '
   WRITE(*,*) 'Number of data sets read = ', nd

   WRITE(*,*) ' '
   WRITE(*,*) '     Grid Size     Quantity'
   WRITE(*,*) ' '

   DO n = 1, nd
     WRITE(*,'(2(F14.6))') x(n), f(n)
   ENDDO

!..Compute the grid refinement ratio, r, between each pair.

   DO n = 1, nd-1
     r(n) = x(n+1) / x(n)
   ENDDO

!..Estimate the order of convergence using the first three data pairs 
!..and assuming that the grid refinement ratio is constant, r(1) = r(2).
!..This is done using Eqn. 5.10.6.1 of (Roache, 1998).

   p = log( ( f(3) - f(2) ) / ( f(2) - f(1) ) )  /  log( r(1) )

   WRITE(*,*) ' '
   WRITE(*,*) 'Order of convergence using first three finest grid '
   WRITE(*,*) 'and assuming constant grid refinement (Eqn. 5.10.6.1) '
   WRITE(*,*) 'Order of Convergence, p = ', p

!..Perform Richardson extrapolation to estimate a zero grid value of f.

   fexact = f(1) + ( f(1) - f(2) ) / ( r(1)**p - 1.0 )

   WRITE(*,*) ' '
   WRITE(*,*) 'Richardson Extrapolation: Use above order of convergence'
   WRITE(*,*) 'and first and second finest grids (Eqn. 5.4.1) '
   WRITE(*,*) 'Estimate to zero grid value, f_exact = ', fexact

!..Compute Grid Convergence Index (GCI) for each fine grid using Eqn. 5.6.1
!..from Roache's book. Use factor of safety as recommended on page 123.

   IF ( nd .ge. 3 ) then
     fsafe = 1.25
   ELSE
     fsafe = 3.0
   ENDIF

   DO n = 1, nd-1
     gcif(n) = fsafe * ( abs( f(n+1) - f(n) ) / f(n) ) / ( r(n)**p - 1.0 )
   ENDDO

   WRITE(*,*) ' '
   WRITE(*,*) 'Grid Convergence Index on fine grids. Uses p from above.'
   WRITE(*,*) 'Factor of Safety = ', fsafe
   WRITE(*,*) ' '
   WRITE(*,*) '  Grid     Refinement            '
   WRITE(*,*) '  Step      Ratio, r       GCI(%)'
   DO n = 1, nd-1
     WRITE(*,'(2x,i2,i3,2(f14.6))') n, n+1, r(n), gcif(n)*100.0
   ENDDO

!..Examine if asymptotic range has been achieved by examining ratio 
!..of Eqn. 5.10.5.2 of Roache's book is one.

   IF ( nd .ge. 3 ) then

     WRITE(*,*) ' '
     WRITE(*,*) 'Checking for asymptotic range using Eqn. 5.10.5.2.'
     WRITE(*,*) 'A ratio of 1.0 indicates asymptotic range.'
     WRITE(*,*) ' '
     WRITE(*,*) ' Grid Range    Ratio'
     DO n = 1, nd-2
       ratio = r(n)**p * gcif(n) / gcif(n+1)
       WRITE(*,'(2x,i1,i1,i2,i1,f14.6)') n, n+1, n+1, n+2, ratio
     ENDDO

   ENDIF


!..Write Trailer.

   WRITE(*,*) ' '
   WRITE(*,*) '--- End of VERIFY ---'
   WRITE(*,*) ' '
 
!--------------------------------------------------------------------------

    END PROGRAM verify