source: trunk/LMDZ.PLUTO/libf/muphypluto/lint_dset.F90 @ 3607

Last change on this file since 3607 was 3560, checked in by debatzbr, 8 weeks ago

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