source: LMDZ5/branches/LMDZ5-dev032011/libf/bibio/pchfe.F @ 1495

Last change on this file since 1495 was 1425, checked in by lguez, 14 years ago

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

File size: 9.9 KB
Line 
1*DECK PCHFE
2      SUBROUTINE PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
3C***BEGIN PROLOGUE  PCHFE
4C***PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
5C            points.  May be used by itself for Hermite interpolation,
6C            or as an evaluator for PCHIM or PCHIC.
7C***LIBRARY   SLATEC (PCHIP)
8C***CATEGORY  E3
9C***TYPE      SINGLE PRECISION (PCHFE-S, DPCHFE-D)
10C***KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
11C             PIECEWISE CUBIC EVALUATION
12C***AUTHOR  Fritsch, F. N., (LLNL)
13C             Lawrence Livermore National Laboratory
14C             P.O. Box 808  (L-316)
15C             Livermore, CA  94550
16C             FTS 532-4275, (510) 422-4275
17C***DESCRIPTION
18C
19C          PCHFE:  Piecewise Cubic Hermite Function Evaluator
20C
21C     Evaluates the cubic Hermite function defined by  N, X, F, D  at
22C     the points  XE(J), J=1(1)NE.
23C
24C     To provide compatibility with PCHIM and PCHIC, includes an
25C     increment between successive values of the F- and D-arrays.
26C
27C ----------------------------------------------------------------------
28C
29C  Calling sequence:
30C
31C        PARAMETER  (INCFD = ...)
32C        INTEGER  N, NE, IERR
33C        REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
34C        LOGICAL  SKIP
35C
36C        CALL  PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
37C
38C   Parameters:
39C
40C     N -- (input) number of data points.  (Error return if N.LT.2 .)
41C
42C     X -- (input) real array of independent variable values.  The
43C           elements of X must be strictly increasing:
44C                X(I-1) .LT. X(I),  I = 2(1)N.
45C           (Error return if not.)
46C
47C     F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
48C           the value corresponding to X(I).
49C
50C     D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
51C           the value corresponding to X(I).
52C
53C     INCFD -- (input) increment between successive values in F and D.
54C           (Error return if  INCFD.LT.1 .)
55C
56C     SKIP -- (input/output) logical variable which should be set to
57C           .TRUE. if the user wishes to skip checks for validity of
58C           preceding parameters, or to .FALSE. otherwise.
59C           This will save time in case these checks have already
60C           been performed (say, in PCHIM or PCHIC).
61C           SKIP will be set to .TRUE. on normal return.
62C
63C     NE -- (input) number of evaluation points.  (Error return if
64C           NE.LT.1 .)
65C
66C     XE -- (input) real array of points at which the function is to be
67C           evaluated.
68C
69C          NOTES:
70C           1. The evaluation will be most efficient if the elements
71C              of XE are increasing relative to X;
72C              that is,   XE(J) .GE. X(I)
73C              implies    XE(K) .GE. X(I),  all K.GE.J .
74C           2. If any of the XE are outside the interval [X(1),X(N)],
75C              values are extrapolated from the nearest extreme cubic,
76C              and a warning error is returned.
77C
78C     FE -- (output) real array of values of the cubic Hermite function
79C           defined by  N, X, F, D  at the points  XE.
80C
81C     IERR -- (output) error flag.
82C           Normal return:
83C              IERR = 0  (no errors).
84C           Warning error:
85C              IERR.GT.0  means that extrapolation was performed at
86C                 IERR points.
87C           "Recoverable" errors:
88C              IERR = -1  if N.LT.2 .
89C              IERR = -2  if INCFD.LT.1 .
90C              IERR = -3  if the X-array is not strictly increasing.
91C              IERR = -4  if NE.LT.1 .
92C             (The FE-array has not been changed in any of these cases.)
93C               NOTE:  The above errors are checked in the order listed,
94C                   and following arguments have **NOT** been validated.
95C
96C***REFERENCES  (NONE)
97C***ROUTINES CALLED  CHFEV, XERMSG
98C***REVISION HISTORY  (YYMMDD)
99C   811020  DATE WRITTEN
100C   820803  Minor cosmetic changes for release 1.
101C   870707  Minor cosmetic changes to prologue.
102C   890531  Changed all specific intrinsics to generic.  (WRB)
103C   890831  Modified array declarations.  (WRB)
104C   890831  REVISION DATE from Version 3.2
105C   891214  Prologue converted to Version 4.0 format.  (BAB)
106C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
107C***END PROLOGUE  PCHFE
108C  Programming notes:
109C
110C     1. To produce a double precision version, simply:
111C        a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they
112C           occur,
113C        b. Change the real declaration to double precision,
114C
115C     2. Most of the coding between the call to CHFEV and the end of
116C        the IR-loop could be eliminated if it were permissible to
117C        assume that XE is ordered relative to X.
118C
119C     3. CHFEV does not assume that X1 is less than X2.  thus, it would
120C        be possible to write a version of PCHFE that assumes a strict-
121C        ly decreasing X-array by simply running the IR-loop backwards
122C        (and reversing the order of appropriate tests).
123C
124C     4. The present code has a minor bug, which I have decided is not
125C        worth the effort that would be required to fix it.
126C        If XE contains points in [X(N-1),X(N)], followed by points .LT.
127C        X(N-1), followed by points .GT.X(N), the extrapolation points
128C        will be counted (at least) twice in the total returned in IERR.
129C
130C  DECLARE ARGUMENTS.
131C
132      INTEGER  N, INCFD, NE, IERR
133      REAL  X(*), F(INCFD,*), D(INCFD,*), XE(*), FE(*)
134      LOGICAL  SKIP
135C
136C  DECLARE LOCAL VARIABLES.
137C
138      INTEGER  I, IERC, IR, J, JFIRST, NEXT(2), NJ
139C
140C  VALIDITY-CHECK ARGUMENTS.
141C
142C***FIRST EXECUTABLE STATEMENT  PCHFE
143      IF (SKIP)  GO TO 5
144C
145      IF ( N.LT.2 )  GO TO 5001
146      IF ( INCFD.LT.1 )  GO TO 5002
147      DO 1  I = 2, N
148         IF ( X(I).LE.X(I-1) )  GO TO 5003
149    1 CONTINUE
150C
151C  FUNCTION DEFINITION IS OK, GO ON.
152C
153    5 CONTINUE
154      IF ( NE.LT.1 )  GO TO 5004
155      IERR = 0
156      SKIP = .TRUE.
157C
158C  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
159C                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
160      JFIRST = 1
161      IR = 2
162   10 CONTINUE
163C
164C     SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
165C
166         IF (JFIRST .GT. NE)  GO TO 5000
167C
168C     LOCATE ALL POINTS IN INTERVAL.
169C
170         DO 20  J = JFIRST, NE
171            IF (XE(J) .GE. X(IR))  GO TO 30
172   20    CONTINUE
173         J = NE + 1
174         GO TO 40
175C
176C     HAVE LOCATED FIRST POINT BEYOND INTERVAL.
177C
178   30    CONTINUE
179         IF (IR .EQ. N)  J = NE + 1
180C
181   40    CONTINUE
182         NJ = J - JFIRST
183C
184C     SKIP EVALUATION IF NO POINTS IN INTERVAL.
185C
186         IF (NJ .EQ. 0)  GO TO 50
187C
188C     EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
189C
190C       ----------------------------------------------------------------
191        CALL CHFEV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR),
192     *              NJ, XE(JFIRST), FE(JFIRST), NEXT, IERC)
193C       ----------------------------------------------------------------
194         IF (IERC .LT. 0)  GO TO 5005
195C
196         IF (NEXT(2) .EQ. 0)  GO TO 42
197C        IF (NEXT(2) .GT. 0)  THEN
198C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
199C           RIGHT OF X(IR).
200C
201            IF (IR .LT. N)  GO TO 41
202C           IF (IR .EQ. N)  THEN
203C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
204               IERR = IERR + NEXT(2)
205               GO TO 42
206   41       CONTINUE
207C           ELSE
208C              WE SHOULD NEVER HAVE GOTTEN HERE.
209               GO TO 5005
210C           ENDIF
211C        ENDIF
212   42    CONTINUE
213C
214         IF (NEXT(1) .EQ. 0)  GO TO 49
215C        IF (NEXT(1) .GT. 0)  THEN
216C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
217C           LEFT OF X(IR-1).
218C
219            IF (IR .GT. 2)  GO TO 43
220C           IF (IR .EQ. 2)  THEN
221C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
222               IERR = IERR + NEXT(1)
223               GO TO 49
224   43       CONTINUE
225C           ELSE
226C              XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
227C              EVALUATION INTERVAL.
228C
229C              FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
230               DO 44  I = JFIRST, J-1
231                  IF (XE(I) .LT. X(IR-1))  GO TO 45
232   44          CONTINUE
233C              NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
234C                     IN CHFEV.
235               GO TO 5005
236C
237   45          CONTINUE
238C              RESET J.  (THIS WILL BE THE NEW JFIRST.)
239               J = I
240C
241C              NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
242               DO 46  I = 1, IR-1
243                  IF (XE(J) .LT. X(I)) GO TO 47
244   46          CONTINUE
245C              NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
246C
247   47          CONTINUE
248C              AT THIS POINT, EITHER  XE(J) .LT. X(1)
249C                 OR      X(I-1) .LE. XE(J) .LT. X(I) .
250C              RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
251C              CYCLING.
252               IR = MAX(1, I-1)
253C           ENDIF
254C        ENDIF
255   49    CONTINUE
256C
257         JFIRST = J
258C
259C     END OF IR-LOOP.
260C
261   50 CONTINUE
262      IR = IR + 1
263      IF (IR .LE. N)  GO TO 10
264C
265C  NORMAL RETURN.
266C
267 5000 CONTINUE
268      RETURN
269C
270C  ERROR RETURNS.
271C
272 5001 CONTINUE
273C     N.LT.2 RETURN.
274      IERR = -1
275      CALL XERMSG ('SLATEC', 'PCHFE',
276     +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
277      RETURN
278C
279 5002 CONTINUE
280C     INCFD.LT.1 RETURN.
281      IERR = -2
282      CALL XERMSG ('SLATEC', 'PCHFE', 'INCREMENT LESS THAN ONE', IERR,
283     +   1)
284      RETURN
285C
286 5003 CONTINUE
287C     X-ARRAY NOT STRICTLY INCREASING.
288      IERR = -3
289      CALL XERMSG ('SLATEC', 'PCHFE', 'X-ARRAY NOT STRICTLY INCREASING'
290     +   , IERR, 1)
291      RETURN
292C
293 5004 CONTINUE
294C     NE.LT.1 RETURN.
295      IERR = -4
296      CALL XERMSG ('SLATEC', 'PCHFE',
297     +   'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
298      RETURN
299C
300 5005 CONTINUE
301C     ERROR RETURN FROM CHFEV.
302C   *** THIS CASE SHOULD NEVER OCCUR ***
303      IERR = -5
304      CALL XERMSG ('SLATEC', 'PCHFE',
305     +   'ERROR RETURN FROM CHFEV -- FATAL', IERR, 2)
306      RETURN
307C------------- LAST LINE OF PCHFE FOLLOWS ------------------------------
308      END
Note: See TracBrowser for help on using the repository browser.