         subroutine  eno_interpolation(qi, xi,     Nxi, 
     &                                 q,  x_imh,  Nx,
     &                                 size_sten,  geom,
     &                                 direction, V_imh )

c------------------------------------------------------------
c Fri Feb 22 09:22:57 PST 2002
c------------------------------------------------------------

         implicit      none

         character*17  cdnm
         parameter   ( cdnm = 'eno_interpolation' ) 

         logical       ltrace

         parameter   ( ltrace = .false. )
        

c------------------------------------------------------------
c  Interpolated values (output) at coordinates xi. (input)
c------------------------------------------------------------
         integer       Nxi

         real*8        qi ( Nxi )
         real*8        xi ( Nxi )

c------------------------------------------------------------
c  Cell averages and cell borders coordinates. (input)
c------------------------------------------------------------

         integer       Nx

         real*8        q     ( Nx )
         real*8        x_imh ( Nx )

c------------------------------------------------------------
c  Mass function. (internal)
c------------------------------------------------------------

         real*8        V_imh ( Nx ) 
         real*8        aux   ( 1000 ) 

c------------------------------------------------------------
c  Size of the eno stencil. (input)
c------------------------------------------------------------
         integer       size_sten

c------------------------------------------------------------
c   geom = 0 for cartesian       (input)
c        = 1 for spherical symmetry.
c------------------------------------------------------------

         integer       geom

c------------------------------------------------------------
c   direction =  1 for right reconstruction. (input)
c             = -1 for left  reconstruction.
c
c               (-1)       |      (1)
c                          |
c              =======>    |    <========
c                          |
c                          |
c  Only afects if xi is a border of a cell.
c------------------------------------------------------------
         integer       direction

         integer       pt,  pt2
         integer       i,   s   
c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,':          xi=',xi
            write(0,*) cdnm,':         Nxi=',Nxi
            write(0,*) cdnm,':          q=',q
            write(0,*) cdnm,':          x=',x_imh
            write(0,*) cdnm,':         Nx=',Nx
            write(0,*) cdnm,':  size_sten=',size_sten
            write(0,*) cdnm,':  direction=',direction
         endif
         if ( (size_sten+1)*(size_sten+2)/2 .gt. 1000 ) then
             write(0,*) cdnm,': Increase allocation for',
     &                  ' auxiliary array aux and recompile.'
             stop
         endif
         if     ( geom  .eq. 0 ) then
             call calc_V(V_imh, q, x_imh, Nx)
         elseif ( geom .eq. 1 ) then
             call calc_sph_V(V_imh, q, x_imh, Nx)
         endif
         pt2 = 1 
         do pt = 1, Nxi 
            do while( xi(pt) .ge. x_imh(pt2))
               pt2 = pt2 + 1 
            enddo
            if ( (direction .eq. -1) .and. 
     &           (xi(pt) .eq. x_imh(pt2-1)) 
     &         ) then
                     i = pt2 - 2 
            else
                     i = pt2 - 1  
            endif
            
            if ( ltrace ) then
               write(0,*) cdnm,': i=',i 
               write(0,*) cdnm,': xi(pt)=',xi(pt) 
            endif 
            pt2 = pt2 - 1
            if (pt2 .le. 0) pt2 = 1  
            call calc_stencil(V_imh, i, size_sten, 
     &                        x_imh, s, Nx, aux)        
            if ( geom .eq. 0 ) then
                call calc_int_value(qi(pt), xi(pt),
     &                          q(i-s), x_imh(i-s),
     &                          size_sten)
            elseif ( geom .eq. 1 ) then
                call calc_sph_int_value(qi(pt), xi(pt),
     &                                  q(i-s), x_imh(i-s),
     &                                  size_sten)
            endif
            if ( ltrace ) then
               write(0,*) cdnm,': qi(pt)=',qi(pt)
               write(0,*) cdnm,':-----------------'
            endif
         enddo
         if ( ltrace ) then
            write(0,*) cdnm,': qi=',qi
            write(0,*) cdnm,': xi=',xi
         endif
         return
         end
c===============================================================
c===============================================================
c===============================================================
         
         subroutine  calc_V(V_imh, q, x_imh, Nx)

c------------------------------------------------------------
c Fri Feb 15 13:20:15 PST 2002
c------------------------------------------------------------

         implicit      none

         character*6   cdnm
         parameter   ( cdnm = 'calc_V' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       Nx

         real*8        V_imh ( Nx )
         real*8        q     ( Nx )  
         real*8        x_imh ( Nx )

         integer       i
 
c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,': Nx = ', Nx
            write(0,*) cdnm,': q  = ', q
            write(0,*) cdnm,': x_imh  = ', x_imh
         endif
         V_imh (1) = 0.0d0 
         do i = 2, Nx 
            V_imh (i)= V_imh(i-1) + q(i-1) * (x_imh(i)-x_imh(i-1)) 
         enddo
         if ( ltrace ) then
             write(0,*) cdnm,': V_imh =',V_imh
         endif
         return
         end

c===============================================================
c===============================================================
c===============================================================
   
         subroutine  calc_sph_V(V_imh, q, x_imh, Nx)

c------------------------------------------------------------
c Mon Feb 25 09:41:47 PST 2002
c------------------------------------------------------------

         implicit      none

         character*10  cdnm
         parameter   ( cdnm = 'calc_sph_V' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       Nx

         real*8        V_imh ( Nx )
         real*8        q     ( Nx )  
         real*8        x_imh ( Nx )

         integer       i
 
c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,': Nx = ', Nx
            write(0,*) cdnm,': q  = ', q
            write(0,*) cdnm,': x_imh  = ', x_imh
         endif
         V_imh (1) = 0.0d0 
         do i = 2, Nx 
            V_imh (i)= V_imh(i-1) + q(i-1) * 
     &                 (x_imh(i)**3-x_imh(i-1)**3) 
         enddo
         if ( ltrace ) then
             write(0,*) cdnm,': V_imh =',V_imh
         endif
         return
         end

c===============================================================
c===============================================================
c===============================================================
 
         subroutine  calc_int_value(qi, xi, q, x_imh, k)

c------------------------------------------------------------
c Mon Feb 25 09:48:35 PST 2002
c------------------------------------------------------------

         implicit      none

         character*18   cdnm
         parameter   ( cdnm = 'calc_int_value' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         real*8        qi
         real*8        xi

         integer       k
        
         real*8        q     ( k )
         real*8        x_imh ( k )

         real*8        prod_deno2 
         real*8        sum_num2

         integer       j,     m

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,':  '
            write(0,*) cdnm,':  '
            write(0,*) cdnm,':  '
         endif
         qi = 0.0d0 
         do m = 0, k 
            do j = 0, m - 1 
               qi = qi + q(j+1)*(x_imh(j+2)-x_imh(j+1)) 
     &                * sum_num2(m,xi,x_imh,k) / prod_deno2(m,x_imh,k)
            enddo
         enddo
         if ( ltrace )  then
            write(0,*) cdnm,': qi=',qi
         endif
         return
         end
c============================================================
c============================================================
c============================================================
c============================================================
         real*8   function sum_num2(m,xi,x_imh,k)

         implicit      none

         character*8   cdnm
         parameter   ( cdnm = 'sum_num2' )

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       k

         real*8        x_imh( k )
         real*8        xi

         integer       m
         integer       l
         real*8        prod_num2

c============================================================
c============================================================
         if ( ltrace ) then
         endif
         sum_num2 = 0.0d0
         do l = 0, k
            if ( l .ne. m ) sum_num2 = sum_num2 + 
     &                      prod_num2(k,m,l,xi,x_imh)
         enddo
         if ( ltrace ) then
            write(0,*) cdnm,': sum_num2=',sum_num2
         endif
         return
         end      
c============================================================
c============================================================
c============================================================
c============================================================
         real*8   function prod_num2(k,m,l,xi, x_imh)

         implicit      none

         character*9   cdnm
         parameter   ( cdnm = 'prod_num2' )

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       k

         real*8        x_imh( k )
         real*8        xi
        
         integer       m,   l
         integer       q

c============================================================
c============================================================
         if ( ltrace ) then
         endif
         prod_num2 = 1.0d0
         do q = 0, k
            if ( q .ne. m .and. q .ne. l )
     &          prod_num2 = prod_num2*(xi-x_imh(q+1))
         enddo
         if ( ltrace ) then
            write(0,*) cdnm,': prod_num2=',prod_num2
         endif
         return
         end
c============================================================
c============================================================
c============================================================
c============================================================
         real*8   function prod_deno2(m,x_imh,k)

         implicit      none

         character*10  cdnm
         parameter   ( cdnm = 'prod_deno2' )

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       k
         real*8        x_imh( k )
         integer       m
         integer       l

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,': k=',k
            write(0,*) cdnm,': m=',m
         endif
         prod_deno2 = 1.0d0
         do l = 0, k
            if ( l .ne. m )
     &          prod_deno2 = prod_deno2*(x_imh(m+1)-x_imh(l+1))
         enddo
         if ( ltrace ) then
            write(0,*) cdnm,': prod_deno2=',prod_deno2
         endif
         return
         end
c============================================================
c============================================================
c============================================================
         subroutine  calc_sph_int_value(qi, xi, q, x_imh, k)

c------------------------------------------------------------
c Mon Feb 25 09:48:35 PST 2002
c------------------------------------------------------------

         implicit      none

         character*18   cdnm
         parameter   ( cdnm = 'calc_sph_int_value' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         real*8        qi
         real*8        xi

         integer       k
        
         real*8        q     ( k )
         real*8        x_imh ( k )

         real*8        prod_deno2 
         real*8        sum_num2 

         integer       j,     m

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,':  '
            write(0,*) cdnm,':  '
            write(0,*) cdnm,':  '
         endif
         qi = 0.0d0 
         do m = 0, k 
            do j = 0, m - 1 
               qi = qi + q(j+1)*(x_imh(j+2)**3-x_imh(j+1)**3)
     &               / (3.0d0*xi**2) 
     &             * sum_num2(m,xi,x_imh,k) / prod_deno2(m,x_imh,k)
            enddo
         enddo
         if ( ltrace )  then
            write(0,*) cdnm,': qi=',qi
         endif
         return
         end
c============================================================
c============================================================
c============================================================

         subroutine  calc_stencil(V_imh, i, k, x_imh, 
     &                            s, Nx, aux)

c------------------------------------------------------------
c    Fri Feb 15 13:40:47 PST 2002
c------------------------------------------------------------

         implicit      none

         character*12  cdnm
         parameter   ( cdnm = 'calc_stencil' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       Nx

c------------------------------------------------------------
c   Primitive function ( Mass function )
c------------------------------------------------------------
         real*8        V_imh( Nx )

c------------------------------------------------------------
c   Coordinates of the cell boundaries.
c------------------------------------------------------------
         real*8        x_imh( Nx )

         logical       uniform
         parameter   ( uniform = .false. )

c------------------------------------------------------------
c
c   i --> We are calculating the stencil for
c         cell i
c
c   k --> Size of the stencil that we want to create.
c 
c   s --> Positions to the left.  (output)
c
c------------------------------------------------------------

         integer       i
         integer       k

         integer       s

         integer       pt
c------------------------------------------------------------
c   Functions.
c------------------------------------------------------------

         real*8        newton_div_dif

         real*8        NDD,   NDD_left,  NDD_right

c------------------------------------------------------------
c Fliping variable 
c to keep the stencil as center as possible.
c------------------------------------------------------------

         logical       flip

         real*8        aux ( (k+1)*(k+2)/2 )

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,': Nx=',Nx
            write(0,*) cdnm,':  i=',i
            write(0,*) cdnm,':  k=',k
            write(0,*) cdnm,':  V_imh=',V_imh
            write(0,*) cdnm,':  x_imh=',x_imh
         endif
         flip = .false.
         pt = 1
         NDD= newton_div_dif(V_imh(i), x_imh(i), pt, 
     &                       uniform, aux)  
         s = 1
         do while (pt .lt. k)
            pt  = pt + 1
            NDD_left = newton_div_dif(V_imh(i-s),  x_imh(i-s), 
     &                               pt, uniform, aux) 
            NDD_right= newton_div_dif(V_imh(i-s+1),x_imh(i-s+1), 
     &                               pt, uniform, aux) 
            if ( abs(NDD_right) .eq. abs(NDD_left))  then
                 flip = .not. flip 
c                write(0,*) flip
                 if ( .not. flip ) s = s + 1
            elseif (abs(NDD_right) .gt. abs(NDD_left))  then
                 s = s + 1 
            endif
         enddo
         s = s - 1
           
         if ( ltrace ) then
                 write(0,*) cdnm,': s=',s 
             do  pt = i-s, i - s + k - 1 
                 write(0,*) cdnm,': i=',pt 
             enddo 
         endif 
         return
         end

c============================================================
c============================================================
c============================================================

         real*8 function newton_div_dif(V, x, j, uniform, aux)

c------------------------------------------------------------
c Fri Feb 15 09:58:56 PST 2002
c------------------------------------------------------------

         implicit      none

         character*14   cdnm
         parameter    ( cdnm = 'newton_div_dif' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         logical       uniform

         integer       j

         real*8        V ( j+1 )
         real*8        x ( j+1 )

         real*8        aux ( (j+1)*(j+2)/2 ) 

         integer       i

         real*8        dx

         real*8        bino

         real*8        fac

         integer       s,     k

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,':  uniform =', uniform
            write(0,*) cdnm,':        j =', j
            write(0,*) cdnm,':        V =', V
            write(0,*) cdnm,':        x =', x
            write(0,*) cdnm,':        int =', (j+1)*(j+2)/2
         endif
         newton_div_dif = 0.0d0
         if ( uniform ) then
            dx     = x(2) - x(1) 
            do i = 0, j  
               newton_div_dif = newton_div_dif + 
     &                          (-1.0d0)**i * bino(j,i) * V(j+1-i) 
            enddo
            newton_div_dif = newton_div_dif / (fac(j)*dx**(j)) 
         else
            do k = 0, j 
               do s = 0, j - k
                  if ( k .eq. 0 ) then 
                     if ( ltrace ) then  
                        write(0,*) cdnm,': int=',
     &                             k*(2*(j+1)-k+1)/2+(s+1)
                        write(0,*) cdnm,': V=',V(s+1)
                     endif
                     aux(k*(2*(j+1)-k+1)/2+(s+1)) = V(s+1)  
                  else
                     aux(k*(2*(j+1)-k+1)/2+(s+1)) = 
     &                 ( aux((k-1)*(2*(j+1)-k+2)/2+(s+2)) -
     &                   aux((k-1)*(2*(j+1)-k+2)/2+(s+1))
     &                 ) /
     &                 (x(s+k+1)-x(s+1))  
                  endif
                  if ( ltrace ) then
                     write(0,*) cdnm,':  aux=',aux
                  endif
               enddo
            enddo
            newton_div_dif =  aux( (j+1)*(j+2)/2 ) 
         endif
         if ( ltrace ) then
            write(0,*) cdnm,':        aux=',aux
            write(0,*) cdnm,':        newton_div_dif =', newton_div_dif
         endif
         return
         end

c============================================================
c============================================================
c============================================================
   
         subroutine  eno_recos_qL(q, k, s, rqL_iph)

c------------------------------------------------------------
c  Feb 15 18:01:23 PST 2002
c------------------------------------------------------------

         implicit      none

         character*12  cdnm
         parameter   ( cdnm = 'eno_recos_qL' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       k
         real*8        q ( k )
         real*8        c ( k )

         integer       s

         integer       j

         real*8        rqL_iph

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,':  q=',q
            write(0,*) cdnm,':  k=',k
            write(0,*) cdnm,':  s=',s
         endif
         call calc_coef(s, k, c)
         rqL_iph = 0.0d0 
         do j = 1, k 
            rqL_iph = rqL_iph  + c(j) * q(j)
         enddo
         if ( ltrace) then
            write(0,*) cdnm,': rqL_iph =', rqL_iph
         endif
         return
         end

c============================================================
c============================================================
c============================================================

         subroutine  eno_recos_qR(q, k, s, rqR_iph)

c------------------------------------------------------------
c  Feb 15 18:01:23 PST 2002
c------------------------------------------------------------

         implicit      none

         character*12  cdnm
         parameter   ( cdnm = 'eno_recos_qR' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       k
         real*8        q ( k )
         real*8        c ( k )

         integer       s

         integer       j

         real*8        rqR_iph

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,':  q=',q
            write(0,*) cdnm,':  k=',k
            write(0,*) cdnm,':  s=',s
         endif
         call calc_coef(s, k, c)
         rqR_iph = 0.0d0 
         do j = 1, k 
            rqR_iph = rqR_iph  + c(j) * q(j)
         enddo
         if ( ltrace ) then
            write(0,*) cdnm,': rqR_iph =', rqR_iph
         endif
         return
         end

c============================================================
c============================================================
c============================================================

         real*8 function fac(n)

c------------------------------------------------------------
c    Fri Feb 15 10:57:36 PST 2002
c------------------------------------------------------------

         implicit      none

         character*3   cdnm
         parameter   ( cdnm = 'fac' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       n  
         integer       i  

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,': n = ', n
         endif
         i   = 1.0d0
         fac = 1.0d0
         if ( n .eq. 0 ) then
             fac = 1.0d0  
         else if ( n .gt. 0 ) then
            do while (i .le. n) 
               fac = fac * i 
               i = i + 1
            enddo
         else
            write(0,*) cdnm,': Cant compute factorial'
            write(0,*) cdnm,': n=',n
            fac = -1.0d0
         endif
         return
         end

c============================================================
c============================================================
c============================================================

         real*8 function bino(n, k)

c------------------------------------------------------------
c  Fri Feb 15 11:05:42 PST 2002
c------------------------------------------------------------

         implicit      none

         character*4   cdnm
         parameter   ( cdnm = 'bino' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       n,   k,   nmk

         real*8        fac

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,': n=  ', n
            write(0,*) cdnm,': k=  ', k
         endif
         nmk =  n - k 
         if ( n .lt. 0.0d0 ) then
            write(0,*) cdnm,': Cant compute binomial coefficient'
            write(0,*) cdnm,': n = ', n 
            bino = -1.0d0
         elseif ( k .lt. 0.0d0 ) then
            write(0,*) cdnm,': Cant compute binomial coefficient'
            write(0,*) cdnm,': k = ', k 
            bino = -1.0d0
         elseif ( nmk .lt. 0.0d0 )  then
            write(0,*) cdnm,': Cant compute binomial coefficient'
            write(0,*) cdnm,': n-k = ', nmk 
            bino = -1.0d0
         else
            bino = fac(n) / (fac(k) * fac(nmk))
         endif 
         return
         end

c============================================================
c============================================================
c============================================================
 
         subroutine  calc_coef(s, k, c )

c------------------------------------------------------------
c  Fri Feb 15 16:04:14 PST 2002
c------------------------------------------------------------
c------------------------------------------------------------
c  Computation of the coeficients for uniform grid.
c  for Lagrange polynomials for the boundaries.
c  See formula (2.21) in
c
c          NASA/CR-97-206253
c          ICASE Report No. 97-65
c          
c------------------------------------------------------------

         implicit      none

         character*8   cdnm
         parameter   ( cdnm = 'calc_coef' ) 

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       k
        
         real*8        c( k )

         integer       s          

         integer       m,   j

         real*8        sum_num
         real*8        prod_deno

c============================================================
c============================================================
         if ( ltrace ) then
            write(0,*) cdnm,': s = ', s 
            write(0,*) cdnm,': k = ', k
         endif
         do j = 0, k-1 
            c(j+1) = 0.0d0  
            do m = j+1 , k
                 c(j+1) = c(j+1) +  sum_num(s,k,m)  / 
     &                              prod_deno(m,k)
            enddo
         enddo
         if ( ltrace ) then
            write(0,*) cdnm,': c = ', c 
         endif
         return
         end
c============================================================
c============================================================
c============================================================
c============================================================
         real*8   function sum_num(r,k,m)
       
         implicit      none  

         character*7   cdnm
         parameter   ( cdnm = 'sum_num' )

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       r,   k,   m
         integer       l 
 
         real*8        prod_num

c============================================================
c============================================================
         if ( ltrace ) then        
            write(0,*) cdnm,': r=',r 
            write(0,*) cdnm,': k=',k 
            write(0,*) cdnm,': m=',m 
         endif
         sum_num = 0.0d0 
         do l = 0, k
            if ( l .ne. m ) sum_num = sum_num + prod_num(k,m,l,r)  
         enddo
         if ( ltrace ) then
            write(0,*) cdnm,': sum_num=',sum_num 
         endif
         return
         end
c============================================================
c============================================================
c============================================================
c============================================================
         real*8   function prod_num(k,m,l,r)
       
         implicit      none  

         character*8   cdnm
         parameter   ( cdnm = 'prod_num' )

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       r,   k,   m,   l
         integer       q 

c============================================================
c============================================================
         if ( ltrace ) then        
            write(0,*) cdnm,': r=',r 
            write(0,*) cdnm,': k=',k 
            write(0,*) cdnm,': m=',m 
            write(0,*) cdnm,': l=',l 
         endif
         prod_num = 1.0d0 
         do q = 0, k
            if ( q .ne. m .and. q .ne. l ) 
     &          prod_num = prod_num*(r-q+1.0d0)
         enddo
         if ( ltrace ) then
            write(0,*) cdnm,': prod_num=',prod_num 
         endif
         return
         end
c============================================================
c============================================================
c============================================================
c============================================================
         real*8   function prod_deno(m,k)
       
         implicit      none  

         character*9   cdnm
         parameter   ( cdnm = 'prod_deno' )

         logical       ltrace
         parameter   ( ltrace = .false. )

         integer       k,   m
         integer       l 

c============================================================
c============================================================
         if ( ltrace ) then        
            write(0,*) cdnm,': k=',k 
            write(0,*) cdnm,': m=',m 
         endif
         prod_deno = 1.0d0 
         do l = 0, k
            if ( l .ne. m ) 
     &          prod_deno = prod_deno*(m-l)
         enddo
         if ( ltrace ) then
            write(0,*) cdnm,': prod_deno=',prod_deno 
         endif
         return
         end
