source: LMDZ4/trunk/libf/bibio/interpolation.F90 @ 2221

Last change on this file since 2221 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

File size: 3.7 KB
Line 
1! $Id$
2module interpolation
3
4  ! From Press et al., 1996, version 2.10a
5  ! B3 Interpolation and Extrapolation
6
7  IMPLICIT NONE
8
9contains
10
11  pure FUNCTION locate(xx,x)
12
13    REAL, DIMENSION(:), INTENT(IN) :: xx
14    REAL, INTENT(IN) :: x
15    INTEGER  locate
16
17    ! Given an array xx(1:N), and given a value x, returns a value j,
18    ! between 0 and N, such that x is between xx(j) and xx(j + 1). xx
19    ! must be monotonic, either increasing or decreasing. j = 0 or j =
20    ! N is returned to indicate that x is out of range. This
21    ! procedure should not be called with a zero-sized array argument.
22    ! See notes.
23
24    INTEGER  n,jl,jm,ju
25    LOGICAL  ascnd
26
27    !----------------------------
28
29    n=size(xx)
30    ascnd = (xx(n) >= xx(1))
31    ! (True if ascending order of table, false otherwise.)
32    ! Initialize lower and upper limits:
33    jl=0
34    ju=n+1
35    do while (ju-jl > 1)
36       jm=(ju+jl)/2 ! Compute a midpoint,
37       if (ascnd .eqv. (x >= xx(jm))) then
38          jl=jm ! and replace either the lower limit
39       else
40          ju=jm ! or the upper limit, as appropriate.
41       end if
42    end do
43    ! {ju == jl + 1}
44
45    ! {(ascnd .and. xx(jl) <= x < xx(jl+1))
46    !  .neqv.
47    !  (.not. ascnd .and. xx(jl+1) <= x < xx(jl))}
48
49    ! Then set the output, being careful with the endpoints:
50    if (x == xx(1)) then
51       locate=1
52    else if (x == xx(n)) then
53       locate=n-1
54    else
55       locate=jl
56    end if
57
58  END FUNCTION locate
59
60  !***************************
61
62  pure SUBROUTINE hunt(xx,x,jlo)
63
64    ! Given an array xx(1:N ), and given a value x, returns a value
65    ! jlo such that x is between xx(jlo) and xx(jlo+1). xx must be
66    ! monotonic, either increasing or decreasing. jlo = 0 or jlo = N is
67    ! returned to indicate that x is out of range. jlo on input is taken as
68    ! the initial guess for jlo on output.
69    ! Modified so that it uses the information "jlo = 0" on input.
70
71    INTEGER, INTENT(INOUT) :: jlo
72    REAL, INTENT(IN) :: x
73    REAL, DIMENSION(:), INTENT(IN) :: xx
74    INTEGER  n,inc,jhi,jm
75    LOGICAL  ascnd, hunt_up
76
77    !-----------------------------------------------------
78
79    n=size(xx)
80    ascnd = (xx(n) >= xx(1))
81    ! (True if ascending order of table, false otherwise.)
82    if (jlo < 0 .or. jlo > n) then
83       ! Input guess not useful. Go immediately to bisection.
84       jlo=0
85       jhi=n+1
86    else
87       inc=1 ! Set the hunting increment.
88       if (jlo == 0) then
89          hunt_up = .true.
90       else
91          hunt_up = x >= xx(jlo) .eqv. ascnd
92       end if
93       if (hunt_up) then ! Hunt up:
94          do
95             jhi=jlo+inc
96             if (jhi > n) then ! Done hunting, since off end of table.
97                jhi=n+1
98                exit
99             else
100                if (x < xx(jhi) .eqv. ascnd) exit
101                jlo=jhi ! Not done hunting,
102                inc=inc+inc ! so double the increment
103             end if
104          end do ! and try again.
105       else ! Hunt down:
106          jhi=jlo
107          do
108             jlo=jhi-inc
109             if (jlo < 1) then ! Done hunting, since off end of table.
110                jlo=0
111                exit
112             else
113                if (x >= xx(jlo) .eqv. ascnd) exit
114                jhi=jlo ! Not done hunting,
115                inc=inc+inc ! so double the increment
116             end if
117          end do ! and try again.
118       end if
119    end if ! Done hunting, value bracketed.
120
121    do ! Hunt is done, so begin the final bisection phase:
122       if (jhi-jlo <= 1) then
123          if (x == xx(n)) jlo=n-1
124          if (x == xx(1)) jlo=1
125          exit
126       else
127          jm=(jhi+jlo)/2
128          if (x >= xx(jm) .eqv. ascnd) then
129             jlo=jm
130          else
131             jhi=jm
132          end if
133       end if
134    end do
135
136  END SUBROUTINE hunt
137
138end module interpolation
Note: See TracBrowser for help on using the repository browser.