source: LMDZ4/trunk/libf/cosp/array_lib.F90 @ 4934

Last change on this file since 4934 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.9 KB
Line 
1! ARRAY_LIB: Array procedures for F90
2! Compiled/Modified:
3!   07/01/06  John Haynes (haynes@atmos.colostate.edu)
4!
5! infind (function)
6! lin_interpolate (function)
7 
8  module array_lib
9  implicit none
10
11  contains
12
13! ----------------------------------------------------------------------------
14! function INFIND
15! ----------------------------------------------------------------------------
16  function infind(list,val,sort,dist)
17  use m_mrgrnk
18  implicit none
19!
20! Purpose:
21!   Finds the index of an array that is closest to a value, plus the
22!   difference between the value found and the value specified
23!
24! Inputs:
25!   [list]   an array of sequential values
26!   [val]    a value to locate
27! Optional input:
28!   [sort]   set to 1 if [list] is in unknown/non-sequential order
29!
30! Returns:
31!   index of [list] that is closest to [val]
32!
33! Optional output:
34!   [dist]   set to variable containing [list([result])] - [val]
35!
36! Requires:
37!   mrgrnk library
38!
39! Created:
40!   10/16/03  John Haynes (haynes@atmos.colostate.edu)
41! Modified:
42!   01/31/06  IDL to Fortran 90
43 
44! ----- INPUTS -----
45  real*8, dimension(:), intent(in) :: list
46  real*8, intent(in) :: val 
47  integer, intent(in), optional :: sort
48 
49! ----- OUTPUTS -----
50  integer*4 :: infind
51  real*8, intent(out), optional :: dist
52
53! ----- INTERNAL -----
54  real*8, dimension(size(list)) :: lists
55  integer*4 :: nlist, result, tmp(1), sort_list
56  integer*4, dimension(size(list)) :: mask, idx
57
58  if (present(sort)) then
59    sort_list = sort
60  else
61    sort_list = 0
62  endif 
63
64  nlist = size(list)
65  if (sort_list == 1) then
66    call mrgrnk(list,idx)
67    lists = list(idx)
68  else
69    lists = list
70  endif
71
72  if (val >= lists(nlist)) then
73    result = nlist
74  else if (val <= lists(1)) then
75    result = 1
76  else
77    mask(:) = 0
78    where (lists < val) mask = 1
79      tmp = minloc(mask,1)
80      if (abs(lists(tmp(1)-1)-val) < abs(lists(tmp(1))-val)) then
81        result = tmp(1) - 1
82      else
83        result = tmp(1)
84      endif
85  endif
86  if (present(dist)) dist = lists(result)-val
87  if (sort_list == 1) then
88    infind = idx(result)
89  else
90    infind = result
91  endif
92
93  end function infind
94
95! ----------------------------------------------------------------------------
96! function LIN_INTERPOLATE
97! ---------------------------------------------------------------------------- 
98  subroutine lin_interpolate(yarr,xarr,yyarr,xxarr,tol)
99  use m_mrgrnk
100  implicit none
101!
102! Purpose:
103!   linearly interpolate a set of y2 values given a set of y1,x1,x2
104!
105! Inputs:
106!   [yarr]    an array of y1 values
107!   [xarr]    an array of x1 values
108!   [xxarr]   an array of x2 values
109!   [tol]     maximum distance for a match
110!
111! Output:
112!   [yyarr]   interpolated array of y2 values
113!
114! Requires:
115!   mrgrnk library
116!
117! Created:
118!   06/07/06  John Haynes (haynes@atmos.colostate.edu)
119
120! ----- INPUTS -----
121  real*8, dimension(:), intent(in) :: yarr, xarr, xxarr
122  real*8, intent(in) :: tol
123
124! ----- OUTPUTS -----
125  real*8, dimension(size(xxarr)), intent(out) :: yyarr
126
127! ----- INTERNAL -----
128  real*8, dimension(size(xarr)) :: ysort, xsort
129  integer*4, dimension(size(xarr)) :: ist
130  integer*4 :: nx, nxx, i, iloc
131  real*8 :: d, m
132
133  nx = size(xarr)
134  nxx = size(xxarr)
135
136! // xsort, ysort are sorted versions of xarr, yarr 
137  call mrgrnk(xarr,ist)
138  ysort = yarr(ist)
139  xsort = xarr(ist)
140 
141  do i=1,nxx
142    iloc = infind(xsort,xxarr(i),dist=d)
143    if (d > tol) then
144      print *, 'interpolation error'
145      stop
146    endif
147    if (iloc == nx) then
148!     :: set to the last value
149      yyarr(i) = ysort(nx)
150    else
151!     :: is there another closeby value?
152      if (abs(xxarr(i)-xsort(iloc+1)) < 2*tol) then
153!       :: yes, do a linear interpolation     
154        m = (ysort(iloc+1)-ysort(iloc))/(xsort(iloc+1)-xsort(iloc))
155        yyarr(i) = ysort(iloc) + m*(xxarr(i)-xsort(iloc))
156      else
157!       :: no, set to the only nearby value
158        yyarr(i) = ysort(iloc)
159      endif
160    endif
161  enddo
162 
163  end subroutine lin_interpolate
164
165  end module array_lib
Note: See TracBrowser for help on using the repository browser.