source: LMDZ5/branches/testing/libf/misc/pchsp_95_m.F90 @ 2298

Last change on this file since 2298 was 2298, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes -r2237:2291 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 4.2 KB
RevLine 
[1687]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.