! 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 (XXMA) 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