! This module encapsulates misc. math functions

MODULE lmdz_libmath_pch
  IMPLICIT NONE; PRIVATE
  PUBLIC pchfe_95, pchsp_95
CONTAINS
  REAL FUNCTION cbrt(x)
    ! Return the cubic root for positive or negative x
    REAL :: x

    cbrt = sign(1., x) * (abs(x)**(1. / 3.))
  END FUNCTION cbrt

  SUBROUTINE CHFEV(X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR)
    !***BEGIN PROLOGUE  CHFEV
    !***PURPOSE  Evaluate a cubic polynomial given in Hermite form at an
    ! array of points.  While designed for use by PCHFE, it may
    ! be useful directly as an evaluator for a piecewise cubic
    ! Hermite function in applications, such as graphing, where
    ! the interval is known in advance.
    !***LIBRARY   SLATEC (PCHIP)
    !***CATEGORY  E3
    !***TYPE      SINGLE PRECISION (CHFEV-S, DCHFEV-D)
    !***KEYWORDS  CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION,
    !  PCHIP
    !***AUTHOR  Fritsch, F. N., (LLNL)
    !  Lawrence Livermore National Laboratory
    !  P.O. Box 808  (L-316)
    !  Livermore, CA  94550
    !  FTS 532-4275, (510) 422-4275
    !***DESCRIPTION
    !
    !      CHFEV:  Cubic Hermite Function EValuator
    !
    ! Evaluates the cubic polynomial determined by function values
    ! F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points
    ! XE(J), J=1(1)NE.
    !
    ! ----------------------------------------------------------------------
    !
    !  Calling sequence:
    !
    !    INTEGER  NE, NEXT(2), IERR
    !    REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
    !
    !    CALL  CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)
    !
    !   Parameters:
    !
    ! X1,X2 -- (input) endpoints of interval of definition of cubic.
    !       (Error return if  X1.EQ.X2 .)
    !
    ! F1,F2 -- (input) values of function at X1 and X2, respectively.
    !
    ! D1,D2 -- (input) values of derivative at X1 and X2, respectively.
    !
    ! NE -- (input) number of evaluation points.  (Error return if
    !       NE.LT.1 .)
    !
    ! XE -- (input) real array of points at which the function is to be
    !       evaluated.  If any of the XE are outside the interval
    !       [X1,X2], a warning error is returned in NEXT.
    !
    ! FE -- (output) real array of values of the cubic function defined
    !       by  X1,X2, F1,F2, D1,D2  at the points  XE.
    !
    ! NEXT -- (output) integer array indicating number of extrapolation
    !       points:
    !        NEXT(1) = number of evaluation points to left of interval.
    !        NEXT(2) = number of evaluation points to right of interval.
    !
    ! IERR -- (output) error flag.
    !       Normal return:
    !          IERR = 0  (no errors).
    !       "Recoverable" errors:
    !          IERR = -1  if NE.LT.1 .
    !          IERR = -2  if X1.EQ.X2 .
    !            (The FE-array has not been changed in either case.)
    !
    !***REFERENCES  (NONE)
    !***ROUTINES CALLED  XERMSG
    !***REVISION HISTORY  (YYMMDD)
    !   811019  DATE WRITTEN
    !   820803  Minor cosmetic changes for release 1.
    !   890411  Added SAVE statements (Vers. 3.2).
    !   890531  Changed all specific intrinsics to generic.  (WRB)
    !   890703  Corrected category record.  (WRB)
    !   890703  REVISION DATE from Version 3.2
    !   891214  Prologue converted to Version 4.0 format.  (BAB)
    !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    !***END PROLOGUE  CHFEV
    !  Programming notes:
    !
    ! To produce a double precision version, simply:
    !    a. Change CHFEV to DCHFEV wherever it occurs,
    !    b. Change the real declaration to double precision, and
    !    c. Change the constant ZERO to double precision.
    !
    !  DECLARE ARGUMENTS.
    !
    INTEGER :: NE, NEXT(2), IERR
    REAL :: X1, X2, F1, F2, D1, D2, XE(*), FE(*)
    !
    !  DECLARE LOCAL VARIABLES.
    !
    INTEGER :: I
    REAL :: C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO
    SAVE ZERO
    DATA  ZERO /0./
    !
    !  VALIDITY-CHECK ARGUMENTS.
    !
    !***FIRST EXECUTABLE STATEMENT  CHFEV
    IF (NE < 1)  GO TO 5001
    H = X2 - X1
    IF (H == ZERO)  GO TO 5002
    !
    !  INITIALIZE.
    !
    IERR = 0
    NEXT(1) = 0
    NEXT(2) = 0
    XMI = MIN(ZERO, H)
    XMA = MAX(ZERO, H)
    !
    !  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
    !
    DELTA = (F2 - F1) / H
    DEL1 = (D1 - DELTA) / H
    DEL2 = (D2 - DELTA) / H
    ! (DELTA IS NO LONGER NEEDED.)
    C2 = -(DEL1 + DEL1 + DEL2)
    C3 = (DEL1 + DEL2) / H
    ! (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
    !
    !  EVALUATION LOOP.
    !
    DO I = 1, NE
      X = XE(I) - X1
      FE(I) = F1 + X * (D1 + X * (C2 + X * C3))
      ! COUNT EXTRAPOLATION POINTS.
      IF (X<XMI)  NEXT(1) = NEXT(1) + 1
      IF (X>XMA)  NEXT(2) = NEXT(2) + 1
      ! (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
    END DO
    !
    !  NORMAL RETURN.
    !
    RETURN
    !
    !  ERROR RETURNS.
    !
    5001   CONTINUE
    ! NE.LT.1 RETURN.
    IERR = -1
    CALL XERMSG ('SLATEC', 'CHFEV', &
            'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
    RETURN
    !
    5002   CONTINUE
    ! X1.EQ.X2 RETURN.
    IERR = -2
    CALL XERMSG ('SLATEC', 'CHFEV', 'INTERVAL ENDPOINTS EQUAL', IERR, &
            1)
    RETURN
    !------------- LAST LINE OF CHFEV FOLLOWS ------------------------------
  END SUBROUTINE CHFEV

  SUBROUTINE PCHFE(N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
    !***BEGIN PROLOGUE  PCHFE
    !***PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
    ! points.  May be used by itself for Hermite interpolation,
    ! or as an evaluator for PCHIM or PCHIC.
    !***LIBRARY   SLATEC (PCHIP)
    !***CATEGORY  E3
    !***TYPE      SINGLE PRECISION (PCHFE-S, DPCHFE-D)
    !***KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
    !  PIECEWISE CUBIC EVALUATION
    !***AUTHOR  Fritsch, F. N., (LLNL)
    !  Lawrence Livermore National Laboratory
    !  P.O. Box 808  (L-316)
    !  Livermore, CA  94550
    !  FTS 532-4275, (510) 422-4275
    !***DESCRIPTION
    !
    !      PCHFE:  Piecewise Cubic Hermite Function Evaluator
    !
    ! Evaluates the cubic Hermite function defined by  N, X, F, D  at
    ! the points  XE(J), J=1(1)NE.
    !
    ! To provide compatibility with PCHIM and PCHIC, includes an
    ! increment between successive values of the F- and D-arrays.
    !
    ! ----------------------------------------------------------------------
    !
    !  Calling sequence:
    !
    !    PARAMETER  (INCFD = ...)
    !    INTEGER  N, NE, IERR
    !    REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
    !    LOGICAL  SKIP
    !
    !    CALL  PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
    !
    !   Parameters:
    !
    ! N -- (input) number of data points.  (Error return if N.LT.2 .)
    !
    ! X -- (input) real array of independent variable values.  The
    !       elements of X must be strictly increasing:
    !            X(I-1) .LT. X(I),  I = 2(1)N.
    !       (Error return if not.)
    !
    ! F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
    !       the value corresponding to X(I).
    !
    ! D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
    !       the value corresponding to X(I).
    !
    ! INCFD -- (input) increment between successive values in F and D.
    !       (Error return if  INCFD.LT.1 .)
    !
    ! SKIP -- (input/output) logical variable which should be set to
    !       .TRUE. if the user wishes to skip checks for validity of
    !       preceding parameters, or to .FALSE. otherwise.
    !       This will save time in case these checks have already
    !       been performed (say, in PCHIM or PCHIC).
    !       SKIP will be set to .TRUE. on normal return.
    !
    ! NE -- (input) number of evaluation points.  (Error return if
    !       NE.LT.1 .)
    !
    ! XE -- (input) real array of points at which the function is to be
    !       evaluated.
    !
    !      NOTES:
    !       1. The evaluation will be most efficient if the elements
    !          of XE are increasing relative to X;
    !          that is,   XE(J) .GE. X(I)
    !          implies    XE(K) .GE. X(I),  all K.GE.J .
    !       2. If any of the XE are outside the interval [X(1),X(N)],
    !          values are extrapolated from the nearest extreme cubic,
    !          and a warning error is returned.
    !
    ! FE -- (output) real array of values of the cubic Hermite function
    !       defined by  N, X, F, D  at the points  XE.
    !
    ! IERR -- (output) error flag.
    !       Normal return:
    !          IERR = 0  (no errors).
    !       Warning error:
    !          IERR.GT.0  means that extrapolation was performed at
    !             IERR points.
    !       "Recoverable" errors:
    !          IERR = -1  if N.LT.2 .
    !          IERR = -2  if INCFD.LT.1 .
    !          IERR = -3  if the X-array is not strictly increasing.
    !          IERR = -4  if NE.LT.1 .
    !         (The FE-array has not been changed in any of these cases.)
    !           NOTE:  The above errors are checked in the order listed,
    !               and following arguments have **NOT** been validated.
    !
    !***REFERENCES  (NONE)
    !***ROUTINES CALLED  CHFEV, XERMSG
    !***REVISION HISTORY  (YYMMDD)
    !   811020  DATE WRITTEN
    !   820803  Minor cosmetic changes for release 1.
    !   870707  Minor cosmetic changes to prologue.
    !   890531  Changed all specific intrinsics to generic.  (WRB)
    !   890831  Modified array declarations.  (WRB)
    !   890831  REVISION DATE from Version 3.2
    !   891214  Prologue converted to Version 4.0 format.  (BAB)
    !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    !***END PROLOGUE  PCHFE
    !  Programming notes:
    !
    ! 1. To produce a double precision version, simply:
    !    a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they
    !       occur,
    !    b. Change the real declaration to double precision,
    !
    ! 2. Most of the coding between the CALL to CHFEV and the end of
    !    the IR-loop could be eliminated if it were permissible to
    !    assume that XE is ordered relative to X.
    !
    ! 3. CHFEV does not assume that X1 is less than X2.  thus, it would
    !    be possible to write a version of PCHFE that assumes a strict-
    !    ly decreasing X-array by simply running the IR-loop backwards
    !    (and reversing the order of appropriate tests).
    !
    ! 4. The present code has a minor bug, which I have decided is not
    !    worth the effort that would be required to fix it.
    !    If XE contains points in [X(N-1),X(N)], followed by points .LT.
    !    X(N-1), followed by points .GT.X(N), the extrapolation points
    !    will be counted (at least) twice in the total returned in IERR.
    !
    !  DECLARE ARGUMENTS.
    !
    INTEGER :: N, INCFD, NE, IERR
    REAL :: X(*), F(INCFD, *), D(INCFD, *), XE(*), FE(*)
    LOGICAL :: SKIP
    !
    !  DECLARE LOCAL VARIABLES.
    !
    INTEGER :: I, IERC, IR, J, JFIRST, NEXT(2), NJ
    !
    !  VALIDITY-CHECK ARGUMENTS.
    !
    !***FIRST EXECUTABLE STATEMENT  PCHFE
    IF (SKIP)  GO TO 5
    !
    IF (N<2)  GO TO 5001
    IF (INCFD<1)  GO TO 5002
    DO I = 2, N
      IF (X(I)<=X(I - 1))  GO TO 5003
    END DO
    !
    !  FUNCTION DEFINITION IS OK, GO ON.
    !
    5   CONTINUE
    IF (NE<1)  GO TO 5004
    IERR = 0
    SKIP = .TRUE.
    !
    !  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
    !                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
    JFIRST = 1
    IR = 2
    10   CONTINUE
    !
    ! SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
    !
    IF (JFIRST > NE)  GO TO 5000
    !
    ! LOCATE ALL POINTS IN INTERVAL.
    !
    DO J = JFIRST, NE
      IF (XE(J) >= X(IR))  GO TO 30
    END DO
    J = NE + 1
    GO TO 40
    !
    ! HAVE LOCATED FIRST POINT BEYOND INTERVAL.
    !
    30   CONTINUE
    IF (IR == N)  J = NE + 1
    !
    40   CONTINUE
    NJ = J - JFIRST
    !
    ! SKIP EVALUATION IF NO POINTS IN INTERVAL.
    !
    IF (NJ == 0)  GO TO 50
    !
    ! EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
    !
    !   ----------------------------------------------------------------
    CALL CHFEV (X(IR - 1), X(IR), F(1, IR - 1), F(1, IR), D(1, IR - 1), D(1, IR), &
            NJ, XE(JFIRST), FE(JFIRST), NEXT, IERC)
    ! ----------------------------------------------------------------
    IF (IERC < 0)  GO TO 5005
    !
    IF (NEXT(2) == 0)  GO TO 42
    ! IF (NEXT(2) .GT. 0)  THEN
    !    IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
    !    RIGHT OF X(IR).
    !
    IF (IR < N)  GO TO 41
    ! IF (IR .EQ. N)  THEN
    !    THESE ARE ACTUALLY EXTRAPOLATION POINTS.
    IERR = IERR + NEXT(2)
    GO TO 42
    41    CONTINUE
    ! ELSE
    !    WE SHOULD NEVER HAVE GOTTEN HERE.
    GO TO 5005
    ! ENDIF
    ! ENDIF
    42   CONTINUE
    !
    IF (NEXT(1) == 0)  GO TO 49
    ! IF (NEXT(1) .GT. 0)  THEN
    !    IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
    !    LEFT OF X(IR-1).
    !
    IF (IR > 2)  GO TO 43
    ! IF (IR .EQ. 2)  THEN
    !    THESE ARE ACTUALLY EXTRAPOLATION POINTS.
    IERR = IERR + NEXT(1)
    GO TO 49
    43    CONTINUE
    ! ELSE
    !    XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
    !    EVALUATION INTERVAL.
    !
    !          FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
    DO I = JFIRST, J - 1
      IF (XE(I) < X(IR - 1))  GO TO 45
    END DO
    ! NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
    !        IN CHFEV.
    GO TO 5005
    !
    45       CONTINUE
    ! RESET J.  (THIS WILL BE THE NEW JFIRST.)
    J = I
    !
    !          NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
    DO I = 1, IR - 1
      IF (XE(J) < X(I)) GO TO 47
    END DO
    ! NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
    !
    47       CONTINUE
    ! AT THIS POINT, EITHER  XE(J) .LT. X(1)
    !    OR      X(I-1) .LE. XE(J) .LT. X(I) .
    ! RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
    ! CYCLING.
    IR = MAX(1, I - 1)
    ! ENDIF
    ! ENDIF
    49   CONTINUE
    !
    JFIRST = J
    !
    ! END OF IR-LOOP.
    !
    50   CONTINUE
    IR = IR + 1
    IF (IR <= N)  GO TO 10
    !
    !  NORMAL RETURN.
    !
    5000   CONTINUE
    RETURN
    !
    !  ERROR RETURNS.
    !
    5001   CONTINUE
    ! N.LT.2 RETURN.
    IERR = -1
    CALL XERMSG ('SLATEC', 'PCHFE', &
            'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
    RETURN
    !
    5002   CONTINUE
    ! INCFD.LT.1 RETURN.
    IERR = -2
    CALL XERMSG ('SLATEC', 'PCHFE', 'INCREMENT LESS THAN ONE', IERR, &
            1)
    RETURN
    !
    5003   CONTINUE
    ! X-ARRAY NOT STRICTLY INCREASING.
    IERR = -3
    CALL XERMSG ('SLATEC', 'PCHFE', 'X-ARRAY NOT STRICTLY INCREASING' &
            , IERR, 1)
    RETURN
    !
    5004   CONTINUE
    ! NE.LT.1 RETURN.
    IERR = -4
    CALL XERMSG ('SLATEC', 'PCHFE', &
            'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
    RETURN
    !
    5005   CONTINUE
    ! ERROR RETURN FROM CHFEV.
    !   *** THIS CASE SHOULD NEVER OCCUR ***
    IERR = -5
    CALL XERMSG ('SLATEC', 'PCHFE', &
            'ERROR RETURN FROM CHFEV -- FATAL', IERR, 2)
    RETURN
    !------------- LAST LINE OF PCHFE FOLLOWS ------------------------------
  END SUBROUTINE PCHFE

  SUBROUTINE PCHFE_95(X, F, D, SKIP, XE, FE, IERR)

    ! PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
    !            points.  May be used by itself for Hermite interpolation,
    !            or as an evaluator for PCHIM or PCHIC.
    ! CATEGORY  E3
    ! KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
    !             PIECEWISE CUBIC EVALUATION

    !          PCHFE:  Piecewise Cubic Hermite Function Evaluator
    ! Evaluates the cubic Hermite function defined by  X, F, D  at
    ! the points  XE.

    use lmdz_assert_eq, only: assert_eq

    REAL, intent(in) :: X(:) ! real array of independent variable values
    ! The elements of X must be strictly increasing.

    REAL, intent(in) :: F(:) ! real array of function values
    ! F(I) is the value corresponding to X(I).

    REAL, intent(in) :: D(:) ! real array of derivative values
    ! D(I) is the value corresponding to X(I).

    LOGICAL, intent(inout) :: SKIP
    ! request to skip checks for validity of "x"
    ! If "skip" is false then "pchfe" will check that size(x) >= 2 and
    ! "x" is in strictly ascending order.
    ! Setting "skip" to true will save time in case these checks have
    ! already been performed (say, in "PCHIM" or "PCHIC").
    ! "SKIP" will be set to TRUE on normal return.

    real, intent(in) :: XE(:) ! points at which the function is to be evaluated
    ! NOTES:
    ! 1. The evaluation will be most efficient if the elements of XE
    ! are increasing relative to X.
    ! That is,   XE(J) .GE. X(I)
    ! implies    XE(K) .GE. X(I),  all K.GE.J
    ! 2. If any of the XE are outside the interval [X(1),X(N)], values
    ! are extrapolated from the nearest extreme cubic, and a warning
    ! error is returned.

    real, intent(out) :: FE(:) ! values of the cubic Hermite function
    ! defined by X, F, D at the points XE

    integer, intent(out) :: IERR ! error flag
    ! Normal return:
    ! IERR = 0  no error
    ! Warning error:
    ! IERR > 0  means that extrapolation was performed at IERR points
    ! "Recoverable" errors:
    !              IERR = -1  if N < 2
    !              IERR = -3  if the X-array is not strictly increasing
    !              IERR = -4  if NE < 1
    ! NOTE: The above errors are checked in the order listed, and
    ! following arguments have **NOT** been validated.

    ! Variables local to the procedure:

    INTEGER  N, NE

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

    n = assert_eq(size(x), size(f), size(d), "PCHFE_95 n")
    ne = assert_eq(size(xe), size(fe), "PCHFE_95 ne")
    CALL PCHFE(N, X, F, D, 1, SKIP, NE, XE, FE, IERR)

  END SUBROUTINE  PCHFE_95

  FUNCTION pchsp_95(x, f, ibeg, iend, vc_beg, vc_end)

    ! PURPOSE: Set derivatives needed to determine the Hermite
    ! representation of the cubic spline interpolant to given data,
    ! with specified boundary conditions.

    ! Part of the "pchip" package.

    ! CATEGORY: E1A

    ! KEYWORDS: cubic hermite interpolation, piecewise cubic
    ! interpolation, spline interpolation

    ! DESCRIPTION: "pchsp" stands for "Piecewise Cubic Hermite Spline"
    ! Computes the Hermite representation of the cubic spline
    ! interpolant to the data given in X and F satisfying the boundary
    ! conditions specified by Ibeg, iend, vc_beg and VC_end.

    ! The resulting piecewise cubic Hermite function may be evaluated
    ! by "pchfe" or "pchfd".

    ! NOTE: This is a modified version of C. de Boor's cubic spline
    ! routine "cubspl".

    ! REFERENCE: Carl de Boor, A Practical Guide to Splines, Springer,
    ! 2001, pages 43-47

    use lmdz_assert_eq, only: assert_eq

    real, intent(in) :: x(:)
    ! independent variable values
    ! The elements of X must be strictly increasing:
    !                X(I-1) < X(I),  I = 2...N.
    !           (Error return if not.)
    ! (error if size(x) < 2)

    real, intent(in) :: f(:)
    !     dependent variable values to be interpolated
    !  F(I) is value corresponding to X(I).

    INTEGER, intent(in) :: ibeg
    !     desired boundary condition at beginning of data

    !        IBEG = 0  to set pchsp_95(1) so that the third derivative is con-
    !              tinuous at X(2).  This is the "not a knot" condition
    !              provided by de Boor's cubic spline routine CUBSPL.
    !              This is the default boundary condition.
    !        IBEG = 1  if first derivative at X(1) is given in VC_BEG.
    !        IBEG = 2  if second derivative at X(1) is given in VC_BEG.
    !        IBEG = 3  to use the 3-point difference formula for pchsp_95(1).
    !              (Reverts to the default boundary condition if size(x) < 3 .)
    !        IBEG = 4  to use the 4-point difference formula for pchsp_95(1).
    !              (Reverts to the default boundary condition if size(x) < 4 .)

    !          NOTES:
    !           1. An error return is taken if IBEG is out of range.
    !           2. For the "natural" boundary condition, use IBEG=2 and
    !              VC_BEG=0.

    INTEGER, intent(in) :: iend
    !           IC(2) = IEND, desired condition at end of data.
    !  IEND may take on the same values as IBEG, but applied to
    !  derivative at X(N). In case IEND = 1 or 2, The value is given in VC_END.

    !          NOTES:
    !           1. An error return is taken if IEND is out of range.
    !           2. For the "natural" boundary condition, use IEND=2 and
    !              VC_END=0.

    REAL, intent(in), optional :: vc_beg
    ! desired boundary value, as indicated above.
    !           VC_BEG need be set only if IBEG = 1 or 2 .

    REAL, intent(in), optional :: vc_end
    ! desired boundary value, as indicated above.
    !           VC_END need be set only if Iend = 1 or 2 .

    real pchsp_95(size(x))
    ! derivative values at the data points
    !           These values will determine the cubic spline interpolant
    !           with the requested boundary conditions.
    !           The value corresponding to X(I) is stored in
    !                PCHSP_95(I),  I=1...N.

    ! LOCAL VARIABLES:
    real wk(2, size(x)) ! real array of working storage
    INTEGER n ! number of data points
    integer ierr, ic(2)
    real vc(2)

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

    n = assert_eq(size(x), size(f), "pchsp_95 n")
    if ((ibeg == 1 .or. ibeg == 2) .and. .not. present(vc_beg)) then
      print *, "vc_beg required for IBEG = 1 or 2"
      stop 1
    end if
    if ((iend == 1 .or. iend == 2) .and. .not. present(vc_end)) then
      print *, "vc_end required for IEND = 1 or 2"
      stop 1
    end if
    ic = (/ibeg, iend/)
    if (present(vc_beg)) vc(1) = vc_beg
    if (present(vc_end)) vc(2) = vc_end
    CALL PCHSP(IC, VC, N, X, F, pchsp_95, 1, WK, size(WK), IERR)
    if (ierr /= 0) stop 1

  END FUNCTION pchsp_95

  REAL FUNCTION PCHDF(K, X, S, IERR)
    !***BEGIN PROLOGUE  PCHDF
    !***SUBSIDIARY
    !***PURPOSE  Computes divided differences for PCHCE and PCHSP
    !***LIBRARY   SLATEC (PCHIP)
    !***TYPE      SINGLE PRECISION (PCHDF-S, DPCHDF-D)
    !***AUTHOR  Fritsch, F. N., (LLNL)
    !***DESCRIPTION
    !
    !      PCHDF:   PCHIP Finite Difference Formula
    !
    ! Uses a divided difference formulation to compute a K-point approx-
    ! imation to the derivative at X(K) based on the data in X and S.
    !
    ! Called by  PCHCE  and  PCHSP  to compute 3- and 4-point boundary
    ! derivative approximations.
    !
    ! ----------------------------------------------------------------------
    !
    ! On input:
    !    K      is the order of the desired derivative approximation.
    !           K must be at least 3 (error return if not).
    !    X      contains the K values of the independent variable.
    !           X need not be ordered, but the values **MUST** be
    !           distinct.  (Not checked here.)
    !    S      contains the associated slope values:
    !              S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1.
    !           (Note that S need only be of length K-1.)
    !
    ! On return:
    !    S      will be destroyed.
    !    IERR   will be set to -1 if K.LT.2 .
    !    PCHDF  will be set to the desired derivative approximation if
    !           IERR=0 or to zero if IERR=-1.
    !
    ! ----------------------------------------------------------------------
    !
    !***SEE ALSO  PCHCE, PCHSP
    !***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
    !             Verlag, New York, 1978, pp. 10-16.
    !***ROUTINES CALLED  XERMSG
    !***REVISION HISTORY  (YYMMDD)
    !   820503  DATE WRITTEN
    !   820805  Converted to SLATEC library version.
    !   870813  Minor cosmetic changes.
    !   890411  Added SAVE statements (Vers. 3.2).
    !   890411  REVISION DATE from Version 3.2
    !   891214  Prologue converted to Version 4.0 format.  (BAB)
    !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    !   900328  Added TYPE section.  (WRB)
    !   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
    !   920429  Revised format and order of references.  (WRB,FNF)
    !   930503  Improved purpose.  (FNF)
    !***END PROLOGUE  PCHDF
    !
    !**End
    !
    !  DECLARE ARGUMENTS.
    !
    INTEGER :: K, IERR
    REAL :: X(K), S(K)
    !
    !  DECLARE LOCAL VARIABLES.
    !
    INTEGER :: I, J
    REAL :: VALUE, ZERO
    SAVE ZERO
    DATA  ZERO /0./
    !
    !  CHECK FOR LEGAL VALUE OF K.
    !
    !***FIRST EXECUTABLE STATEMENT  PCHDF
    IF (K < 3)  GO TO 5001
    !
    !  COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL.
    !
    DO J = 2, K - 1
      DO I = 1, K - J
        S(I) = (S(I + 1) - S(I)) / (X(I + J) - X(I))
      END DO
    END DO
    !
    !  EVALUATE DERIVATIVE AT X(K).
    !
    VALUE = S(1)
    DO I = 2, K - 1
      VALUE = S(I) + VALUE * (X(K) - X(I))
    END DO
    !
    !  NORMAL RETURN.
    !
    IERR = 0
    PCHDF = VALUE
    RETURN
    !
    !  ERROR RETURN.
    !
    5001   CONTINUE
    ! K.LT.3 RETURN.
    IERR = -1
    CALL XERMSG ('SLATEC', 'PCHDF', 'K LESS THAN THREE', IERR, 1)
    PCHDF = ZERO
    RETURN
    !------------- LAST LINE OF PCHDF FOLLOWS ------------------------------
  END FUNCTION PCHDF

  SUBROUTINE PCHSP(IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
    !***BEGIN PROLOGUE  PCHSP
    !***PURPOSE  Set derivatives needed to determine the Hermite represen-
    ! tation of the cubic spline interpolant to given data, with
    ! specified boundary conditions.
    !***LIBRARY   SLATEC (PCHIP)
    !***CATEGORY  E1A
    !***TYPE      SINGLE PRECISION (PCHSP-S, DPCHSP-D)
    !***KEYWORDS  CUBIC HERMITE INTERPOLATION, PCHIP,
    !  PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION
    !***AUTHOR  Fritsch, F. N., (LLNL)
    !  Lawrence Livermore National Laboratory
    !  P.O. Box 808  (L-316)
    !  Livermore, CA  94550
    !  FTS 532-4275, (510) 422-4275
    !***DESCRIPTION
    !
    !      PCHSP:   Piecewise Cubic Hermite Spline
    !
    ! Computes the Hermite representation of the cubic spline inter-
    ! polant to the data given in X and F satisfying the boundary
    ! conditions specified by IC and VC.
    !
    ! To facilitate two-dimensional applications, includes an increment
    ! between successive values of the F- and D-arrays.
    !
    ! The resulting piecewise cubic Hermite function may be evaluated
    ! by PCHFE or PCHFD.
    !
    ! NOTE:  This is a modified version of C. de Boor's cubic spline
    !        routine CUBSPL.
    !
    ! ----------------------------------------------------------------------
    !
    !  Calling sequence:
    !
    !    PARAMETER  (INCFD = ...)
    !    INTEGER  IC(2), N, NWK, IERR
    !    REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
    !
    !    CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
    !
    !   Parameters:
    !
    ! IC -- (input) integer array of length 2 specifying desired
    !       boundary conditions:
    !       IC(1) = IBEG, desired condition at beginning of data.
    !       IC(2) = IEND, desired condition at end of data.
    !
    !       IBEG = 0  to set D(1) so that the third derivative is con-
    !          tinuous at X(2).  This is the "not a knot" condition
    !          provided by de Boor's cubic spline routine CUBSPL.
    !          < This is the default boundary condition. >
    !       IBEG = 1  if first derivative at X(1) is given in VC(1).
    !       IBEG = 2  if second derivative at X(1) is given in VC(1).
    !       IBEG = 3  to use the 3-point difference formula for D(1).
    !                 (Reverts to the default b.c. if N.LT.3 .)
    !       IBEG = 4  to use the 4-point difference formula for D(1).
    !                 (Reverts to the default b.c. if N.LT.4 .)
    !      NOTES:
    !       1. An error return is taken if IBEG is out of range.
    !       2. For the "natural" boundary condition, use IBEG=2 and
    !          VC(1)=0.
    !
    !       IEND may take on the same values as IBEG, but applied to
    !       derivative at X(N).  In case IEND = 1 or 2, the value is
    !       given in VC(2).
    !
    !      NOTES:
    !       1. An error return is taken if IEND is out of range.
    !       2. For the "natural" boundary condition, use IEND=2 and
    !          VC(2)=0.
    !
    ! VC -- (input) real array of length 2 specifying desired boundary
    !       values, as indicated above.
    !       VC(1) need be set only if IC(1) = 1 or 2 .
    !       VC(2) need be set only if IC(2) = 1 or 2 .
    !
    ! N -- (input) number of data points.  (Error return if N.LT.2 .)
    !
    ! X -- (input) real array of independent variable values.  The
    !       elements of X must be strictly increasing:
    !            X(I-1) .LT. X(I),  I = 2(1)N.
    !       (Error return if not.)
    !
    ! F -- (input) real array of dependent variable values to be inter-
    !       polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
    !
    ! D -- (output) real array of derivative values at the data points.
    !       These values will determine the cubic spline interpolant
    !       with the requested boundary conditions.
    !       The value corresponding to X(I) is stored in
    !            D(1+(I-1)*INCFD),  I=1(1)N.
    !       No other entries in D are changed.
    !
    ! INCFD -- (input) increment between successive values in F and D.
    !       This argument is provided primarily for 2-D applications.
    !       (Error return if  INCFD.LT.1 .)
    !
    ! WK -- (scratch) real array of working storage.
    !
    ! NWK -- (input) length of work array.
    !       (Error return if NWK.LT.2*N .)
    !
    ! IERR -- (output) error flag.
    !       Normal return:
    !          IERR = 0  (no errors).
    !       "Recoverable" errors:
    !          IERR = -1  if N.LT.2 .
    !          IERR = -2  if INCFD.LT.1 .
    !          IERR = -3  if the X-array is not strictly increasing.
    !          IERR = -4  if IBEG.LT.0 or IBEG.GT.4 .
    !          IERR = -5  if IEND.LT.0 of IEND.GT.4 .
    !          IERR = -6  if both of the above are true.
    !          IERR = -7  if NWK is too small.
    !           NOTE:  The above errors are checked in the order listed,
    !               and following arguments have **NOT** been validated.
    !         (The D-array has not been changed in any of these cases.)
    !          IERR = -8  in case of trouble solving the linear system
    !                     for the interior derivative values.
    !         (The D-array may have been changed in this case.)
    !         (             Do **NOT** use it!                )
    !
    !***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
    !             Verlag, New York, 1978, pp. 53-59.
    !***ROUTINES CALLED  PCHDF, XERMSG
    !***REVISION HISTORY  (YYMMDD)
    !   820503  DATE WRITTEN
    !   820804  Converted to SLATEC library version.
    !   870707  Minor cosmetic changes to prologue.
    !   890411  Added SAVE statements (Vers. 3.2).
    !   890703  Corrected category record.  (WRB)
    !   890831  Modified array declarations.  (WRB)
    !   890831  REVISION DATE from Version 3.2
    !   891214  Prologue converted to Version 4.0 format.  (BAB)
    !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    !   920429  Revised format and order of references.  (WRB,FNF)
    !***END PROLOGUE  PCHSP
    !  Programming notes:
    !
    ! To produce a double precision version, simply:
    !    a. Change PCHSP to DPCHSP wherever it occurs,
    !    b. Change the real declarations to double precision, and
    !    c. Change the constants ZERO, HALF, ... to double precision.
    !
    !  DECLARE ARGUMENTS.
    !
    INTEGER :: IC(2), N, INCFD, NWK, IERR
    REAL :: VC(2), X(*), F(INCFD, *), D(INCFD, *), WK(2, *)
    !
    !  DECLARE LOCAL VARIABLES.
    !
    INTEGER :: IBEG, IEND, INDEX, J, NM1
    REAL :: G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO
    SAVE ZERO, HALF, ONE, TWO, THREE
    REAL :: PCHDF
    !
    DATA  ZERO /0./, HALF /0.5/, ONE /1./, TWO /2./, THREE /3./
    !
    !  VALIDITY-CHECK ARGUMENTS.
    !
    !***FIRST EXECUTABLE STATEMENT  PCHSP
    IF (N<2)  GO TO 5001
    IF (INCFD<1)  GO TO 5002
    DO J = 2, N
      IF (X(J)<=X(J - 1))  GO TO 5003
    END DO
    !
    IBEG = IC(1)
    IEND = IC(2)
    IERR = 0
    IF ((IBEG<0).OR.(IBEG>4))  IERR = IERR - 1
    IF ((IEND<0).OR.(IEND>4))  IERR = IERR - 2
    IF (IERR<0)  GO TO 5004
    !
    !  FUNCTION DEFINITION IS OK -- GO ON.
    !
    IF (NWK < 2 * N)  GO TO 5007
    !
    !  COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO,
    !  COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.).
    DO J = 2, N
      WK(1, J) = X(J) - X(J - 1)
      WK(2, J) = (F(1, J) - F(1, J - 1)) / WK(1, J)
    END DO
    !
    !  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
    !
    IF (IBEG>N)  IBEG = 0
    IF (IEND>N)  IEND = 0
    !
    !  SET UP FOR BOUNDARY CONDITIONS.
    !
    IF ((IBEG==1).OR.(IBEG==2))  THEN
      D(1, 1) = VC(1)
    ELSE IF (IBEG > 2)  THEN
      ! PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
      DO J = 1, IBEG
        INDEX = IBEG - J + 1
        ! INDEX RUNS FROM IBEG DOWN TO 1.
        XTEMP(J) = X(INDEX)
        IF (J < IBEG)  STEMP(J) = WK(2, INDEX)
      END DO
      ! --------------------------------
      D(1, 1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
      ! --------------------------------
      IF (IERR /= 0)  GO TO 5009
      IBEG = 1
    ENDIF
    !
    IF ((IEND==1).OR.(IEND==2))  THEN
      D(1, N) = VC(2)
    ELSE IF (IEND > 2)  THEN
      ! PICK UP LAST IEND POINTS.
      DO J = 1, IEND
        INDEX = N - IEND + J
        ! INDEX RUNS FROM N+1-IEND UP TO N.
        XTEMP(J) = X(INDEX)
        IF (J < IEND)  STEMP(J) = WK(2, INDEX + 1)
      END DO
      ! --------------------------------
      D(1, N) = PCHDF (IEND, XTEMP, STEMP, IERR)
      ! --------------------------------
      IF (IERR /= 0)  GO TO 5009
      IEND = 1
    ENDIF
    !
    ! --------------------( BEGIN CODING FROM CUBSPL )--------------------
    !
    !  **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF
    !  F  AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM-
    !  INATION, WITH S(J) ENDING UP IN D(1,J), ALL J.
    ! WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE.
    !
    !  CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM
    !         WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1)
    !
    IF (IBEG == 0)  THEN
      IF (N == 2)  THEN
        ! NO CONDITION AT LEFT END AND N = 2.
        WK(2, 1) = ONE
        WK(1, 1) = ONE
        D(1, 1) = TWO * WK(2, 2)
      ELSE
        ! NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2.
        WK(2, 1) = WK(1, 3)
        WK(1, 1) = WK(1, 2) + WK(1, 3)
        D(1, 1) = ((WK(1, 2) + TWO * WK(1, 1)) * WK(2, 2) * WK(1, 3) &
                + WK(1, 2)**2 * WK(2, 3)) / WK(1, 1)
      ENDIF
    ELSE IF (IBEG == 1)  THEN
      ! SLOPE PRESCRIBED AT LEFT END.
      WK(2, 1) = ONE
      WK(1, 1) = ZERO
    ELSE
      ! SECOND DERIVATIVE PRESCRIBED AT LEFT END.
      WK(2, 1) = TWO
      WK(1, 1) = ONE
      D(1, 1) = THREE * WK(2, 2) - HALF * WK(1, 2) * D(1, 1)
    ENDIF
    !
    !  IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND
    !  CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH
    !  EQUATION READS    WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J).
    !
    NM1 = N - 1
    IF (NM1 > 1)  THEN
      DO J = 2, NM1
        IF (WK(2, J - 1) == ZERO)  GO TO 5008
        G = -WK(1, J + 1) / WK(2, J - 1)
        D(1, J) = G * D(1, J - 1) &
                + THREE * (WK(1, J) * WK(2, J + 1) + WK(1, J + 1) * WK(2, J))
        WK(2, J) = G * WK(1, J - 1) + TWO * (WK(1, J) + WK(1, J + 1))
      END DO
    ENDIF
    !
    !  CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM
    !       (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N)
    !
    ! IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK-
    ! SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT
    ! AT THIS POINT.
    IF (IEND == 1)  GO TO 30
    !
    IF (IEND == 0)  THEN
      IF (N==2 .AND. IBEG==0)  THEN
        ! NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2.
        D(1, 2) = WK(2, 2)
        GO TO 30
      ELSE IF ((N==2) .OR. (N==3 .AND. IBEG==0))  THEN
        ! EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*
        ! NOT-A-KNOT AT LEFT END POINT).
        D(1, N) = TWO * WK(2, N)
        WK(2, N) = ONE
        IF (WK(2, N - 1) == ZERO)  GO TO 5008
        G = -ONE / WK(2, N - 1)
      ELSE
        ! NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR  ALSO NOT-A-
        ! KNOT AT LEFT END POINT.
        G = WK(1, N - 1) + WK(1, N)
        ! DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES).
        D(1, N) = ((WK(1, N) + TWO * G) * WK(2, N) * WK(1, N - 1) &
                + WK(1, N)**2 * (F(1, N - 1) - F(1, N - 2)) / WK(1, N - 1)) / G
        IF (WK(2, N - 1) == ZERO)  GO TO 5008
        G = -G / WK(2, N - 1)
        WK(2, N) = WK(1, N - 1)
      ENDIF
    ELSE
      ! SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT.
      D(1, N) = THREE * WK(2, N) + HALF * WK(1, N) * D(1, N)
      WK(2, N) = TWO
      IF (WK(2, N - 1) == ZERO)  GO TO 5008
      G = -ONE / WK(2, N - 1)
    ENDIF
    !
    !  COMPLETE FORWARD PASS OF GAUSS ELIMINATION.
    !
    WK(2, N) = G * WK(1, N - 1) + WK(2, N)
    IF (WK(2, N) == ZERO)   GO TO 5008
    D(1, N) = (G * D(1, N - 1) + D(1, N)) / WK(2, N)
    !
    !  CARRY OUT BACK SUBSTITUTION
    !
    30   CONTINUE
    DO J = NM1, 1, -1
      IF (WK(2, J) == ZERO)  GO TO 5008
      D(1, J) = (D(1, J) - WK(1, J) * D(1, J + 1)) / WK(2, J)
    END DO
    ! --------------------(  END  CODING FROM CUBSPL )--------------------
    !
    !  NORMAL RETURN.
    !
    RETURN
    !
    !  ERROR RETURNS.
    !
    5001   CONTINUE
    ! N.LT.2 RETURN.
    IERR = -1
    CALL XERMSG ('SLATEC', 'PCHSP', &
            'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
    RETURN
    !
    5002   CONTINUE
    ! INCFD.LT.1 RETURN.
    IERR = -2
    CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR, &
            1)
    RETURN
    !
    5003   CONTINUE
    ! X-ARRAY NOT STRICTLY INCREASING.
    IERR = -3
    CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING' &
            , IERR, 1)
    RETURN
    !
    5004   CONTINUE
    ! IC OUT OF RANGE RETURN.
    IERR = IERR - 3
    CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1)
    RETURN
    !
    5007   CONTINUE
    ! NWK TOO SMALL RETURN.
    IERR = -7
    CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1)
    RETURN
    !
    5008   CONTINUE
    ! SINGULAR SYSTEM.
    !   *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES   ***
    !   *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). ***
    IERR = -8
    CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR, &
            1)
    RETURN
    !
    5009   CONTINUE
    ! ERROR RETURN FROM PCHDF.
    !   *** THIS CASE SHOULD NEVER OCCUR ***
    IERR = -9
    CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR, &
            1)
    RETURN
    !------------- LAST LINE OF PCHSP FOLLOWS ------------------------------
  END SUBROUTINE PCHSP


END MODULE lmdz_libmath_pch