source: LMDZ4/trunk/libf/bibio/pchsp_95.F90

Last change on this file 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: 4.2 KB
Line 
1module pchsp_95_m
2
3  implicit none
4
5contains
6
7  function pchsp_95(x, f, ibeg, iend, vc_beg, vc_end)
8
9    ! PURPOSE: Set derivatives needed to determine the Hermite
10    ! representation of the cubic spline interpolant to given data,
11    ! with specified boundary conditions.
12
13    ! Part of the "pchip" package.
14
15    ! CATEGORY: E1A
16
17    ! KEYWORDS: cubic hermite interpolation, piecewise cubic
18    ! interpolation, spline interpolation
19
20    ! DESCRIPTION: "pchsp" stands for "Piecewise Cubic Hermite Spline"
21    ! Computes the Hermite representation of the cubic spline
22    ! interpolant to the data given in X and F satisfying the boundary
23    ! conditions specified by Ibeg, iend, vc_beg and VC_end.
24
25    ! The resulting piecewise cubic Hermite function may be evaluated
26    ! by "pchfe" or "pchfd".
27
28    ! NOTE: This is a modified version of C. de Boor's cubic spline
29    ! routine "cubspl".
30
31    ! REFERENCE: Carl de Boor, A Practical Guide to Splines, Springer,
32    ! 2001, pages 43-47
33
34    use assert_eq_m, only: assert_eq
35
36    real, intent(in):: x(:)
37    ! independent variable values
38    ! The elements of X must be strictly increasing:
39    !                X(I-1) < X(I),  I = 2...N.
40    !           (Error return if not.)
41    ! (error if size(x) < 2)
42
43    real, intent(in):: f(:)
44    !     dependent variable values to be interpolated
45    !  F(I) is value corresponding to X(I).
46
47    INTEGER, intent(in):: ibeg
48    !     desired boundary condition at beginning of data
49
50    !        IBEG = 0  to set pchsp_95(1) so that the third derivative is con-
51    !              tinuous at X(2).  This is the "not a knot" condition
52    !              provided by de Boor's cubic spline routine CUBSPL.
53    !              This is the default boundary condition.
54    !        IBEG = 1  if first derivative at X(1) is given in VC_BEG.
55    !        IBEG = 2  if second derivative at X(1) is given in VC_BEG.
56    !        IBEG = 3  to use the 3-point difference formula for pchsp_95(1).
57    !              (Reverts to the default boundary condition if size(x) < 3 .)
58    !        IBEG = 4  to use the 4-point difference formula for pchsp_95(1).
59    !              (Reverts to the default boundary condition if size(x) < 4 .)
60
61    !          NOTES:
62    !           1. An error return is taken if IBEG is out of range.
63    !           2. For the "natural" boundary condition, use IBEG=2 and
64    !              VC_BEG=0.
65
66    INTEGER, intent(in):: iend
67    !           IC(2) = IEND, desired condition at end of data.
68    !  IEND may take on the same values as IBEG, but applied to
69    !  derivative at X(N). In case IEND = 1 or 2, The value is given in VC_END.
70
71    !          NOTES:
72    !           1. An error return is taken if IEND is out of range.
73    !           2. For the "natural" boundary condition, use IEND=2 and
74    !              VC_END=0.
75
76    REAL, intent(in), optional:: vc_beg
77    ! desired boundary value, as indicated above.
78    !           VC_BEG need be set only if IBEG = 1 or 2 .
79
80    REAL, intent(in), optional:: vc_end
81    ! desired boundary value, as indicated above.
82    !           VC_END need be set only if Iend = 1 or 2 .
83
84    real pchsp_95(size(x))
85    ! derivative values at the data points
86    !           These values will determine the cubic spline interpolant
87    !           with the requested boundary conditions.
88    !           The value corresponding to X(I) is stored in
89    !                PCHSP_95(I),  I=1...N.
90
91    ! LOCAL VARIABLES:
92    real wk(2, size(x)) ! real array of working storage
93    INTEGER n ! number of data points
94    integer ierr, ic(2)
95    real vc(2)
96
97    !-------------------------------------------------------------------
98
99    n = assert_eq(size(x), size(f), "pchsp_95 n")
100    if ((ibeg == 1 .or. ibeg == 2) .and. .not. present(vc_beg)) then
101       print *, "vc_beg required for IBEG = 1 or 2"
102       stop 1
103    end if
104    if ((iend == 1 .or. iend == 2) .and. .not. present(vc_end)) then
105       print *, "vc_end required for IEND = 1 or 2"
106       stop 1
107    end if
108    ic = (/ibeg, iend/)
109    if (present(vc_beg)) vc(1) = vc_beg
110    if (present(vc_end)) vc(2) = vc_end
111    call PCHSP(IC, VC, N, X, F, pchsp_95, 1, WK, size(WK), IERR)
112    if (ierr /= 0) stop 1
113
114  END function pchsp_95
115
116end module pchsp_95_m
Note: See TracBrowser for help on using the repository browser.