Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lintdset.f90 @ 3583

Last change on this file since 3583 was 1793, checked in by jvatant, 8 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

Line
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
41	MODULE 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)+f1yd
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)+f10yd
485	f1 = f00 (1d0-yd)+f10yd
486
487	res = f0(1d0-zd)+f1zd
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)+f100yd
531	f10 = f010(1d0-yd)+f110yd
532	f01 = f001(1d0-yd)+f101yd
533	f11 = f011(1d0-yd)+f111yd
534
535	f0 = f00 (1d0-zd)+f10zd
536	f1 = f01 (1d0-zd)+f11zd
537
538	res = f0(1d0-td)+f1td
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) + f1000yd
595	f100 = f0100(1d0-yd) + f1100yd
596	f010 = f0010(1d0-yd) + f1010yd
597	f110 = f0110(1d0-yd) + f1110yd
598	f101 = f0101(1d0-yd) + f1101yd
599	f011 = f0011(1d0-yd) + f1011yd
600	f111 = f0111(1d0-yd) + f1111yd
601
602	f00 = f000(1d0-zd)+f100zd
603	f10 = f010(1d0-zd)+f110zd
604	f01 = f001(1d0-zd)+f101zd
605	f11 = f011(1d0-zd)+f111zd
606
607	f0 = f00 (1d0-td)+f10td
608	f1 = f01 (1d0-td)+f11td
609
610	res = f0(1d0-wd)+f1wd
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)+f1yd
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)+f10yd
709	f1 = f00 (1d0-yd)+f10yd
710
711	res(j) = f0(1d0-zd)+f1zd
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)+f100yd
761	f10 = f010(1d0-yd)+f110yd
762	f01 = f001(1d0-yd)+f101yd
763	f11 = f011(1d0-yd)+f111yd
764
765	f0 = f00 (1d0-zd)+f10zd
766	f1 = f01 (1d0-zd)+f11zd
767
768	res(j) = f0(1d0-td)+f1td
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) + f1000yd
831	f100 = f0100(1d0-yd) + f1100yd
832	f010 = f0010(1d0-yd) + f1010yd
833	f110 = f0110(1d0-yd) + f1110yd
834	f101 = f0101(1d0-yd) + f1101yd
835	f011 = f0011(1d0-yd) + f1011yd
836	f111 = f0111(1d0-yd) + f1111yd
837
838	f00 = f000(1d0-zd)+f100zd
839	f10 = f010(1d0-zd)+f110zd
840	f01 = f001(1d0-zd)+f101zd
841	f11 = f011(1d0-zd)+f111zd
842
843	f0 = f00 (1d0-td)+f10td
844	f1 = f01 (1d0-td)+f11td
845
846	res(j) = f0(1d0-wd)+f1wd
847	ENDDO
848	ret = .true.
849	END FUNCTION hdcdl5d_dset_ve
850
851	END MODULE LINTDSET
852

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format