Index: LMDZ6/branches/Amaury_dev/libf/misc/lmdz_libmath_pch.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/misc/lmdz_libmath_pch.f90	(revision 5195)
+++ 	(revision )
@@ -1,1150 +1,0 @@
-! This module encapsulates misc. math functions
-
-MODULE lmdz_libmath_pch
-  USE lmdz_xer, ONLY: xermsg
-  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
-
-    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
Index: LMDZ6/branches/Amaury_dev/libf/misc/lmdz_xer.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/misc/lmdz_xer.f90	(revision 5195)
+++ 	(revision )
@@ -1,1123 +1,0 @@
-! Contains "legacy" xer* functions required by some other very old external code
-
-MODULE lmdz_xer
-  IMPLICIT NONE; PRIVATE
-  PUBLIC xermsg
-CONTAINS
-
-  !DECK XERMSG
-  SUBROUTINE XERMSG(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
-    IMPLICIT NONE
-    !***BEGIN PROLOGUE  XERMSG
-    !***PURPOSE  Process error messages for SLATEC and other libraries.
-    !***LIBRARY   SLATEC (XERROR)
-    !***CATEGORY  R3C
-    !***TYPE      ALL (XERMSG-A)
-    !***KEYWORDS  ERROR MESSAGE, XERROR
-    !***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
-    !***DESCRIPTION
-
-    !   XERMSG processes a diagnostic message in a manner determined by the
-    !   value of LEVEL and the current value of the library error control
-    !   flag, KONTRL.  See SUBROUTINE XSETF for details.
-
-    !    LIBRAR   A character constant (or character variable) with the name
-    !         of the library.  This will be 'SLATEC' for the SLATEC
-    !         Common Math Library.  The error handling package is
-    !         general enough to be used by many libraries
-    !         simultaneously, so it is desirable for the routine that
-    !         detects and reports an error to identify the library name
-    !         as well as the routine name.
-
-    !    SUBROU   A character constant (or character variable) with the name
-    !             of the routine that detected the error.  Usually it is the
-    !         name of the routine that is calling XERMSG.  There are
-    !         some instances where a user callable library routine calls
-    !         lower level subsidiary routines where the error is
-    !         detected.  In such cases it may be more informative to
-    !         supply the name of the routine the user called rather than
-    !         the name of the subsidiary routine that detected the
-    !         error.
-
-    !    MESSG    A character constant (or character variable) with the text
-    !         of the error or warning message.  In the example below,
-    !         the message is a character constant that contains a
-    !         generic message.
-
-    !               CALL XERMSG ('SLATEC', 'MMPY',
-    !              *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION',
-    !              *3, 1)
-
-    !         It is possible (and is sometimes desirable) to generate a
-    !         specific message--e.g., one that contains actual numeric
-    !         values.  Specific numeric values can be converted into
-    !         character strings using formatted WRITE statements into
-    !         character variables.  This is called standard Fortran
-    !         internal file I/O and is exemplified in the first three
-    !         lines of the following example.  You can also catenate
-    !         substrings of characters to construct the error message.
-    !         Here is an example showing the use of both writing to
-    !         an internal file and catenating character strings.
-
-    !               CHARACTER*5 CHARN, CHARL
-    !               WRITE (CHARN,10) N
-    !               WRITE (CHARL,10) LDA
-    !            10 FORMAT(I5)
-    !               CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN//
-    !              *   ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'//
-    !              *   CHARL, 3, 1)
-
-    !         There are two subtleties worth mentioning.  One is that
-    !         the // for character catenation is used to construct the
-    !         error message so that no single character constant is
-    !         continued to the next line.  This avoids confusion as to
-    !         whether there are trailing blanks at the end of the line.
-    !         The second is that by catenating the parts of the message
-    !             as an actual argument rather than encoding the entire
-    !         message into one large character variable, we avoid
-    !         having to know how long the message will be in order to
-    !         declare an adequate length for that large character
-    !         variable.  XERMSG calls XERPRN to print the message using
-    !         multiple lines if necessary.  If the message is very long,
-    !         XERPRN will break it into pieces of 72 characters (as
-    !         requested by XERMSG) for printing on multiple lines.
-    !         Also, XERMSG asks XERPRN to prefix each line with ' *  '
-    !         so that the total line length could be 76 characters.
-    !         Note also that XERPRN scans the error message backwards
-    !         to ignore trailing blanks.  Another feature is that
-    !         the substring '$$' is treated as a new line sentinel
-    !         by XERPRN.  If you want to construct a multiline
-    !         message without having to count out multiples of 72
-    !         characters, just use '$$' as a separator.  '$$'
-    !         obviously must occur within 72 characters of the
-    !         start of each line to have its intended effect since
-    !         XERPRN is asked to wrap around at 72 characters in
-    !         addition to looking for '$$'.
-
-    !    NERR     An integer value that is chosen by the library routine's
-    !         author.  It must be in the range -99 to 999 (three
-    !         printable digits).  Each distinct error should have its
-    !         own error number.  These error numbers should be described
-    !         in the machine readable documentation for the routine.
-    !         The error numbers need be unique only within each routine,
-    !         so it is reasonable for each routine to start enumerating
-    !         errors from 1 and proceeding to the next integer.
-
-    !    LEVEL    An integer value in the range 0 to 2 that indicates the
-    !         level (severity) of the error.  Their meanings are
-
-    !        -1  A warning message.  This is used if it is not clear
-    !            that there really is an error, but the user's attention
-    !            may be needed.  An attempt is made to only print this
-    !            message once.
-
-    !         0  A warning message.  This is used if it is not clear
-    !            that there really is an error, but the user's attention
-    !            may be needed.
-
-    !         1  A recoverable error.  This is used even if the error is
-    !            so serious that the routine cannot return any useful
-    !            answer.  If the user has told the error package to
-    !            return after recoverable errors, then XERMSG will
-    !            return to the Library routine which can then return to
-    !            the user's routine.  The user may also permit the error
-    !            package to terminate the program upon encountering a
-    !            recoverable error.
-
-    !         2  A fatal error.  XERMSG will not return to its caller
-    !            after it receives a fatal error.  This level should
-    !            hardly ever be used; it is much better to allow the
-    !            user a chance to recover.  An example of one of the few
-    !            cases in which it is permissible to declare a level 2
-    !            error is a reverse communication Library routine that
-    !            is likely to be called repeatedly until it integrates
-    !            across some interval.  If there is a serious error in
-    !            the input such that another step cannot be taken and
-    !            the Library routine is called again without the input
-    !            error having been corrected by the caller, the Library
-    !            routine will probably be called forever with improper
-    !            input.  In this case, it is reasonable to declare the
-    !            error to be fatal.
-
-    !    Each of the arguments to XERMSG is input; none will be modified by
-    !    XERMSG.  A routine may make multiple calls to XERMSG with warning
-    !    level messages; however, after a CALL to XERMSG with a recoverable
-    !    error, the routine should return to the user.  Do not try to call
-    !    XERMSG with a second recoverable error after the first recoverable
-    !    error because the error package saves the error number.  The user
-    !    can retrieve this error number by calling another entry point in
-    !    the error handling package and then clear the error number when
-    !    recovering from the error.  Calling XERMSG in succession causes the
-    !    old error number to be overwritten by the latest error number.
-    !    This is considered harmless for error numbers associated with
-    !    warning messages but must not be done for error numbers of serious
-    !    errors.  After a CALL to XERMSG with a recoverable error, the user
-    !    must be given a chance to CALL NUMXER or XERCLR to retrieve or
-    !    clear the error number.
-    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
-    !             Error-handling Package, SAND82-0800, Sandia
-    !             Laboratories, 1982.
-    !***ROUTINES CALLED  FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE
-    !***REVISION HISTORY  (YYMMDD)
-    !   880101  DATE WRITTEN
-    !   880621  REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988.
-    !       THERE ARE TWO BASIC CHANGES.
-    !       1.  A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO
-    !           PRINT MESSAGES.  THIS ROUTINE WILL BREAK LONG MESSAGES
-    !           INTO PIECES FOR PRINTING ON MULTIPLE LINES.  '$$' IS
-    !           ACCEPTED AS A NEW LINE SENTINEL.  A PREFIX CAN BE
-    !           ADDED TO EACH LINE TO BE PRINTED.  XERMSG USES EITHER
-    !           ' ***' OR ' *  ' AND LONG MESSAGES ARE BROKEN EVERY
-    !           72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE
-    !           LENGTH OUTPUT CAN NOW BE AS GREAT AS 76.
-    !       2.  THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE
-    !           FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE
-    !           OF LOWER CASE.
-    !   880708  REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30.
-    !       THE PRINCIPAL CHANGES ARE
-    !       1.  CLARIFY COMMENTS IN THE PROLOGUES
-    !       2.  RENAME XRPRNT TO XERPRN
-    !       3.  REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES
-    !           SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE /
-    !           CHARACTER FOR NEW RECORDS.
-    !   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
-    !       CLEAN UP THE CODING.
-    !   890721  REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN
-    !       PREFIX.
-    !   891013  REVISED TO CORRECT COMMENTS.
-    !   891214  Prologue converted to Version 4.0 format.  (WRB)
-    !   900510  Changed test on NERR to be -9999999 < NERR < 99999999, but
-    !       NERR .NE. 0, and on LEVEL to be -2 < LEVEL < 3.  Added
-    !       LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and
-    !       XERCTL to XERCNT.  (RWC)
-    !   920501  Reformatted the REFERENCES section.  (WRB)
-    !***END PROLOGUE  XERMSG
-    CHARACTER(len = *) :: LIBRAR, SUBROU, MESSG
-    CHARACTER(len = 8) :: XLIBR, XSUBR
-    CHARACTER(len = 72) :: TEMP
-    CHARACTER(len = 20) :: LFIRST
-    INTEGER :: NERR, LEVEL, LKNTRL
-    INTEGER :: MAXMES, KDUMMY, I, KOUNT, LERR, LLEVEL
-    INTEGER :: MKNTRL, LTEMP
-    !***FIRST EXECUTABLE STATEMENT  XERMSG
-    LKNTRL = J4SAVE (2, 0, .FALSE.)
-    MAXMES = J4SAVE (4, 0, .FALSE.)
-
-    !   LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL.
-    !   MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE
-    !      SHOULD BE PRINTED.
-
-    !   WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN
-    !      CALLING XERMSG.  THE ERROR NUMBER SHOULD BE POSITIVE,
-    !      AND THE LEVEL SHOULD BE BETWEEN 0 AND 2.
-
-    IF (NERR<-9999999 .OR. NERR>99999999 .OR. NERR==0 .OR. &
-            LEVEL<-1 .OR. LEVEL>2) THEN
-      CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' // &
-              'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ ' // &
-              'JOB ABORT DUE TO FATAL ERROR.', 72)
-      CALL XERSVE (' ', ' ', ' ', 0, 0, 0, KDUMMY)
-      CALL XERHLT (' ***XERMSG -- INVALID INPUT')
-      RETURN
-    ENDIF
-
-    !   RECORD THE MESSAGE.
-
-    I = J4SAVE (1, NERR, .TRUE.)
-    CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT)
-
-    !   HANDLE PRINT-ONCE WARNING MESSAGES.
-
-    IF (LEVEL==-1 .AND. KOUNT>1) RETURN
-
-    !   ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG.
-
-    XLIBR = LIBRAR
-    XSUBR = SUBROU
-    LFIRST = MESSG
-    LERR = NERR
-    LLEVEL = LEVEL
-    CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL)
-
-    LKNTRL = MAX(-2, MIN(2, LKNTRL))
-    MKNTRL = ABS(LKNTRL)
-
-    !   SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS
-    !   ZERO AND THE ERROR IS NOT FATAL.
-
-    IF (LEVEL<2 .AND. LKNTRL==0) GO TO 30
-    IF (LEVEL==0 .AND. KOUNT>MAXMES) GO TO 30
-    IF (LEVEL==1 .AND. KOUNT>MAXMES .AND. MKNTRL==1) GO TO 30
-    IF (LEVEL==2 .AND. KOUNT>MAX(1, MAXMES)) GO TO 30
-
-    !   ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A
-    !   MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS)
-    !   AND SENDING IT OUT VIA XERPRN.  PRINT ONLY IF CONTROL FLAG
-    !   IS NOT ZERO.
-
-    IF (LKNTRL /= 0) THEN
-      TEMP(1:21) = 'MESSAGE FROM ROUTINE '
-      I = MIN(LEN(SUBROU), 16)
-      TEMP(22:21 + I) = SUBROU(1:I)
-      TEMP(22 + I:33 + I) = ' IN LIBRARY '
-      LTEMP = 33 + I
-      I = MIN(LEN(LIBRAR), 16)
-      TEMP(LTEMP + 1:LTEMP + I) = LIBRAR (1:I)
-      TEMP(LTEMP + I + 1:LTEMP + I + 1) = '.'
-      LTEMP = LTEMP + I + 1
-      CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
-    ENDIF
-
-    !   IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE
-    !   PRINTING THE MESSAGE.  THE INTRODUCTORY LINE TELLS THE CHOICE
-    !   FROM EACH OF THE FOLLOWING THREE OPTIONS.
-    !   1.  LEVEL OF THE MESSAGE
-    !          'INFORMATIVE MESSAGE'
-    !          'POTENTIALLY RECOVERABLE ERROR'
-    !          'FATAL ERROR'
-    !   2.  WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE
-    !          'PROG CONTINUES'
-    !          'PROG ABORTED'
-    !   3.  WHETHER OR NOT A TRACEBACK WAS REQUESTED.  (THE TRACEBACK
-    !       MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS
-    !       WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.)
-    !          'TRACEBACK REQUESTED'
-    !          'TRACEBACK NOT REQUESTED'
-    !   NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT
-    !   EXCEED 74 CHARACTERS.
-    !   WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED.
-
-    IF (LKNTRL > 0) THEN
-
-      !   THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
-
-      IF (LEVEL <= 0) THEN
-        TEMP(1:20) = 'INFORMATIVE MESSAGE,'
-        LTEMP = 20
-      ELSEIF (LEVEL == 1) THEN
-        TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,'
-        LTEMP = 30
-      ELSE
-        TEMP(1:12) = 'FATAL ERROR,'
-        LTEMP = 12
-      ENDIF
-
-      !   THEN WHETHER THE PROGRAM WILL CONTINUE.
-
-      IF ((MKNTRL==2 .AND. LEVEL>=1) .OR. &
-              (MKNTRL==1 .AND. LEVEL==2)) THEN
-        TEMP(LTEMP + 1:LTEMP + 14) = ' PROG ABORTED,'
-        LTEMP = LTEMP + 14
-      ELSE
-        TEMP(LTEMP + 1:LTEMP + 16) = ' PROG CONTINUES,'
-        LTEMP = LTEMP + 16
-      ENDIF
-
-      !   FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK.
-
-      IF (LKNTRL > 0) THEN
-        TEMP(LTEMP + 1:LTEMP + 20) = ' TRACEBACK REQUESTED'
-        LTEMP = LTEMP + 20
-      ELSE
-        TEMP(LTEMP + 1:LTEMP + 24) = ' TRACEBACK NOT REQUESTED'
-        LTEMP = LTEMP + 24
-      ENDIF
-      CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
-    ENDIF
-
-    !   NOW SEND OUT THE MESSAGE.
-
-    CALL XERPRN (' *  ', -1, MESSG, 72)
-
-    !   IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A
-    !      TRACEBACK.
-
-    IF (LKNTRL > 0) THEN
-      WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
-      DO I = 16, 22
-        IF (TEMP(I:I) /= ' ') GO TO 20
-      END DO
-
-      20   CALL XERPRN (' *  ', -1, TEMP(1:15) // TEMP(I:23), 72)
-    ENDIF
-
-    !   IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE.
-
-    IF (LKNTRL /= 0) THEN
-      CALL XERPRN (' *  ', -1, ' ', 72)
-      CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72)
-      CALL XERPRN ('    ', 0, ' ', 72)
-    ENDIF
-
-    !   IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE
-    !   CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN.
-
-    30   IF (LEVEL<=0 .OR. (LEVEL==1 .AND. MKNTRL<=1)) RETURN
-
-    !   THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A
-    !   FATAL ERROR.  PRINT THE REASON FOR THE ABORT AND THE ERROR
-    !   SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT.
-
-    IF (LKNTRL>0 .AND. KOUNT<MAX(1, MAXMES)) THEN
-      IF (LEVEL == 1) THEN
-        CALL XERPRN &
-                (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72)
-      ELSE
-        CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72)
-      ENDIF
-      CALL XERSVE (' ', ' ', ' ', -1, 0, 0, KDUMMY)
-      CALL XERHLT (' ')
-    ELSE
-      CALL XERHLT (MESSG)
-    ENDIF
-
-  END SUBROUTINE XERMSG
-
-  !DECK I1MACH
-  INTEGER FUNCTION I1MACH(I)
-    IMPLICIT NONE
-    !***BEGIN PROLOGUE  I1MACH
-    !***PURPOSE  Return integer machine dependent constants.
-    !***LIBRARY   SLATEC
-    !***CATEGORY  R1
-    !***TYPE      INTEGER (I1MACH-I)
-    !***KEYWORDS  MACHINE CONSTANTS
-    !***AUTHOR  Fox, P. A., (Bell Labs)
-    ! Hall, A. D., (Bell Labs)
-    ! Schryer, N. L., (Bell Labs)
-    !***DESCRIPTION
-
-    !   I1MACH can be used to obtain machine-dependent parameters for the
-    !   local machine environment.  It is a function subprogram with one
-    !   (input) argument and can be referenced as follows:
-
-    !    K = I1MACH(I)
-
-    !   where I=1,...,16.  The (output) value of K above is determined by
-    !   the (input) value of I.  The results for various values of I are
-    !   discussed below.
-
-    !   I/O unit numbers:
-    ! I1MACH( 1) = the standard input unit.
-    ! I1MACH( 2) = the standard output unit.
-    ! I1MACH( 3) = the standard punch unit.
-    ! I1MACH( 4) = the standard error message unit.
-
-    !   Words:
-    ! I1MACH( 5) = the number of bits per integer storage unit.
-    ! I1MACH( 6) = the number of characters per integer storage unit.
-
-    !   Integers:
-    ! assume integers are represented in the S-digit, base-A form
-
-    !            sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
-
-    !            where 0 .LE. X(I) .LT. A for I=0,...,S-1.
-    ! I1MACH( 7) = A, the base.
-    ! I1MACH( 8) = S, the number of base-A digits.
-    ! I1MACH( 9) = A**S - 1, the largest magnitude.
-
-    !   Floating-Point Numbers:
-    ! Assume floating-point numbers are represented in the T-digit,
-    ! base-B form
-    !            sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
-
-    !            where 0 .LE. X(I) .LT. B for I=1,...,T,
-    !            0 .LT. X(1), and EMIN .LE. E .LE. EMAX.
-    ! I1MACH(10) = B, the base.
-
-    !   Single-Precision:
-    ! I1MACH(11) = T, the number of base-B digits.
-    ! I1MACH(12) = EMIN, the smallest exponent E.
-    ! I1MACH(13) = EMAX, the largest exponent E.
-
-    !   Double-Precision:
-    ! I1MACH(14) = T, the number of base-B digits.
-    ! I1MACH(15) = EMIN, the smallest exponent E.
-    ! I1MACH(16) = EMAX, the largest exponent E.
-
-    !   To alter this function for a particular environment, the desired
-    !   set of DATA statements should be activated by removing the C from
-    !   column 1.  Also, the values of I1MACH(1) - I1MACH(4) should be
-    !   checked for consistency with the local operating system.
-
-    !***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
-    !             a portable library, ACM Transactions on Mathematical
-    !             Software 4, 2 (June 1978), pp. 177-188.
-    !***ROUTINES CALLED  (NONE)
-    !***REVISION HISTORY  (YYMMDD)
-    !   750101  DATE WRITTEN
-    !   891012  Added VAX G-floating constants.  (WRB)
-    !   891012  REVISION DATE from Version 3.2
-    !   891214  Prologue converted to Version 4.0 format.  (BAB)
-    !   900618  Added DEC RISC constants.  (WRB)
-    !   900723  Added IBM RS 6000 constants.  (WRB)
-    !   901009  Correct I1MACH(7) for IBM Mainframes. Should be 2 not 16.
-    !       (RWC)
-    !   910710  Added HP 730 constants.  (SMR)
-    !   911114  Added Convex IEEE constants.  (WRB)
-    !   920121  Added SUN -r8 compiler option constants.  (WRB)
-    !   920229  Added Touchstone Delta i860 constants.  (WRB)
-    !   920501  Reformatted the REFERENCES section.  (WRB)
-    !   920625  Added Convex -p8 and -pd8 compiler option constants.
-    !       (BKS, WRB)
-    !   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)
-    !   930618  Corrected I1MACH(5) for Convex -p8 and -pd8 compiler
-    !       options.  (DWL, RWC and WRB).
-    !   100623  Use Fortran 95 intrinsic functions (Lionel GUEZ)
-    !***END PROLOGUE  I1MACH
-
-    INTEGER :: IMACH(16), OUTPUT
-    SAVE IMACH
-    EQUIVALENCE (IMACH(4), OUTPUT)
-    INTEGER :: I
-    !***FIRST EXECUTABLE STATEMENT  I1MACH
-    IMACH(1) = 5
-    IMACH(2) = 6
-    IMACH(3) = 6
-    IMACH(4) = 6
-    IMACH(5) = bit_size(0)
-    IMACH(6) = IMACH(5) / 8
-    IMACH(7) = radix(0)
-    IMACH(8) = digits(0)
-    IMACH(9) = huge(0)
-    IMACH(10) = radix(0.)
-    IMACH(11) = digits(0.)
-    IMACH(12) = minexponent(0.)
-    IMACH(13) = maxexponent(0.)
-    IMACH(14) = digits(0d0)
-    IMACH(15) = minexponent(0d0)
-    IMACH(16) = maxexponent(0d0)
-    IF (I < 1  .OR.  I > 16) GO TO 10
-
-    I1MACH = IMACH(I)
-    RETURN
-
-    10   CONTINUE
-    WRITE (UNIT = OUTPUT, FMT = 9000)
-    9000   FORMAT ('1ERROR    1 IN I1MACH - I OUT OF BOUNDS')
-
-    ! CALL FDUMP
-
-    STOP
-  END FUNCTION I1MACH
-
-  !DECK XGETUA
-  SUBROUTINE XGETUA(IUNITA, N)
-    IMPLICIT NONE
-    !***BEGIN PROLOGUE  XGETUA
-    !***PURPOSE  Return unit number(s) to which error messages are being
-    ! sent.
-    !***LIBRARY   SLATEC (XERROR)
-    !***CATEGORY  R3C
-    !***TYPE      ALL (XGETUA-A)
-    !***KEYWORDS  ERROR, XERROR
-    !***AUTHOR  Jones, R. E., (SNLA)
-    !***DESCRIPTION
-
-    ! Abstract
-    !    XGETUA may be called to determine the unit number or numbers
-    !    to which error messages are being sent.
-    !    These unit numbers may have been set by a CALL to XSETUN,
-    !    or a CALL to XSETUA, or may be a default value.
-
-    ! Description of Parameters
-    !  --Output--
-    !    IUNIT - an array of one to five unit numbers, depending
-    !            on the value of N.  A value of zero refers to the
-    !            default unit, as defined by the I1MACH machine
-    !            constant routine.  Only IUNIT(1),...,IUNIT(N) are
-    !            defined by XGETUA.  The values of IUNIT(N+1),...,
-    !            IUNIT(5) are not defined (for N .LT. 5) or altered
-    !            in any way by XGETUA.
-    !    N     - the number of units to which copies of the
-    !            error messages are being sent.  N will be in the
-    !            range from 1 to 5.
-
-    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
-    !             Error-handling Package, SAND82-0800, Sandia
-    !             Laboratories, 1982.
-    !***ROUTINES CALLED  J4SAVE
-    !***REVISION HISTORY  (YYMMDD)
-    !   790801  DATE WRITTEN
-    !   861211  REVISION DATE from Version 3.2
-    !   891214  Prologue converted to Version 4.0 format.  (BAB)
-    !   920501  Reformatted the REFERENCES section.  (WRB)
-    !***END PROLOGUE  XGETUA
-    DIMENSION IUNITA(5)
-    INTEGER :: IUNITA, N, INDEX, I
-    !***FIRST EXECUTABLE STATEMENT  XGETUA
-    N = J4SAVE(5, 0, .FALSE.)
-    DO I = 1, N
-      INDEX = I + 4
-      IF (I==1) INDEX = 3
-      IUNITA(I) = J4SAVE(INDEX, 0, .FALSE.)
-    END DO
-
-  END SUBROUTINE XGETUA
-
-  !DECK XERSVE
-  SUBROUTINE XERSVE(LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, &
-          ICOUNT)
-    IMPLICIT NONE
-    !***BEGIN PROLOGUE  XERSVE
-    !***SUBSIDIARY
-    !***PURPOSE  Record that an error has occurred.
-    !***LIBRARY   SLATEC (XERROR)
-    !***CATEGORY  R3
-    !***TYPE      ALL (XERSVE-A)
-    !***KEYWORDS  ERROR, XERROR
-    !***AUTHOR  Jones, R. E., (SNLA)
-    !***DESCRIPTION
-
-    ! *Usage:
-
-    !    INTEGER  KFLAG, NERR, LEVEL, ICOUNT
-    !    CHARACTER * (len) LIBRAR, SUBROU, MESSG
-
-    !    CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT)
-
-    ! *Arguments:
-
-    !    LIBRAR :IN    is the library that the message is from.
-    !    SUBROU :IN    is the SUBROUTINE that the message is from.
-    !    MESSG  :IN    is the message to be saved.
-    !    KFLAG  :IN    indicates the action to be performed.
-    !                  when KFLAG > 0, the message in MESSG is saved.
-    !                  when KFLAG=0 the tables will be dumped and
-    !                  cleared.
-    !                  when KFLAG < 0, the tables will be dumped and
-    !                  not cleared.
-    !    NERR   :IN    is the error number.
-    !    LEVEL  :IN    is the error severity.
-    !    ICOUNT :OUT   the number of times this message has been seen,
-    !                  or zero if the table has overflowed and does not
-    !                  contain this message specifically.  When KFLAG=0,
-    !                  ICOUNT will not be altered.
-
-    ! *Description:
-
-    !   Record that this error occurred and possibly dump and clear the
-    !   tables.
-
-    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
-    !             Error-handling Package, SAND82-0800, Sandia
-    !             Laboratories, 1982.
-    !***ROUTINES CALLED  I1MACH, XGETUA
-    !***REVISION HISTORY  (YYMMDD)
-    !   800319  DATE WRITTEN
-    !   861211  REVISION DATE from Version 3.2
-    !   891214  Prologue converted to Version 4.0 format.  (BAB)
-    !   900413  Routine modified to remove reference to KFLAG.  (WRB)
-    !   900510  Changed to add LIBRARY NAME and SUBROUTINE to calling
-    !       sequence, use IF-THEN-ELSE, make number of saved entries
-    !       easily changeable, changed routine name from XERSAV to
-    !       XERSVE.  (RWC)
-    !   910626  Added LIBTAB and SUBTAB to SAVE statement.  (BKS)
-    !   920501  Reformatted the REFERENCES section.  (WRB)
-    !***END PROLOGUE  XERSVE
-    INTEGER, PARAMETER :: LENTAB = 10
-    INTEGER :: LUN(5)
-    CHARACTER(len = *) :: LIBRAR, SUBROU, MESSG
-    CHARACTER(len = 8) :: LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB
-    CHARACTER(len = 20) :: MESTAB(LENTAB), MES
-    DIMENSION NERTAB(LENTAB), LEVTAB(LENTAB), KOUNT(LENTAB)
-    SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
-    DATA KOUNTX/0/, NMSG/0/
-    INTEGER :: NERR, LEVEL, KONTRL
-    INTEGER :: NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
-    INTEGER :: KFLAG, ICOUNT, NUNIT, KUNIT, IUNIT, I
-    !***FIRST EXECUTABLE STATEMENT  XERSVE
-
-    IF (KFLAG<=0) THEN
-
-      !    Dump the table.
-
-      IF (NMSG==0) RETURN
-
-      !    Print to each unit.
-
-      CALL XGETUA (LUN, NUNIT)
-      DO KUNIT = 1, NUNIT
-        IUNIT = LUN(KUNIT)
-        IF (IUNIT==0) IUNIT = I1MACH(4)
-
-        !       Print the table header.
-
-        WRITE (IUNIT, 9000)
-
-        !       Print body of table.
-
-        DO I = 1, NMSG
-          WRITE (IUNIT, 9010) LIBTAB(I), SUBTAB(I), MESTAB(I), &
-                  NERTAB(I), LEVTAB(I), KOUNT(I)
-        END DO
-
-        !       Print number of other errors.
-
-        IF (KOUNTX/=0) WRITE (IUNIT, 9020) KOUNTX
-        WRITE (IUNIT, 9030)
-      END DO
-
-      !    Clear the error tables.
-
-      IF (KFLAG==0) THEN
-        NMSG = 0
-        KOUNTX = 0
-      ENDIF
-    ELSE
-
-      !    PROCESS A MESSAGE...
-      !    SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,
-      !    OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.
-
-      LIB = LIBRAR
-      SUB = SUBROU
-      MES = MESSG
-      DO I = 1, NMSG
-        IF (LIB==LIBTAB(I) .AND. SUB==SUBTAB(I) .AND. &
-                MES==MESTAB(I) .AND. NERR==NERTAB(I) .AND. &
-                LEVEL==LEVTAB(I)) THEN
-          KOUNT(I) = KOUNT(I) + 1
-          ICOUNT = KOUNT(I)
-          RETURN
-        ENDIF
-      END DO
-
-      IF (NMSG<LENTAB) THEN
-
-        !       Empty slot found for new message.
-
-        NMSG = NMSG + 1
-        LIBTAB(I) = LIB
-        SUBTAB(I) = SUB
-        MESTAB(I) = MES
-        NERTAB(I) = NERR
-        LEVTAB(I) = LEVEL
-        KOUNT (I) = 1
-        ICOUNT = 1
-      ELSE
-
-        !       Table is full.
-
-        KOUNTX = KOUNTX + 1
-        ICOUNT = 0
-      ENDIF
-    ENDIF
-    RETURN
-
-    ! Formats.
-
-    9000   FORMAT ('0          ERROR MESSAGE SUMMARY' / &
-            ' LIBRARY    SUBROUTINE MESSAGE START             NERR', &
-            '     LEVEL     COUNT')
-    9010   FORMAT (1X, A, 3X, A, 3X, A, 3I10)
-    9020   FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED = ', I10)
-    9030   FORMAT (1X)
-  END SUBROUTINE XERSVE
-
-  !DECK XERPRN
-  SUBROUTINE XERPRN(PREFIX, NPREF, MESSG, NWRAP)
-    IMPLICIT NONE
-    !***BEGIN PROLOGUE  XERPRN
-    !***SUBSIDIARY
-    !***PURPOSE  Print error messages processed by XERMSG.
-    !***LIBRARY   SLATEC (XERROR)
-    !***CATEGORY  R3C
-    !***TYPE      ALL (XERPRN-A)
-    !***KEYWORDS  ERROR MESSAGES, PRINTING, XERROR
-    !***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
-    !***DESCRIPTION
-
-    ! This routine sends one or more lines to each of the (up to five)
-    ! LOGICAL units to which error messages are to be sent.  This routine
-    ! is called several times by XERMSG, sometimes with a single line to
-    ! print and sometimes with a (potentially very long) message that may
-    ! wrap around into multiple lines.
-
-    ! PREFIX  Input argument of type CHARACTER.  This argument contains
-    !     characters to be put at the beginning of each line before
-    !     the body of the message.  No more than 16 characters of
-    !     PREFIX will be used.
-
-    ! NPREF   Input argument of type INTEGER.  This argument is the number
-    !     of characters to use from PREFIX.  If it is negative, the
-    !     intrinsic function LEN is used to determine its length.  If
-    !         it is zero, PREFIX is not used.  If it exceeds 16 or if
-    !     LEN(PREFIX) exceeds 16, ONLY the first 16 characters will be
-    !     used.  If NPREF is positive and the length of PREFIX is less
-    !     than NPREF, a copy of PREFIX extended with blanks to length
-    !     NPREF will be used.
-
-    ! MESSG   Input argument of type CHARACTER.  This is the text of a
-    !     message to be printed.  If it is a long message, it will be
-    !     broken into pieces for printing on multiple lines.  Each line
-    !     will start with the appropriate prefix and be followed by a
-    !     piece of the message.  NWRAP is the number of characters per
-    !     piece; that is, after each NWRAP characters, we break and
-    !     start a new line.  In addition the characters '$$' embedded
-    !     in MESSG are a sentinel for a new line.  The counting of
-    !     characters up to NWRAP starts over for each new line.  The
-    !     value of NWRAP typically used by XERMSG is 72 since many
-    !     older error messages in the SLATEC Library are laid out to
-    !     rely on wrap-around every 72 characters.
-
-    ! NWRAP   Input argument of type INTEGER.  This gives the maximum size
-    !     piece into which to break MESSG for printing on multiple
-    !     lines.  An embedded '$$' ends a line, and the count restarts
-    !     at the following character.  If a line break does not occur
-    !     on a blank (it would split a word) that word is moved to the
-    !     next line.  Values of NWRAP less than 16 will be treated as
-    !     16.  Values of NWRAP greater than 132 will be treated as 132.
-    !     The actual line length will be NPREF + NWRAP after NPREF has
-    !     been adjusted to fall between 0 and 16 and NWRAP has been
-    !     adjusted to fall between 16 and 132.
-
-    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
-    !             Error-handling Package, SAND82-0800, Sandia
-    !             Laboratories, 1982.
-    !***ROUTINES CALLED  I1MACH, XGETUA
-    !***REVISION HISTORY  (YYMMDD)
-    !   880621  DATE WRITTEN
-    !   880708  REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF
-    !       JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK
-    !       THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE
-    !       SLASH CHARACTER IN FORMAT STATEMENTS.
-    !   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
-    !       STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK
-    !       LINES TO BE PRINTED.
-    !   890721  REVISED TO ADD A NEW FEATURE.  A NEGATIVE VALUE OF NPREF
-    !       CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH.
-    !   891013  REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH.
-    !   891214  Prologue converted to Version 4.0 format.  (WRB)
-    !   900510  Added code to break messages between words.  (RWC)
-    !   920501  Reformatted the REFERENCES section.  (WRB)
-    !***END PROLOGUE  XERPRN
-    CHARACTER(len = *) :: PREFIX, MESSG
-    INTEGER :: NPREF, NWRAP
-    CHARACTER(len = 148) :: CBUFF
-    INTEGER :: IU(5), NUNIT
-    CHARACTER(len = 2) :: NEWLIN
-    PARAMETER (NEWLIN = '$$')
-    INTEGER :: N, I, LPREF, LWRAP, LENMSG, NEXTC
-    INTEGER :: LPIECE, IDELTA
-    !***FIRST EXECUTABLE STATEMENT  XERPRN
-    CALL XGETUA(IU, NUNIT)
-
-    !   A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD
-    !   ERROR MESSAGE UNIT INSTEAD.  I1MACH(4) RETRIEVES THE STANDARD
-    !   ERROR MESSAGE UNIT.
-
-    N = I1MACH(4)
-    DO I = 1, NUNIT
-      IF (IU(I) == 0) IU(I) = N
-    END DO
-
-    !   LPREF IS THE LENGTH OF THE PREFIX.  THE PREFIX IS PLACED AT THE
-    !   BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING
-    !   THE REST OF THIS ROUTINE.
-
-    IF (NPREF < 0) THEN
-      LPREF = LEN(PREFIX)
-    ELSE
-      LPREF = NPREF
-    ENDIF
-    LPREF = MIN(16, LPREF)
-    IF (LPREF /= 0) CBUFF(1:LPREF) = PREFIX
-
-    !   LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE
-    !   TIME FROM MESSG TO PRINT ON ONE LINE.
-
-    LWRAP = MAX(16, MIN(132, NWRAP))
-
-    !   SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS.
-
-    LENMSG = LEN(MESSG)
-    N = LENMSG
-    DO I = 1, N
-      IF (MESSG(LENMSG:LENMSG) /= ' ') GO TO 30
-      LENMSG = LENMSG - 1
-    END DO
-    30   CONTINUE
-
-    !   IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
-
-    IF (LENMSG == 0) THEN
-      CBUFF(LPREF + 1:LPREF + 1) = ' '
-      DO I = 1, NUNIT
-        WRITE(IU(I), '(A)') CBUFF(1:LPREF + 1)
-      END DO
-      RETURN
-    ENDIF
-
-    !   SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING
-    !   STARTS.  FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL.
-    !   WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT.
-    !   WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED.
-
-    !   WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL.  THE
-    !   INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE
-    !   OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH
-    !   OF THE SECOND ARGUMENT.
-
-    !   THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE
-    !   FOLLOWING ORDER.  WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER
-    !   OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT
-    !   POSITION NEXTC.
-
-    !   LPIECE .EQ. 0   THE NEW LINE SENTINEL DOES NOT OCCUR IN THE
-    !                   REMAINDER OF THE CHARACTER STRING.  LPIECE
-    !                   SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC,
-    !                   WHICHEVER IS LESS.
-
-    !   LPIECE .EQ. 1   THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC:
-    !                   NEXTC).  LPIECE IS EFFECTIVELY ZERO, AND WE
-    !                   PRINT NOTHING TO AVOID PRODUCING UNNECESSARY
-    !                   BLANK LINES.  THIS TAKES CARE OF THE SITUATION
-    !                   WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF
-    !                   EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE
-    !                   SENTINEL FOLLOWED BY MORE CHARACTERS.  NEXTC
-    !                   SHOULD BE INCREMENTED BY 2.
-
-    !   LPIECE .GT. LWRAP+1  REDUCE LPIECE TO LWRAP.
-
-    !   ELSE            THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1
-    !                   RESET LPIECE = LPIECE-1.  NOTE THAT THIS
-    !                   PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ.
-    !                   LWRAP+1.  THAT IS, THE SENTINEL FALLS EXACTLY
-    !                   AT THE END OF A LINE.
-
-    NEXTC = 1
-    50   LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN)
-    IF (LPIECE == 0) THEN
-
-      !   THERE WAS NO NEW LINE SENTINEL FOUND.
-
-      IDELTA = 0
-      LPIECE = MIN(LWRAP, LENMSG + 1 - NEXTC)
-      IF (LPIECE < LENMSG + 1 - NEXTC) THEN
-        DO I = LPIECE + 1, 2, -1
-          IF (MESSG(NEXTC + I - 1:NEXTC + I - 1) == ' ') THEN
-            LPIECE = I - 1
-            IDELTA = 1
-            GOTO 54
-          ENDIF
-        END DO
-      ENDIF
-      54   CBUFF(LPREF + 1:LPREF + LPIECE) = MESSG(NEXTC:NEXTC + LPIECE - 1)
-      NEXTC = NEXTC + LPIECE + IDELTA
-    ELSEIF (LPIECE == 1) THEN
-
-      !   WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1).
-      !   DON'T PRINT A BLANK LINE.
-
-      NEXTC = NEXTC + 2
-      GO TO 50
-    ELSEIF (LPIECE > LWRAP + 1) THEN
-
-      !   LPIECE SHOULD BE SET DOWN TO LWRAP.
-
-      IDELTA = 0
-      LPIECE = LWRAP
-      DO I = LPIECE + 1, 2, -1
-        IF (MESSG(NEXTC + I - 1:NEXTC + I - 1) == ' ') THEN
-          LPIECE = I - 1
-          IDELTA = 1
-          GOTO 58
-        ENDIF
-      END DO
-      58   CBUFF(LPREF + 1:LPREF + LPIECE) = MESSG(NEXTC:NEXTC + LPIECE - 1)
-      NEXTC = NEXTC + LPIECE + IDELTA
-    ELSE
-
-      !   IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1.
-      !   WE SHOULD DECREMENT LPIECE BY ONE.
-
-      LPIECE = LPIECE - 1
-      CBUFF(LPREF + 1:LPREF + LPIECE) = MESSG(NEXTC:NEXTC + LPIECE - 1)
-      NEXTC = NEXTC + LPIECE + 2
-    ENDIF
-
-    !   PRINT
-
-    DO I = 1, NUNIT
-      WRITE(IU(I), '(A)') CBUFF(1:LPREF + LPIECE)
-    END DO
-
-    IF (NEXTC <= LENMSG) GO TO 50
-
-  END SUBROUTINE XERPRN
-
-  !DECK XERHLT
-  SUBROUTINE XERHLT(MESSG)
-    !***BEGIN PROLOGUE  XERHLT
-    !***SUBSIDIARY
-    !***PURPOSE  Abort program execution and print error message.
-    !***LIBRARY   SLATEC (XERROR)
-    !***CATEGORY  R3C
-    !***TYPE      ALL (XERHLT-A)
-    !***KEYWORDS  ABORT PROGRAM EXECUTION, ERROR, XERROR
-    !***AUTHOR  Jones, R. E., (SNLA)
-    !***DESCRIPTION
-
-    ! Abstract
-    !    ***Note*** machine dependent routine
-    !    XERHLT aborts the execution of the program.
-    !    The error message causing the abort is given in the calling
-    !    sequence, in case one needs it for printing on a dayfile,
-    !    for example.
-
-    ! Description of Parameters
-    !    MESSG is as in XERMSG.
-
-    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
-    !             Error-handling Package, SAND82-0800, Sandia
-    !             Laboratories, 1982.
-    !***ROUTINES CALLED  (NONE)
-    !***REVISION HISTORY  (YYMMDD)
-    !   790801  DATE WRITTEN
-    !   861211  REVISION DATE from Version 3.2
-    !   891214  Prologue converted to Version 4.0 format.  (BAB)
-    !   900206  Routine changed from user-callable to subsidiary.  (WRB)
-    !   900510  Changed calling sequence to delete length of character
-    !       and changed routine name from XERABT to XERHLT.  (RWC)
-    !   920501  Reformatted the REFERENCES section.  (WRB)
-    !***END PROLOGUE  XERHLT
-    CHARACTER(len = *) :: MESSG
-    !***FIRST EXECUTABLE STATEMENT  XERHLT
-    STOP
-  END SUBROUTINE XERHLT
-
-  !DECK XERCNT
-  SUBROUTINE XERCNT(LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
-    IMPLICIT NONE
-    !***BEGIN PROLOGUE  XERCNT
-    !***SUBSIDIARY
-    !***PURPOSE  Allow user control over handling of errors.
-    !***LIBRARY   SLATEC (XERROR)
-    !***CATEGORY  R3C
-    !***TYPE      ALL (XERCNT-A)
-    !***KEYWORDS  ERROR, XERROR
-    !***AUTHOR  Jones, R. E., (SNLA)
-    !***DESCRIPTION
-
-    ! Abstract
-    !    Allows user control over handling of individual errors.
-    !    Just after each message is recorded, but before it is
-    !    processed any further (i.e., before it is printed or
-    !    a decision to abort is made), a CALL is made to XERCNT.
-    !    If the user has provided his own version of XERCNT, he
-    !    can then override the value of KONTROL used in processing
-    !    this message by redefining its value.
-    !    KONTRL may be set to any value from -2 to 2.
-    !    The meanings for KONTRL are the same as in XSETF, except
-    !    that the value of KONTRL changes only for this message.
-    !    If KONTRL is set to a value outside the range from -2 to 2,
-    !    it will be moved back into that range.
-
-    ! Description of Parameters
-
-    !  --Input--
-    !    LIBRAR - the library that the routine is in.
-    !    SUBROU - the SUBROUTINE that XERMSG is being called from
-    !    MESSG  - the first 20 characters of the error message.
-    !    NERR   - same as in the CALL to XERMSG.
-    !    LEVEL  - same as in the CALL to XERMSG.
-    !    KONTRL - the current value of the control flag as set
-    !             by a CALL to XSETF.
-
-    !  --Output--
-    !    KONTRL - the new value of KONTRL.  If KONTRL is not
-    !             defined, it will remain at its original value.
-    !             This changed value of control affects only
-    !             the current occurrence of the current message.
-
-    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
-    !             Error-handling Package, SAND82-0800, Sandia
-    !             Laboratories, 1982.
-    !***ROUTINES CALLED  (NONE)
-    !***REVISION HISTORY  (YYMMDD)
-    !   790801  DATE WRITTEN
-    !   861211  REVISION DATE from Version 3.2
-    !   891214  Prologue converted to Version 4.0 format.  (BAB)
-    !   900206  Routine changed from user-callable to subsidiary.  (WRB)
-    !   900510  Changed calling sequence to INCLUDE LIBRARY and SUBROUTINE
-    !       names, changed routine name from XERCTL to XERCNT.  (RWC)
-    !   920501  Reformatted the REFERENCES section.  (WRB)
-    !***END PROLOGUE  XERCNT
-    CHARACTER(len = *) :: LIBRAR, SUBROU, MESSG
-    INTEGER :: NERR, LEVEL, KONTRL
-    !***FIRST EXECUTABLE STATEMENT  XERCNT
-
-  END SUBROUTINE XERCNT
-
-  !DECK J4SAVE
-  INTEGER FUNCTION J4SAVE(IWHICH, IVALUE, ISET)
-    IMPLICIT NONE
-    !***BEGIN PROLOGUE  J4SAVE
-    !***SUBSIDIARY
-    !***PURPOSE  Save or reCALL global variables needed by error
-    ! handling routines.
-    !***LIBRARY   SLATEC (XERROR)
-    !***TYPE      INTEGER (J4SAVE-I)
-    !***KEYWORDS  ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR
-    !***AUTHOR  Jones, R. E., (SNLA)
-    !***DESCRIPTION
-
-    ! Abstract
-    !    J4SAVE saves and recalls several global variables needed
-    !    by the library error handling routines.
-
-    ! Description of Parameters
-    !  --Input--
-    !    IWHICH - Index of item desired.
-    !            = 1 Refers to current error number.
-    !            = 2 Refers to current error control flag.
-    !            = 3 Refers to current unit number to which error
-    !                messages are to be sent.  (0 means use standard.)
-    !            = 4 Refers to the maximum number of times any
-    !                 message is to be printed (as set by XERMAX).
-    !            = 5 Refers to the total number of units to which
-    !                 each error message is to be written.
-    !            = 6 Refers to the 2nd unit for error messages
-    !            = 7 Refers to the 3rd unit for error messages
-    !            = 8 Refers to the 4th unit for error messages
-    !            = 9 Refers to the 5th unit for error messages
-    !    IVALUE - The value to be set for the IWHICH-th parameter,
-    !             if ISET is .TRUE. .
-    !    ISET   - If ISET=.TRUE., the IWHICH-th parameter will BE
-    !             given the value, IVALUE.  If ISET=.FALSE., the
-    !             IWHICH-th parameter will be unchanged, and IVALUE
-    !             is a dummy parameter.
-    !  --Output--
-    !    The (old) value of the IWHICH-th parameter will be returned
-    !    in the function value, J4SAVE.
-
-    !***SEE ALSO  XERMSG
-    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
-    !             Error-handling Package, SAND82-0800, Sandia
-    !             Laboratories, 1982.
-    !***ROUTINES CALLED  (NONE)
-    !***REVISION HISTORY  (YYMMDD)
-    !   790801  DATE WRITTEN
-    !   891214  Prologue converted to Version 4.0 format.  (BAB)
-    !   900205  Minor modifications to prologue.  (WRB)
-    !   900402  Added TYPE section.  (WRB)
-    !   910411  Added KEYWORDS section.  (WRB)
-    !   920501  Reformatted the REFERENCES section.  (WRB)
-    !***END PROLOGUE  J4SAVE
-    LOGICAL :: ISET
-    INTEGER :: IPARAM(9)
-    SAVE IPARAM
-    DATA IPARAM(1), IPARAM(2), IPARAM(3), IPARAM(4)/0, 2, 0, 10/
-    DATA IPARAM(5)/1/
-    DATA IPARAM(6), IPARAM(7), IPARAM(8), IPARAM(9)/0, 0, 0, 0/
-    INTEGER :: IWHICH, IVALUE
-    !***FIRST EXECUTABLE STATEMENT  J4SAVE
-    J4SAVE = IPARAM(IWHICH)
-    IF (ISET) IPARAM(IWHICH) = IVALUE
-
-  END FUNCTION J4SAVE
-
-
-END MODULE lmdz_xer
Index: LMDZ6/branches/Amaury_dev/libf/misc/slatec_libmath_pch.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/misc/slatec_libmath_pch.f90	(revision 5218)
+++ LMDZ6/branches/Amaury_dev/libf/misc/slatec_libmath_pch.f90	(revision 5218)
@@ -0,0 +1,1150 @@
+! This module encapsulates misc. math functions
+
+MODULE lmdz_libmath_pch
+  USE slatec_xer, ONLY: xermsg
+  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
+
+    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
Index: LMDZ6/branches/Amaury_dev/libf/misc/slatec_xer.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/misc/slatec_xer.f90	(revision 5218)
+++ LMDZ6/branches/Amaury_dev/libf/misc/slatec_xer.f90	(revision 5218)
@@ -0,0 +1,1123 @@
+! Contains "legacy" xer* functions required by some other very old external code
+
+MODULE lmdz_xer
+  IMPLICIT NONE; PRIVATE
+  PUBLIC xermsg
+CONTAINS
+
+  !DECK XERMSG
+  SUBROUTINE XERMSG(LIBRAR, SUBROU, MESSG, NERR, LEVEL)
+    IMPLICIT NONE
+    !***BEGIN PROLOGUE  XERMSG
+    !***PURPOSE  Process error messages for SLATEC and other libraries.
+    !***LIBRARY   SLATEC (XERROR)
+    !***CATEGORY  R3C
+    !***TYPE      ALL (XERMSG-A)
+    !***KEYWORDS  ERROR MESSAGE, XERROR
+    !***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
+    !***DESCRIPTION
+
+    !   XERMSG processes a diagnostic message in a manner determined by the
+    !   value of LEVEL and the current value of the library error control
+    !   flag, KONTRL.  See SUBROUTINE XSETF for details.
+
+    !    LIBRAR   A character constant (or character variable) with the name
+    !         of the library.  This will be 'SLATEC' for the SLATEC
+    !         Common Math Library.  The error handling package is
+    !         general enough to be used by many libraries
+    !         simultaneously, so it is desirable for the routine that
+    !         detects and reports an error to identify the library name
+    !         as well as the routine name.
+
+    !    SUBROU   A character constant (or character variable) with the name
+    !             of the routine that detected the error.  Usually it is the
+    !         name of the routine that is calling XERMSG.  There are
+    !         some instances where a user callable library routine calls
+    !         lower level subsidiary routines where the error is
+    !         detected.  In such cases it may be more informative to
+    !         supply the name of the routine the user called rather than
+    !         the name of the subsidiary routine that detected the
+    !         error.
+
+    !    MESSG    A character constant (or character variable) with the text
+    !         of the error or warning message.  In the example below,
+    !         the message is a character constant that contains a
+    !         generic message.
+
+    !               CALL XERMSG ('SLATEC', 'MMPY',
+    !              *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION',
+    !              *3, 1)
+
+    !         It is possible (and is sometimes desirable) to generate a
+    !         specific message--e.g., one that contains actual numeric
+    !         values.  Specific numeric values can be converted into
+    !         character strings using formatted WRITE statements into
+    !         character variables.  This is called standard Fortran
+    !         internal file I/O and is exemplified in the first three
+    !         lines of the following example.  You can also catenate
+    !         substrings of characters to construct the error message.
+    !         Here is an example showing the use of both writing to
+    !         an internal file and catenating character strings.
+
+    !               CHARACTER*5 CHARN, CHARL
+    !               WRITE (CHARN,10) N
+    !               WRITE (CHARL,10) LDA
+    !            10 FORMAT(I5)
+    !               CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN//
+    !              *   ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'//
+    !              *   CHARL, 3, 1)
+
+    !         There are two subtleties worth mentioning.  One is that
+    !         the // for character catenation is used to construct the
+    !         error message so that no single character constant is
+    !         continued to the next line.  This avoids confusion as to
+    !         whether there are trailing blanks at the end of the line.
+    !         The second is that by catenating the parts of the message
+    !             as an actual argument rather than encoding the entire
+    !         message into one large character variable, we avoid
+    !         having to know how long the message will be in order to
+    !         declare an adequate length for that large character
+    !         variable.  XERMSG calls XERPRN to print the message using
+    !         multiple lines if necessary.  If the message is very long,
+    !         XERPRN will break it into pieces of 72 characters (as
+    !         requested by XERMSG) for printing on multiple lines.
+    !         Also, XERMSG asks XERPRN to prefix each line with ' *  '
+    !         so that the total line length could be 76 characters.
+    !         Note also that XERPRN scans the error message backwards
+    !         to ignore trailing blanks.  Another feature is that
+    !         the substring '$$' is treated as a new line sentinel
+    !         by XERPRN.  If you want to construct a multiline
+    !         message without having to count out multiples of 72
+    !         characters, just use '$$' as a separator.  '$$'
+    !         obviously must occur within 72 characters of the
+    !         start of each line to have its intended effect since
+    !         XERPRN is asked to wrap around at 72 characters in
+    !         addition to looking for '$$'.
+
+    !    NERR     An integer value that is chosen by the library routine's
+    !         author.  It must be in the range -99 to 999 (three
+    !         printable digits).  Each distinct error should have its
+    !         own error number.  These error numbers should be described
+    !         in the machine readable documentation for the routine.
+    !         The error numbers need be unique only within each routine,
+    !         so it is reasonable for each routine to start enumerating
+    !         errors from 1 and proceeding to the next integer.
+
+    !    LEVEL    An integer value in the range 0 to 2 that indicates the
+    !         level (severity) of the error.  Their meanings are
+
+    !        -1  A warning message.  This is used if it is not clear
+    !            that there really is an error, but the user's attention
+    !            may be needed.  An attempt is made to only print this
+    !            message once.
+
+    !         0  A warning message.  This is used if it is not clear
+    !            that there really is an error, but the user's attention
+    !            may be needed.
+
+    !         1  A recoverable error.  This is used even if the error is
+    !            so serious that the routine cannot return any useful
+    !            answer.  If the user has told the error package to
+    !            return after recoverable errors, then XERMSG will
+    !            return to the Library routine which can then return to
+    !            the user's routine.  The user may also permit the error
+    !            package to terminate the program upon encountering a
+    !            recoverable error.
+
+    !         2  A fatal error.  XERMSG will not return to its caller
+    !            after it receives a fatal error.  This level should
+    !            hardly ever be used; it is much better to allow the
+    !            user a chance to recover.  An example of one of the few
+    !            cases in which it is permissible to declare a level 2
+    !            error is a reverse communication Library routine that
+    !            is likely to be called repeatedly until it integrates
+    !            across some interval.  If there is a serious error in
+    !            the input such that another step cannot be taken and
+    !            the Library routine is called again without the input
+    !            error having been corrected by the caller, the Library
+    !            routine will probably be called forever with improper
+    !            input.  In this case, it is reasonable to declare the
+    !            error to be fatal.
+
+    !    Each of the arguments to XERMSG is input; none will be modified by
+    !    XERMSG.  A routine may make multiple calls to XERMSG with warning
+    !    level messages; however, after a CALL to XERMSG with a recoverable
+    !    error, the routine should return to the user.  Do not try to call
+    !    XERMSG with a second recoverable error after the first recoverable
+    !    error because the error package saves the error number.  The user
+    !    can retrieve this error number by calling another entry point in
+    !    the error handling package and then clear the error number when
+    !    recovering from the error.  Calling XERMSG in succession causes the
+    !    old error number to be overwritten by the latest error number.
+    !    This is considered harmless for error numbers associated with
+    !    warning messages but must not be done for error numbers of serious
+    !    errors.  After a CALL to XERMSG with a recoverable error, the user
+    !    must be given a chance to CALL NUMXER or XERCLR to retrieve or
+    !    clear the error number.
+    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+    !             Error-handling Package, SAND82-0800, Sandia
+    !             Laboratories, 1982.
+    !***ROUTINES CALLED  FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE
+    !***REVISION HISTORY  (YYMMDD)
+    !   880101  DATE WRITTEN
+    !   880621  REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988.
+    !       THERE ARE TWO BASIC CHANGES.
+    !       1.  A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO
+    !           PRINT MESSAGES.  THIS ROUTINE WILL BREAK LONG MESSAGES
+    !           INTO PIECES FOR PRINTING ON MULTIPLE LINES.  '$$' IS
+    !           ACCEPTED AS A NEW LINE SENTINEL.  A PREFIX CAN BE
+    !           ADDED TO EACH LINE TO BE PRINTED.  XERMSG USES EITHER
+    !           ' ***' OR ' *  ' AND LONG MESSAGES ARE BROKEN EVERY
+    !           72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE
+    !           LENGTH OUTPUT CAN NOW BE AS GREAT AS 76.
+    !       2.  THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE
+    !           FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE
+    !           OF LOWER CASE.
+    !   880708  REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30.
+    !       THE PRINCIPAL CHANGES ARE
+    !       1.  CLARIFY COMMENTS IN THE PROLOGUES
+    !       2.  RENAME XRPRNT TO XERPRN
+    !       3.  REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES
+    !           SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE /
+    !           CHARACTER FOR NEW RECORDS.
+    !   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
+    !       CLEAN UP THE CODING.
+    !   890721  REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN
+    !       PREFIX.
+    !   891013  REVISED TO CORRECT COMMENTS.
+    !   891214  Prologue converted to Version 4.0 format.  (WRB)
+    !   900510  Changed test on NERR to be -9999999 < NERR < 99999999, but
+    !       NERR .NE. 0, and on LEVEL to be -2 < LEVEL < 3.  Added
+    !       LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and
+    !       XERCTL to XERCNT.  (RWC)
+    !   920501  Reformatted the REFERENCES section.  (WRB)
+    !***END PROLOGUE  XERMSG
+    CHARACTER(len = *) :: LIBRAR, SUBROU, MESSG
+    CHARACTER(len = 8) :: XLIBR, XSUBR
+    CHARACTER(len = 72) :: TEMP
+    CHARACTER(len = 20) :: LFIRST
+    INTEGER :: NERR, LEVEL, LKNTRL
+    INTEGER :: MAXMES, KDUMMY, I, KOUNT, LERR, LLEVEL
+    INTEGER :: MKNTRL, LTEMP
+    !***FIRST EXECUTABLE STATEMENT  XERMSG
+    LKNTRL = J4SAVE (2, 0, .FALSE.)
+    MAXMES = J4SAVE (4, 0, .FALSE.)
+
+    !   LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL.
+    !   MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE
+    !      SHOULD BE PRINTED.
+
+    !   WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN
+    !      CALLING XERMSG.  THE ERROR NUMBER SHOULD BE POSITIVE,
+    !      AND THE LEVEL SHOULD BE BETWEEN 0 AND 2.
+
+    IF (NERR<-9999999 .OR. NERR>99999999 .OR. NERR==0 .OR. &
+            LEVEL<-1 .OR. LEVEL>2) THEN
+      CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' // &
+              'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ ' // &
+              'JOB ABORT DUE TO FATAL ERROR.', 72)
+      CALL XERSVE (' ', ' ', ' ', 0, 0, 0, KDUMMY)
+      CALL XERHLT (' ***XERMSG -- INVALID INPUT')
+      RETURN
+    ENDIF
+
+    !   RECORD THE MESSAGE.
+
+    I = J4SAVE (1, NERR, .TRUE.)
+    CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT)
+
+    !   HANDLE PRINT-ONCE WARNING MESSAGES.
+
+    IF (LEVEL==-1 .AND. KOUNT>1) RETURN
+
+    !   ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG.
+
+    XLIBR = LIBRAR
+    XSUBR = SUBROU
+    LFIRST = MESSG
+    LERR = NERR
+    LLEVEL = LEVEL
+    CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL)
+
+    LKNTRL = MAX(-2, MIN(2, LKNTRL))
+    MKNTRL = ABS(LKNTRL)
+
+    !   SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS
+    !   ZERO AND THE ERROR IS NOT FATAL.
+
+    IF (LEVEL<2 .AND. LKNTRL==0) GO TO 30
+    IF (LEVEL==0 .AND. KOUNT>MAXMES) GO TO 30
+    IF (LEVEL==1 .AND. KOUNT>MAXMES .AND. MKNTRL==1) GO TO 30
+    IF (LEVEL==2 .AND. KOUNT>MAX(1, MAXMES)) GO TO 30
+
+    !   ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A
+    !   MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS)
+    !   AND SENDING IT OUT VIA XERPRN.  PRINT ONLY IF CONTROL FLAG
+    !   IS NOT ZERO.
+
+    IF (LKNTRL /= 0) THEN
+      TEMP(1:21) = 'MESSAGE FROM ROUTINE '
+      I = MIN(LEN(SUBROU), 16)
+      TEMP(22:21 + I) = SUBROU(1:I)
+      TEMP(22 + I:33 + I) = ' IN LIBRARY '
+      LTEMP = 33 + I
+      I = MIN(LEN(LIBRAR), 16)
+      TEMP(LTEMP + 1:LTEMP + I) = LIBRAR (1:I)
+      TEMP(LTEMP + I + 1:LTEMP + I + 1) = '.'
+      LTEMP = LTEMP + I + 1
+      CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
+    ENDIF
+
+    !   IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE
+    !   PRINTING THE MESSAGE.  THE INTRODUCTORY LINE TELLS THE CHOICE
+    !   FROM EACH OF THE FOLLOWING THREE OPTIONS.
+    !   1.  LEVEL OF THE MESSAGE
+    !          'INFORMATIVE MESSAGE'
+    !          'POTENTIALLY RECOVERABLE ERROR'
+    !          'FATAL ERROR'
+    !   2.  WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE
+    !          'PROG CONTINUES'
+    !          'PROG ABORTED'
+    !   3.  WHETHER OR NOT A TRACEBACK WAS REQUESTED.  (THE TRACEBACK
+    !       MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS
+    !       WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.)
+    !          'TRACEBACK REQUESTED'
+    !          'TRACEBACK NOT REQUESTED'
+    !   NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT
+    !   EXCEED 74 CHARACTERS.
+    !   WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED.
+
+    IF (LKNTRL > 0) THEN
+
+      !   THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
+
+      IF (LEVEL <= 0) THEN
+        TEMP(1:20) = 'INFORMATIVE MESSAGE,'
+        LTEMP = 20
+      ELSEIF (LEVEL == 1) THEN
+        TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,'
+        LTEMP = 30
+      ELSE
+        TEMP(1:12) = 'FATAL ERROR,'
+        LTEMP = 12
+      ENDIF
+
+      !   THEN WHETHER THE PROGRAM WILL CONTINUE.
+
+      IF ((MKNTRL==2 .AND. LEVEL>=1) .OR. &
+              (MKNTRL==1 .AND. LEVEL==2)) THEN
+        TEMP(LTEMP + 1:LTEMP + 14) = ' PROG ABORTED,'
+        LTEMP = LTEMP + 14
+      ELSE
+        TEMP(LTEMP + 1:LTEMP + 16) = ' PROG CONTINUES,'
+        LTEMP = LTEMP + 16
+      ENDIF
+
+      !   FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK.
+
+      IF (LKNTRL > 0) THEN
+        TEMP(LTEMP + 1:LTEMP + 20) = ' TRACEBACK REQUESTED'
+        LTEMP = LTEMP + 20
+      ELSE
+        TEMP(LTEMP + 1:LTEMP + 24) = ' TRACEBACK NOT REQUESTED'
+        LTEMP = LTEMP + 24
+      ENDIF
+      CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
+    ENDIF
+
+    !   NOW SEND OUT THE MESSAGE.
+
+    CALL XERPRN (' *  ', -1, MESSG, 72)
+
+    !   IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A
+    !      TRACEBACK.
+
+    IF (LKNTRL > 0) THEN
+      WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
+      DO I = 16, 22
+        IF (TEMP(I:I) /= ' ') GO TO 20
+      END DO
+
+      20   CALL XERPRN (' *  ', -1, TEMP(1:15) // TEMP(I:23), 72)
+    ENDIF
+
+    !   IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE.
+
+    IF (LKNTRL /= 0) THEN
+      CALL XERPRN (' *  ', -1, ' ', 72)
+      CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72)
+      CALL XERPRN ('    ', 0, ' ', 72)
+    ENDIF
+
+    !   IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE
+    !   CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN.
+
+    30   IF (LEVEL<=0 .OR. (LEVEL==1 .AND. MKNTRL<=1)) RETURN
+
+    !   THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A
+    !   FATAL ERROR.  PRINT THE REASON FOR THE ABORT AND THE ERROR
+    !   SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT.
+
+    IF (LKNTRL>0 .AND. KOUNT<MAX(1, MAXMES)) THEN
+      IF (LEVEL == 1) THEN
+        CALL XERPRN &
+                (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72)
+      ELSE
+        CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72)
+      ENDIF
+      CALL XERSVE (' ', ' ', ' ', -1, 0, 0, KDUMMY)
+      CALL XERHLT (' ')
+    ELSE
+      CALL XERHLT (MESSG)
+    ENDIF
+
+  END SUBROUTINE XERMSG
+
+  !DECK I1MACH
+  INTEGER FUNCTION I1MACH(I)
+    IMPLICIT NONE
+    !***BEGIN PROLOGUE  I1MACH
+    !***PURPOSE  Return integer machine dependent constants.
+    !***LIBRARY   SLATEC
+    !***CATEGORY  R1
+    !***TYPE      INTEGER (I1MACH-I)
+    !***KEYWORDS  MACHINE CONSTANTS
+    !***AUTHOR  Fox, P. A., (Bell Labs)
+    ! Hall, A. D., (Bell Labs)
+    ! Schryer, N. L., (Bell Labs)
+    !***DESCRIPTION
+
+    !   I1MACH can be used to obtain machine-dependent parameters for the
+    !   local machine environment.  It is a function subprogram with one
+    !   (input) argument and can be referenced as follows:
+
+    !    K = I1MACH(I)
+
+    !   where I=1,...,16.  The (output) value of K above is determined by
+    !   the (input) value of I.  The results for various values of I are
+    !   discussed below.
+
+    !   I/O unit numbers:
+    ! I1MACH( 1) = the standard input unit.
+    ! I1MACH( 2) = the standard output unit.
+    ! I1MACH( 3) = the standard punch unit.
+    ! I1MACH( 4) = the standard error message unit.
+
+    !   Words:
+    ! I1MACH( 5) = the number of bits per integer storage unit.
+    ! I1MACH( 6) = the number of characters per integer storage unit.
+
+    !   Integers:
+    ! assume integers are represented in the S-digit, base-A form
+
+    !            sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
+
+    !            where 0 .LE. X(I) .LT. A for I=0,...,S-1.
+    ! I1MACH( 7) = A, the base.
+    ! I1MACH( 8) = S, the number of base-A digits.
+    ! I1MACH( 9) = A**S - 1, the largest magnitude.
+
+    !   Floating-Point Numbers:
+    ! Assume floating-point numbers are represented in the T-digit,
+    ! base-B form
+    !            sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
+
+    !            where 0 .LE. X(I) .LT. B for I=1,...,T,
+    !            0 .LT. X(1), and EMIN .LE. E .LE. EMAX.
+    ! I1MACH(10) = B, the base.
+
+    !   Single-Precision:
+    ! I1MACH(11) = T, the number of base-B digits.
+    ! I1MACH(12) = EMIN, the smallest exponent E.
+    ! I1MACH(13) = EMAX, the largest exponent E.
+
+    !   Double-Precision:
+    ! I1MACH(14) = T, the number of base-B digits.
+    ! I1MACH(15) = EMIN, the smallest exponent E.
+    ! I1MACH(16) = EMAX, the largest exponent E.
+
+    !   To alter this function for a particular environment, the desired
+    !   set of DATA statements should be activated by removing the C from
+    !   column 1.  Also, the values of I1MACH(1) - I1MACH(4) should be
+    !   checked for consistency with the local operating system.
+
+    !***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
+    !             a portable library, ACM Transactions on Mathematical
+    !             Software 4, 2 (June 1978), pp. 177-188.
+    !***ROUTINES CALLED  (NONE)
+    !***REVISION HISTORY  (YYMMDD)
+    !   750101  DATE WRITTEN
+    !   891012  Added VAX G-floating constants.  (WRB)
+    !   891012  REVISION DATE from Version 3.2
+    !   891214  Prologue converted to Version 4.0 format.  (BAB)
+    !   900618  Added DEC RISC constants.  (WRB)
+    !   900723  Added IBM RS 6000 constants.  (WRB)
+    !   901009  Correct I1MACH(7) for IBM Mainframes. Should be 2 not 16.
+    !       (RWC)
+    !   910710  Added HP 730 constants.  (SMR)
+    !   911114  Added Convex IEEE constants.  (WRB)
+    !   920121  Added SUN -r8 compiler option constants.  (WRB)
+    !   920229  Added Touchstone Delta i860 constants.  (WRB)
+    !   920501  Reformatted the REFERENCES section.  (WRB)
+    !   920625  Added Convex -p8 and -pd8 compiler option constants.
+    !       (BKS, WRB)
+    !   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)
+    !   930618  Corrected I1MACH(5) for Convex -p8 and -pd8 compiler
+    !       options.  (DWL, RWC and WRB).
+    !   100623  Use Fortran 95 intrinsic functions (Lionel GUEZ)
+    !***END PROLOGUE  I1MACH
+
+    INTEGER :: IMACH(16), OUTPUT
+    SAVE IMACH
+    EQUIVALENCE (IMACH(4), OUTPUT)
+    INTEGER :: I
+    !***FIRST EXECUTABLE STATEMENT  I1MACH
+    IMACH(1) = 5
+    IMACH(2) = 6
+    IMACH(3) = 6
+    IMACH(4) = 6
+    IMACH(5) = bit_size(0)
+    IMACH(6) = IMACH(5) / 8
+    IMACH(7) = radix(0)
+    IMACH(8) = digits(0)
+    IMACH(9) = huge(0)
+    IMACH(10) = radix(0.)
+    IMACH(11) = digits(0.)
+    IMACH(12) = minexponent(0.)
+    IMACH(13) = maxexponent(0.)
+    IMACH(14) = digits(0d0)
+    IMACH(15) = minexponent(0d0)
+    IMACH(16) = maxexponent(0d0)
+    IF (I < 1  .OR.  I > 16) GO TO 10
+
+    I1MACH = IMACH(I)
+    RETURN
+
+    10   CONTINUE
+    WRITE (UNIT = OUTPUT, FMT = 9000)
+    9000   FORMAT ('1ERROR    1 IN I1MACH - I OUT OF BOUNDS')
+
+    ! CALL FDUMP
+
+    STOP
+  END FUNCTION I1MACH
+
+  !DECK XGETUA
+  SUBROUTINE XGETUA(IUNITA, N)
+    IMPLICIT NONE
+    !***BEGIN PROLOGUE  XGETUA
+    !***PURPOSE  Return unit number(s) to which error messages are being
+    ! sent.
+    !***LIBRARY   SLATEC (XERROR)
+    !***CATEGORY  R3C
+    !***TYPE      ALL (XGETUA-A)
+    !***KEYWORDS  ERROR, XERROR
+    !***AUTHOR  Jones, R. E., (SNLA)
+    !***DESCRIPTION
+
+    ! Abstract
+    !    XGETUA may be called to determine the unit number or numbers
+    !    to which error messages are being sent.
+    !    These unit numbers may have been set by a CALL to XSETUN,
+    !    or a CALL to XSETUA, or may be a default value.
+
+    ! Description of Parameters
+    !  --Output--
+    !    IUNIT - an array of one to five unit numbers, depending
+    !            on the value of N.  A value of zero refers to the
+    !            default unit, as defined by the I1MACH machine
+    !            constant routine.  Only IUNIT(1),...,IUNIT(N) are
+    !            defined by XGETUA.  The values of IUNIT(N+1),...,
+    !            IUNIT(5) are not defined (for N .LT. 5) or altered
+    !            in any way by XGETUA.
+    !    N     - the number of units to which copies of the
+    !            error messages are being sent.  N will be in the
+    !            range from 1 to 5.
+
+    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+    !             Error-handling Package, SAND82-0800, Sandia
+    !             Laboratories, 1982.
+    !***ROUTINES CALLED  J4SAVE
+    !***REVISION HISTORY  (YYMMDD)
+    !   790801  DATE WRITTEN
+    !   861211  REVISION DATE from Version 3.2
+    !   891214  Prologue converted to Version 4.0 format.  (BAB)
+    !   920501  Reformatted the REFERENCES section.  (WRB)
+    !***END PROLOGUE  XGETUA
+    DIMENSION IUNITA(5)
+    INTEGER :: IUNITA, N, INDEX, I
+    !***FIRST EXECUTABLE STATEMENT  XGETUA
+    N = J4SAVE(5, 0, .FALSE.)
+    DO I = 1, N
+      INDEX = I + 4
+      IF (I==1) INDEX = 3
+      IUNITA(I) = J4SAVE(INDEX, 0, .FALSE.)
+    END DO
+
+  END SUBROUTINE XGETUA
+
+  !DECK XERSVE
+  SUBROUTINE XERSVE(LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, &
+          ICOUNT)
+    IMPLICIT NONE
+    !***BEGIN PROLOGUE  XERSVE
+    !***SUBSIDIARY
+    !***PURPOSE  Record that an error has occurred.
+    !***LIBRARY   SLATEC (XERROR)
+    !***CATEGORY  R3
+    !***TYPE      ALL (XERSVE-A)
+    !***KEYWORDS  ERROR, XERROR
+    !***AUTHOR  Jones, R. E., (SNLA)
+    !***DESCRIPTION
+
+    ! *Usage:
+
+    !    INTEGER  KFLAG, NERR, LEVEL, ICOUNT
+    !    CHARACTER * (len) LIBRAR, SUBROU, MESSG
+
+    !    CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT)
+
+    ! *Arguments:
+
+    !    LIBRAR :IN    is the library that the message is from.
+    !    SUBROU :IN    is the SUBROUTINE that the message is from.
+    !    MESSG  :IN    is the message to be saved.
+    !    KFLAG  :IN    indicates the action to be performed.
+    !                  when KFLAG > 0, the message in MESSG is saved.
+    !                  when KFLAG=0 the tables will be dumped and
+    !                  cleared.
+    !                  when KFLAG < 0, the tables will be dumped and
+    !                  not cleared.
+    !    NERR   :IN    is the error number.
+    !    LEVEL  :IN    is the error severity.
+    !    ICOUNT :OUT   the number of times this message has been seen,
+    !                  or zero if the table has overflowed and does not
+    !                  contain this message specifically.  When KFLAG=0,
+    !                  ICOUNT will not be altered.
+
+    ! *Description:
+
+    !   Record that this error occurred and possibly dump and clear the
+    !   tables.
+
+    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+    !             Error-handling Package, SAND82-0800, Sandia
+    !             Laboratories, 1982.
+    !***ROUTINES CALLED  I1MACH, XGETUA
+    !***REVISION HISTORY  (YYMMDD)
+    !   800319  DATE WRITTEN
+    !   861211  REVISION DATE from Version 3.2
+    !   891214  Prologue converted to Version 4.0 format.  (BAB)
+    !   900413  Routine modified to remove reference to KFLAG.  (WRB)
+    !   900510  Changed to add LIBRARY NAME and SUBROUTINE to calling
+    !       sequence, use IF-THEN-ELSE, make number of saved entries
+    !       easily changeable, changed routine name from XERSAV to
+    !       XERSVE.  (RWC)
+    !   910626  Added LIBTAB and SUBTAB to SAVE statement.  (BKS)
+    !   920501  Reformatted the REFERENCES section.  (WRB)
+    !***END PROLOGUE  XERSVE
+    INTEGER, PARAMETER :: LENTAB = 10
+    INTEGER :: LUN(5)
+    CHARACTER(len = *) :: LIBRAR, SUBROU, MESSG
+    CHARACTER(len = 8) :: LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB
+    CHARACTER(len = 20) :: MESTAB(LENTAB), MES
+    DIMENSION NERTAB(LENTAB), LEVTAB(LENTAB), KOUNT(LENTAB)
+    SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
+    DATA KOUNTX/0/, NMSG/0/
+    INTEGER :: NERR, LEVEL, KONTRL
+    INTEGER :: NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
+    INTEGER :: KFLAG, ICOUNT, NUNIT, KUNIT, IUNIT, I
+    !***FIRST EXECUTABLE STATEMENT  XERSVE
+
+    IF (KFLAG<=0) THEN
+
+      !    Dump the table.
+
+      IF (NMSG==0) RETURN
+
+      !    Print to each unit.
+
+      CALL XGETUA (LUN, NUNIT)
+      DO KUNIT = 1, NUNIT
+        IUNIT = LUN(KUNIT)
+        IF (IUNIT==0) IUNIT = I1MACH(4)
+
+        !       Print the table header.
+
+        WRITE (IUNIT, 9000)
+
+        !       Print body of table.
+
+        DO I = 1, NMSG
+          WRITE (IUNIT, 9010) LIBTAB(I), SUBTAB(I), MESTAB(I), &
+                  NERTAB(I), LEVTAB(I), KOUNT(I)
+        END DO
+
+        !       Print number of other errors.
+
+        IF (KOUNTX/=0) WRITE (IUNIT, 9020) KOUNTX
+        WRITE (IUNIT, 9030)
+      END DO
+
+      !    Clear the error tables.
+
+      IF (KFLAG==0) THEN
+        NMSG = 0
+        KOUNTX = 0
+      ENDIF
+    ELSE
+
+      !    PROCESS A MESSAGE...
+      !    SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,
+      !    OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.
+
+      LIB = LIBRAR
+      SUB = SUBROU
+      MES = MESSG
+      DO I = 1, NMSG
+        IF (LIB==LIBTAB(I) .AND. SUB==SUBTAB(I) .AND. &
+                MES==MESTAB(I) .AND. NERR==NERTAB(I) .AND. &
+                LEVEL==LEVTAB(I)) THEN
+          KOUNT(I) = KOUNT(I) + 1
+          ICOUNT = KOUNT(I)
+          RETURN
+        ENDIF
+      END DO
+
+      IF (NMSG<LENTAB) THEN
+
+        !       Empty slot found for new message.
+
+        NMSG = NMSG + 1
+        LIBTAB(I) = LIB
+        SUBTAB(I) = SUB
+        MESTAB(I) = MES
+        NERTAB(I) = NERR
+        LEVTAB(I) = LEVEL
+        KOUNT (I) = 1
+        ICOUNT = 1
+      ELSE
+
+        !       Table is full.
+
+        KOUNTX = KOUNTX + 1
+        ICOUNT = 0
+      ENDIF
+    ENDIF
+    RETURN
+
+    ! Formats.
+
+    9000   FORMAT ('0          ERROR MESSAGE SUMMARY' / &
+            ' LIBRARY    SUBROUTINE MESSAGE START             NERR', &
+            '     LEVEL     COUNT')
+    9010   FORMAT (1X, A, 3X, A, 3X, A, 3I10)
+    9020   FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED = ', I10)
+    9030   FORMAT (1X)
+  END SUBROUTINE XERSVE
+
+  !DECK XERPRN
+  SUBROUTINE XERPRN(PREFIX, NPREF, MESSG, NWRAP)
+    IMPLICIT NONE
+    !***BEGIN PROLOGUE  XERPRN
+    !***SUBSIDIARY
+    !***PURPOSE  Print error messages processed by XERMSG.
+    !***LIBRARY   SLATEC (XERROR)
+    !***CATEGORY  R3C
+    !***TYPE      ALL (XERPRN-A)
+    !***KEYWORDS  ERROR MESSAGES, PRINTING, XERROR
+    !***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
+    !***DESCRIPTION
+
+    ! This routine sends one or more lines to each of the (up to five)
+    ! LOGICAL units to which error messages are to be sent.  This routine
+    ! is called several times by XERMSG, sometimes with a single line to
+    ! print and sometimes with a (potentially very long) message that may
+    ! wrap around into multiple lines.
+
+    ! PREFIX  Input argument of type CHARACTER.  This argument contains
+    !     characters to be put at the beginning of each line before
+    !     the body of the message.  No more than 16 characters of
+    !     PREFIX will be used.
+
+    ! NPREF   Input argument of type INTEGER.  This argument is the number
+    !     of characters to use from PREFIX.  If it is negative, the
+    !     intrinsic function LEN is used to determine its length.  If
+    !         it is zero, PREFIX is not used.  If it exceeds 16 or if
+    !     LEN(PREFIX) exceeds 16, ONLY the first 16 characters will be
+    !     used.  If NPREF is positive and the length of PREFIX is less
+    !     than NPREF, a copy of PREFIX extended with blanks to length
+    !     NPREF will be used.
+
+    ! MESSG   Input argument of type CHARACTER.  This is the text of a
+    !     message to be printed.  If it is a long message, it will be
+    !     broken into pieces for printing on multiple lines.  Each line
+    !     will start with the appropriate prefix and be followed by a
+    !     piece of the message.  NWRAP is the number of characters per
+    !     piece; that is, after each NWRAP characters, we break and
+    !     start a new line.  In addition the characters '$$' embedded
+    !     in MESSG are a sentinel for a new line.  The counting of
+    !     characters up to NWRAP starts over for each new line.  The
+    !     value of NWRAP typically used by XERMSG is 72 since many
+    !     older error messages in the SLATEC Library are laid out to
+    !     rely on wrap-around every 72 characters.
+
+    ! NWRAP   Input argument of type INTEGER.  This gives the maximum size
+    !     piece into which to break MESSG for printing on multiple
+    !     lines.  An embedded '$$' ends a line, and the count restarts
+    !     at the following character.  If a line break does not occur
+    !     on a blank (it would split a word) that word is moved to the
+    !     next line.  Values of NWRAP less than 16 will be treated as
+    !     16.  Values of NWRAP greater than 132 will be treated as 132.
+    !     The actual line length will be NPREF + NWRAP after NPREF has
+    !     been adjusted to fall between 0 and 16 and NWRAP has been
+    !     adjusted to fall between 16 and 132.
+
+    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+    !             Error-handling Package, SAND82-0800, Sandia
+    !             Laboratories, 1982.
+    !***ROUTINES CALLED  I1MACH, XGETUA
+    !***REVISION HISTORY  (YYMMDD)
+    !   880621  DATE WRITTEN
+    !   880708  REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF
+    !       JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK
+    !       THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE
+    !       SLASH CHARACTER IN FORMAT STATEMENTS.
+    !   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
+    !       STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK
+    !       LINES TO BE PRINTED.
+    !   890721  REVISED TO ADD A NEW FEATURE.  A NEGATIVE VALUE OF NPREF
+    !       CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH.
+    !   891013  REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH.
+    !   891214  Prologue converted to Version 4.0 format.  (WRB)
+    !   900510  Added code to break messages between words.  (RWC)
+    !   920501  Reformatted the REFERENCES section.  (WRB)
+    !***END PROLOGUE  XERPRN
+    CHARACTER(len = *) :: PREFIX, MESSG
+    INTEGER :: NPREF, NWRAP
+    CHARACTER(len = 148) :: CBUFF
+    INTEGER :: IU(5), NUNIT
+    CHARACTER(len = 2) :: NEWLIN
+    PARAMETER (NEWLIN = '$$')
+    INTEGER :: N, I, LPREF, LWRAP, LENMSG, NEXTC
+    INTEGER :: LPIECE, IDELTA
+    !***FIRST EXECUTABLE STATEMENT  XERPRN
+    CALL XGETUA(IU, NUNIT)
+
+    !   A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD
+    !   ERROR MESSAGE UNIT INSTEAD.  I1MACH(4) RETRIEVES THE STANDARD
+    !   ERROR MESSAGE UNIT.
+
+    N = I1MACH(4)
+    DO I = 1, NUNIT
+      IF (IU(I) == 0) IU(I) = N
+    END DO
+
+    !   LPREF IS THE LENGTH OF THE PREFIX.  THE PREFIX IS PLACED AT THE
+    !   BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING
+    !   THE REST OF THIS ROUTINE.
+
+    IF (NPREF < 0) THEN
+      LPREF = LEN(PREFIX)
+    ELSE
+      LPREF = NPREF
+    ENDIF
+    LPREF = MIN(16, LPREF)
+    IF (LPREF /= 0) CBUFF(1:LPREF) = PREFIX
+
+    !   LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE
+    !   TIME FROM MESSG TO PRINT ON ONE LINE.
+
+    LWRAP = MAX(16, MIN(132, NWRAP))
+
+    !   SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS.
+
+    LENMSG = LEN(MESSG)
+    N = LENMSG
+    DO I = 1, N
+      IF (MESSG(LENMSG:LENMSG) /= ' ') GO TO 30
+      LENMSG = LENMSG - 1
+    END DO
+    30   CONTINUE
+
+    !   IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
+
+    IF (LENMSG == 0) THEN
+      CBUFF(LPREF + 1:LPREF + 1) = ' '
+      DO I = 1, NUNIT
+        WRITE(IU(I), '(A)') CBUFF(1:LPREF + 1)
+      END DO
+      RETURN
+    ENDIF
+
+    !   SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING
+    !   STARTS.  FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL.
+    !   WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT.
+    !   WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED.
+
+    !   WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL.  THE
+    !   INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE
+    !   OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH
+    !   OF THE SECOND ARGUMENT.
+
+    !   THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE
+    !   FOLLOWING ORDER.  WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER
+    !   OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT
+    !   POSITION NEXTC.
+
+    !   LPIECE .EQ. 0   THE NEW LINE SENTINEL DOES NOT OCCUR IN THE
+    !                   REMAINDER OF THE CHARACTER STRING.  LPIECE
+    !                   SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC,
+    !                   WHICHEVER IS LESS.
+
+    !   LPIECE .EQ. 1   THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC:
+    !                   NEXTC).  LPIECE IS EFFECTIVELY ZERO, AND WE
+    !                   PRINT NOTHING TO AVOID PRODUCING UNNECESSARY
+    !                   BLANK LINES.  THIS TAKES CARE OF THE SITUATION
+    !                   WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF
+    !                   EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE
+    !                   SENTINEL FOLLOWED BY MORE CHARACTERS.  NEXTC
+    !                   SHOULD BE INCREMENTED BY 2.
+
+    !   LPIECE .GT. LWRAP+1  REDUCE LPIECE TO LWRAP.
+
+    !   ELSE            THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1
+    !                   RESET LPIECE = LPIECE-1.  NOTE THAT THIS
+    !                   PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ.
+    !                   LWRAP+1.  THAT IS, THE SENTINEL FALLS EXACTLY
+    !                   AT THE END OF A LINE.
+
+    NEXTC = 1
+    50   LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN)
+    IF (LPIECE == 0) THEN
+
+      !   THERE WAS NO NEW LINE SENTINEL FOUND.
+
+      IDELTA = 0
+      LPIECE = MIN(LWRAP, LENMSG + 1 - NEXTC)
+      IF (LPIECE < LENMSG + 1 - NEXTC) THEN
+        DO I = LPIECE + 1, 2, -1
+          IF (MESSG(NEXTC + I - 1:NEXTC + I - 1) == ' ') THEN
+            LPIECE = I - 1
+            IDELTA = 1
+            GOTO 54
+          ENDIF
+        END DO
+      ENDIF
+      54   CBUFF(LPREF + 1:LPREF + LPIECE) = MESSG(NEXTC:NEXTC + LPIECE - 1)
+      NEXTC = NEXTC + LPIECE + IDELTA
+    ELSEIF (LPIECE == 1) THEN
+
+      !   WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1).
+      !   DON'T PRINT A BLANK LINE.
+
+      NEXTC = NEXTC + 2
+      GO TO 50
+    ELSEIF (LPIECE > LWRAP + 1) THEN
+
+      !   LPIECE SHOULD BE SET DOWN TO LWRAP.
+
+      IDELTA = 0
+      LPIECE = LWRAP
+      DO I = LPIECE + 1, 2, -1
+        IF (MESSG(NEXTC + I - 1:NEXTC + I - 1) == ' ') THEN
+          LPIECE = I - 1
+          IDELTA = 1
+          GOTO 58
+        ENDIF
+      END DO
+      58   CBUFF(LPREF + 1:LPREF + LPIECE) = MESSG(NEXTC:NEXTC + LPIECE - 1)
+      NEXTC = NEXTC + LPIECE + IDELTA
+    ELSE
+
+      !   IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1.
+      !   WE SHOULD DECREMENT LPIECE BY ONE.
+
+      LPIECE = LPIECE - 1
+      CBUFF(LPREF + 1:LPREF + LPIECE) = MESSG(NEXTC:NEXTC + LPIECE - 1)
+      NEXTC = NEXTC + LPIECE + 2
+    ENDIF
+
+    !   PRINT
+
+    DO I = 1, NUNIT
+      WRITE(IU(I), '(A)') CBUFF(1:LPREF + LPIECE)
+    END DO
+
+    IF (NEXTC <= LENMSG) GO TO 50
+
+  END SUBROUTINE XERPRN
+
+  !DECK XERHLT
+  SUBROUTINE XERHLT(MESSG)
+    !***BEGIN PROLOGUE  XERHLT
+    !***SUBSIDIARY
+    !***PURPOSE  Abort program execution and print error message.
+    !***LIBRARY   SLATEC (XERROR)
+    !***CATEGORY  R3C
+    !***TYPE      ALL (XERHLT-A)
+    !***KEYWORDS  ABORT PROGRAM EXECUTION, ERROR, XERROR
+    !***AUTHOR  Jones, R. E., (SNLA)
+    !***DESCRIPTION
+
+    ! Abstract
+    !    ***Note*** machine dependent routine
+    !    XERHLT aborts the execution of the program.
+    !    The error message causing the abort is given in the calling
+    !    sequence, in case one needs it for printing on a dayfile,
+    !    for example.
+
+    ! Description of Parameters
+    !    MESSG is as in XERMSG.
+
+    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+    !             Error-handling Package, SAND82-0800, Sandia
+    !             Laboratories, 1982.
+    !***ROUTINES CALLED  (NONE)
+    !***REVISION HISTORY  (YYMMDD)
+    !   790801  DATE WRITTEN
+    !   861211  REVISION DATE from Version 3.2
+    !   891214  Prologue converted to Version 4.0 format.  (BAB)
+    !   900206  Routine changed from user-callable to subsidiary.  (WRB)
+    !   900510  Changed calling sequence to delete length of character
+    !       and changed routine name from XERABT to XERHLT.  (RWC)
+    !   920501  Reformatted the REFERENCES section.  (WRB)
+    !***END PROLOGUE  XERHLT
+    CHARACTER(len = *) :: MESSG
+    !***FIRST EXECUTABLE STATEMENT  XERHLT
+    STOP
+  END SUBROUTINE XERHLT
+
+  !DECK XERCNT
+  SUBROUTINE XERCNT(LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
+    IMPLICIT NONE
+    !***BEGIN PROLOGUE  XERCNT
+    !***SUBSIDIARY
+    !***PURPOSE  Allow user control over handling of errors.
+    !***LIBRARY   SLATEC (XERROR)
+    !***CATEGORY  R3C
+    !***TYPE      ALL (XERCNT-A)
+    !***KEYWORDS  ERROR, XERROR
+    !***AUTHOR  Jones, R. E., (SNLA)
+    !***DESCRIPTION
+
+    ! Abstract
+    !    Allows user control over handling of individual errors.
+    !    Just after each message is recorded, but before it is
+    !    processed any further (i.e., before it is printed or
+    !    a decision to abort is made), a CALL is made to XERCNT.
+    !    If the user has provided his own version of XERCNT, he
+    !    can then override the value of KONTROL used in processing
+    !    this message by redefining its value.
+    !    KONTRL may be set to any value from -2 to 2.
+    !    The meanings for KONTRL are the same as in XSETF, except
+    !    that the value of KONTRL changes only for this message.
+    !    If KONTRL is set to a value outside the range from -2 to 2,
+    !    it will be moved back into that range.
+
+    ! Description of Parameters
+
+    !  --Input--
+    !    LIBRAR - the library that the routine is in.
+    !    SUBROU - the SUBROUTINE that XERMSG is being called from
+    !    MESSG  - the first 20 characters of the error message.
+    !    NERR   - same as in the CALL to XERMSG.
+    !    LEVEL  - same as in the CALL to XERMSG.
+    !    KONTRL - the current value of the control flag as set
+    !             by a CALL to XSETF.
+
+    !  --Output--
+    !    KONTRL - the new value of KONTRL.  If KONTRL is not
+    !             defined, it will remain at its original value.
+    !             This changed value of control affects only
+    !             the current occurrence of the current message.
+
+    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+    !             Error-handling Package, SAND82-0800, Sandia
+    !             Laboratories, 1982.
+    !***ROUTINES CALLED  (NONE)
+    !***REVISION HISTORY  (YYMMDD)
+    !   790801  DATE WRITTEN
+    !   861211  REVISION DATE from Version 3.2
+    !   891214  Prologue converted to Version 4.0 format.  (BAB)
+    !   900206  Routine changed from user-callable to subsidiary.  (WRB)
+    !   900510  Changed calling sequence to INCLUDE LIBRARY and SUBROUTINE
+    !       names, changed routine name from XERCTL to XERCNT.  (RWC)
+    !   920501  Reformatted the REFERENCES section.  (WRB)
+    !***END PROLOGUE  XERCNT
+    CHARACTER(len = *) :: LIBRAR, SUBROU, MESSG
+    INTEGER :: NERR, LEVEL, KONTRL
+    !***FIRST EXECUTABLE STATEMENT  XERCNT
+
+  END SUBROUTINE XERCNT
+
+  !DECK J4SAVE
+  INTEGER FUNCTION J4SAVE(IWHICH, IVALUE, ISET)
+    IMPLICIT NONE
+    !***BEGIN PROLOGUE  J4SAVE
+    !***SUBSIDIARY
+    !***PURPOSE  Save or reCALL global variables needed by error
+    ! handling routines.
+    !***LIBRARY   SLATEC (XERROR)
+    !***TYPE      INTEGER (J4SAVE-I)
+    !***KEYWORDS  ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR
+    !***AUTHOR  Jones, R. E., (SNLA)
+    !***DESCRIPTION
+
+    ! Abstract
+    !    J4SAVE saves and recalls several global variables needed
+    !    by the library error handling routines.
+
+    ! Description of Parameters
+    !  --Input--
+    !    IWHICH - Index of item desired.
+    !            = 1 Refers to current error number.
+    !            = 2 Refers to current error control flag.
+    !            = 3 Refers to current unit number to which error
+    !                messages are to be sent.  (0 means use standard.)
+    !            = 4 Refers to the maximum number of times any
+    !                 message is to be printed (as set by XERMAX).
+    !            = 5 Refers to the total number of units to which
+    !                 each error message is to be written.
+    !            = 6 Refers to the 2nd unit for error messages
+    !            = 7 Refers to the 3rd unit for error messages
+    !            = 8 Refers to the 4th unit for error messages
+    !            = 9 Refers to the 5th unit for error messages
+    !    IVALUE - The value to be set for the IWHICH-th parameter,
+    !             if ISET is .TRUE. .
+    !    ISET   - If ISET=.TRUE., the IWHICH-th parameter will BE
+    !             given the value, IVALUE.  If ISET=.FALSE., the
+    !             IWHICH-th parameter will be unchanged, and IVALUE
+    !             is a dummy parameter.
+    !  --Output--
+    !    The (old) value of the IWHICH-th parameter will be returned
+    !    in the function value, J4SAVE.
+
+    !***SEE ALSO  XERMSG
+    !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
+    !             Error-handling Package, SAND82-0800, Sandia
+    !             Laboratories, 1982.
+    !***ROUTINES CALLED  (NONE)
+    !***REVISION HISTORY  (YYMMDD)
+    !   790801  DATE WRITTEN
+    !   891214  Prologue converted to Version 4.0 format.  (BAB)
+    !   900205  Minor modifications to prologue.  (WRB)
+    !   900402  Added TYPE section.  (WRB)
+    !   910411  Added KEYWORDS section.  (WRB)
+    !   920501  Reformatted the REFERENCES section.  (WRB)
+    !***END PROLOGUE  J4SAVE
+    LOGICAL :: ISET
+    INTEGER :: IPARAM(9)
+    SAVE IPARAM
+    DATA IPARAM(1), IPARAM(2), IPARAM(3), IPARAM(4)/0, 2, 0, 10/
+    DATA IPARAM(5)/1/
+    DATA IPARAM(6), IPARAM(7), IPARAM(8), IPARAM(9)/0, 0, 0, 0/
+    INTEGER :: IWHICH, IVALUE
+    !***FIRST EXECUTABLE STATEMENT  J4SAVE
+    J4SAVE = IPARAM(IWHICH)
+    IF (ISET) IPARAM(IWHICH) = IVALUE
+
+  END FUNCTION J4SAVE
+
+
+END MODULE lmdz_xer
Index: LMDZ6/branches/Amaury_dev/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmd/phys_output_write_mod.F90	(revision 5195)
+++ LMDZ6/branches/Amaury_dev/libf/phylmd/phys_output_write_mod.F90	(revision 5218)
@@ -2921,10 +2921,7 @@
             i = INDEX(outiso, '_', .TRUE.)
             outiso = outiso(1:i-1)//outiso(i+1:LEN_TRIM(outiso))
-            ! dayly
-            CALL xios_set_fieldgroup_attr("dd_2D_"//TRIM(outiso), enabled=.TRUE.)
-            CALL xios_set_fieldgroup_attr("dd_3D_"//TRIM(outiso), enabled=.TRUE.)
-            ! monthly
-            CALL xios_set_fieldgroup_attr("mm_2D_"//TRIM(outiso), enabled=.TRUE.)
-            CALL xios_set_fieldgroup_attr("mm_3D_"//TRIM(outiso), enabled=.TRUE.)
+
+            CALL xios_set_fieldgroup_attr("iso2D_"//TRIM(outiso), enabled=.TRUE.)
+            CALL xios_set_fieldgroup_attr("iso3D_"//TRIM(outiso), enabled=.TRUE.)
 
           ENDDO
Index: LMDZ6/branches/Amaury_dev/libf/phylmdiso/reevap.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/phylmdiso/reevap.F90	(revision 5195)
+++ LMDZ6/branches/Amaury_dev/libf/phylmdiso/reevap.F90	(revision 5218)
@@ -28,5 +28,5 @@
     REAL, DIMENSION(klon,klev,nqtot), INTENT(IN) ::     qx
     REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t_eva
-    REAL, DIMENSION(klon,klev,nqtot), INTENT(OUT) ::        d_qx_eva
+    REAL, DIMENSION(klon,klev,nqtot), INTENT(OUT) ::    d_qx_eva
 
     REAL za,zb,zdelta,zlvdcp,zlsdcp
@@ -35,5 +35,5 @@
     !--------Stochastic Boundary Layer Triggering: ALE_BL--------
 
-DO ixt=1,1+ntiso
+DO ixt = 1, 1+ntiso
     ! Re-evaporer l'eau liquide nuageuse
 
@@ -43,7 +43,7 @@
 !PRINT *,'rrevap ; fl_cor_ebil:',fl_cor_ebil,' iflag_ice_thermo:',iflag_ice_thermo,' RVTMP2',RVTMP2
     DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
-       DO i = 1, klon
+      DO i = 1, klon
 
-        IF (ixt==1) THEN
+        IF (ixt == 1) THEN ! water
          IF (fl_cor_ebil > 0) THEN
           !zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(q_seri(i,k)+ql_seri(i,k)+qs_seri(i,k)))
@@ -51,5 +51,5 @@
           zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*(qx(i,k,ivapcur)+qx(i,k,iliqcur)+qx(i,k,isolcur)))
           zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*(qx(i,k,ivapcur)+qx(i,k,iliqcur)+qx(i,k,isolcur)))
-         else
+         ELSE
            zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*qx(i,k,ivapcur))
            !jyg<
@@ -57,9 +57,27 @@
            !                  A verifier !!!
            zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*qx(i,k,ivapcur))
-         end if
+         END IF
          IF (iflag_ice_thermo == 0) THEN
              zlsdcp=zlvdcp
-         ENDIF
+         END IF
           !>jyg
+        END IF
+        IF (iflag_ice_thermo == 0) THEN
+          !pas necessaire a priori
+
+          zdelta = MAX(0., SIGN(1., RTT - t_seri(i, k)))
+          zdelta = 0.
+          zb = MAX(0.0, qx(i, k, iliqcur))
+          IF (ixt == 1) THEN
+            za = - MAX(0.0, qx(i, k, iliqcur)) &
+                    * (zlvdcp * (1. - zdelta) + zlsdcp * zdelta)
+            d_t_eva(i, k) = za
+          END IF
+          d_qx_eva(i, k, ivapcur) = zb
+          d_qx_eva(i, k, iliqcur) = -qx(i, k, iliqcur)
+          d_qx_eva(i, k, isolcur) = 0.
+
+        ELSE
+
              !CR: on r\'e-\'evapore eau liquide et glace
 
@@ -68,9 +86,10 @@
              !        za = - MAX(0.0,ql_seri(i,k)) &
              !             * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
-             za = - MAX(0.0,qx(i,k,iliqcur))*zlvdcp &
-                  - MAX(0.0,qx(i,k,iliqcur))*zlsdcp
-             d_t_eva(i,k) = za
+          IF (ixt == 1) THEN
+            za = - MAX(0.0, qx(i, k, iliqcur)) * zlvdcp &
+                    - MAX(0.0, qx(i, k, iliqcur)) * zlsdcp
+            d_t_eva(i, k) = za
 
-          endif !if (ixt.EQ.1) THEN
+          END IF
             !zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
             !d_q_eva(i,k) = zb
@@ -82,11 +101,11 @@
             d_qx_eva(i,k,iliqcur) = -qx(i,k,iliqcur)
             d_qx_eva(i,k,isolcur) = -qx(i,k,isolcur)
+        END IF
 
 
-       ENDDO
-    ENDDO
+      END DO
+    END DO
 
-    enddo ! do ixt=1,1+niso*(nzone +1)
-
+    END DO ! DO ixt = 1, 1+niso*(nzone +1)
 
 
