Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

lint_dset.F90 @ 3608

Last change on this file since 3608 was 3560, checked in by debatzbr, 8 weeks ago
Addition of the microphysics model in moments.
File size: 40.4 KB

Rev	Line
[3560]	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
	38	MODULE 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)+f1yd
	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)+f10yd
	482	f1 = f00 (1d0-yd)+f10yd
	483
	484	res = f0(1d0-zd)+f1zd
	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)+f100yd
	528	f10 = f010(1d0-yd)+f110yd
	529	f01 = f001(1d0-yd)+f101yd
	530	f11 = f011(1d0-yd)+f111yd
	531
	532	f0 = f00 (1d0-zd)+f10zd
	533	f1 = f01 (1d0-zd)+f11zd
	534
	535	res = f0(1d0-td)+f1td
	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) + f1000yd
	592	f100 = f0100(1d0-yd) + f1100yd
	593	f010 = f0010(1d0-yd) + f1010yd
	594	f110 = f0110(1d0-yd) + f1110yd
	595	f101 = f0101(1d0-yd) + f1101yd
	596	f011 = f0011(1d0-yd) + f1011yd
	597	f111 = f0111(1d0-yd) + f1111yd
	598
	599	f00 = f000(1d0-zd)+f100zd
	600	f10 = f010(1d0-zd)+f110zd
	601	f01 = f001(1d0-zd)+f101zd
	602	f11 = f011(1d0-zd)+f111zd
	603
	604	f0 = f00 (1d0-td)+f10td
	605	f1 = f01 (1d0-td)+f11td
	606
	607	res = f0(1d0-wd)+f1wd
	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)+f1yd
	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)+f10yd
	706	f1 = f00 (1d0-yd)+f10yd
	707
	708	res(j) = f0(1d0-zd)+f1zd
	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)+f100yd
	758	f10 = f010(1d0-yd)+f110yd
	759	f01 = f001(1d0-yd)+f101yd
	760	f11 = f011(1d0-yd)+f111yd
	761
	762	f0 = f00 (1d0-zd)+f10zd
	763	f1 = f01 (1d0-zd)+f11zd
	764
	765	res(j) = f0(1d0-td)+f1td
	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) + f1000yd
	828	f100 = f0100(1d0-yd) + f1100yd
	829	f010 = f0010(1d0-yd) + f1010yd
	830	f110 = f0110(1d0-yd) + f1110yd
	831	f101 = f0101(1d0-yd) + f1101yd
	832	f011 = f0011(1d0-yd) + f1011yd
	833	f111 = f0111(1d0-yd) + f1111yd
	834
	835	f00 = f000(1d0-zd)+f100zd
	836	f10 = f010(1d0-zd)+f110zd
	837	f01 = f001(1d0-zd)+f101zd
	838	f11 = f011(1d0-zd)+f111zd
	839
	840	f0 = f00 (1d0-td)+f10td
	841	f1 = f01 (1d0-td)+f11td
	842
	843	res(j) = f0(1d0-wd)+f1wd
	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.

Download in other formats:

Original Format