Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lintgen.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: 44.4 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: lintgen.f90
36	!! summary: Linear interpolation generic interfaces
37	!! author: burgalat
38	!! date: 2010-2014
39
40	MODULE LINTGEN
41	!! Generic linear interpolation module
42	USE LINT_PREC
43	USE LINTCORE
44	IMPLICIT NONE
45
46	PUBLIC :: lint_gen, hdcd_lint_gen
47	PRIVATE
48
49	!> Interface to multi-variate linear interpolation
50	!!
51	!! For __n__, the number of function's variables, user must provide:
52	!!
53	!! - __n__ scalar(s) with the coordinate(s) of the point to interpolate.
54	!! - __n__ vector(s) with the tabulated values of each coordinate.
55	!! - a single __n__-D array with the value of tabulated function at each
56	!! tabulated coordinate.
57	!! - A function, satisfying [[lintcore(module):locate(interface)]] interface,
58	!! that locates the neighbourhood of the point to interpolate.
59	!! - A scalar that stores the output value.
60	INTERFACE lint_gen
61	MODULE PROCEDURE l1d_gen_sc,l2d_gen_sc,l3d_gen_sc,l4d_gen_sc,l5d_gen_sc
62	MODULE PROCEDURE l1d_gen_ve,l2d_gen_ve,l3d_gen_ve,l4d_gen_ve,l5d_gen_ve
63	END INTERFACE
64
65	!> Interface to multi-variate hard-coded linear interpolation
66	!!
67	!! Same remarks as in [[lintgen(module):lint_gen(interface)]] apply here.
68	INTERFACE hdcd_lint_gen
69	MODULE PROCEDURE hdcdl1d_gen_sc, hdcdl2d_gen_sc, hdcdl3d_gen_sc, &
70	hdcdl4d_gen_sc, hdcdl5d_gen_sc
71	MODULE PROCEDURE hdcdl1d_gen_ve, hdcdl2d_gen_ve, hdcdl3d_gen_ve, &
72	hdcdl4d_gen_ve, hdcdl5d_gen_ve
73	END INTERFACE
74
75	CONTAINS
76
77	! GENERIC versions
78	! ----------------
79
80	FUNCTION l1d_gen_sc(x,cx,cv, locator,res) RESULT(ret)
81	!! Wrapper interface for 1D linear interpolation (scalar)
82	!!
83	!! @warning
84	!! On error, __res__ output value is undefined.
85	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate
86	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
87	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cv !! Tabulated function values
88	REAL(kind=wp), INTENT(out) :: res !! Interpolated value
89	PROCEDURE(locate) :: locator !! Locator function
90	LOGICAL :: ret !! .true. on success, .false. otherwise
91	REAL(kind=wp), DIMENSION(2,2) :: g1d
92	INTEGER :: i,ix
93	ret = .false.
94	ix = locator(x,cx) ; IF (ix == 0) RETURN
95	DO i=0,1
96	g1d(i+1,1) = cx(ix+i)
97	g1d(i+1,2) = cv(ix+i)
98	ENDDO
99	res = lintc_((/x/),g1d)
100	ret = .true.
101	END FUNCTION l1d_gen_sc
102
103	FUNCTION l2d_gen_sc(x,y,cx,cy,cv,locator,res) RESULT(ret)
104	!! Wrapper interface for 2D linear interpolation (scalar)
105	!!
106	!! @warning
107	!! On error, __res__ output value is undefined.
108	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate
109	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate
110	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
111	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
112	REAL(kind=wp), INTENT(in), DIMENSION(:,:) :: cv !! Tabulated function values
113	REAL(kind=wp), INTENT(out) :: res !! Interpolated value
114	PROCEDURE(locate) :: locator !! Locator function
115	LOGICAL :: ret !! .true. on success, .false. otherwise
116	REAL(kind=wp), DIMENSION(4,3) :: g2d
117	INTEGER :: ix0,iy0,i,a,b
118	ret = .false.
119	ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN
120	iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN
121	DO i=1,4
122	a=ix0+MOD((i-1),2) ; g2d(i,1) = cx(a)
123	b=iy0+MOD((i-1)/2,2) ; g2d(i,2) = cy(b)
124	g2d(i,3) = cv(a,b)
125	ENDDO
126	res = lintc_((/x,y/),g2d)
127	ret = .true.
128	END FUNCTION l2d_gen_sc
129
130	FUNCTION l3d_gen_sc(x,y,z,cx,cy,cz,cv,locator,res) RESULT(ret)
131	!! Wrapper interface for 3D linear interpolation (scalar)
132	!!
133	!! @warning
134	!! On error, __res__ output value is undefined.
135	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate
136	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate
137	REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate
138	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
139	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
140	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
141	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:) :: cv !! Tabulated function values
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 :: ix0,iy0,iz0,i,a,b,c
147	ret = .false.
148	ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN
149	iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN
150	iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN
151	DO i=1,8
152	a=ix0+MOD((i-1),2) ; g3d(i,1) = cx(a)
153	b=iy0+MOD((i-1)/2,2) ; g3d(i,2) = cy(b)
154	c=iz0+MOD((i-1)/4,2) ; g3d(i,3) = cz(c)
155	g3d(i,4) = cv(a,b,c)
156	ENDDO
157	res = lintc_((/x,y,z/),g3d)
158	ret = .true.
159	END FUNCTION l3d_gen_sc
160
161	FUNCTION l4d_gen_sc(x,y,z,t,cx,cy,cz,ct,cv,locator,res) RESULT(ret)
162	!! Wrapper interface for 4D linear interpolation (scalar)
163	!!
164	!! @warning
165	!! On error, __res__ output value is undefined.
166	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate
167	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate
168	REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate
169	REAL(kind=wp), INTENT(in) :: t !! T coordinate of the point to interpolate
170	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
171	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
172	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
173	REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates
174	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:) :: cv !! Tabulated function values
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 :: ix0,iy0,iz0,it0
180	INTEGER :: i,a,b,c,d
181	ret = .false.
182	ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN
183	iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN
184	iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN
185	it0 = locator(t,ct) ; IF (it0 == 0) RETURN
186	DO i=1,16
187	a=ix0+MOD((i-1),2) ; g4d(i,1) = cx(a)
188	b=iy0+MOD((i-1)/2,2) ; g4d(i,2) = cy(b)
189	c=iz0+MOD((i-1)/4,2) ; g4d(i,3) = cz(c)
190	d=it0+MOD((i-1)/8,2) ; g4d(i,4) = ct(d)
191	g4d(i,5) = cv(a,b,c,d)
192	ENDDO
193	res = lintc_((/x,y,z,t/),g4d)
194	ret = .true.
195	END FUNCTION l4d_gen_sc
196
197	FUNCTION l5d_gen_sc(x,y,z,t,w,cx,cy,cz,ct,cw,cv,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	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
208	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
209	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
210	REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates
211	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cw !! Tabulated function W coordinates
212	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:,:) :: cv !! Tabulated function values
213	REAL(kind=wp), INTENT(out) :: res !! Interpolated value
214	PROCEDURE(locate) :: locator !! Locator function
215	LOGICAL :: ret !! .true. on success, .false. otherwise
216	REAL(kind=wp), DIMENSION(32,6) :: g5d
217	INTEGER :: ix0,iy0,iz0,it0,iw0
218	INTEGER :: i,a,b,c,d,e
219	ret = .false.
220	ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN
221	iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN
222	iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN
223	it0 = locator(t,ct) ; IF (it0 == 0) RETURN
224	iw0 = locator(w,cw) ; IF (iw0 == 0) RETURN
225	DO i=1,32
226	a=ix0+MOD((i-1),2) ; g5d(i,1) = cx(a)
227	b=iy0+MOD((i-1)/2,2) ; g5d(i,2) = cy(b)
228	c=iz0+MOD((i-1)/4,2) ; g5d(i,3) = cz(c)
229	d=it0+MOD((i-1)/8,2) ; g5d(i,4) = ct(d)
230	e=iw0+MOD((i-1)/16,2) ; g5d(i,5) = cw(e)
231	g5d(i,6) = cv(a,b,c,d,e)
232	ENDDO
233	res = lintc_((/x,y,z,t,w/),g5d)
234	ret = .true.
235	END FUNCTION l5d_gen_sc
236
237	! HARD-CODED versions
238	! -------------------
239
240	FUNCTION hdcdl1d_gen_sc(x,cx,cv,locator,res) RESULT(ret)
241	!! Hard-coded 1D linear interpolation (scalar)
242	!!
243	!! @warning
244	!! On error, __res__ output value is undefined.
245	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate
246	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
247	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cv !! Tabulated function values
248	REAL(kind=wp), INTENT(out) :: res !! Interpolated value
249	PROCEDURE(locate) :: locator !! Locator function
250	LOGICAL :: ret !! .true. on success, .false. otherwise
251	REAL(kind=wp) :: xd
252	INTEGER :: ix0,ix1
253	ret = .false.
254	ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN
255	ix1 = ix0+1
256	xd = (x-cx(ix0))/(cx(ix1)-cx(ix0))
257	res = cv(ix0)(1d0-xd)+cv(ix1)xd
258	ret = .true.
259	END FUNCTION hdcdl1d_gen_sc
260
261	FUNCTION hdcdl2d_gen_sc(x,y,cx,cy,cv,locator,res) RESULT(ret)
262	!! Hard-coded 2D linear interpolation (scalar)
263	!!
264	!! @warning
265	!! On error, __res__ output value is undefined.
266	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate
267	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate
268	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
269	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
270	REAL(kind=wp), INTENT(in), DIMENSION(:,:) :: cv !! Tabulated function values
271	REAL(kind=wp), INTENT(out) :: res !! Interpolated value
272	PROCEDURE(locate) :: locator !! Locator function
273	LOGICAL :: ret !! .true. on success, .false. otherwise
274	REAL(kind=wp) :: xd,yd
275	REAL(kind=wp) :: f0,f1
276	INTEGER :: ix0,iy0,ix1,iy1
277	ret = .false.
278	ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN
279	iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN
280	ix1 = ix0 + 1 ; iy1 = iy0 + 1
281	xd = (x-cx(ix0))/(cx(ix1)-cx(ix0))
282	yd = (y-cy(iy0))/(cy(iy1)-cy(iy0))
283	f0 = cv(ix0,iy0)(1d0-xd)+cv(ix1,iy0)xd
284	f1 = cv(ix0,iy1)(1d0-xd)+cv(ix1,iy1)xd
285	res = f0(1d0-yd)+f1yd
286	ret = .true.
287	END FUNCTION hdcdl2d_gen_sc
288
289	FUNCTION hdcdl3d_gen_sc(x,y,z,cx,cy,cz,cv,locator,res) RESULT(ret)
290	!! Hard-coded 3D linear interpolation (scalar)
291	!!
292	!! @warning
293	!! On error, __res__ output value is undefined.
294	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate
295	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate
296	REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate
297	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
298	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
299	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
300	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:) :: cv !! Tabulated function values
301	REAL(kind=wp), INTENT(out) :: res !! Interpolated value
302	PROCEDURE(locate) :: locator !! Locator function
303	LOGICAL :: ret !! .true. on success, .false. otherwise
304	REAL(kind=wp) :: xd,yd,zd
305	REAL(kind=wp) :: f00,f10,f01,f11,f0,f1
306	INTEGER :: ix0,iy0,iz0,ix1,iy1,iz1
307	ret = .false.
308	ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN
309	iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN
310	iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN
311	ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1
312	xd = (x-cx(ix0))/(cx(ix1)-cx(ix0))
313	yd = (y-cy(iy0))/(cy(iy1)-cy(iy0))
314	zd = (z-cz(iz0))/(cz(iz1)-cz(iz0))
315
316	f00 = cv(ix0,iy0,iz0)(1d0-xd)+cv(ix1,iy0,iz0)xd
317	f10 = cv(ix0,iy1,iz0)(1d0-xd)+cv(ix1,iy1,iz0)xd
318	f01 = cv(ix0,iy0,iz1)(1d0-xd)+cv(ix1,iy0,iz1)xd
319	f11 = cv(ix0,iy1,iz1)(1d0-xd)+cv(ix1,iy1,iz1)xd
320
321	f0 = f00 (1d0-yd)+f10yd
322	f1 = f00 (1d0-yd)+f10yd
323
324	res = f0(1d0-zd)+f1zd
325	ret = .true.
326	END FUNCTION hdcdl3d_gen_sc
327
328	FUNCTION hdcdl4d_gen_sc(x,y,z,t,cx,cy,cz,ct,cv,locator,res) RESULT(ret)
329	!! Hard-coded 4D linear interpolation (scalar)
330	!!
331	!! @warning
332	!! On error, __res__ output value is undefined.
333	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate
334	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate
335	REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate
336	REAL(kind=wp), INTENT(in) :: t !! T coordinate of the point to interpolate
337	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
338	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
339	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
340	REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates
341	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:) :: cv !! Tabulated function values
342	REAL(kind=wp), INTENT(out) :: res !! Interpolated value
343	PROCEDURE(locate) :: locator !! Locator function
344	LOGICAL :: ret !! .true. on success, .false. otherwise
345	REAL(kind=wp) :: xd,yd,zd,td
346	REAL(kind=wp) :: f000,f100,f010,f110,f001,f101,f011,f111, &
347	f00,f10,f01,f11,f0,f1
348	INTEGER :: ix0,iy0,iz0,it0,ix1,iy1,iz1,it1
349	ret = .false.
350	ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN
351	iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN
352	iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN
353	it0 = locator(t,ct) ; IF (it0 == 0) RETURN
354
355	ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ; it1=it0+1
356	xd = (x-cx(ix0))/(cx(ix1)-cx(ix0))
357	yd = (y-cy(iy0))/(cy(iy1)-cy(iy0))
358	zd = (z-cz(iz0))/(cz(iz1)-cz(iz0))
359	td = (t-ct(it0))/(ct(it1)-ct(it0))
360
361	f000 = cv(ix0,iy0,iz0,it0)(1d0-xd)+cv(ix1,iy0,iz0,it0)xd
362	f100 = cv(ix0,iy1,iz0,it0)(1d0-xd)+cv(ix1,iy1,iz0,it0)xd
363	f010 = cv(ix0,iy0,iz1,it0)(1d0-xd)+cv(ix1,iy0,iz1,it0)xd
364	f110 = cv(ix0,iy1,iz1,it0)(1d0-xd)+cv(ix1,iy1,iz1,it0)xd
365	f001 = cv(ix0,iy0,iz0,it1)(1d0-xd)+cv(ix1,iy0,iz0,it1)xd
366	f101 = cv(ix0,iy1,iz0,it1)(1d0-xd)+cv(ix1,iy1,iz0,it1)xd
367	f011 = cv(ix0,iy0,iz1,it1)(1d0-xd)+cv(ix1,iy0,iz1,it1)xd
368	f111 = cv(ix0,iy1,iz1,it1)(1d0-xd)+cv(ix1,iy1,iz1,it1)xd
369
370	f00 = f000(1d0-yd)+f100yd
371	f10 = f010(1d0-yd)+f110yd
372	f01 = f001(1d0-yd)+f101yd
373	f11 = f011(1d0-yd)+f111yd
374
375	f0 = f00 (1d0-zd)+f10zd
376	f1 = f01 (1d0-zd)+f11zd
377
378	res = f0(1d0-td)+f1td
379	ret = .true.
380	END FUNCTION hdcdl4d_gen_sc
381
382	FUNCTION hdcdl5d_gen_sc(x,y,z,t,w,cx,cy,cz,ct,cw,cv,locator,res) RESULT(ret)
383	!! Hard-coded 5D linear interpolation (scalar)
384	!!
385	!! @warning
386	!! On error, __res__ output value is undefined.
387	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to interpolate
388	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to interpolate
389	REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to interpolate
390	REAL(kind=wp), INTENT(in) :: t !! T coordinate of the point to interpolate
391	REAL(kind=wp), INTENT(in) :: w !! W coordinate of the point to interpolate
392	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
393	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
394	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
395	REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates
396	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cw !! Tabulated function W coordinates
397	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:,:) :: cv !! Tabulated function values
398	REAL(kind=wp), INTENT(out) :: res !! Interpolated value
399	PROCEDURE(locate) :: locator !! Locator function
400	LOGICAL :: ret !! .true. on success, .false. otherwise
401	REAL(kind=wp) :: xd,yd,zd,td,wd
402	REAL(kind=wp) :: f0000,f1000,f0100,f1100,f0010,f1010,f0110,f1110, &
403	f0001,f1001,f0101,f1101,f0011,f1011,f0111,f1111, &
404	f000,f100,f010,f110,f001,f101,f011,f111, &
405	f00,f10,f01,f11,f0,f1
406	INTEGER :: ix0,iy0,iz0,it0,iw0,ix1,iy1,iz1,it1,iw1
407	ret = .false.
408	ix0 = locator(x,cx) ; IF (ix0 == 0) RETURN
409	iy0 = locator(y,cy) ; IF (iy0 == 0) RETURN
410	iz0 = locator(z,cz) ; IF (iz0 == 0) RETURN
411	it0 = locator(t,ct) ; IF (it0 == 0) RETURN
412	iw0 = locator(w,cw) ; IF (iw0 == 0) RETURN
413
414	ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ; it1=it0+1 ; iw1=iw0+1
415	xd = (x-cx(ix0))/(cx(ix1)-cx(ix0))
416	yd = (y-cy(iy0))/(cy(iy1)-cy(iy0))
417	zd = (z-cz(iz0))/(cz(iz1)-cz(iz0))
418	td = (t-ct(it0))/(ct(it1)-ct(it0))
419	wd = (w-cw(iw0))/(cw(iw1)-cw(iw0))
420
421	f0000 = cv(ix0,iy0,iz0,it0,iw0)(1d0-xd)+cv(ix1,iy0,iz0,it0,iw0)xd
422	f1000 = cv(ix0,iy1,iz0,it0,iw0)(1d0-xd)+cv(ix0,iy1,iz0,it0,iw0)xd
423	f0100 = cv(ix0,iy0,iz1,it0,iw0)(1d0-xd)+cv(ix1,iy0,iz1,it0,iw0)xd
424	f1100 = cv(ix0,iy1,iz1,it0,iw0)(1d0-xd)+cv(ix1,iy1,iz1,it0,iw0)xd
425	f0010 = cv(ix0,iy0,iz0,it1,iw0)(1d0-xd)+cv(ix1,iy0,iz0,it1,iw0)xd
426	f1010 = cv(ix0,iy1,iz0,it1,iw0)(1d0-xd)+cv(ix1,iy1,iz0,it1,iw0)xd
427	f0110 = cv(ix0,iy0,iz1,it1,iw0)(1d0-xd)+cv(ix1,iy0,iz1,it1,iw0)xd
428	f1110 = cv(ix0,iy1,iz1,it1,iw0)(1d0-xd)+cv(ix1,iy1,iz1,it1,iw0)xd
429	f0001 = cv(ix0,iy0,iz0,it0,iw1)(1d0-xd)+cv(ix1,iy0,iz0,it0,iw1)xd
430	f1001 = cv(ix0,iy1,iz0,it0,iw1)(1d0-xd)+cv(ix0,iy1,iz0,it0,iw1)xd
431	f0101 = cv(ix0,iy0,iz1,it0,iw1)(1d0-xd)+cv(ix1,iy0,iz1,it0,iw1)xd
432	f1101 = cv(ix0,iy1,iz1,it0,iw1)(1d0-xd)+cv(ix1,iy1,iz1,it0,iw1)xd
433	f0011 = cv(ix0,iy0,iz0,it1,iw1)(1d0-xd)+cv(ix1,iy0,iz0,it1,iw1)xd
434	f1011 = cv(ix0,iy1,iz0,it1,iw1)(1d0-xd)+cv(ix1,iy1,iz0,it1,iw1)xd
435	f0111 = cv(ix0,iy0,iz1,it1,iw1)(1d0-xd)+cv(ix1,iy0,iz1,it1,iw1)xd
436	f1111 = cv(ix0,iy1,iz1,it1,iw1)(1d0-xd)+cv(ix1,iy1,iz1,it1,iw1)xd
437
438	f000 = f0000(1d0-yd) + f1000yd
439	f100 = f0100(1d0-yd) + f1100yd
440	f010 = f0010(1d0-yd) + f1010yd
441	f110 = f0110(1d0-yd) + f1110yd
442	f101 = f0101(1d0-yd) + f1101yd
443	f011 = f0011(1d0-yd) + f1011yd
444	f111 = f0111(1d0-yd) + f1111yd
445
446	f00 = f000(1d0-zd)+f100zd
447	f10 = f010(1d0-zd)+f110zd
448	f01 = f001(1d0-zd)+f101zd
449	f11 = f011(1d0-zd)+f111zd
450
451	f0 = f00 (1d0-td)+f10td
452	f1 = f01 (1d0-td)+f11td
453
454	res = f0(1d0-wd)+f1wd
455	ret = .true.
456	END FUNCTION hdcdl5d_gen_sc
457
458	!--------------------
459	! VECTORIZED VERSIONS
460	!--------------------
461
462	! GENERIC versions
463	! ----------------
464
465	FUNCTION l1d_gen_ve(x,cx,cv,locator,res) RESULT(ret)
466	!! Wrapper interface for 1D linear interpolation (vector)
467	!!
468	!! @warning
469	!! On error, __res__ output vector is undefined (i.e. not allocated).
470	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
471	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
472	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cv !! Tabulated function values
473	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
474	PROCEDURE(locate) :: locator !! Locator function
475	LOGICAL :: ret !! .true. on success, .false. otherwise
476	REAL(kind=wp), DIMENSION(2,2) :: g1d
477	INTEGER :: n,j,i,ix
478	ret = .false.
479	n = SIZE(x) ; ALLOCATE(res(n))
480	DO j=1,n
481	ix = locator(x(j),cx)
482	IF (ix == 0) THEN
483	DEALLOCATE(res) ; RETURN
484	ENDIF
485	DO i=0,1
486	g1d(i+1,1) = cx(ix+i)
487	g1d(i+1,2) = cv(ix+i)
488	ENDDO
489	res(j) = lintc_((/x(j)/),g1d)
490	ENDDO
491	ret = .true.
492	END FUNCTION l1d_gen_ve
493
494	FUNCTION l2d_gen_ve(x,y,cx,cy,cv,locator,res) RESULT(ret)
495	!! Wrapper interface for 2D linear interpolation (vector)
496	!!
497	!! @warning
498	!! On error, __res__ output vector is undefined (i.e. not allocated).
499	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
500	REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate
501	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
502	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
503	REAL(kind=wp), INTENT(in), DIMENSION(:,:) :: cv !! Tabulated function values
504	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
505	PROCEDURE(locate) :: locator !! Locator function
506	LOGICAL :: ret !! .true. on success, .false. otherwise
507	REAL(kind=wp), DIMENSION(4,3) :: g2d
508	INTEGER :: n,j,ix0,iy0,i,a,b
509	ret = .false.
510	n = SIZE(x) ; ALLOCATE(res(n))
511	DO j=1,n
512	ix0 = locator(x(j),cx)
513	iy0 = locator(y(j),cy)
514	IF (ix0 == 0 .OR. iy0 == 0) THEN
515	DEALLOCATE(res) ; RETURN
516	ENDIF
517	DO i=1,4
518	a=ix0+MOD((i-1),2) ; g2d(i,1) = cx(a)
519	b=iy0+MOD((i-1)/2,2) ; g2d(i,2) = cy(b)
520	g2d(i,3) = cv(a,b)
521	ENDDO
522	res(j) = lintc_((/x(j),y(j)/),g2d)
523	ENDDO
524	ret = .true.
525	END FUNCTION l2d_gen_ve
526
527	FUNCTION l3d_gen_ve(x,y,z,cx,cy,cz,cv,locator,res) RESULT(ret)
528	!! Wrapper interface for 3D linear interpolation (vector)
529	!!
530	!! @warning
531	!! On error, __res__ output vector is undefined (i.e. not allocated).
532	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
533	REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate
534	REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate
535	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
536	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
537	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
538	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:) :: cv !! Tabulated function values
539	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
540	PROCEDURE(locate) :: locator !! Locator function
541	LOGICAL :: ret !! .true. on success, .false. otherwise
542	REAL(kind=wp), DIMENSION(8,4) :: g3d
543	INTEGER :: ix0,iy0,iz0
544	INTEGER :: n,j,i,a,b,c
545	ret = .false.
546	n = SIZE(x) ; ALLOCATE(res(n))
547	DO j=1,n
548	ix0 = locator(x(j),cx)
549	iy0 = locator(y(j),cy)
550	iz0 = locator(z(j),cz)
551	IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0) THEN
552	DEALLOCATE(res) ; RETURN
553	ENDIF
554	DO i=1,8
555	a=ix0+MOD((i-1),2) ; g3d(i,1) = cx(a)
556	b=iy0+MOD((i-1)/2,2) ; g3d(i,2) = cy(b)
557	c=iz0+MOD((i-1)/4,2) ; g3d(i,3) = cz(c)
558	g3d(i,4) = cv(a,b,c)
559	ENDDO
560	res(j) = lintc_((/x(j),y(j),z(j)/),g3d)
561	ENDDO
562	ret = .true.
563	END FUNCTION l3d_gen_ve
564
565	FUNCTION l4d_gen_ve(x,y,z,t,cx,cy,cz,ct,cv,locator,res) RESULT(ret)
566	!! Wrapper interface for 4D linear interpolation (vector)
567	!!
568	!! @warning
569	!! On error, __res__ output vector is undefined (i.e. not allocated).
570	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
571	REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate
572	REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate
573	REAL(kind=wp), INTENT(in), DIMENSION(:) :: t !! T coordinate of the points to interpolate
574	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
575	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
576	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
577	REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates
578	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:) :: cv !! Tabulated function values
579	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
580	PROCEDURE(locate) :: locator !! Locator function
581	LOGICAL :: ret !! .true. on success, .false. otherwise
582	REAL(kind=wp), DIMENSION(16,5) :: g4d
583	INTEGER :: ix0,iy0,iz0,it0
584	INTEGER :: n,j,i,a,b,c,d
585	ret = .false.
586	n = SIZE(x) ; ALLOCATE(res(n))
587	DO j=1,n
588	ix0 = locator(x(j),cx)
589	iy0 = locator(y(j),cy)
590	iz0 = locator(z(j),cz)
591	it0 = locator(t(j),ct)
592	IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0 .OR. it0 == 0) THEN
593	DEALLOCATE(res) ; RETURN
594	ENDIF
595	DO i=1,16
596	a=ix0+MOD((i-1),2) ; g4d(i,1) = cx(a)
597	b=iy0+MOD((i-1)/2,2) ; g4d(i,2) = cy(b)
598	c=iz0+MOD((i-1)/4,2) ; g4d(i,3) = cz(c)
599	d=it0+MOD((i-1)/8,2) ; g4d(i,4) = ct(d)
600	g4d(i,5) = cv(a,b,c,d)
601	ENDDO
602	res(j) = lintc_((/x(j),y(j),z(j),t(j)/),g4d)
603	ENDDO
604	ret = .true.
605	END FUNCTION l4d_gen_ve
606
607	FUNCTION l5d_gen_ve(x,y,z,t,w,cx,cy,cz,ct,cw,cv,locator,res) RESULT(ret)
608	!! Wrapper interface for 5D linear interpolation (vector)
609	!!
610	!! @warning
611	!! On error, __res__ output vector is undefined (i.e. not allocated).
612	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
613	REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate
614	REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate
615	REAL(kind=wp), INTENT(in), DIMENSION(:) :: t !! T coordinate of the points to interpolate
616	REAL(kind=wp), INTENT(in), DIMENSION(:) :: w !! W coordinate of the points to interpolate
617	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
618	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
619	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
620	REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates
621	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cw !! Tabulated function W coordinates
622	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:,:) :: cv !! Tabulated function values
623	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
624	PROCEDURE(locate) :: locator !! Locator function
625	LOGICAL :: ret !! .true. on success, .false. otherwise
626	REAL(kind=wp), DIMENSION(32,6) :: g5d
627	INTEGER :: ix0,iy0,iz0,it0,iw0
628	INTEGER :: n,j,i,a,b,c,d,e
629	ret = .false.
630	n = SIZE(x) ; ALLOCATE(res(n))
631	DO j=1,n
632	ix0 = locator(x(j),cx)
633	iy0 = locator(y(j),cy)
634	iz0 = locator(z(j),cz)
635	it0 = locator(t(j),ct)
636	iw0 = locator(w(j),cw)
637	IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0 .OR. it0 == 0 .OR. iw0 == 0) THEN
638	DEALLOCATE(res) ; RETURN
639	ENDIF
640	DO i=1,32
641	a=ix0+MOD((i-1),2) ; g5d(i,1) = cx(a)
642	b=iy0+MOD((i-1)/2,2) ; g5d(i,2) = cy(b)
643	c=iz0+MOD((i-1)/4,2) ; g5d(i,3) = cz(c)
644	d=it0+MOD((i-1)/8,2) ; g5d(i,4) = ct(d)
645	e=iw0+MOD((i-1)/16,2) ; g5d(i,5) = cw(e)
646	g5d(i,6) = cv(a,b,c,d,e)
647	ENDDO
648	res(j) = lintc_((/x,y,z,t,w/),g5d)
649	ENDDO
650	ret = .true.
651	END FUNCTION l5d_gen_ve
652
653	! HARD-CODED versions
654	! -------------------
655
656	FUNCTION hdcdl1d_gen_ve(x,cx,cv,locator,res) RESULT(ret)
657	!! Hard-coded 1D linear interpolation (vector)
658	!!
659	!! @warning
660	!! On error, __res__ output vector is undefined (i.e. not allocated).
661	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
662	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
663	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cv !! Tabulated function values
664	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
665	PROCEDURE(locate) :: locator !! Locator function
666	LOGICAL :: ret !! .true. on success, .false. otherwise
667	REAL(kind=wp) :: xd
668	INTEGER :: ix0,ix1
669	INTEGER :: n,j
670	ret = .false.
671	n = SIZE(x) ; ALLOCATE(res(n))
672	DO j=1,n
673	ix0 = locator(x(j),cx)
674	IF (ix0 == 0) THEN
675	DEALLOCATE(res) ; RETURN
676	ENDIF
677	ix1 = ix0+1
678	xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0))
679	res(j) = cv(ix0)(1d0-xd)+cv(ix1)xd
680	ENDDO
681	ret = .true.
682	END FUNCTION hdcdl1d_gen_ve
683
684	FUNCTION hdcdl2d_gen_ve(x,y,cx,cy,cv,locator,res) RESULT(ret)
685	!! Hard-coded 2D linear interpolation (vector)
686	!!
687	!! @warning
688	!! On error, __res__ output vector is undefined (i.e. not allocated).
689	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
690	REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate
691	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
692	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
693	REAL(kind=wp), INTENT(in), DIMENSION(:,:) :: cv !! Tabulated function values
694	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
695	PROCEDURE(locate) :: locator !! Locator function
696	LOGICAL :: ret !! .true. on success, .false. otherwise
697	REAL(kind=wp) :: xd,yd
698	REAL(kind=wp) :: f0,f1
699	INTEGER :: ix0,iy0,ix1,iy1
700	INTEGER :: n,j
701	ret = .false.
702	n = SIZE(x) ; ALLOCATE(res(n))
703	DO j=1,n
704	ix0 = locator(x(j),cx)
705	iy0 = locator(y(j),cy)
706	IF (ix0 == 0 .OR. iy0 == 0) THEN
707	DEALLOCATE(res) ; RETURN
708	ENDIF
709	ix1 = ix0 + 1 ; iy1 = iy0 + 1
710	xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0))
711	yd = (y(j)-cy(iy0))/(cy(iy1)-cy(iy0))
712	f0 = cv(ix0,iy0)(1d0-xd)+cv(ix1,iy0)xd
713	f1 = cv(ix0,iy1)(1d0-xd)+cv(ix1,iy1)xd
714	res(j) = f0(1d0-yd)+f1yd
715	ENDDO
716	ret = .true.
717	END FUNCTION hdcdl2d_gen_ve
718
719	FUNCTION hdcdl3d_gen_ve(x,y,z,cx,cy,cz,cv,locator,res) RESULT(ret)
720	!! Hard-coded 3D linear interpolation (vector)
721	!!
722	!! @warning
723	!! On error, __res__ output vector is undefined (i.e. not allocated).
724	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
725	REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate
726	REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate
727	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
728	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
729	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
730	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:) :: cv !! Tabulated function values
731	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
732	PROCEDURE(locate) :: locator !! Locator function
733	LOGICAL :: ret !! .true. on success, .false. otherwise
734	REAL(kind=wp) :: xd,yd,zd
735	REAL(kind=wp) :: f00,f10,f01,f11,f0,f1
736	INTEGER :: ix0,iy0,iz0,ix1,iy1,iz1
737	INTEGER :: n,j
738	ret = .false.
739	n = SIZE(x) ; ALLOCATE(res(n))
740	DO j=1,n
741	ix0 = locator(x(j),cx)
742	iy0 = locator(y(j),cy)
743	iz0 = locator(z(j),cz)
744	IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0) THEN
745	DEALLOCATE(res) ; RETURN
746	ENDIF
747	ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1
748	xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0))
749	yd = (y(j)-cy(iy0))/(cy(iy1)-cy(iy0))
750	zd = (z(j)-cz(iz0))/(cz(iz1)-cz(iz0))
751
752	f00 = cv(ix0,iy0,iz0)(1d0-xd)+cv(ix1,iy0,iz0)xd
753	f10 = cv(ix0,iy1,iz0)(1d0-xd)+cv(ix1,iy1,iz0)xd
754	f01 = cv(ix0,iy0,iz1)(1d0-xd)+cv(ix1,iy0,iz1)xd
755	f11 = cv(ix0,iy1,iz1)(1d0-xd)+cv(ix1,iy1,iz1)xd
756
757	f0 = f00 (1d0-yd)+f10yd
758	f1 = f00 (1d0-yd)+f10yd
759
760	res(j) = f0(1d0-zd)+f1zd
761	ENDDO
762	ret = .true.
763	END FUNCTION hdcdl3d_gen_ve
764
765	FUNCTION hdcdl4d_gen_ve(x,y,z,t,cx,cy,cz,ct,cv,locator,res) RESULT(ret)
766	!! Hard-coded 4D linear interpolation (vector)
767	!!
768	!! @warning
769	!! On error, __res__ output vector is undefined (i.e. not allocated).
770	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
771	REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate
772	REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate
773	REAL(kind=wp), INTENT(in), DIMENSION(:) :: t !! T coordinate of the points to interpolate
774	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
775	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
776	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
777	REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates
778	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:) :: cv !! Tabulated function values
779	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
780	PROCEDURE(locate) :: locator !! Locator function
781	LOGICAL :: ret !! .true. on success, .false. otherwise
782	REAL(kind=wp) :: xd,yd,zd,td
783	REAL(kind=wp) :: f000,f100,f010,f110,f001,f101,f011,f111, &
784	f00,f10,f01,f11,f0,f1
785	INTEGER :: ix0,iy0,iz0,it0,ix1,iy1,iz1,it1
786	INTEGER :: n,j
787	ret = .false.
788	n = SIZE(x) ; ALLOCATE(res(n))
789	DO j=1,n
790	ix0 = locator(x(j),cx)
791	iy0 = locator(y(j),cy)
792	iz0 = locator(z(j),cz)
793	it0 = locator(t(j),ct)
794	IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0 .OR. it0 == 0) THEN
795	DEALLOCATE(res) ; RETURN
796	ENDIF
797	ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ; it1=it0+1
798	xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0))
799	yd = (y(j)-cy(iy0))/(cy(iy1)-cy(iy0))
800	zd = (z(j)-cz(iz0))/(cz(iz1)-cz(iz0))
801	td = (t(j)-ct(it0))/(ct(it1)-ct(it0))
802
803	f000 = cv(ix0,iy0,iz0,it0)(1d0-xd)+cv(ix1,iy0,iz0,it0)xd
804	f100 = cv(ix0,iy1,iz0,it0)(1d0-xd)+cv(ix1,iy1,iz0,it0)xd
805	f010 = cv(ix0,iy0,iz1,it0)(1d0-xd)+cv(ix1,iy0,iz1,it0)xd
806	f110 = cv(ix0,iy1,iz1,it0)(1d0-xd)+cv(ix1,iy1,iz1,it0)xd
807	f001 = cv(ix0,iy0,iz0,it1)(1d0-xd)+cv(ix1,iy0,iz0,it1)xd
808	f101 = cv(ix0,iy1,iz0,it1)(1d0-xd)+cv(ix1,iy1,iz0,it1)xd
809	f011 = cv(ix0,iy0,iz1,it1)(1d0-xd)+cv(ix1,iy0,iz1,it1)xd
810	f111 = cv(ix0,iy1,iz1,it1)(1d0-xd)+cv(ix1,iy1,iz1,it1)xd
811
812	f00 = f000(1d0-yd)+f100yd
813	f10 = f010(1d0-yd)+f110yd
814	f01 = f001(1d0-yd)+f101yd
815	f11 = f011(1d0-yd)+f111yd
816
817	f0 = f00 (1d0-zd)+f10zd
818	f1 = f01 (1d0-zd)+f11zd
819
820	res(j) = f0(1d0-td)+f1td
821	ENDDO
822	ret = .true.
823	END FUNCTION hdcdl4d_gen_ve
824
825	FUNCTION hdcdl5d_gen_ve(x,y,z,t,w,cx,cy,cz,ct,cw,cv,locator,res) RESULT(ret)
826	!! Hard-coded 5D linear interpolation (vector)
827	!!
828	!! @warning
829	!! On error, __res__ output vector is undefined (i.e. not allocated).
830	REAL(kind=wp), INTENT(in), DIMENSION(:) :: x !! X coordinate of the points to interpolate
831	REAL(kind=wp), INTENT(in), DIMENSION(:) :: y !! Y coordinate of the points to interpolate
832	REAL(kind=wp), INTENT(in), DIMENSION(:) :: z !! Z coordinate of the points to interpolate
833	REAL(kind=wp), INTENT(in), DIMENSION(:) :: t !! T coordinate of the points to interpolate
834	REAL(kind=wp), INTENT(in), DIMENSION(:) :: w !! W coordinate of the points to interpolate
835	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cx !! Tabulated function X coordinates
836	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cy !! Tabulated function Y coordinates
837	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cz !! Tabulated function Z coordinates
838	REAL(kind=wp), INTENT(in), DIMENSION(:) :: ct !! Tabulated function T coordinates
839	REAL(kind=wp), INTENT(in), DIMENSION(:) :: cw !! Tabulated function W coordinates
840	REAL(kind=wp), INTENT(in), DIMENSION(:,:,:,:,:) :: cv !! Tabulated function values
841	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: res !! Interpolated values
842	PROCEDURE(locate) :: locator !! Locator function
843	LOGICAL :: ret !! .true. on success, .false. otherwise
844	REAL(kind=wp) :: xd,yd,zd,td,wd
845	REAL(kind=wp) :: f0000,f1000,f0100,f1100,f0010,f1010,f0110,f1110, &
846	f0001,f1001,f0101,f1101,f0011,f1011,f0111,f1111, &
847	f000,f100,f010,f110,f001,f101,f011,f111, &
848	f00,f10,f01,f11,f0,f1
849	INTEGER :: ix0,iy0,iz0,it0,iw0,ix1,iy1,iz1,it1,iw1
850	INTEGER :: n,j
851	ret = .false.
852	n = SIZE(x) ; ALLOCATE(res(n))
853	DO j=1,n
854	ix0 = locator(x(j),cx)
855	iy0 = locator(y(j),cy)
856	iz0 = locator(z(j),cz)
857	it0 = locator(t(j),ct)
858	iw0 = locator(w(j),cw)
859	IF (ix0 == 0 .OR. iy0 == 0 .OR. iz0 == 0 .OR. it0 == 0 .OR. iw0 == 0) THEN
860	DEALLOCATE(res) ; RETURN
861	ENDIF
862	ix1=ix0+1 ; iy1=iy0+1 ; iz1=iz0+1 ;
863	xd = (x(j)-cx(ix0))/(cx(ix1)-cx(ix0))
864	yd = (y(j)-cy(iy0))/(cy(iy1)-cy(iy0))
865	zd = (z(j)-cz(iz0))/(cz(iz1)-cz(iz0))
866	td = (t(j)-ct(it0))/(ct(it1)-ct(it0))
867	wd = (w(j)-cw(it0))/(cw(it1)-cw(it0))
868
869	f0000 = cv(ix0,iy0,iz0,it0,iw0)(1d0-xd)+cv(ix1,iy0,iz0,it0,iw0)xd
870	f1000 = cv(ix0,iy1,iz0,it0,iw0)(1d0-xd)+cv(ix0,iy1,iz0,it0,iw0)xd
871	f0100 = cv(ix0,iy0,iz1,it0,iw0)(1d0-xd)+cv(ix1,iy0,iz1,it0,iw0)xd
872	f1100 = cv(ix0,iy1,iz1,it0,iw0)(1d0-xd)+cv(ix1,iy1,iz1,it0,iw0)xd
873	f0010 = cv(ix0,iy0,iz0,it1,iw0)(1d0-xd)+cv(ix1,iy0,iz0,it1,iw0)xd
874	f1010 = cv(ix0,iy1,iz0,it1,iw0)(1d0-xd)+cv(ix1,iy1,iz0,it1,iw0)xd
875	f0110 = cv(ix0,iy0,iz1,it1,iw0)(1d0-xd)+cv(ix1,iy0,iz1,it1,iw0)xd
876	f1110 = cv(ix0,iy1,iz1,it1,iw0)(1d0-xd)+cv(ix1,iy1,iz1,it1,iw0)xd
877	f0001 = cv(ix0,iy0,iz0,it0,iw1)(1d0-xd)+cv(ix1,iy0,iz0,it0,iw1)xd
878	f1001 = cv(ix0,iy1,iz0,it0,iw1)(1d0-xd)+cv(ix0,iy1,iz0,it0,iw1)xd
879	f0101 = cv(ix0,iy0,iz1,it0,iw1)(1d0-xd)+cv(ix1,iy0,iz1,it0,iw1)xd
880	f1101 = cv(ix0,iy1,iz1,it0,iw1)(1d0-xd)+cv(ix1,iy1,iz1,it0,iw1)xd
881	f0011 = cv(ix0,iy0,iz0,it1,iw1)(1d0-xd)+cv(ix1,iy0,iz0,it1,iw1)xd
882	f1011 = cv(ix0,iy1,iz0,it1,iw1)(1d0-xd)+cv(ix1,iy1,iz0,it1,iw1)xd
883	f0111 = cv(ix0,iy0,iz1,it1,iw1)(1d0-xd)+cv(ix1,iy0,iz1,it1,iw1)xd
884	f1111 = cv(ix0,iy1,iz1,it1,iw1)(1d0-xd)+cv(ix1,iy1,iz1,it1,iw1)xd
885
886	f000 = f0000(1d0-yd) + f1000yd
887	f100 = f0100(1d0-yd) + f1100yd
888	f010 = f0010(1d0-yd) + f1010yd
889	f110 = f0110(1d0-yd) + f1110yd
890	f101 = f0101(1d0-yd) + f1101yd
891	f011 = f0011(1d0-yd) + f1011yd
892	f111 = f0111(1d0-yd) + f1111yd
893
894	f00 = f000(1d0-zd)+f100zd
895	f10 = f010(1d0-zd)+f110zd
896	f01 = f001(1d0-zd)+f101zd
897	f11 = f011(1d0-zd)+f111zd
898
899	f0 = f00 (1d0-td)+f10td
900	f1 = f01 (1d0-td)+f11td
901
902	res(j) = f0(1d0-wd)+f1wd
903	ENDDO
904	ret = .true.
905	END FUNCTION hdcdl5d_gen_ve
906
907	END MODULE LINTGEN

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

Download in other formats:

Original Format