source: trunk/AGRIF/AGRIF_FILES/modinterpbasic.F @ 396

!
! $Id$
!
C     AGRIF (Adaptive Grid Refinement In Fortran)
C
C     Copyright (C) 2003 Laurent Debreu (Laurent.Debreu@imag.fr)
C                        Christophe Vouland (Christophe.Vouland@imag.fr)    
C
C     This program is free software; you can redistribute it and/or modify
C     it under the terms of the GNU General Public License as published by
C     the Free Software Foundation; either version 2 of the License, or
C     (at your option) any later version.
C
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License
C     along with this program; if not, write to the Free Software
C     Foundation, Inc., 59 Temple Place- Suite 330, Boston, MA 02111-1307, USA.
C
C
C
CCC   Module Agrif_Interpbasic
C
      Module Agrif_Interpbasic
C
CCC   Description:
CCC   Module containing different procedures of interpolation (linear,lagrange,
CCC   spline,...) used in the Agrif_Interpolation module.
C
C     Modules used:
      USE Agrif_types
C
      IMPLICIT NONE
C             
      CONTAINS
C     Define procedures contained in this module
C 
C
C     **************************************************************************  
CCC   Subroutine Linear1d  
C     ************************************************************************** 
C 
      Subroutine Linear1d(x,y,np,nc,
     &                    s_parent,s_child,ds_parent,ds_child) 
C
CCC   Description:
CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from
CCC   its parent grid (vector x).  
C
CC    Method:
C
C     Declarations:
C
      
C        
C     Arguments
      INTEGER             :: np,nc      
      REAL, DIMENSION(np) :: x      
      REAL, DIMENSION(nc) :: y  
      REAL                :: s_parent,s_child,ds_parent,ds_child
C
C     Local scalars
      INTEGER :: i,coeffraf,locind_parent_left
      REAL    :: ypos,globind_parent_left
C
C

      coeffraf = nint(ds_parent/ds_child)
C
      if (coeffraf == 1) then
C
          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
C        
          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
C
          return
C
      endif                          
C
      ypos = s_child      
C
      do i = 1,nc-1
C
        locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
C
        globind_parent_left = s_parent 
     &                        + (locind_parent_left - 1)*ds_parent
C        
        y(i) = ((globind_parent_left + ds_parent - ypos)
     &          *x(locind_parent_left)
     &        + (ypos - globind_parent_left)
     &          *x(locind_parent_left+1))
     &         / ds_parent                
C
        ypos = ypos + ds_child
C
      enddo
C
      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
C
      if (locind_parent_left == np) then
C
          y(nc) = x(np)
C
        else
C
          globind_parent_left = s_parent 
     &                        + (locind_parent_left - 1)*ds_parent  
C      
          y(nc) = ((globind_parent_left + ds_parent - ypos)
     &            *x(locind_parent_left)
     &          + (ypos - globind_parent_left)
     &            *x(locind_parent_left+1))
     &           / ds_parent     
C
      endif                                          
C           
      Return
C
C       
      End Subroutine Linear1d   
C
C
C
C     **************************************************************************  
CCC   Subroutine Lagrange1d  
C     **************************************************************************
C
      Subroutine Lagrange1d(x,y,np,nc,
     &                      s_parent,s_child,ds_parent,ds_child)
C
CCC   Description:
CCC   Subroutine to do a lagrange 1D interpolation on a child grid (vector y) 
CCC   from its parent grid (vector x).  
C
CC    Method:
C
C     Declarations:
C
      
C             
C     Arguments
      INTEGER             :: np,nc      
      REAL, DIMENSION(np) :: x      
      REAL, DIMENSION(nc) :: y  
      REAL                :: s_parent,s_child,ds_parent,ds_child 
C
C     Local scalars
      INTEGER :: i,coeffraf,locind_parent_left
      REAL    :: ypos,globind_parent_left
      REAL    :: X1,X2,X3 
C
C 
      if (np <= 2) then
C      
          Call Linear1D(x,y,np,nc,
     &                  s_parent,s_child,ds_parent,ds_child)
C         
         Return
C 
      endif
C
      coeffraf = nint(ds_parent/ds_child)
C
      if (coeffraf == 1) then
C
          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
C        
          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
C
          return
C
      endif
C
      ypos = s_child      
C
      do i = 1,nc
C
        locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)
C
        globind_parent_left = s_parent 
     &                        + (locind_parent_left - 1)*ds_parent 
C
        if (locind_parent_left+2 <= np) then            
C
            X1 = (x(locind_parent_left+1)-x(locind_parent_left))
     &           /ds_parent
C 
            X2 = (x(locind_parent_left+2)-x(locind_parent_left+1))
     &           /ds_parent
C
            X3 = (X2 - X1)/(2.*ds_parent)                 
C
            y(i) = x(locind_parent_left) + 
     &             (ypos - globind_parent_left)*X1 +
     &             (ypos - globind_parent_left)*
     &             (ypos - globind_parent_left - ds_parent)*X3 
C
          elseif (locind_parent_left+1 <= np) then 
C
            X1 = (x(locind_parent_left)-x(locind_parent_left-1))
     &           /ds_parent
C     
            X2 = (x(locind_parent_left+1)-x(locind_parent_left))
     &           /ds_parent
C 
            X3 = (X2 - X1)/(2.*ds_parent)                 
C
            y(i) = x(locind_parent_left-1) + 
     &             (ypos - globind_parent_left - ds_parent)*X1 +
     &             (ypos - globind_parent_left - ds_parent)*
     &             (ypos - globind_parent_left)*X3
C
          else
C
            X1 = (x(locind_parent_left-1)-x(locind_parent_left-2))
     &           /ds_parent
C     
            X2 = (x(locind_parent_left)-x(locind_parent_left-1))
     &           /ds_parent
C 
            X3 = (X2 - X1)/(2.*ds_parent)                 
C
            y(i) = x(locind_parent_left-2) + 
     &             (ypos - globind_parent_left - 2.*ds_parent)*X1 +
     &             (ypos - globind_parent_left - 2.*ds_parent)*
     &             (ypos - globind_parent_left - ds_parent)*X3
C
        endif   
C       
        ypos = ypos + ds_child 
C
      enddo
C
      return
C
C       
      End Subroutine Lagrange1d 
C
C
C     **************************************************************************  
CCC   Subroutine Constant1d  
C     ************************************************************************** 
C 
      Subroutine constant1d(x,y,np,nc,
     &                    s_parent,s_child,ds_parent,ds_child) 
C
CCC   Description:
CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from
CCC   its parent grid (vector x).  
C
CC    Method:
C
C     Declarations:
C
      
C        
C     Arguments
      INTEGER             :: np,nc      
      REAL, DIMENSION(np) :: x      
      REAL, DIMENSION(nc) :: y  
      REAL                :: s_parent,s_child,ds_parent,ds_child
C
C     Local scalars
      INTEGER :: i,coeffraf,locind_parent
      REAL    :: ypos
C
C

      coeffraf = nint(ds_parent/ds_child)
C
      if (coeffraf == 1) then
C
          locind_parent = 1 + nint((s_child - s_parent)/ds_parent)
C        
          y(1:nc) = x(locind_parent:locind_parent+nc-1)
C
          return
C
      endif                          
C
      ypos = s_child      
C
      do i = 1,nc
C
        locind_parent = 1 + nint((ypos - s_parent)/ds_parent)
C
        y(i) = x(locind_parent)
C
        ypos = ypos + ds_child
C
      enddo
C           
      Return
C
C       
      End Subroutine constant1d   
C
C     **************************************************************************  
CCC   Subroutine Linear1dconserv
C     ************************************************************************** 
C 
      Subroutine Linear1dconserv(x,y,np,nc,
     &                    s_parent,s_child,ds_parent,ds_child) 
C
CCC   Description:
CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from
CCC   its parent grid (vector x).  
C
CC    Method:
C
C     Declarations:
C
      Implicit none
C        
C     Arguments
      Integer             :: np,nc      
      Real, Dimension(np) :: x      
      Real, Dimension(nc) :: y
      Real, Dimension(:),Allocatable :: ytemp
      Real                :: s_parent,s_child,ds_parent,ds_child
C
C     Local scalars
      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
      Real    :: ypos
      integer :: i1,i2,ii
      real :: xpmin,xpmax,slope
      INTEGER :: diffmod
      REAL :: xdiffmod

C
C

      coeffraf = nint(ds_parent/ds_child)
C
      If (coeffraf == 1) Then
C
          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
C        
          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
C
          return
C
      End If
C              
      diffmod = 0
      IF (mod(coeffraf,2) == 0) diffmod = 1  

      xdiffmod = real(diffmod)/2.
                         
      allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
C
      ypos = s_child  
C      
      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 

      locind_parent_last = 1 +
     &  agrif_ceiling((ypos +(nc - 1) *ds_child - s_parent)/ds_parent)
     
      xpmin = s_parent + (locind_parent_left-1)*ds_parent 
      xpmax = s_parent + (locind_parent_last-1)*ds_parent          
      
      i1 = 1+agrif_int((xpmin-s_child)/ds_child)
      i2 = 1+agrif_int((xpmax-s_child)/ds_child)

      i = i1
      
      if (locind_parent_left == 1) then
        slope=
     &   (x(locind_parent_left+1)-x(locind_parent_left))/(coeffraf)
      else
         slope=
     &   (x(locind_parent_left+1)-x(locind_parent_left-1))/(2.*coeffraf)
      endif
      
        do ii=i-coeffraf/2+diffmod,i+coeffraf/2
          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
        enddo 

        locind_parent_left = locind_parent_left + 1
                      
      do i=i1 +  coeffraf, i2 - coeffraf,coeffraf
        slope=
     &   (x(locind_parent_left+1)-x(locind_parent_left-1))/(2.*coeffraf)
        do ii=i-coeffraf/2+diffmod,i+coeffraf/2
          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
        enddo
        locind_parent_left = locind_parent_left + 1
      enddo
      
      i = i2
      
      if (locind_parent_left == np) then
        slope=
     &   (x(locind_parent_left)-x(locind_parent_left-1))/(coeffraf)
      else
         slope=
     &   (x(locind_parent_left+1)-x(locind_parent_left-1))/(2.*coeffraf)
      endif
      
        do ii=i-coeffraf/2+diffmod,nc
          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
        enddo      
C
      y(1:nc)=ytemp(1:nc)                                   
C           
      deallocate(ytemp)
      Return
C       
      End Subroutine Linear1dconserv
      
C
C     **************************************************************************  
CCC   Subroutine Linear1dconservlim
C     ************************************************************************** 
C 
      Subroutine Linear1dconservlim(x,y,np,nc,
     &                    s_parent,s_child,ds_parent,ds_child) 
C
CCC   Description:
CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from
CCC   its parent grid (vector x).  
C
CC    Method:
C
C     Declarations:
C
      Implicit none
C        
C     Arguments
      Integer             :: np,nc      
      Real, Dimension(np) :: x      
      Real, Dimension(nc) :: y
      Real, Dimension(:),Allocatable :: ytemp
      Real                :: s_parent,s_child,ds_parent,ds_child
C
C     Local scalars
      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
      Real    :: ypos
      integer :: i1,i2,ii
      real :: xpmin,xpmax,slope
      INTEGER :: diffmod
      real :: xdiffmod
C
C

      coeffraf = nint(ds_parent/ds_child)
C
      If (coeffraf == 1) Then
C
          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
C        
          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
C
          return
C
      End If
C
      IF (coeffraf .NE.3) THEN
      print *,'LINEARCONSERVLIM not ready for refinement ratio = ',
     &   coeffraf
      stop
      ENDIF   
      
      diffmod = 0
      IF (mod(coeffraf,2) == 0) diffmod = 1        

      xdiffmod = real(diffmod)/2.
                         
      allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
C
      ypos = s_child  
C      
      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 

      locind_parent_last = 1 +
     &  agrif_ceiling((ypos +(nc - 1) *ds_child - s_parent)/ds_parent)
     
      xpmin = s_parent + (locind_parent_left-1)*ds_parent 
      xpmax = s_parent + (locind_parent_last-1)*ds_parent          
      
      i1 = 1+agrif_int((xpmin-s_child)/ds_child)
      i2 = 1+agrif_int((xpmax-s_child)/ds_child)

      i = i1
      
      if (locind_parent_left == 1) then
        slope=0.       
      else
        slope = vanleer(x(locind_parent_left-1:locind_parent_left+1))
        slope = slope / coeffraf    
      endif
      
        do ii=i-coeffraf/2+diffmod,i+coeffraf/2
          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
        enddo 

        locind_parent_left = locind_parent_left + 1
                      
      do i=i1 +  coeffraf, i2 - coeffraf,coeffraf     
        slope = vanleer(x(locind_parent_left-1:locind_parent_left+1))
        slope = slope / coeffraf

        do ii=i-coeffraf/2+diffmod,i+coeffraf/2
          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
        enddo
        locind_parent_left = locind_parent_left + 1
      enddo
      
      i = i2
      
      if (locind_parent_left == np) then
        slope=0.     
      else
        slope = vanleer(x(locind_parent_left-1:locind_parent_left+1))
        slope = slope / coeffraf      
      endif
      
        do ii=i-coeffraf/2+diffmod,nc
          ytemp(ii) = x(locind_parent_left)+(ii-i-xdiffmod/2.)*slope
        enddo      
C
      y(1:nc)=ytemp(1:nc)                                   
C           
      deallocate(ytemp)
      Return
C       
      End Subroutine Linear1dconservlim      
C         

C     **************************************************************************  
CCC   Subroutine ppm1d
C     ************************************************************************** 
C 
      Subroutine ppm1d(x,y,np,nc,
     &                    s_parent,s_child,ds_parent,ds_child) 
C
CCC   Description:
CCC   Subroutine to do a 1D interpolation and apply monotonicity constraints
CCC   using piecewise parabolic method  
CCC   on a child grid (vector y) from its parent grid (vector x).
CC    Method:
C
C     Declarations:
C
      Implicit none
C        
C     Arguments
      Integer             :: np,nc      
      Real, Dimension(np) :: x      
      Real, Dimension(nc) :: y
      Real, Dimension(:),Allocatable :: ytemp
      Real                :: s_parent,s_child,ds_parent,ds_child
C
C     Local scalars
      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
      Integer :: iparent,ipos,pos,nmin,nmax
      Real    :: ypos
      integer :: i1,jj
      Real :: xpmin,cavg,a,b
C      
      Real :: xrmin,xrmax,am3,s2,s1  
      Real, Dimension(np) :: dela,xr,xl,delta,a6,slope,slope2
      Real, Dimension(:),Allocatable  :: diff,diff2,diff3    
      INTEGER :: diffmod
C      
      coeffraf = nint(ds_parent/ds_child)
C
      If (coeffraf == 1) Then
          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
          return
      End If
C      
      Allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
      ypos = s_child  
C
      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
      locind_parent_last = 1 +
     &      agrif_ceiling((ypos +(nc - 1) 
     &      *ds_child - s_parent)/ds_parent)  
C
      xpmin = s_parent + (locind_parent_left-1)*ds_parent       
      i1 = 1+agrif_int((xpmin-s_child)/ds_child)        
C     
      Allocate( diff(coeffraf),diff2(coeffraf),diff3(coeffraf) )
C      
         diff(:) = ds_child/ds_parent
C      
      Do i=1,coeffraf
         a = real(i-1)*ds_child/ds_parent
         b = real(i)*ds_child/ds_parent
         diff2(i) = 0.5*(b*b - a*a)  
         diff3(i) = (1./3.)*(b*b*b - a*a*a)
      End do
C
      if( locind_parent_last+2 <= np ) then
           nmax = locind_parent_last+2    
      else if( locind_parent_last+1 <= np ) then
           nmax = locind_parent_last+1
      else
           nmax = locind_parent_last 
      endif     
C      
      if(locind_parent_left-1 >= 1) then
          nmin = locind_parent_left-1
      else 
          nmin = locind_parent_left
      endif    
C 
      Do i = nmin,nmax
         slope(i) = x(i) - x(i-1)
         slope2(i) = 2.*abs(slope(i))
      Enddo
C
      Do i = nmin,nmax-1
         dela(i) = 0.5 * ( slope(i) + slope(i+1) )
C Van Leer slope limiter
         dela(i) = min( abs(dela(i)),slope2(i),
     &                  slope2(i+1) )*sign(1.,dela(i))
         IF( slope(i)*slope(i+1) <= 0. ) dela(i) = 0.
      Enddo
C
      Do i = nmin,nmax-2
         xr(i) = x(i) + (1./2.)*slope(i+1) + (-1./6.)*dela(i+1)
     &                                     + ( 1./6. )*dela(i)
      Enddo
C
      Do i = nmin,nmax-2
         xrmin = min(x(i),x(i+1))
         xrmax = max(x(i),x(i+1))
         xr(i) = min(xr(i),xrmax)
         xr(i) = max(xr(i),xrmin)
         xl(i+1) = xr(i)         
      Enddo
C apply parabolic monotonicity
       Do i = locind_parent_left,locind_parent_last
          If( ( (xr(i)-x(i))* (x(i)-xl(i)) ) .le. 0. ) then
             xl(i) = x(i) 
             xr(i) = x(i)
          Endif          
          delta(i) = xr(i) - xl(i)
          am3 = 3. * x(i)
          s1  = am3 - 2. * xr(i)
          s2  = am3 - 2. * xl(i)
          IF( delta(i) * (xl(i) - s1) .le. 0. ) xl(i) = s1
          IF( delta(i) * (s2 - xr(i)) .le. 0. ) xr(i) = s2
          delta(i) = xr(i) - xl(i)
          a6(i) = 6.*x(i)-3.*(xl(i) +xr(i))
C
       End do   
C
        diffmod = 0
   IF (mod(coeffraf,2) == 0) diffmod = 1           
C
        ipos = i1
C               
        Do iparent = locind_parent_left,locind_parent_last       
             pos=1
             cavg = 0.
             Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
C
               ytemp(jj) = (diff(pos)*xl(iparent)   
     &             + diff2(pos)
     &             *  (delta(iparent)+a6(iparent))
     &             - diff3(pos)*a6(iparent))*coeffraf
                              
               cavg = cavg + ytemp(jj)
               pos = pos+1 
             End do 
             ipos = ipos + coeffraf
C
        End do     
C
C
        y(1:nc)=ytemp(1:nc)                                 
        deallocate(ytemp)                
        deallocate(diff, diff2, diff3)
      Return
      End Subroutine ppm1d
C                              
C     **************************************************************************  
CCC   Subroutine eno1d
C     ************************************************************************** 
C 
      Subroutine eno1d(x,y,np,nc,
     &                    s_parent,s_child,ds_parent,ds_child) 
C
CCC   Description:
CCC   ---- p 163-164 Computational gasdynamics ----
CCC   Subroutine to do a 1D interpolation 
CCC   using piecewise polynomial ENO reconstruction technique  
CCC   on a child grid (vector y) from its parent grid (vector x).
CC    Method:
C
C     Declarations:
C
      Implicit none
C        
C     Arguments
      Integer             :: np,nc      
      Real, Dimension(np) :: x      
      Real, Dimension(nc) :: y
      Real, Dimension(:),Allocatable :: ytemp
      Real                :: s_parent,s_child,ds_parent,ds_child
C
C     Local scalars
      Integer :: i,coeffraf,locind_parent_left,locind_parent_last
      Integer :: ipos,pos
      Real    :: ypos,xi
      integer :: i1,jj
      Real :: xpmin,cavg
C 
      Real, Dimension(3,np) :: dd,c
      Integer :: left
C           
      Real, DImension(1:np+1) :: xhalf
      Real, Dimension(:,:),Allocatable  :: Xbar 
      INTEGER :: diffmod     
C
      coeffraf = nint(ds_parent/ds_child)
C      
      If (coeffraf == 1) Then
          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent)
          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1)
          return
      End If     
      
      diffmod = 0
      IF (mod(coeffraf,2) == 0) diffmod = 1  
C      
      Allocate(ytemp(-2*coeffraf:nc+2*coeffraf))
      ypos = s_child  
      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
      locind_parent_last = 1 +
     &      agrif_ceiling((ypos +(nc - 1) *ds_child - 
     &      s_parent)/ds_parent)       
      xpmin = s_parent + (locind_parent_left-1)*ds_parent      
      i1 = 1+agrif_int((xpmin-s_child)/ds_child)           
C      
      xhalf(np+1) = np + 0.5
      Do i = 1,np
          xhalf(i) = i - 0.5
      Enddo
C
C compute divided differences
C
      dd(1,1:np) = x(1:np)
      dd(2,1:np-1) = 0.5*( dd(1,2:np) - dd(1,1:np-1) )
      dd(3,1:np-2) = (1./3.)*( dd(2,2:np-1) - dd(2,1:np-2) )
C
      Allocate( Xbar( coeffraf,2 ) )
      xi = 0.5
      Do i = 1,coeffraf
        Xbar(i,1) = (i-1)*ds_child/ds_parent - xi
        Xbar(i,2) = i*ds_child/ds_parent - xi
      Enddo
C
      ipos = i1
C           
      DO i = locind_parent_left,locind_parent_last           
         left = i           
         do jj = 2,3
             If(abs(dd(jj,left)) .gt. abs(dd(jj,left-1)))
     &               left = left-1            
         enddo
C           
C  convert to Taylor series form
C
         Call Taylor(i,xhalf(left:left+2),dd(1:3,left),c(1:3,i))
      ENDDO     
C
C evaluate the reconstruction on each cell
C
       DO i = locind_parent_left,locind_parent_last  
C
         cavg = 0.
         pos = 1.            
C         
         Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2
           ytemp(jj) =(c(1,i)*(Xbar(pos,2)-Xbar(pos,1))
     &                +c(2,i)*(Xbar(pos,2)*Xbar(pos,2)-
     &                         Xbar(pos,1)*Xbar(pos,1))
     &                +c(3,i)*(Xbar(pos,2)*Xbar(pos,2)*Xbar(pos,2)-
     &                         Xbar(pos,1)*Xbar(pos,1)*Xbar(pos,1)))
     &                         *coeffraf
           cavg = cavg + ytemp(jj)
           pos = pos+1
         Enddo                
         ipos = ipos + coeffraf                 
      ENDDO
C
      y(1:nc)=ytemp(1:nc)                                            
      deallocate(ytemp,Xbar)                 
C      
      Return       
      End Subroutine eno1d
C 
C      
C     **************************************************************************  
CCC   Subroutine taylor
C     ************************************************************************** 
C       
      subroutine taylor(ind,xhalf,dd,c)      
C      
      Integer :: ind
      real,dimension(3) :: dd,c     
      real,dimension(0:3,0:3) :: d 
      real,dimension(3) :: xhalf    
      integer ::i,j
C      
C
      d(0,0:3)=1.
      do i = 1,3
         d(i,0)=(ind-xhalf(i))*d(i-1,0)
      enddo  
C      
      do i = 1,3
         do j = 1,3-i
           d(i,j) = d(i,j-1) + (ind-xhalf(i+j))*d(i-1,j)
         enddo
      enddo
C         
      do j = 1,3        
         c(j) = 0.
         do i=0,3-j
            c(j) = c(j) + d(i,j)*dd(i+j)         
         enddo
      enddo
C      
      end subroutine taylor    
      
      
      REAL FUNCTION vanleer(tab)
      REAL, DIMENSION(3) :: tab
       real res1
       real p1,p2,p3
       
       p1=(tab(3)-tab(1))/2.
       p2=2.*(tab(2)-tab(1))
       p3=2.*(tab(3)-tab(2))
       
       if ((p1>0.).AND.(p2>0.).AND.(p3>0)) then
          res1=minval((/p1,p2,p3/))
       elseif ((p1<0.).AND.(p2<0.).AND.(p3<0)) then
          res1=maxval((/p1,p2,p3/))
       else
          res1=0.
       endif
          
          vanleer = res1   
      
      
      END FUNCTION vanleer 

C      
      End Module Agrif_Interpbasic