source: LMDZ6/trunk/libf/misc/pchfe.f90 @ 5425

!DECK PCHFE
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.LT.2 )  GO TO 5001
  IF ( INCFD.LT.1 )  GO TO 5002
  DO  I = 2, N
     IF ( X(I).LE.X(I-1) )  GO TO 5003
  END DO
  !
  !  FUNCTION DEFINITION IS OK, GO ON.
  !
    5   CONTINUE
  IF ( NE.LT.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 .GT. NE)  GO TO 5000
  !
  ! LOCATE ALL POINTS IN INTERVAL.
  !
     DO  J = JFIRST, NE
        IF (XE(J) .GE. X(IR))  GO TO 30
     END DO
     J = NE + 1
     GO TO 40
  !
  ! HAVE LOCATED FIRST POINT BEYOND INTERVAL.
  !
   30   CONTINUE
     IF (IR .EQ. N)  J = NE + 1
  !
   40   CONTINUE
     NJ = J - JFIRST
  !
  ! SKIP EVALUATION IF NO POINTS IN INTERVAL.
  !
     IF (NJ .EQ. 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 .LT. 0)  GO TO 5005
  !
     IF (NEXT(2) .EQ. 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 .LT. 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) .EQ. 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 .GT. 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) .LT. 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) .LT. 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 .LE. 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