c**************************************************************************** 
c  Reference for this technique can be found in
c 
c          AKIMA, H., 1970:  A NEW METHOD OF INTERPOLATION AND SMOOTH 
c             CURVE FITTING BASED ON LOCAL PROCEDURES.  J. ASSOC. 
c             COMPUT. MACH., 17, 582-602. 
c 
c  Subroutine akima computes a smooth curve through x(i), y(i),
c  i=1,nobs where nobs is greater than 2. The interpolated values yint
c  are computed at xint with an array size of nint. 
c
c  No extrapolation is done in this routine. Any xint value smaller than
c  x(1) or larger than x(nobs) will not have a valid yint value.
c
c  Written by Dr. Pat Fitzpatrick 8/1/97
c 
c  Affiliations:
c
c  Dr. Pat Fitzpatrick, Asst Professor of Meteorology
c  Dept. of Physics, Atmospheric Sciences, & General Science
c  Jackson State University
c
c****************************************************************************
      subroutine akima (nobs,x,y,nint,xint,yint)

      implicit none
      integer nobs, nint, i, j, L
      real a(4), b(4)
      real x(nobs),y(nobs),xint(nint),yint(nint)
      real w2, w3, m2, m3, a0, b0, t1, t2, p0, p1, p2, p3

c*********************************************************************
c  SET UP INITIAL INTERPOLATION COEFFICIENTS
c*********************************************************************
      a(1) = 3.*(x(2)-x(1)) - 2.*(x(3)-x(2))
      b(1) = 3.*(y(2)-y(1)) - 2.*(y(3)-y(2))
      a(2) = 2.*(x(2)-x(1)) -    (x(3)-x(2))
      b(2) = 2.*(y(2)-y(1)) -    (y(3)-y(2))
      a(3) =    (x(2)-x(1)) 
      b(3) =    (y(2)-y(1))

      if (xint(1).eq.x(1)) then
         yint(1)=y(1)
      endif

c*********************************************************************
c  LOOP THROUGH TOTAL POINTS NPTS IN THE INPUT DATA.
c*********************************************************************
      do i=1,nobs

         if (i.lt.(nobs-1)) then
            a(4) =    (x( i+2  )-x( i+1  ))
            b(4) =    (y( i+2  )-y( i+1  ))
         else if (i.eq.(nobs-1)) then
            a(4) = 2.*(x(nobs  )-x(nobs-1))
     &             -  (x(nobs-1)-x(nobs-2))
            b(4) = 2.*(y(nobs  )-y(nobs-1))
     &             -  (y(nobs-1)-y(nobs-2))
         else
            a(4) = 3.*(x(nobs  )-x(nobs-1))
     &            -2.*(x(nobs-1)-x(nobs-2))
            b(4) = 3.*(y(nobs  )-y(nobs-1))
     &            -2.*(y(nobs-1)-y(nobs-2))
         endif

         w2 = abs(a(3)*b(4)-a(4)*b(3)) 
         w3 = abs(a(1)*b(2)-a(2)*b(1)) 
      
         a0 = w2*a(2)+w3*a(3)
         b0 = w2*b(2)+w3*b(3)

c*********************************************************************
c If a0 is equal to 0, then set it equal to average of slope
c m2 and m3 (see confusing discussion on pg 591 in Akima (1970).
c
c My experience has shown this can occur for a line where y is
c constant for several adjacent values of x.
c
c NOTE!!!! This part has not been debugged very carefully. It may
c require further debugging. However, early tests show it seems to
c work.
c
c*********************************************************************
         if (a0.eq.0.0) then
            m3=(y(3)-y(2))/(x(3)-x(2))
            m2=(y(2)-y(1))/(x(2)-x(1))
            t2=0.5*(m3+m2)
         else
            t2=b0/a0
         endif

         if (i.ne.1) then
            p0=y(i-1)
            p1=t1
            p2=(3.*(y(i)-y(i-1))/(x(i)-x(i-1))-2.*t1-t2)/(x(i)-x(i-1))
            p3=(t1+t2-2.*(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i)-x(i-1))**2

c*********************************************************************
c  Compute yint if xint gt x(i-1) and le x(i)
c
c  See equation 2 in Akima (1970).
c*********************************************************************
            do L=1,nint
                if (xint(L).gt.x(i-1).and.xint(L).le.x(i)) then
                  yint(L)=p0+p1*(xint(L)-x(i-1))+p2*(xint(L)-x(i-1))**2
     &                                     +p3*(xint(L)-x(i-1))**3
               endif
            enddo
         endif

         t1=t2

c*********************************************************************
c  Set coefficients 2,3,4 to be equal to 1,2,3 for next loop
c*********************************************************************
         do j=1,3
            a(j) = a(j+1)
            b(j) = b(j+1)
         enddo

      enddo

      return
      end