! $Id: splinx1.f90 2022-11-02 09:51:28Z ychen $

!****s* fsi/splinx1 *
!
! NAME
!   splinx1  - natural cubic spline interpolation
!
! SYNOPSIS
!    call splinx1(x, y, xx, yy, der)
!
! DESCRIPTION
!    This subroutine interpolates a grid function onto another
!    grid by means of natural cubic spline. Input function y(x);
!    output grid xx; output function yy; its derivative der.
!
! INPUTS
!    REAL(wp), DIMENSION(:)   :: x        array of input argument data
!    REAL(wp), DIMENSION(:)   :: y        array of input function data
!    REAL(wp),                :: xx       output argument data
!
! OUTPUT
!    REAL(wp)                 :: yy       output function data
!    REAL(wp)                 :: der      derivative of output function
!
! AUTHOR
!   S.V.Sokolovskiy, UCAR
!   Update: Yong Chen, NOAA/NESDIS/STAR, yong.chen@noaa.gov
!
! COPYRIGHT
!   Copyright (c) 2022-2023 Yong Chen <yong.chen@noaa.gov>
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

SUBROUTINE splinx1(x, y, xx, yy, der)

      USE typesizes, ONLY: wp => EightByteReal

      IMPLICIT NONE
      REAL(wp), DIMENSION(:), INTENT(in)     :: x       ! array of input argument data
      REAL(wp), DIMENSION(:), INTENT(in)     :: y       ! array of input function data
      REAL(wp),               INTENT(in)     :: xx      ! output argument data
      REAL(wp),               INTENT(out)    :: yy      ! output function data
      REAL(wp),               INTENT(out)    :: der     ! derivative of output function

      REAL(wp), DIMENSION(:,:), ALLOCATABLE  :: c       ! Matrix for regression
      INTEGER                                :: n       ! array size of x

      REAL(wp)  :: dxr, dyr, dxl, dyl, dx, dy, dyx, d
      INTEGER   :: i, m

      n = SIZE(x)

      ALLOCATE(c(3, n))

      c(1, 1) = 0.0_wp
      c(2, 1) = 0.0_wp
      c(2, n) = 0.0_wp
      dxr = x(2) - x(1)
      dyr = (y(2) - y(1) ) / dxr

      DO  i = 2, n - 1
         dxl = dxr
         dxr = x(i + 1) - x(i)
         dyl = dyr
         dyr = (y(i + 1) - y(i) ) / dxr
         dx = dxr + dxl
         dy = (dyr - dyl) / dx
         c(1, i) = - dxr / (2.0_wp * dx + dxl * c(1, i - 1) )
         c(2, i) = (6.0_wp * dx * dy - dxl * c(2, i - 1) ) / (2.0_wp * dx +  dxl * c(1, i - 1) )
      END DO

      DO  i = n - 1, 2, - 1
         c(2, i) = c(1, i) * c(2, i + 1) + c(2, i)
      END DO

      DO i = 1, n - 1
         dx = x(i + 1) - x(i)
         dy = y(i + 1) - y(i)
         dyx = dy / dx
         c(1, i) = dyx - dx * (c(2, i) / 3.0_wp + c(2, i + 1) / 6.0_wp)
         c(2, i) = c(2, i) / 2.
         c(3, i) = (dy - c(1, i) * dx - c(2, i) * dx**2) / dx**3
      END DO

      IF (xx <= x(1) ) then
         der = c(1, 1)
         yy = y(1) + der * (xx - x(1) )
      ELSEIF (xx >= x(n) ) then
         d = x(n) - x(n - 1)
         der = c(1, n - 1) + 2.0_wp * c(2, n - 1) * d + 3.0_wp * c(3, n - 1) * d**2
         yy = y(n) + der * (xx - x(n) )
      ELSE
         DO i = 1, n
           IF (x(i) <= xx) m = i
         END DO
         d = xx - x(m)
         yy = y(m) + c(1, m) * d+ c(2, m) * d**2 + c(3, m) * d**3
         der = c(1, m) + 2.0_wp * c(2, m) * d+3.0_wp * c(3, m) * d**2
      ENDIF
      DEALLOCATE(c)

END SUBROUTINE splinx1