source: trunk/LMDZ.TITAN/libf/muphytitan/lintdset.f90 @ 3094

Last change on this file since 3094 was 1793, checked in by jvatant, 7 years ago

Making Titan's hazy again, part I
+ Added the source folder libf/muphytitan which contains

YAMMS ( Titan's microphysical model ) from J. Burgalat

+ Modif. compilation files linked to this change
JVO
File size: 38.8 KB
RevLine 
[1793]1! Copyright Université Reims Champagnne-Ardenne (2010-2015)
2! contributor: Jérémie Burgalat
3!
4! jeremie.burgalat@univ-reims.fr
5!
6! This software is a computer program whose purpose is to compute multi-variate
7! linear interpolation.
8!
9! This software is governed by the CeCILL-B license under French law and
10! abiding by the rules of distribution of free software.  You can  use,
11! modify and/ or redistribute the software under the terms of the CeCILL-B
12! license as circulated by CEA, CNRS and INRIA at the following URL
13! "http://www.cecill.info".
14!
15! As a counterpart to the access to the source code and  rights to copy,
16! modify and redistribute granted by the license, users are provided only
17! with a limited warranty  and the software's author,  the holder of the
18! economic rights,  and the successive licensors  have only  limited
19! liability.
20!
21! In this respect, the user's attention is drawn to the risks associated
22! with loading,  using,  modifying and/or developing or reproducing the
23! software by the user in light of its specific status of free software,
24! that may mean  that it is complicated to manipulate,  and  that  also
25! therefore means  that it is reserved for developers  and  experienced
26! professionals having in-depth computer knowledge. Users are therefore
27! encouraged to load and test the software's suitability as regards their
28! requirements in conditions enabling the security of their systems and/or
29! data to be ensured and,  more generally, to use and operate it in the
30! same conditions as regards security.
31!
32! The fact that you are presently reading this means that you have had
33! knowledge of the CeCILL-B license and that you accept its terms.
34
35!! file: lintdset.f90
36!! summary: Linear interpolation generic interfaces
37!! author: burgalat
38!! date: 2010-2014
39
40
41MODULE LINTDSET
42  !! Generic linear interpolation using simple [[datasets(module)]]
43  USE LINT_PREC
44  USE LINTCORE
45  USE DATASETS
46  IMPLICIT NONE
47
48  PUBLIC  :: lint_dset, hdcd_lint_dset,              &
49             dset1d, dset2d, dset3d, dset4d, dset5d, &
50             read_dset, clear_dset, is_in
51  PRIVATE
52
53  !> Interface to multi-variate linear interpolation using data set
54  !!
55  !! For __n__, the number of variables, user must provide:
56  !!
57  !! - __n__ scalar(s)/vector(s) with the coordinate(s) of the point(s) to
58  !!   interpolate.
59  !! - A single __n__-D data set with the tabulated function.
60  !! - A function, satisfying the [[lintcore(module):locate(interface)]] interface, that
61  !!   should locate the neighbourhood of the point to interpolate in the given data
62  !!   set.
63  !! - Either a scalar or vector which stores the interpolated values.
64  !!
65  !! This interface also provides vectorized version of the linear
66  !! interpolation. In this case, input coordinates should be vector of the same
67  !! size as well as the output argument.
68  INTERFACE lint_dset
69    MODULE PROCEDURE l1d_dset_sc, l2d_dset_sc, l3d_dset_sc, l4d_dset_sc, &
70                     l5d_dset_sc
71    MODULE PROCEDURE l1d_dset_ve, l2d_dset_ve, l3d_dset_ve, l4d_dset_ve, &
72                     l5d_dset_ve
73  END INTERFACE
74
75  !> Interface to multi-variate hard-coded linear interpolation using data set
76  !!
77  !! Same remarks as in [[lintdset(module):lint_dset(interface)]] apply here.
78  INTERFACE hdcd_lint_dset
79    MODULE PROCEDURE hdcdl1d_dset_sc, hdcdl2d_dset_sc, hdcdl3d_dset_sc, &
80                     hdcdl4d_dset_sc, hdcdl5d_dset_sc
81    MODULE PROCEDURE hdcdl1d_dset_ve, hdcdl2d_dset_ve, hdcdl3d_dset_ve, &
82                     hdcdl4d_dset_ve, hdcdl5d_dset_ve
83  END INTERFACE
84
85  CONTAINS
86
87  FUNCTION l1d_dset_sc(x,set,locator,res) RESULT(ret)
88    !! Wrapper interface for 1D linear interpolation (scalar)
89    !!
90    !! @warning
91    !! On error, __res__ output value is undefined.
92    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
93    TYPE(dset1d), INTENT(in)   :: set     !! Dataset with the tabulated function
94    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
95    PROCEDURE(locate)          :: locator !! Locator function
96    LOGICAL :: ret                        !! .true. on success, .false. otherwise
97    REAL(kind=wp), DIMENSION(2,2) :: g1d
98    INTEGER                       :: i,ix
99    ret = .false.
100    res = HUGE(res)
101    ix = locator(x,set%x) ; IF (ix == 0) RETURN
102    DO i=0,1
103      g1d(i+1,1) = set%x(ix+i)
104      g1d(i+1,2) = set%data(ix+i)
105    ENDDO
106    res = lintc_((/x/),g1d)
107    ret = .true.
108  END FUNCTION l1d_dset_sc
109
110  FUNCTION l2d_dset_sc(x,y,set,locator,res) RESULT(ret)
111    !! Wrapper interface for 2D linear interpolation (scalar)
112    !!
113    !! @warning
114    !! On error, __res__ output value is undefined.
115    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
116    REAL(kind=wp), INTENT(in)  :: y       !! Y coordinate of the point to interpolate
117    TYPE(dset2d), INTENT(in)   :: set     !! Dataset with the tabulated function
118    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
119    PROCEDURE(locate)          :: locator !! Locator function
120    LOGICAL :: ret                        !! .true. on success, .false. otherwise
121    REAL(kind=wp), DIMENSION(4,3) :: g2d
122    INTEGER                       :: i,ix0,iy0,a,b
123    ret = .false.
124    res = HUGE(res)
125    ix0 = locator(x,set%x) ; IF (ix0 == 0) RETURN
126    iy0 = locator(y,set%y) ; IF (iy0 == 0) RETURN
127    DO i=1,4
128      a=ix0+MOD((i-1),2)   ; g2d(i,1) = set%x(a)
129      b=iy0+MOD((i-1)/2,2) ; g2d(i,2) = set%y(b)
130      g2d(i,3) = set%data(a,b)
131    ENDDO
132    res = lintc_((/x,y/),g2d)
133    ret = .true.
134  END FUNCTION l2d_dset_sc
135
136  FUNCTION l3d_dset_sc(x,y,z,set,locator,res) RESULT(ret)
137    !! Wrapper interface for 3D linear interpolation (scalar)
138    !!
139    !! @warning
140    !! On error, __res__ output value is undefined.
141    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
142    REAL(kind=wp), INTENT(in)  :: y       !! Y coordinate of the point to interpolate
143    REAL(kind=wp), INTENT(in)  :: z       !! Z coordinate of the point to interpolate
144    TYPE(dset3d), INTENT(in)   :: set     !! Dataset with the tabulated function
145    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
146    PROCEDURE(locate)          :: locator !! Locator function
147    LOGICAL :: ret                        !! .true. on success, .false. otherwise
148    REAL(kind=wp), DIMENSION(8,4) :: g3d
149    INTEGER                       :: i,ix0,iy0,iz0,a,b,c
150    ret = .false.
151    res = HUGE(res)
152    ix0 = locator(x,set%x) ; IF (ix0 == 0) RETURN
153    iy0 = locator(y,set%y) ; IF (iy0 == 0) RETURN
154    iz0 = locator(z,set%z) ; IF (iz0 == 0) RETURN
155    DO i=1,8
156      a=ix0+MOD((i-1),2)   ; g3d(i,1) = set%x(a)
157      b=iy0+MOD((i-1)/2,2) ; g3d(i,2) = set%y(b)
158      c=iz0+MOD((i-1)/4,2) ; g3d(i,3) = set%z(c)
159      g3d(i,4) = set%data(a,b,c)
160    ENDDO
161    res = lintc_((/x,y,z/),g3d)
162    ret = .true.
163  END FUNCTION l3d_dset_sc
164
165  FUNCTION l4d_dset_sc(x,y,z,t,set,locator,res) RESULT(ret)
166    !! Wrapper interface for 4D linear interpolation (scalar)
167    !!
168    !! @warning
169    !! On error, __res__ output value is undefined.
170    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
171    REAL(kind=wp), INTENT(in)  :: y       !! Y coordinate of the point to interpolate
172    REAL(kind=wp), INTENT(in)  :: z       !! Z coordinate of the point to interpolate
173    REAL(kind=wp), INTENT(in)  :: t       !! T coordinate of the point to interpolate
174    TYPE(dset4d), INTENT(in)   :: set     !! Dataset with the tabulated function
175    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
176    PROCEDURE(locate)          :: locator !! Locator function
177    LOGICAL :: ret                        !! .true. on success, .false. otherwise
178    REAL(kind=wp), DIMENSION(16,5) :: g4d
179    INTEGER                        :: i,ix0,iy0,iz0,it0,a,b,c,d
180    ret = .false.
181    res = HUGE(res)
182    ix0 = locator(x,set%x) ; IF (ix0 == 0) RETURN
183    iy0 = locator(y,set%y) ; IF (iy0 == 0) RETURN
184    iz0 = locator(z,set%z) ; IF (iz0 == 0) RETURN
185    it0 = locator(t,set%t) ; IF (it0 == 0) RETURN
186    DO i=1,16
187      a=ix0+MOD((i-1),2)   ; g4d(i,1) = set%x(a)
188      b=iy0+MOD((i-1)/2,2) ; g4d(i,2) = set%y(b)
189      c=iz0+MOD((i-1)/4,2) ; g4d(i,3) = set%z(c)
190      d=it0+MOD((i-1)/8,2) ; g4d(i,4) = set%t(d)
191      g4d(i,5) = set%data(a,b,c,d)
192    ENDDO
193    res = lintc_((/x,y,z,t/),g4d)
194    ret = .true.
195  END FUNCTION l4d_dset_sc
196
197  FUNCTION l5d_dset_sc(x,y,z,t,w,set,locator,res) RESULT(ret)
198    !! Wrapper interface for 5D linear interpolation (scalar)
199    !!
200    !! @warning
201    !! On error, __res__ output value is undefined.
202    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
203    REAL(kind=wp), INTENT(in)  :: y       !! Y coordinate of the point to interpolate
204    REAL(kind=wp), INTENT(in)  :: z       !! Z coordinate of the point to interpolate
205    REAL(kind=wp), INTENT(in)  :: t       !! T coordinate of the point to interpolate
206    REAL(kind=wp), INTENT(in)  :: w       !! W coordinate of the point to interpolate
207    TYPE(dset5d), INTENT(in)   :: set     !! Dataset with the tabulated function
208    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
209    PROCEDURE(locate)          :: locator !! Locator function
210    LOGICAL :: ret                        !! .true. on success, .false. otherwise
211    REAL(kind=wp), DIMENSION(32,6) :: g5d
212    INTEGER                        :: i,ix0,iy0,iz0,it0,iw0,a,b,c,d,e
213    ret = .false.
214    res = HUGE(res)
215    ix0 = locator(x,set%x) ; IF (ix0 == 0) RETURN
216    iy0 = locator(y,set%y) ; IF (iy0 == 0) RETURN
217    iz0 = locator(z,set%z) ; IF (iz0 == 0) RETURN
218    it0 = locator(t,set%t) ; IF (it0 == 0) RETURN
219    iw0 = locator(w,set%w) ; IF (iw0 == 0) RETURN
220    DO i=1,32
221      a=ix0+MOD((i-1),2)    ; g5d(i,1) = set%x(a)
222      b=iy0+MOD((i-1)/2,2)  ; g5d(i,2) = set%y(b)
223      c=iz0+MOD((i-1)/4,2)  ; g5d(i,3) = set%z(c)
224      d=it0+MOD((i-1)/8,2)  ; g5d(i,4) = set%t(d)
225      e=iw0+MOD((i-1)/16,2) ; g5d(i,5) = set%w(e)
226      g5d(i,6) = set%data(a,b,c,d,e)
227    ENDDO
228    res = lintc_((/x,y,z,t,w/),g5d)
229    ret = .true.
230  END FUNCTION l5d_dset_sc
231
232  FUNCTION l1d_dset_ve(x,set,locator,res) RESULT(ret)
233    !! Wrapper interface for 1D linear interpolation (vector)
234    !!
235    !! @warning
236    !! On error, __res__ output vector is undefined (i.e. not allocated).
237    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
238    TYPE(dset1d), INTENT(in)                              :: set     !! Dataset with the tabulated function
239    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
240    PROCEDURE(locate)                                     :: locator !! Locator function
241    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
242    REAL(kind=wp), DIMENSION(2,2) :: g1d
243    INTEGER                       :: nv,i,j,ix0
244    ret = .false.
245    nv = SIZE(x) ; ALLOCATE(res(nv))
246    DO j=1,nv
247      ix0 = locator(x(j),set%x)
248      IF (ix0 == 0) THEN
249        DEALLOCATE(res) ; RETURN
250      ENDIF
251      DO i=0,1
252        g1d(i+1,1) = set%x(ix0+i)
253        g1d(i+1,2) = set%data(ix0+i)
254      ENDDO
255      res(j) = lintc_((/x(j)/),g1d)
256    ENDDO
257    ret = .true.
258  END FUNCTION l1d_dset_ve
259
260  FUNCTION l2d_dset_ve(x,y,set,locator,res) RESULT(ret)
261    !! Wrapper interface for 2D linear interpolation (vector)
262    !!
263    !! @warning
264    !! On error, __res__ output vector is undefined (i.e. not allocated).
265    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
266    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: y       !! Y coordinate of the points to interpolate
267    TYPE(dset2d), INTENT(in)                              :: set     !! Dataset with the tabulated function
268    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
269    PROCEDURE(locate)                                     :: locator !! Locator function
270    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
271    REAL(kind=wp), DIMENSION(4,3) :: g2d
272    INTEGER                       :: nv,i,j,ix0,iy0,a,b
273    ret = .false.
274    nv = SIZE(x) ; ALLOCATE(res(nv))
275    DO j=1,nv
276      ix0 = locator(x(j),set%x)
277      iy0 = locator(y(j),set%y)
278      IF (ix0 == 0 .OR. iy0 == 0) THEN
279        DEALLOCATE(res) ; RETURN
280      ENDIF
281      DO i=1,4
282        a=ix0+MOD((i-1),2)   ; g2d(i,1) = set%x(a)
283        b=iy0+MOD((i-1)/2,2) ; g2d(i,2) = set%y(b)
284        g2d(i,3) = set%data(a,b)
285      ENDDO
286      res(j) = lintc_((/x(j),y(j)/),g2d)
287    ENDDO
288    ret = .true.
289  END FUNCTION l2d_dset_ve
290
291  FUNCTION l3d_dset_ve(x,y,z,set,locator,res) RESULT(ret)
292    !! Wrapper interface for 3D linear interpolation (vector)
293    !!
294    !! @warning
295    !! On error, __res__ output vector is undefined (i.e. not allocated).
296    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
297    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: y       !! Y coordinate of the points to interpolate
298    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: z       !! Z coordinate of the points to interpolate
299    TYPE(dset3d), INTENT(in)                              :: set     !! Dataset with the tabulated function
300    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
301    PROCEDURE(locate)                                     :: locator !! Locator function
302    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
303    REAL(kind=wp), DIMENSION(8,4) :: g3d
304    INTEGER                       :: nv,i,j,ix0,iy0,iz0,a,b,c
305    ret = .false.
306    nv = SIZE(x) ; ALLOCATE(res(nv))
307    DO j=1,nv
308      ix0 = locator(x(j),set%x)
309      iy0 = locator(y(j),set%y)
310      iz0 = locator(z(j),set%z)
311      IF (ix0==0 .OR. iy0==0 .OR. iz0==0) THEN
312        DEALLOCATE(res) ; RETURN
313      ENDIF
314      DO i=1,8
315        a=ix0+MOD((i-1),2)   ; g3d(i,1) = set%x(a)
316        b=iy0+MOD((i-1)/2,2) ; g3d(i,2) = set%y(b)
317        c=iz0+MOD((i-1)/4,2) ; g3d(i,3) = set%z(c)
318        g3d(i,4) = set%data(a,b,c)
319      ENDDO
320      res(j) = lintc_((/x(j),y(j),z(j)/),g3d)
321    ENDDO
322    ret = .true.
323  END FUNCTION l3d_dset_ve
324
325  FUNCTION l4d_dset_ve(x,y,z,t,set,locator,res) RESULT(ret)
326    !! Wrapper interface for 4D linear interpolation (vector)
327    !!
328    !! @warning
329    !! On error, __res__ output vector is undefined (i.e. not allocated).
330    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
331    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: y       !! Y coordinate of the points to interpolate
332    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: z       !! Z coordinate of the points to interpolate
333    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: t       !! T coordinate of the points to interpolate
334    TYPE(dset4d), INTENT(in)                              :: set     !! Dataset with the tabulated function
335    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
336    PROCEDURE(locate)                                     :: locator !! Locator function
337    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
338    REAL(kind=wp), DIMENSION(16,5) :: g4d
339    INTEGER                        :: nv,i,j,ix0,iy0,iz0,it0,a,b,c,d
340    ret = .false.
341    nv = SIZE(x) ; ALLOCATE(res(nv))
342    DO j=1,nv
343      ix0 = locator(x(j),set%x)
344      iy0 = locator(y(j),set%y)
345      iz0 = locator(z(j),set%z)
346      it0 = locator(t(j),set%t)
347      IF (ix0==0 .OR. iy0==0 .OR. iz0==0 .OR. it0==0) THEN
348        DEALLOCATE(res) ; RETURN
349      ENDIF
350      DO i=1,16
351        a=ix0+MOD((i-1),2)   ; g4d(i,1) = set%x(a)
352        b=iy0+MOD((i-1)/2,2) ; g4d(i,2) = set%y(b)
353        c=iz0+MOD((i-1)/4,2) ; g4d(i,3) = set%z(c)
354        d=it0+MOD((i-1)/8,2) ; g4d(i,4) = set%t(d)
355        g4d(i,5) = set%data(a,b,c,d)
356      ENDDO
357      res(j) = lintc_((/x(j),y(j),z(j),t(j)/),g4d)
358    ENDDO
359    ret = .true.
360  END FUNCTION l4d_dset_ve
361
362  FUNCTION l5d_dset_ve(x,y,z,t,w,set,locator,res) RESULT(ret)
363    !! Wrapper interface for 5D linear interpolation (vector)
364    !!
365    !! @warning
366    !! On error, __res__ output vector is undefined (i.e. not allocated).
367    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
368    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: y       !! Y coordinate of the points to interpolate
369    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: z       !! Z coordinate of the points to interpolate
370    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: t       !! T coordinate of the points to interpolate
371    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: w       !! W coordinate of the points to interpolate
372    TYPE(dset5d), INTENT(in)                              :: set     !! Dataset with the tabulated function
373    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
374    PROCEDURE(locate)                                     :: locator !! Locator function
375    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
376    REAL(kind=wp), DIMENSION(32,6) :: g5d
377    INTEGER                        :: nv,i,j,ix0,iy0,iz0,it0,iw0,a,b,c,d,e
378    ret = .false.
379    nv = SIZE(x) ; ALLOCATE(res(nv))
380    DO j=1,nv
381      ix0 = locator(x(j),set%x)
382      iy0 = locator(y(j),set%y)
383      iz0 = locator(z(j),set%z)
384      it0 = locator(t(j),set%t)
385      iw0 = locator(w(j),set%w)
386      IF (ix0==0 .OR. iy0==0 .OR. iz0==0 .OR. it0==0 .OR. iw0==0) THEN
387        DEALLOCATE(res) ; RETURN
388      ENDIF
389      DO i=1,32
390        a=ix0+MOD((i-1),2)    ; g5d(i,1) = set%x(a)
391        b=iy0+MOD((i-1)/2,2)  ; g5d(i,2) = set%y(b)
392        c=iz0+MOD((i-1)/4,2)  ; g5d(i,3) = set%z(c)
393        d=it0+MOD((i-1)/8,2)  ; g5d(i,4) = set%t(d)
394        e=iw0+MOD((i-1)/16,2) ; g5d(i,5) = set%w(e)
395        g5d(i,6) = set%data(a,b,c,d,e)
396      ENDDO
397      res(j) = lintc_((/x(j),y(j),z(j),t(j),w(j)/),g5d)
398    ENDDO
399    ret = .true.
400  END FUNCTION l5d_dset_ve
401
402  !--------------------!
403  ! HARD CODED VERSION !
404  !--------------------!
405
406  FUNCTION hdcdl1d_dset_sc(x,set, locator,res) RESULT(ret)
407    !! Hard-coded for 1D linear interpolation (scalar)
408    !!
409    !! @warning
410    !! On error, __res__ output value is undefined.
411    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
412    TYPE(dset1d), INTENT(in)   :: set     !! Dataset with the tabulated function
413    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
414    PROCEDURE(locate)          :: locator !! Locator function
415    LOGICAL :: ret                        !! .true. on success, .false. otherwise
416    REAL(kind=wp) :: xd
417    INTEGER       :: ix0,ix1
418    ret = .false.
419    res = HUGE(res)
420    ix0 = locator(x,set%x) ; IF (ix0 == 0) RETURN
421    ix1 = ix0+1
422    xd = (x-set%x(ix0))/(set%x(ix1)-set%x(ix0))
423    res = set%data(ix0)*(1d0-xd)+set%data(ix1)*xd
424    ret = .true.
425  END FUNCTION hdcdl1d_dset_sc
426
427  FUNCTION hdcdl2d_dset_sc(x,y,set,locator,res) RESULT(ret)
428    !! Hard-coded for 2D linear interpolation (scalar)
429    !!
430    !! @warning
431    !! On error, __res__ output value is undefined.
432    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
433    REAL(kind=wp), INTENT(in)  :: y       !! Y coordinate of the point to interpolate
434    TYPE(dset2d), INTENT(in)   :: set     !! Dataset with the tabulated function
435    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
436    PROCEDURE(locate)          :: locator !! Locator function
437    LOGICAL :: ret                        !! .true. on success, .false. otherwise
438    REAL(kind=wp) :: xd,yd
439    REAL(kind=wp) :: f0,f1
440    INTEGER       :: ix0,iy0,ix1,iy1
441    ret = .false.
442    res = HUGE(res)
443    ix0 = locator(x,set%x) ; IF (ix0 == 0) RETURN
444    iy0 = locator(y,set%y) ; IF (iy0 == 0) RETURN
445    ix1 = ix0 + 1 ; iy1 = iy0 + 1
446    xd  = (x-set%x(ix0))/(set%x(ix1)-set%x(ix0))
447    yd  = (y-set%y(iy0))/(set%y(iy1)-set%y(iy0))
448    f0  = set%data(ix0,iy0)*(1d0-xd)+set%data(ix1,iy0)*xd
449    f1  = set%data(ix0,iy1)*(1d0-xd)+set%data(ix1,iy1)*xd
450    res = f0*(1d0-yd)+f1*yd
451    ret = .true.
452  END FUNCTION hdcdl2d_dset_sc
453
454  FUNCTION hdcdl3d_dset_sc(x,y,z,set,locator,res) RESULT(ret)
455    !! Hard-coded for 3D linear interpolation (scalar)
456    !!
457    !! @warning
458    !! On error, __res__ output value is undefined.
459    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
460    REAL(kind=wp), INTENT(in)  :: y       !! Y coordinate of the point to interpolate
461    REAL(kind=wp), INTENT(in)  :: z       !! Z coordinate of the point to interpolate
462    TYPE(dset3d), INTENT(in)   :: set     !! Dataset with the tabulated function
463    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
464    PROCEDURE(locate)          :: locator !! Locator function
465    LOGICAL :: ret                        !! .true. on success, .false. otherwise
466    REAL(kind=wp) :: xd,yd,zd
467    REAL(kind=wp) :: f00,f10,f01,f11,f0,f1
468    INTEGER       :: ix0,iy0,iz0,ix1,iy1,iz1
469    ret = .false.
470    res = HUGE(res)
471    ix0 = locator(x,set%x) ; IF (ix0 == 0) RETURN
472    iy0 = locator(y,set%y) ; IF (iy0 == 0) RETURN
473    iz0 = locator(z,set%z) ; IF (iz0 == 0) RETURN
474    ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1
475    xd = (x-set%x(ix0))/(set%x(ix1)-set%x(ix0))
476    yd = (y-set%y(iy0))/(set%y(iy1)-set%y(iy0))
477    zd = (z-set%z(iz0))/(set%z(iz1)-set%z(iz0))
478
479    f00 = set%data(ix0,iy0,iz0)*(1d0-xd)+set%data(ix1,iy0,iz0)*xd
480    f10 = set%data(ix0,iy1,iz0)*(1d0-xd)+set%data(ix1,iy1,iz0)*xd
481    f01 = set%data(ix0,iy0,iz1)*(1d0-xd)+set%data(ix1,iy0,iz1)*xd
482    f11 = set%data(ix0,iy1,iz1)*(1d0-xd)+set%data(ix1,iy1,iz1)*xd
483
484    f0 = f00 *(1d0-yd)+f10*yd
485    f1 = f00 *(1d0-yd)+f10*yd
486
487    res = f0*(1d0-zd)+f1*zd
488    ret = .true.
489  END FUNCTION hdcdl3d_dset_sc
490
491  FUNCTION hdcdl4d_dset_sc(x,y,z,t,set,locator,res) RESULT(ret)
492    !! Hard-coded for 4D linear interpolation (scalar)
493    !!
494    !! @warning
495    !! On error, __res__ output value is undefined.
496    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
497    REAL(kind=wp), INTENT(in)  :: y       !! Y coordinate of the point to interpolate
498    REAL(kind=wp), INTENT(in)  :: z       !! Z coordinate of the point to interpolate
499    REAL(kind=wp), INTENT(in)  :: t       !! T coordinate of the point to interpolate
500    TYPE(dset4d), INTENT(in)   :: set     !! Dataset with the tabulated function
501    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
502    PROCEDURE(locate)          :: locator !! Locator function
503    LOGICAL :: ret                        !! .true. on success, .false. otherwise
504    REAL(kind=wp) :: xd,yd,zd,td
505    REAL(kind=wp) :: f000,f100,f010,f110,f001,f101,f011,f111, &
506                     f00,f10,f01,f11,f0,f1
507    INTEGER       :: ix0,iy0,iz0,it0,ix1,iy1,iz1,it1
508    ret = .false.
509    res = HUGE(res)
510    ix0 = locator(x,set%x) ; IF (ix0 == 0) RETURN
511    iy0 = locator(y,set%y) ; IF (iy0 == 0) RETURN
512    iz0 = locator(z,set%z) ; IF (iz0 == 0) RETURN
513    it0 = locator(t,set%t) ; IF (it0 == 0) RETURN
514   
515    ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ; it1=it0+1
516    xd = (x-set%x(ix0))/(set%x(ix1)-set%x(ix0))
517    yd = (y-set%y(iy0))/(set%y(iy1)-set%y(iy0))
518    zd = (z-set%z(iz0))/(set%z(iz1)-set%z(iz0))
519    td = (t-set%t(it0))/(set%t(it1)-set%t(it0))
520
521    f000 = set%data(ix0,iy0,iz0,it0)*(1d0-xd)+set%data(ix1,iy0,iz0,it0)*xd
522    f100 = set%data(ix0,iy1,iz0,it0)*(1d0-xd)+set%data(ix1,iy1,iz0,it0)*xd
523    f010 = set%data(ix0,iy0,iz1,it0)*(1d0-xd)+set%data(ix1,iy0,iz1,it0)*xd
524    f110 = set%data(ix0,iy1,iz1,it0)*(1d0-xd)+set%data(ix1,iy1,iz1,it0)*xd
525    f001 = set%data(ix0,iy0,iz0,it1)*(1d0-xd)+set%data(ix1,iy0,iz0,it1)*xd
526    f101 = set%data(ix0,iy1,iz0,it1)*(1d0-xd)+set%data(ix1,iy1,iz0,it1)*xd
527    f011 = set%data(ix0,iy0,iz1,it1)*(1d0-xd)+set%data(ix1,iy0,iz1,it1)*xd
528    f111 = set%data(ix0,iy1,iz1,it1)*(1d0-xd)+set%data(ix1,iy1,iz1,it1)*xd
529
530    f00 = f000*(1d0-yd)+f100*yd
531    f10 = f010*(1d0-yd)+f110*yd
532    f01 = f001*(1d0-yd)+f101*yd
533    f11 = f011*(1d0-yd)+f111*yd
534
535    f0 = f00 *(1d0-zd)+f10*zd
536    f1 = f01 *(1d0-zd)+f11*zd
537
538    res = f0*(1d0-td)+f1*td
539    ret = .true.
540  END FUNCTION hdcdl4d_dset_sc
541
542  FUNCTION hdcdl5d_dset_sc(x,y,z,t,w,set,locator,res) RESULT(ret)
543    !! Hard-coded for 5D linear interpolation (scalar)
544    !!
545    !! @warning
546    !! On error, __res__ output value is undefined.
547    REAL(kind=wp), INTENT(in)  :: x       !! X coordinate of the point to interpolate
548    REAL(kind=wp), INTENT(in)  :: y       !! Y coordinate of the point to interpolate
549    REAL(kind=wp), INTENT(in)  :: z       !! Z coordinate of the point to interpolate
550    REAL(kind=wp), INTENT(in)  :: t       !! T coordinate of the point to interpolate
551    REAL(kind=wp), INTENT(in)  :: w       !! W coordinate of the point to interpolate
552    TYPE(dset5d), INTENT(in)   :: set     !! Dataset with the tabulated function
553    REAL(kind=wp), INTENT(out) :: res     !! Interpolated value
554    PROCEDURE(locate)          :: locator !! Locator function
555    LOGICAL :: ret                        !! .true. on success, .false. otherwise
556    REAL(kind=wp) :: xd,yd,zd,td,wd
557    REAL(kind=wp) :: f0000,f1000,f0100,f1100,f0010,f1010,f0110,f1110, &
558                     f0001,f1001,f0101,f1101,f0011,f1011,f0111,f1111, &
559                     f000,f100,f010,f110,f001,f101,f011,f111,         &
560                     f00,f10,f01,f11,f0,f1
561    INTEGER       :: ix0,iy0,iz0,it0,iw0,ix1,iy1,iz1,it1,iw1
562    ret = .false.
563    res = HUGE(res)
564    ix0 = locator(x,set%x) ; IF (ix0 == 0) RETURN
565    iy0 = locator(y,set%y) ; IF (iy0 == 0) RETURN
566    iz0 = locator(z,set%z) ; IF (iz0 == 0) RETURN
567    it0 = locator(t,set%t) ; IF (it0 == 0) RETURN
568    iw0 = locator(w,set%w) ; IF (iw0 == 0) RETURN
569   
570    ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1
571    xd = (x-set%x(ix0))/(set%x(ix1)-set%x(ix0))
572    yd = (y-set%y(iy0))/(set%y(iy1)-set%y(iy0))
573    zd = (z-set%z(iz0))/(set%z(iz1)-set%z(iz0))
574    td = (t-set%t(it0))/(set%t(it1)-set%t(it0))
575    wd = (w-set%w(it0))/(set%w(it1)-set%w(it0))
576
577    f0000 = set%data(ix0,iy0,iz0,it0,iw0)*(1d0-xd)+set%data(ix1,iy0,iz0,it0,iw0)*xd
578    f1000 = set%data(ix0,iy1,iz0,it0,iw0)*(1d0-xd)+set%data(ix0,iy1,iz0,it0,iw0)*xd
579    f0100 = set%data(ix0,iy0,iz1,it0,iw0)*(1d0-xd)+set%data(ix1,iy0,iz1,it0,iw0)*xd
580    f1100 = set%data(ix0,iy1,iz1,it0,iw0)*(1d0-xd)+set%data(ix1,iy1,iz1,it0,iw0)*xd
581    f0010 = set%data(ix0,iy0,iz0,it1,iw0)*(1d0-xd)+set%data(ix1,iy0,iz0,it1,iw0)*xd
582    f1010 = set%data(ix0,iy1,iz0,it1,iw0)*(1d0-xd)+set%data(ix1,iy1,iz0,it1,iw0)*xd
583    f0110 = set%data(ix0,iy0,iz1,it1,iw0)*(1d0-xd)+set%data(ix1,iy0,iz1,it1,iw0)*xd
584    f1110 = set%data(ix0,iy1,iz1,it1,iw0)*(1d0-xd)+set%data(ix1,iy1,iz1,it1,iw0)*xd
585    f0001 = set%data(ix0,iy0,iz0,it0,iw1)*(1d0-xd)+set%data(ix1,iy0,iz0,it0,iw1)*xd
586    f1001 = set%data(ix0,iy1,iz0,it0,iw1)*(1d0-xd)+set%data(ix0,iy1,iz0,it0,iw1)*xd
587    f0101 = set%data(ix0,iy0,iz1,it0,iw1)*(1d0-xd)+set%data(ix1,iy0,iz1,it0,iw1)*xd
588    f1101 = set%data(ix0,iy1,iz1,it0,iw1)*(1d0-xd)+set%data(ix1,iy1,iz1,it0,iw1)*xd
589    f0011 = set%data(ix0,iy0,iz0,it1,iw1)*(1d0-xd)+set%data(ix1,iy0,iz0,it1,iw1)*xd
590    f1011 = set%data(ix0,iy1,iz0,it1,iw1)*(1d0-xd)+set%data(ix1,iy1,iz0,it1,iw1)*xd
591    f0111 = set%data(ix0,iy0,iz1,it1,iw1)*(1d0-xd)+set%data(ix1,iy0,iz1,it1,iw1)*xd
592    f1111 = set%data(ix0,iy1,iz1,it1,iw1)*(1d0-xd)+set%data(ix1,iy1,iz1,it1,iw1)*xd
593
594    f000 = f0000*(1d0-yd) + f1000*yd
595    f100 = f0100*(1d0-yd) + f1100*yd
596    f010 = f0010*(1d0-yd) + f1010*yd
597    f110 = f0110*(1d0-yd) + f1110*yd
598    f101 = f0101*(1d0-yd) + f1101*yd
599    f011 = f0011*(1d0-yd) + f1011*yd
600    f111 = f0111*(1d0-yd) + f1111*yd
601
602    f00 = f000*(1d0-zd)+f100*zd
603    f10 = f010*(1d0-zd)+f110*zd
604    f01 = f001*(1d0-zd)+f101*zd
605    f11 = f011*(1d0-zd)+f111*zd
606
607    f0 = f00 *(1d0-td)+f10*td
608    f1 = f01 *(1d0-td)+f11*td
609
610    res = f0*(1d0-wd)+f1*wd
611    ret = .true.
612  END FUNCTION hdcdl5d_dset_sc
613
614  FUNCTION hdcdl1d_dset_ve(x,set,locator,res) RESULT(ret)
615    !! Hard-coded for 1D linear interpolation (vector)
616    !!
617    !! @warning
618    !! On error, __res__ output vector is undefined (i.e. not allocated).
619    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
620    TYPE(dset1d), INTENT(in)                              :: set     !! Dataset with the tabulated function
621    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
622    PROCEDURE(locate)                                     :: locator !! Locator function
623    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
624    REAL(kind=wp) :: xd
625    INTEGER       :: nv,j,ix0,ix1
626    ret = .false.
627    nv = SIZE(x) ; ALLOCATE(res(nv))
628    DO j=1,nv
629      ix0 = locator(x(j),set%x)
630      IF (ix0==0) THEN
631        DEALLOCATE(res) ; RETURN
632      ENDIF
633      ix1 = ix0+1
634      xd = (x(j)-set%x(ix0))/(set%x(ix1)-set%x(ix0))
635      res(j) = set%data(ix0)*(1d0-xd)+set%data(ix1)*xd
636    ENDDO
637    ret = .true.
638  END FUNCTION hdcdl1d_dset_ve
639
640  FUNCTION hdcdl2d_dset_ve(x,y,set,locator,res) RESULT(ret)
641    !! Hard-coded for 2D linear interpolation (vector)
642    !!
643    !! @warning
644    !! On error, __res__ output vector is undefined (i.e. not allocated).
645    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
646    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: y       !! Y coordinate of the points to interpolate
647    TYPE(dset2d), INTENT(in)                              :: set     !! Dataset with the tabulated function
648    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
649    PROCEDURE(locate)                                     :: locator !! Locator function
650    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
651    REAL(kind=wp) :: xd,yd
652    REAL(kind=wp) :: f0,f1
653    INTEGER      :: nv,j,ix0,iy0,ix1,iy1
654    ret = .false.
655    nv = SIZE(x) ; ALLOCATE(res(nv))
656    DO j=1,nv
657      ix0 = locator(x(j),set%x)
658      iy0 = locator(y(j),set%y)
659      IF (ix0==0 .OR. iy0==0) THEN
660        DEALLOCATE(res) ; RETURN
661      ENDIF
662      ix1 = ix0 + 1 ; iy1 = iy0 + 1
663      xd  = (x(j)-set%x(ix0))/(set%x(ix1)-set%x(ix0))
664      yd  = (y(j)-set%y(iy0))/(set%y(iy1)-set%y(iy0))
665      f0  = set%data(ix0,iy0)*(1d0-xd)+set%data(ix1,iy0)*xd
666      f1  = set%data(ix0,iy1)*(1d0-xd)+set%data(ix1,iy1)*xd
667      res(j) = f0*(1d0-yd)+f1*yd
668    ENDDO
669    ret = .true.
670  END FUNCTION hdcdl2d_dset_ve
671
672  FUNCTION hdcdl3d_dset_ve(x,y,z,set,locator,res) RESULT(ret)
673    !! Hard-coded for 3D linear interpolation (vector)
674    !!
675    !! @warning
676    !! On error, __res__ output vector is undefined (i.e. not allocated).
677    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
678    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: y       !! Y coordinate of the points to interpolate
679    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: z       !! Z coordinate of the points to interpolate
680    TYPE(dset3d), INTENT(in)                              :: set     !! Dataset with the tabulated function
681    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
682    PROCEDURE(locate)                                     :: locator !! Locator function
683    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
684    REAL(kind=wp) :: xd,yd,zd
685    REAL(kind=wp) :: f00,f10,f01,f11,f0,f1
686    INTEGER       :: nv,j,ix0,iy0,iz0,ix1,iy1,iz1
687    ret = .false.
688    nv = SIZE(x) ; ALLOCATE(res(nv))
689    DO j=1,nv
690      ix0 = locator(x(j),set%x)
691      iy0 = locator(y(j),set%y)
692      iz0 = locator(z(j),set%z)
693
694      IF (ix0==0 .OR. iy0==0 .OR. iz0==0) THEN
695        DEALLOCATE(res) ; RETURN
696      ENDIF
697
698      ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1
699      xd = (x(j)-set%x(ix0))/(set%x(ix1)-set%x(ix0))
700      yd = (y(j)-set%y(iy0))/(set%y(iy1)-set%y(iy0))
701      zd = (z(j)-set%z(iz0))/(set%z(iz1)-set%z(iz0))
702
703      f00 = set%data(ix0,iy0,iz0)*(1d0-xd)+set%data(ix1,iy0,iz0)*xd
704      f10 = set%data(ix0,iy1,iz0)*(1d0-xd)+set%data(ix1,iy1,iz0)*xd
705      f01 = set%data(ix0,iy0,iz1)*(1d0-xd)+set%data(ix1,iy0,iz1)*xd
706      f11 = set%data(ix0,iy1,iz1)*(1d0-xd)+set%data(ix1,iy1,iz1)*xd
707
708      f0 = f00 *(1d0-yd)+f10*yd
709      f1 = f00 *(1d0-yd)+f10*yd
710
711      res(j) = f0*(1d0-zd)+f1*zd
712    ENDDO
713    ret = .true.
714  END FUNCTION hdcdl3d_dset_ve
715
716  FUNCTION hdcdl4d_dset_ve(x,y,z,t,set,locator,res) RESULT(ret)
717    !! Hard-coded for 4D linear interpolation (vector)
718    !!
719    !! @warning
720    !! On error, __res__ output vector is undefined (i.e. not allocated).
721    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
722    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: y       !! Y coordinate of the points to interpolate
723    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: z       !! Z coordinate of the points to interpolate
724    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: t       !! T coordinate of the points to interpolate
725    TYPE(dset4d), INTENT(in)                              :: set     !! Dataset with the tabulated function
726    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
727    PROCEDURE(locate)                                     :: locator !! Locator function
728    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
729    REAL(kind=wp) :: xd,yd,zd,td
730    REAL(kind=wp) :: f000,f100,f010,f110,f001,f101,f011,f111, &
731                     f00,f10,f01,f11,f0,f1
732    INTEGER       :: nv,j,ix0,iy0,iz0,it0,ix1,iy1,iz1,it1
733    ret = .false.
734    nv = SIZE(x) ; ALLOCATE(res(nv))
735    DO j=1,nv
736      ix0 = locator(x(j),set%x)
737      iy0 = locator(y(j),set%y)
738      iz0 = locator(z(j),set%z)
739      it0 = locator(t(j),set%t)
740   
741      IF (ix0==0 .OR. iy0==0 .OR. iz0==0 .OR. it0==0) THEN
742        DEALLOCATE(res) ; RETURN
743      ENDIF
744
745      ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ; it1=it0+1
746      xd = (x(j)-set%x(ix0))/(set%x(ix1)-set%x(ix0))
747      yd = (y(j)-set%y(iy0))/(set%y(iy1)-set%y(iy0))
748      zd = (z(j)-set%z(iz0))/(set%z(iz1)-set%z(iz0))
749      td = (t(j)-set%t(it0))/(set%t(it1)-set%t(it0))
750
751      f000 = set%data(ix0,iy0,iz0,it0)*(1d0-xd)+set%data(ix1,iy0,iz0,it0)*xd
752      f100 = set%data(ix0,iy1,iz0,it0)*(1d0-xd)+set%data(ix1,iy1,iz0,it0)*xd
753      f010 = set%data(ix0,iy0,iz1,it0)*(1d0-xd)+set%data(ix1,iy0,iz1,it0)*xd
754      f110 = set%data(ix0,iy1,iz1,it0)*(1d0-xd)+set%data(ix1,iy1,iz1,it0)*xd
755      f001 = set%data(ix0,iy0,iz0,it1)*(1d0-xd)+set%data(ix1,iy0,iz0,it1)*xd
756      f101 = set%data(ix0,iy1,iz0,it1)*(1d0-xd)+set%data(ix1,iy1,iz0,it1)*xd
757      f011 = set%data(ix0,iy0,iz1,it1)*(1d0-xd)+set%data(ix1,iy0,iz1,it1)*xd
758      f111 = set%data(ix0,iy1,iz1,it1)*(1d0-xd)+set%data(ix1,iy1,iz1,it1)*xd
759
760      f00 = f000*(1d0-yd)+f100*yd
761      f10 = f010*(1d0-yd)+f110*yd
762      f01 = f001*(1d0-yd)+f101*yd
763      f11 = f011*(1d0-yd)+f111*yd
764
765      f0 = f00 *(1d0-zd)+f10*zd
766      f1 = f01 *(1d0-zd)+f11*zd
767
768      res(j) = f0*(1d0-td)+f1*td
769    ENDDO
770    ret = .true.
771  END FUNCTION hdcdl4d_dset_ve
772
773  FUNCTION hdcdl5d_dset_ve(x,y,z,t,w,set,locator,res) RESULT(ret)
774    !! Hard-coded for 5D linear interpolation (vector)
775    !!
776    !! @warning
777    !! On error, __res__ output vector is undefined (i.e. not allocated).
778    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: x       !! X coordinate of the points to interpolate
779    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: y       !! Y coordinate of the points to interpolate
780    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: z       !! Z coordinate of the points to interpolate
781    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: t       !! T coordinate of the points to interpolate
782    REAL(kind=wp), INTENT(in), DIMENSION(:)               :: w       !! W coordinate of the points to interpolate
783    TYPE(dset5d), INTENT(in)                              :: set     !! Dataset with the tabulated function
784    REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res     !! Interpolated values
785    PROCEDURE(locate)                                     :: locator !! Locator function
786    LOGICAL :: ret                                                   !! .true. on success, .false. otherwise
787    REAL(kind=wp) :: xd,yd,zd,td,wd
788    REAL(kind=wp) :: f0000,f1000,f0100,f1100,f0010,f1010,f0110,f1110, &
789                     f0001,f1001,f0101,f1101,f0011,f1011,f0111,f1111, &
790                     f000,f100,f010,f110,f001,f101,f011,f111,         &
791                     f00,f10,f01,f11,f0,f1
792    INTEGER       :: nv,j,ix0,iy0,iz0,it0,iw0,ix1,iy1,iz1,it1,iw1
793    ret = .false.
794    nv = SIZE(x) ; ALLOCATE(res(nv))
795    DO j=1,nv
796      ix0 = locator(x(j),set%x)
797      iy0 = locator(y(j),set%y)
798      iz0 = locator(z(j),set%z)
799      it0 = locator(t(j),set%t)
800      iw0 = locator(w(j),set%w)
801
802      IF (ix0==0 .OR. iy0==0 .OR. iz0==0 .OR. it0==0 .OR. iw0==0) THEN
803        DEALLOCATE(res) ; RETURN
804      ENDIF
805   
806      ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1
807      xd = (x(j)-set%x(ix0))/(set%x(ix1)-set%x(ix0))
808      yd = (y(j)-set%y(iy0))/(set%y(iy1)-set%y(iy0))
809      zd = (z(j)-set%z(iz0))/(set%z(iz1)-set%z(iz0))
810      td = (t(j)-set%t(it0))/(set%t(it1)-set%t(it0))
811      wd = (w(j)-set%w(it0))/(set%w(it1)-set%w(it0))
812
813      f0000 = set%data(ix0,iy0,iz0,it0,iw0)*(1d0-xd)+set%data(ix1,iy0,iz0,it0,iw0)*xd
814      f1000 = set%data(ix0,iy1,iz0,it0,iw0)*(1d0-xd)+set%data(ix0,iy1,iz0,it0,iw0)*xd
815      f0100 = set%data(ix0,iy0,iz1,it0,iw0)*(1d0-xd)+set%data(ix1,iy0,iz1,it0,iw0)*xd
816      f1100 = set%data(ix0,iy1,iz1,it0,iw0)*(1d0-xd)+set%data(ix1,iy1,iz1,it0,iw0)*xd
817      f0010 = set%data(ix0,iy0,iz0,it1,iw0)*(1d0-xd)+set%data(ix1,iy0,iz0,it1,iw0)*xd
818      f1010 = set%data(ix0,iy1,iz0,it1,iw0)*(1d0-xd)+set%data(ix1,iy1,iz0,it1,iw0)*xd
819      f0110 = set%data(ix0,iy0,iz1,it1,iw0)*(1d0-xd)+set%data(ix1,iy0,iz1,it1,iw0)*xd
820      f1110 = set%data(ix0,iy1,iz1,it1,iw0)*(1d0-xd)+set%data(ix1,iy1,iz1,it1,iw0)*xd
821      f0001 = set%data(ix0,iy0,iz0,it0,iw1)*(1d0-xd)+set%data(ix1,iy0,iz0,it0,iw1)*xd
822      f1001 = set%data(ix0,iy1,iz0,it0,iw1)*(1d0-xd)+set%data(ix0,iy1,iz0,it0,iw1)*xd
823      f0101 = set%data(ix0,iy0,iz1,it0,iw1)*(1d0-xd)+set%data(ix1,iy0,iz1,it0,iw1)*xd
824      f1101 = set%data(ix0,iy1,iz1,it0,iw1)*(1d0-xd)+set%data(ix1,iy1,iz1,it0,iw1)*xd
825      f0011 = set%data(ix0,iy0,iz0,it1,iw1)*(1d0-xd)+set%data(ix1,iy0,iz0,it1,iw1)*xd
826      f1011 = set%data(ix0,iy1,iz0,it1,iw1)*(1d0-xd)+set%data(ix1,iy1,iz0,it1,iw1)*xd
827      f0111 = set%data(ix0,iy0,iz1,it1,iw1)*(1d0-xd)+set%data(ix1,iy0,iz1,it1,iw1)*xd
828      f1111 = set%data(ix0,iy1,iz1,it1,iw1)*(1d0-xd)+set%data(ix1,iy1,iz1,it1,iw1)*xd
829
830      f000 = f0000*(1d0-yd) + f1000*yd
831      f100 = f0100*(1d0-yd) + f1100*yd
832      f010 = f0010*(1d0-yd) + f1010*yd
833      f110 = f0110*(1d0-yd) + f1110*yd
834      f101 = f0101*(1d0-yd) + f1101*yd
835      f011 = f0011*(1d0-yd) + f1011*yd
836      f111 = f0111*(1d0-yd) + f1111*yd
837
838      f00 = f000*(1d0-zd)+f100*zd
839      f10 = f010*(1d0-zd)+f110*zd
840      f01 = f001*(1d0-zd)+f101*zd
841      f11 = f011*(1d0-zd)+f111*zd
842
843      f0 = f00 *(1d0-td)+f10*td
844      f1 = f01 *(1d0-td)+f11*td
845
846      res(j) = f0*(1d0-wd)+f1*wd
847    ENDDO
848    ret = .true.
849  END FUNCTION hdcdl5d_dset_ve
850
851END MODULE LINTDSET
852
Note: See TracBrowser for help on using the repository browser.