Context Navigation

datasets.F90 @ 1862

Last change on this file since 1862 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: 48.1 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: datasets.F90
36	!! summary: Dataset module definition file
37	!! author: J. Burgalat
38	!! date: 2014
39
40	#if HAVE_CONFIG_H
41	#include "config.h"
42	#endif
43
44	MODULE DATASETS
45	!! dataset definitions module
46	!!
47	!! This module defines simple derived types that encapsulate data set variables. For a N-dimension
48	!! dataset, N vectors of \(n_{i}\) elements and a single N-dimensional array of
49	!! \(\prod_{i}^{N} n_{i}\) elements are defined.
50	!!
51	!! If the data set is to be used within [[lintdset(module)]] module then each coordinate values of the
52	!! dataset must be sorted either in ascending or descending order with no duplicates. The module does
53	!! not provide functionnality to sort and to check such requirements.
54	!!
55	!! The module also provides two interfaces to initialize data sets from input data file which can be
56	!! a NetCDF file (if NetCDF support is enabled at compilation time) or an ASCII file. In the latter
57	!! case, the file must contain a header which is formatted as follows:
58	!!
59	!! - The first line must contain only one value which is the number of coordinates (__N__)
60	!! - The second line must contain all the dimension sizes (that is __N__ values)
61	!! - Each other lines should contain __N+1__ columns with, for the __N__ first columns, the
62	!! coordinates values and finally the point value in the last column.
63	!!
64	!! @note
65	!! Note that for ASCII file, data must be ordered so first dimensions vary first. Same requirement
66	!! is needed for NetCDF file but in most cases, it is implicitly done (if dimensions are ordered).
67	USE LINT_PREC
68	#if HAVE_NC_FTN
69	USE NETCDF
70	#endif
71	IMPLICIT NONE
72
73	PRIVATE
74
75	PUBLIC :: read_dset, clear_dset, is_in, debug
76
77	LOGICAL :: debug = .false. !! A control flag to enable verbose mode
78
79	#if HAVE_NC_FTN
80	LOGICAL, PUBLIC, PARAMETER ::nc_supported = .true. !! NetCDF input files are supported
81	#else
82	LOGICAL, PUBLIC, PARAMETER ::nc_supported = .false. !! NetCDF input files are not supported
83	#endif
84
85	!> Initialize a data set from either an ASCII or a NetCDF file
86	!!
87	!! Whatever the kind of file, this interface assumes that you know the dimensionnality of the
88	!! extracted variable. If file content is not compatible with the kind of data set given as
89	!! output argument, or if some I/O error occured, the dataset is cleared.
90	!! @note
91	!! Netcdf reader interface is available only if the library has been compiled with NetCDF support.
92	INTERFACE read_dset
93	#if HAVE_NC_FTN
94	MODULE PROCEDURE ncdf_rd_1d,ncdf_rd_2d,ncdf_rd_3d,ncdf_rd_4d,ncdf_rd_5d
95	#if HAVE_NC4_FTN
96	MODULE PROCEDURE ncdf4_rd_1d,ncdf4_rd_2d,ncdf4_rd_3d,ncdf4_rd_4d,ncdf4_rd_5d
97	#endif
98	#endif
99	MODULE PROCEDURE ascii_rd_1d,ascii_rd_2d,ascii_rd_3d,ascii_rd_4d,ascii_rd_5d
100	END INTERFACE
101
102	!> Clear the given data set
103	INTERFACE clear_dset
104	MODULE PROCEDURE clr_1d_set,clr_2d_set,clr_3d_set,clr_4d_set,clr_5d_set
105	END INTERFACE
106
107	!> Check if given point is within the data set
108	INTERFACE is_in
109	MODULE PROCEDURE is_in_1d,is_in_2d,is_in_3d,is_in_4d,is_in_5d
110	END INTERFACE
111
112	!> Private interface to netcdf informations getters
113	#if HAVE_NC_FTN
114	INTERFACE get_nc_info
115	MODULE PROCEDURE get_nc3_info
116	#if HAVE_NC4_FTN
117	MODULE PROCEDURE get_nc4_info
118	#endif
119	END INTERFACE
120	#endif
121
122	TYPE, PUBLIC :: DSET1D
123	!! A 1D data set
124	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
125	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
126	END TYPE DSET1D
127
128	TYPE, PUBLIC :: DSET2D
129	!! A 2D data set
130	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
131	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: y !! Y coordinate tabulated values
132	REAL(kind=wp), DIMENSION(:,:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
133	END TYPE DSET2D
134
135	TYPE, PUBLIC :: DSET3D
136	!! A 3D data set
137	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
138	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: y !! Y coordinate tabulated values
139	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: z !! Z coordinate tabulated values
140	REAL(kind=wp), DIMENSION(:,:,:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
141	END TYPE DSET3D
142
143	TYPE, PUBLIC :: DSET4D
144	!! A 4D data set
145	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
146	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: y !! Y coordinate tabulated values
147	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: z !! Z coordinate tabulated values
148	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: t !! T coordinate tabulated values
149	REAL(kind=wp), DIMENSION(:,:,:,:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
150	END TYPE DSET4D
151
152	TYPE, PUBLIC :: DSET5D
153	!! A 5D data set
154	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
155	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: y !! Y coordinate tabulated values
156	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: z !! Z coordinate tabulated values
157	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: t !! T coordinate tabulated values
158	REAL(kind=wp), DIMENSION(:), ALLOCATABLE :: w !! W coordinate tabulated values
159	REAL(kind=wp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
160	END TYPE DSET5D
161
162	CONTAINS
163
164
165	FUNCTION is_in_1d(set,x) RESULT(ret)
166	!! Check if point is in the 1D data set
167	TYPE(DSET1D), INTENT(in) :: set !! Dataset object to search in
168	REAL(kind=wp), INTENT(in) :: x !! coordinate of the point to check
169	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
170	REAL(kind=wp) :: l,u
171	ret=.true.
172	l = set%x(1) ; u= set%x(size(set%x))
173	IF ((x>=l.EQV.x<=u).OR.(x<=l.EQV.x>=u)) ret=.false.
174	RETURN
175	END FUNCTION is_in_1d
176
177	FUNCTION is_in_2d(set,x,y) RESULT(ret)
178	!! Check if point is in the 2D data set
179	TYPE(DSET2D), INTENT(in) :: set !! Dataset object to search in
180	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to check
181	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to check
182	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
183	REAL(kind=wp) :: l,u
184	ret=.false.
185	l = set%x(1) ; u= set%x(size(set%x))
186	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
187	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
188	l = set%y(1) ; u= set%y(size(set%y))
189	IF ((y>l.AND.y>u).OR.(y<l.AND.y<u)) RETURN
190	ret=.true.
191	RETURN
192	END FUNCTION is_in_2d
193
194	FUNCTION is_in_3d(set,x,y,z) RESULT(ret)
195	!! Check if point is in the 3D data set
196	TYPE(DSET3D), INTENT(in) :: set !! Dataset object to search in
197	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to check
198	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to check
199	REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to check
200	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
201	REAL(kind=wp) :: l,u
202	ret=.false.
203	l = set%x(1) ; u= set%x(size(set%x))
204	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
205	l = set%y(1) ; u= set%y(size(set%y))
206	IF ((y>l.AND.y>u).OR.(y<l.AND.y<u)) RETURN
207	l = set%z(1) ; u= set%z(size(set%z))
208	IF ((z>l.AND.z>u).OR.(z<l.AND.z<u)) RETURN
209	ret=.true.
210	RETURN
211	END FUNCTION is_in_3d
212
213	FUNCTION is_in_4d(set,x,y,z,t) RESULT(ret)
214	!! Check if point is in the 4D data set
215	TYPE(DSET4D), INTENT(in) :: set !! Dataset object to search in
216	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to check
217	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to check
218	REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to check
219	REAL(kind=wp), INTENT(in) :: t !! T coordinate of the point to check
220	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
221	REAL(kind=wp) :: l,u
222	ret=.false.
223	l = set%x(1) ; u= set%x(size(set%x))
224	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
225	l = set%y(1) ; u= set%y(size(set%y))
226	IF ((y>l.AND.y>u).OR.(y<l.AND.y<u)) RETURN
227	l = set%z(1) ; u= set%z(size(set%z))
228	IF ((z>l.AND.z>u).OR.(z<l.AND.z<u)) RETURN
229	l = set%t(1) ; u= set%t(size(set%t))
230	IF ((t>l.AND.t>u).OR.(t<l.AND.t<u)) RETURN
231	ret=.true.
232	RETURN
233	END FUNCTION is_in_4d
234
235	FUNCTION is_in_5d(set,x,y,z,t,w) RESULT(ret)
236	!! Check if point is in the 4D data set
237	TYPE(DSET5D), INTENT(in) :: set !! Dataset object to search in
238	REAL(kind=wp), INTENT(in) :: x !! X coordinate of the point to check
239	REAL(kind=wp), INTENT(in) :: y !! Y coordinate of the point to check
240	REAL(kind=wp), INTENT(in) :: z !! Z coordinate of the point to check
241	REAL(kind=wp), INTENT(in) :: t !! T coordinate of the point to check
242	REAL(kind=wp), INTENT(in) :: w !! W coordinate of the point to check
243	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
244	REAL(kind=wp) :: l,u
245	ret=.false.
246	l = set%x(1) ; u= set%x(size(set%x))
247	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
248	l = set%y(1) ; u= set%y(size(set%y))
249	IF ((y>l.AND.y>u).OR.(y<l.AND.y<u)) RETURN
250	l = set%z(1) ; u= set%z(size(set%z))
251	IF ((z>l.AND.z>u).OR.(z<l.AND.z<u)) RETURN
252	l = set%t(1) ; u= set%t(size(set%t))
253	IF ((t>l.AND.t>u).OR.(t<l.AND.t<u)) RETURN
254	l = set%w(1) ; u= set%w(size(set%w))
255	IF ((w>l.AND.w>u).OR.(w<l.AND.w<u)) RETURN
256	ret=.true.
257	RETURN
258	END FUNCTION is_in_5d
259
260	FUNCTION ascii_header(path,ds,dp) RESULT(ret)
261	!! Read ASCII file header
262	!!
263	!! The method assumes the header is on two lines :
264	!!
265	!! - the first line must contain a single value which is the number of
266	!! dimensions.
267	!! - the second must contain N values with the size of each dimensions (N
268	!! referring to the first line number).
269	!!
270	!! The file remains open in the logical unit 666 unless some error occured.
271	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
272	INTEGER, INTENT(out), DIMENSION(:) :: ds !! Size of each dimensions
273	INTEGER, INTENT(out), DIMENSION(:) :: dp !! Product of each lower dimension size (1 for 1st)
274	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
275	INTEGER :: nd,i,e
276	CHARACTER(len=15) :: i2s
277	INQUIRE(666,OPENED=ret)
278	IF (ret) THEN
279	WRITE(,) 'ERROR: LUN 666 already used...'
280	ret = .false. ; RETURN
281	ENDIF
282	ret = .false.
283	OPEN(666,file=TRIM(path))
284	READ(666,*) nd
285	IF (nd /= SIZE(ds)) THEN
286	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
287	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be "//TRIM(i2s)//"D"
288	RETURN
289	ENDIF
290	READ(666,*,iostat=e) ds
291	IF (e /= 0) THEN
292	WRITE(*,'(a)') 'ERROR: Cannot get dimensions size'
293	CLOSE(666) ; RETURN
294	ENDIF
295	dp(1) = 1
296	DO i=2,nd
297	dp(i)=PRODUCT(ds(:i-1))
298	ENDDO
299	ret = .true.
300	END FUNCTION ascii_header
301
302	#if HAVE_NC_FTN
303	FUNCTION nc_get_dim(fid,did,var,values,ds) RESULT(ret)
304	!! Get dimensions information of a NetCDF variable
305	!!
306	!! The method gets the values and size of a given dimension.
307	INTEGER, INTENT(in) :: fid
308	!! Id of the NetCDF file
309	INTEGER, INTENT(in) :: did
310	!! Id of the NetCDF dimension
311	CHARACTER(len=*), INTENT(in) :: var
312	!! Name of the related NetCDF variable
313	REAL(kind=wp), INTENT(out), DIMENSION(:), ALLOCATABLE :: values
314	!! Values of the dimension
315	INTEGER, INTENT(out) :: ds
316	!! Size of __values__
317	LOGICAL :: ret
318	!! .true. if no error(s) occured, .false. otherwise
319	INTEGER :: vid,err
320	CHARACTER(len=NF90_MAX_NAME) :: dn
321	CHARACTER(len=15) :: i2s
322	ret = .false.
323	! --- Get dimension informations
324	IF (NF90_INQUIRE_DIMENSION(fid,did,dn,ds) /= NF90_NOERR) THEN
325	IF (debug) THEN
326	WRITE(i2s,*) did ; i2s=TRIM(i2s)
327	WRITE(*,'(a)') "ERROR:"//TRIM(var)//": Cannot find "// &
328	"dimension #"//TRIM(i2s)//" name"
329	ENDIF
330	err = NF90_CLOSE(fid) ; RETURN
331	ENDIF
332	IF (ds < 2) THEN
333	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(dn)//": Invalid dimension "// &
334	"size (<2)"
335	err = NF90_CLOSE(fid) ; RETURN
336	ENDIF
337	IF (NF90_INQ_VARID(fid,TRIM(dn),vid) /= NF90_NOERR) THEN
338	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(dn)//": Cannot get ID"
339	err = NF90_CLOSE(fid) ; RETURN
340	ENDIF
341	ALLOCATE(values(ds))
342	IF (NF90_GET_VAR(fid,vid,values) /= NF90_NOERR) THEN
343	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(dn)//": Cannot get values"
344	err = NF90_CLOSE(fid) ; RETURN
345	ENDIF
346	ret = .true.
347	END FUNCTION nc_get_dim
348
349	FUNCTION get_nc3_info(path,variable,ncid,vid,dimids) RESULT(ret)
350	!! Get variable informations from NetCDF file
351	!!
352	!! The method attempts to read the given NetCDF file header and retrieves
353	!! variable's informations from it:
354	!!
355	!! - parent file ID
356	!! - variable ID
357	!! - variable's dimensions ID
358	!! - variable's dimensions sizes
359	!!
360	!! If no error occured (i.e. the function has returned .true.) then the
361	!! NetCDF file remains open. In other case, the file is closed.
362	CHARACTER(len=*), INTENT(in) :: path
363	!! Path of the NetCDF file
364	CHARACTER(len=*), INTENT(in) :: variable
365	!! Name of the NetCDF variable to extract
366	INTEGER, INTENT(out) :: ncid
367	!! NetCDF file ID
368	INTEGER, INTENT(out) :: vid
369	!! NetCDF Variable ID
370	INTEGER, INTENT(out), DIMENSION(:), ALLOCATABLE :: dimids
371	!! Id of the variable's dimensions. __dimids__ is not allocated on error.
372	LOGICAL :: ret
373	!! .true. if no errors occured, .false. otherwise
374	INTEGER :: fid,ty,nc,err
375	ret = .false.
376	! Opens file
377	IF (NF90_OPEN(TRIM(path),NF90_NOWRITE,fid) /= NF90_NOERR) THEN
378	WRITE(*,'(a)') 'ERROR: Cannot open '//trim(path)
379	RETURN
380	ENDIF
381	ncid = fid
382	! Searches for variable
383	IF (NF90_INQ_VARID(ncid,TRIM(variable),vid) /= NF90_NOERR) THEN
384	WRITE(*,'(a)') 'Cannot find '//TRIM(variable)//' in '//trim(path)
385	nc = NF90_CLOSE(fid)
386	RETURN
387	ENDIF
388	! Get variable infos
389	! 1st call to get type and number of dimensions)
390	IF (NF90_INQUIRE_VARIABLE(ncid,vid,xtype=ty,ndims=nc) /= NF90_NOERR) THEN
391	WRITE(*,'(a)') 'Cannot access to '//TRIM(variable)//' informations'
392	nc = NF90_CLOSE(fid)
393	RETURN
394	ELSE
395	! Checks type
396	IF (ty == NF90_CHAR) THEN
397	WRITE(*,'(a)') 'Inconsistent variable type (should be numeric)'
398	nc = NF90_CLOSE(fid)
399	RETURN
400	ENDIF
401	ALLOCATE(dimids(nc))
402	ENDIF
403	! Gets variable's dimensions informations
404	! first get dimensions id
405	IF (NF90_INQUIRE_VARIABLE(ncid,vid,dimids=dimids) /= NF90_NOERR) THEN
406	WRITE(*,'(a)') 'Cannot access to '//TRIM(variable)//' informations'
407	nc = NF90_CLOSE(fid)
408	DEALLOCATE(dimids)
409	RETURN
410	ENDIF
411	ret = .true.
412	END FUNCTION get_nc3_info
413
414	#if HAVE_NC4_FTN
415	FUNCTION get_nc4_info(path,variable,group,ncid,vid,dimids) RESULT(ret)
416	!! Get variable informations from NetCDF4 file
417	!!
418	!! The method attempts to read the given NetCDF file header and retrieves
419	!! variable's informations from it :
420	!!
421	!! - parent group/file ID
422	!! - variable ID
423	!! - variable's dimensions ID
424	!! - variable's dimensions sizes
425	!!
426	!! If no error occured (i.e. the function has returned .true.) then the
427	!! NetCDF file remains open. In other case, the file is closed.
428	CHARACTER(len=*), INTENT(in) :: path
429	!! Path of the NetCDF file
430	CHARACTER(len=*), INTENT(in) :: variable
431	!! Name of the NetCDF variable to extract
432	CHARACTER(len=*), INTENT(in) :: group
433	!! Fullname of the variable's parent group (should be set to empty string for root group)
434	INTEGER, INTENT(out) :: ncid
435	!! NetCDF group ID
436	INTEGER, INTENT(out) :: vid
437	!! NetCDF Variable ID
438	INTEGER, INTENT(out), DIMENSION(:), ALLOCATABLE :: dimids
439	!! Id of the variable's dimensions. On error, __dimids__ is not allocated.
440	LOGICAL :: ret
441	!! .true. if no errors occured, .false. otherwise
442	INTEGER :: fid,ty,nc,err
443	ret = .false.
444	! Opens file
445	IF (NF90_OPEN(TRIM(path),NF90_NOWRITE,fid) /= NF90_NOERR) THEN
446	WRITE(*,'(a)') 'ERROR: Cannot open '//trim(path)
447	RETURN
448	ENDIF
449	! Searches for group based on its fullname
450	IF (LEN_TRIM(group) == 0) THEN
451	nc = NF90_NOERR ; ncid = fid
452	ELSE
453	nc = NF90_INQ_GRP_FULL_NCID(fid,TRIM(group),ncid)
454	ENDIF
455	SELECT CASE(nc)
456	! NF90_ENOGRP is missing from Netcdf-Fortran4.2 : its value=-125
457	CASE(-125)
458	WRITE(*,'(a)') TRIM(group)//' does not exist in '//TRIM(path)
459	nc = NF90_CLOSE(fid) ; RETURN
460	CASE(NF90_ENOTNC4,NF90_ESTRICTNC3)
461	WRITE(*,'(a)') TRIM(path)//' is not a NetCDF-4 file with "group" feature'
462	nc = NF90_CLOSE(fid) ; RETURN
463	CASE(NF90_EHDFERR)
464	WRITE(*,'(a)') "Too bad, an HDF5 error has been reported..."
465	nc = NF90_CLOSE(fid) ; RETURN
466	END SELECT
467	! Searches for variable
468	IF (NF90_INQ_VARID(ncid,TRIM(variable),vid) /= NF90_NOERR) THEN
469	WRITE(*,'(a)') 'Cannot find '//TRIM(variable)//' in '//trim(path)
470	nc = NF90_CLOSE(fid)
471	RETURN
472	ENDIF
473	! Get variable infos
474	! 1st call to get type and number of dimensions)
475	IF (NF90_INQUIRE_VARIABLE(ncid,vid,xtype=ty,ndims=nc) /= NF90_NOERR) THEN
476	WRITE(*,'(a)') 'Cannot access to '//TRIM(variable)//' informations'
477	nc = NF90_CLOSE(fid)
478	RETURN
479	ELSE
480	! Checks type
481	IF (ty == NF90_CHAR) THEN
482	WRITE(*,'(a)') 'Inconsistent variable type (should be numeric)'
483	nc = NF90_CLOSE(fid)
484	RETURN
485	ENDIF
486	ALLOCATE(dimids(nc))
487	ENDIF
488	! Gets variable's dimensions informations
489	! first get dimensions id
490	IF (NF90_INQUIRE_VARIABLE(ncid,vid,dimids=dimids) /= NF90_NOERR) THEN
491	WRITE(*,'(a)') 'Cannot access to '//TRIM(variable)//' informations'
492	nc = NF90_CLOSE(fid)
493	DEALLOCATE(dimids)
494	RETURN
495	ENDIF
496	ret = .true.
497	END FUNCTION get_nc4_info
498	#endif
499	#endif
500
501	!-------------------------
502	! NetCDF data file readers
503	!-------------------------
504
505	#if HAVE_NC_FTN
506
507
508	FUNCTION ncdf_rd_1d(path,variable,set) RESULT(ret)
509	!! Read a 1D data set from a NetCDF file
510	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
511	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
512	TYPE(DSET1D), INTENT(out) :: set !! Output dataset object
513	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
514	INTEGER :: fi,vi,nd
515	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
516	CHARACTER(len=NF90_MAX_NAME) :: dn
517	CHARACTER(len=15) :: i2s
518	ret = .false.
519	! --- Reset the data set
520	CALL clear_dset(set)
521	! --- Check NetCDF file info
522	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
523	! --- Check dimension size
524	nd = SIZE(di)
525	IF (nd /= 1) THEN
526	IF (debug) THEN
527	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
528	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
529	ENDIF
530	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
531	ENDIF
532	! --- Get coordinates values
533	ALLOCATE(ds(nd))
534	! ------ X coordinate
535	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
536	CALL clear_dset(set) ; RETURN
537	ENDIF
538	! --- Read data
539	ALLOCATE(set%data(ds(1)))
540	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
541	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
542	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
543	ENDIF
544	nd = NF90_CLOSE(fi)
545	ret = .true.
546	RETURN
547	END FUNCTION ncdf_rd_1d
548
549	FUNCTION ncdf_rd_2d(path,variable,set) RESULT(ret)
550	!! Read a 2D data set from a NetCDF file
551	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
552	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
553	TYPE(DSET2D), INTENT(out) :: set !! Output dataset object
554	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
555	INTEGER :: fi,vi,nd
556	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
557	CHARACTER(len=15) :: i2s
558	ret = .false.
559	! --- Reset the data set
560	CALL clear_dset(set)
561	! --- Check NetCDF file info
562	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
563	! --- Check dimension size
564	nd = SIZE(di)
565	IF (nd /= 2) THEN
566	IF (debug) THEN
567	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
568	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
569	ENDIF
570	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
571	ENDIF
572	! --- Get dimensions informations
573	ALLOCATE(ds(nd))
574	! ------ X coordinate
575	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
576	CALL clear_dset(set) ; RETURN
577	ENDIF
578	! ------ Y coordinate
579	IF (.NOT.nc_get_dim(fi,di(2),variable,set%y,ds(2))) THEN
580	CALL clear_dset(set) ; RETURN
581	ENDIF
582	! --- Read data
583	ALLOCATE(set%data(ds(1),ds(2)))
584	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
585	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
586	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
587	ENDIF
588	nd = NF90_CLOSE(fi)
589	ret = .true.
590	RETURN
591	END FUNCTION ncdf_rd_2d
592
593	FUNCTION ncdf_rd_3d(path,variable,set) RESULT(ret)
594	!! Read a 3D data set from a NetCDF file
595	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
596	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
597	TYPE(DSET3D), INTENT(out) :: set !! Output dataset object
598	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
599	INTEGER :: fi,vi,nd
600	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
601	CHARACTER(len=15) :: i2s
602	ret = .false.
603	! --- Reset the data set
604	CALL clear_dset(set)
605	! --- Check NetCDF file info
606	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
607	! --- Check dimension size
608	nd = SIZE(di)
609	IF (nd /= 3) THEN
610	IF (debug) THEN
611	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
612	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
613	ENDIF
614	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
615	ENDIF
616	! --- Get dimensions informations
617	ALLOCATE(ds(nd))
618	! ------ X coordinate
619	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
620	CALL clear_dset(set) ; RETURN
621	ENDIF
622	! ------ Y coordinate
623	IF (.NOT.nc_get_dim(fi,di(2),variable,set%y,ds(2))) THEN
624	CALL clear_dset(set) ; RETURN
625	ENDIF
626	! ------ Z coordinate
627	IF (.NOT.nc_get_dim(fi,di(3),variable,set%z,ds(3))) THEN
628	CALL clear_dset(set) ; RETURN
629	ENDIF
630	! --- Read data
631	ALLOCATE(set%data(ds(1),ds(2),ds(3)))
632	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
633	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
634	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
635	ENDIF
636	nd = NF90_CLOSE(fi)
637	ret = .true.
638	RETURN
639	END FUNCTION ncdf_rd_3d
640
641	FUNCTION ncdf_rd_4d(path,variable,set) RESULT(ret)
642	!! Read a 4D data set from a NetCDF file
643	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
644	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
645	TYPE(DSET4D), INTENT(out) :: set !! Output dataset object
646	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
647	INTEGER :: fi,vi,nd
648	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
649	CHARACTER(len=15) :: i2s
650	ret = .false.
651	! --- Reset the data set
652	CALL clear_dset(set)
653	! --- Check NetCDF file info
654	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
655	! --- Check dimension size
656	nd = SIZE(di)
657	IF (nd /= 4) THEN
658	IF (debug) THEN
659	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
660	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
661	ENDIF
662	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
663	ENDIF
664	! --- Get dimensions informations
665	ALLOCATE(ds(nd))
666	! ------ X coordinate
667	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
668	CALL clear_dset(set) ; RETURN
669	ENDIF
670	! ------ Y coordinate
671	IF (.NOT.nc_get_dim(fi,di(2),variable,set%y,ds(2))) THEN
672	CALL clear_dset(set) ; RETURN
673	ENDIF
674	! ------ Z coordinate
675	IF (.NOT.nc_get_dim(fi,di(3),variable,set%z,ds(3))) THEN
676	CALL clear_dset(set) ; RETURN
677	ENDIF
678	! ------ T coordinate
679	IF (.NOT.nc_get_dim(fi,di(4),variable,set%t,ds(4))) THEN
680	CALL clear_dset(set) ; RETURN
681	ENDIF
682	! --- Read data
683	ALLOCATE(set%data(ds(1),ds(2),ds(3),ds(4)))
684	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
685	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
686	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
687	ENDIF
688	nd = NF90_CLOSE(fi)
689	ret = .true.
690	RETURN
691	END FUNCTION ncdf_rd_4d
692
693	FUNCTION ncdf_rd_5d(path,variable,set) RESULT(ret)
694	!! Read a 5D data set from a NetCDF file
695	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
696	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
697	TYPE(DSET5D), INTENT(out) :: set !! Output dataset object
698	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
699	INTEGER :: fi,vi,nd
700	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
701	CHARACTER(len=15) :: i2s
702	ret = .false.
703	! --- Reset the data set
704	CALL clear_dset(set)
705	! --- Check NetCDF file info
706	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
707	! --- Check dimension size
708	nd = SIZE(di)
709	IF (nd /= 5) THEN
710	IF (debug) THEN
711	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
712	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
713	ENDIF
714	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
715	ENDIF
716	! --- Get dimensions informations
717	ALLOCATE(ds(nd))
718	! ------ X coordinate
719	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
720	CALL clear_dset(set) ; RETURN
721	ENDIF
722	! ------ Y coordinate
723	IF (.NOT.nc_get_dim(fi,di(2),variable,set%y,ds(2))) THEN
724	CALL clear_dset(set) ; RETURN
725	ENDIF
726	! ------ Z coordinate
727	IF (.NOT.nc_get_dim(fi,di(3),variable,set%z,ds(3))) THEN
728	CALL clear_dset(set) ; RETURN
729	ENDIF
730	! ------ T coordinate
731	IF (.NOT.nc_get_dim(fi,di(4),variable,set%t,ds(4))) THEN
732	CALL clear_dset(set) ; RETURN
733	ENDIF
734	! ------ W coordinate
735	IF (.NOT.nc_get_dim(fi,di(5),variable,set%w,ds(5))) THEN
736	CALL clear_dset(set) ; RETURN
737	ENDIF
738	! --- Read data
739	ALLOCATE(set%data(ds(1),ds(2),ds(3),ds(4),ds(5)))
740	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
741	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
742	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
743	ENDIF
744	nd = NF90_CLOSE(fi)
745	ret = .true.
746	RETURN
747	END FUNCTION ncdf_rd_5d
748
749	! NetCDF4 with groups !
750	#if HAVE_NC4_FTN
751
752	FUNCTION ncdf4_rd_1d(path,group,variable,set) RESULT(ret)
753	!! Read a 1D data set from a NetCDF4 file
754	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
755	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
756	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
757	TYPE(DSET1D), INTENT(out) :: set !! Output dataset
758	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
759	INTEGER :: fi,vi,nd
760	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
761	CHARACTER(len=NF90_MAX_NAME) :: dn
762	CHARACTER(len=15) :: i2s
763	ret = .false.
764	! --- Reset the data set
765	CALL clear_dset(set)
766	! --- Check NetCDF file info
767	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) RETURN
768	! --- Check dimension size
769	nd = SIZE(di)
770	IF (nd /= 1) THEN
771	IF (debug) THEN
772	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
773	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
774	ENDIF
775	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
776	ENDIF
777	! --- Get coordinates values
778	ALLOCATE(ds(nd))
779	! ------ X coordinate
780	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
781	CALL clear_dset(set) ; RETURN
782	ENDIF
783	! --- Read data
784	ALLOCATE(set%data(ds(1)))
785	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
786	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
787	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
788	ENDIF
789	nd = NF90_CLOSE(fi)
790	ret = .true.
791	RETURN
792	END FUNCTION ncdf4_rd_1d
793
794	FUNCTION ncdf4_rd_2d(path,group,variable,set) RESULT(ret)
795	!! Read a 2D data set from a NetCDF4 file
796	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
797	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
798	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
799	TYPE(DSET2D), INTENT(out) :: set !! Output dataset
800	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
801	INTEGER :: fi,vi,nd
802	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
803	CHARACTER(len=15) :: i2s
804	ret = .false.
805	! --- Reset the data set
806	CALL clear_dset(set)
807	! --- Check NetCDF file info
808	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) rETURN
809	! --- Check dimension size
810	nd = SIZE(di)
811	IF (nd /= 2) THEN
812	IF (debug) THEN
813	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
814	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
815	ENDIF
816	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
817	ENDIF
818	! --- Get dimensions informations
819	ALLOCATE(ds(nd))
820	! ------ X coordinate
821	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
822	CALL clear_dset(set) ; RETURN
823	ENDIF
824	! ------ Y coordinate
825	IF (.NOT.nc_get_dim(fi,di(2),variable,set%y,ds(2))) THEN
826	CALL clear_dset(set) ; RETURN
827	ENDIF
828	! --- Read data
829	ALLOCATE(set%data(ds(1),ds(2)))
830	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
831	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
832	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
833	ENDIF
834	nd = NF90_CLOSE(fi)
835	ret = .true.
836	RETURN
837	END FUNCTION ncdf4_rd_2d
838
839	FUNCTION ncdf4_rd_3d(path,group,variable,set) RESULT(ret)
840	!! Read a 3D data set from a NetCDF4 file
841	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
842	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
843	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
844	TYPE(DSET3D), INTENT(out) :: set !! Output dataset
845	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
846	INTEGER :: fi,vi,nd
847	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
848	CHARACTER(len=15) :: i2s
849	ret = .false.
850	! --- Reset the data set
851	CALL clear_dset(set)
852	! --- Check NetCDF file info
853	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) RETURN
854	! --- Check dimension size
855	nd = SIZE(di)
856	IF (nd /= 3) THEN
857	IF (debug) THEN
858	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
859	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
860	ENDIF
861	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
862	ENDIF
863	! --- Get dimensions informations
864	ALLOCATE(ds(nd))
865	! ------ X coordinate
866	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
867	CALL clear_dset(set) ; RETURN
868	ENDIF
869	! ------ Y coordinate
870	IF (.NOT.nc_get_dim(fi,di(2),variable,set%y,ds(2))) THEN
871	CALL clear_dset(set) ; RETURN
872	ENDIF
873	! ------ Z coordinate
874	IF (.NOT.nc_get_dim(fi,di(3),variable,set%z,ds(3))) THEN
875	CALL clear_dset(set) ; RETURN
876	ENDIF
877	! --- Read data
878	ALLOCATE(set%data(ds(1),ds(2),ds(3)))
879	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
880	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
881	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
882	ENDIF
883	nd = NF90_CLOSE(fi)
884	ret = .true.
885	RETURN
886	END FUNCTION ncdf4_rd_3d
887
888	FUNCTION ncdf4_rd_4d(path,group,variable,set) RESULT(ret)
889	!! Read a 4D data set from a NetCDF4 file
890	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
891	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
892	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
893	TYPE(DSET4D), INTENT(out) :: set !! Output dataset
894	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
895	INTEGER :: fi,vi,nd
896	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
897	CHARACTER(len=15) :: i2s
898	ret = .false.
899	! --- Reset the data set
900	CALL clear_dset(set)
901	! --- Check NetCDF file info
902	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) RETURN
903	! --- Check dimension size
904	nd = SIZE(di)
905	IF (nd /= 4) THEN
906	IF (debug) THEN
907	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
908	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
909	ENDIF
910	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
911	ENDIF
912	! --- Get dimensions informations
913	ALLOCATE(ds(nd))
914	! ------ X coordinate
915	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
916	CALL clear_dset(set) ; RETURN
917	ENDIF
918	! ------ Y coordinate
919	IF (.NOT.nc_get_dim(fi,di(2),variable,set%y,ds(2))) THEN
920	CALL clear_dset(set) ; RETURN
921	ENDIF
922	! ------ Z coordinate
923	IF (.NOT.nc_get_dim(fi,di(3),variable,set%z,ds(3))) THEN
924	CALL clear_dset(set) ; RETURN
925	ENDIF
926	! ------ T coordinate
927	IF (.NOT.nc_get_dim(fi,di(4),variable,set%t,ds(4))) THEN
928	CALL clear_dset(set) ; RETURN
929	ENDIF
930	! --- Read data
931	ALLOCATE(set%data(ds(1),ds(2),ds(3),ds(4)))
932	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
933	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
934	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
935	ENDIF
936	nd = NF90_CLOSE(fi)
937	ret = .true.
938	RETURN
939	END FUNCTION ncdf4_rd_4d
940
941	FUNCTION ncdf4_rd_5d(path,group,variable,set) RESULT(ret)
942	!! Read a 5D data set from a NetCDF4 file
943	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
944	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
945	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
946	TYPE(DSET5D), INTENT(out) :: set !! Output dataset
947	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
948	INTEGER :: fi,vi,nd
949	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
950	CHARACTER(len=15) :: i2s
951	ret = .false.
952	! --- Reset the data set
953	CALL clear_dset(set)
954	! --- Check NetCDF file info
955	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) RETURN
956	! --- Check dimension size
957	nd = SIZE(di)
958	IF (nd /= 5) THEN
959	IF (debug) THEN
960	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
961	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
962	ENDIF
963	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
964	ENDIF
965	! --- Get dimensions informations
966	ALLOCATE(ds(nd))
967	! ------ X coordinate
968	IF (.NOT.nc_get_dim(fi,di(1),variable,set%x,ds(1))) THEN
969	CALL clear_dset(set) ; RETURN
970	ENDIF
971	! ------ Y coordinate
972	IF (.NOT.nc_get_dim(fi,di(2),variable,set%y,ds(2))) THEN
973	CALL clear_dset(set) ; RETURN
974	ENDIF
975	! ------ Z coordinate
976	IF (.NOT.nc_get_dim(fi,di(3),variable,set%z,ds(3))) THEN
977	CALL clear_dset(set) ; RETURN
978	ENDIF
979	! ------ T coordinate
980	IF (.NOT.nc_get_dim(fi,di(4),variable,set%t,ds(4))) THEN
981	CALL clear_dset(set) ; RETURN
982	ENDIF
983	! ------ W coordinate
984	IF (.NOT.nc_get_dim(fi,di(5),variable,set%w,ds(5))) THEN
985	CALL clear_dset(set) ; RETURN
986	ENDIF
987	! --- Read data
988	ALLOCATE(set%data(ds(1),ds(2),ds(3),ds(4),ds(5)))
989	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
990	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
991	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
992	ENDIF
993	nd = NF90_CLOSE(fi)
994	ret = .true.
995	RETURN
996	END FUNCTION ncdf4_rd_5d
997	#endif
998	#endif
999
1000	!------------------------
1001	! ASCII data file readers
1002	!------------------------
1003
1004	FUNCTION ascii_rd_1d(path,set) RESULT(ret)
1005	!! Read a 1D data set from a ASCII file
1006	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1007	TYPE(dset1d), INTENT(out) :: set !! output 1D dataset
1008	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1009	INTEGER :: i,e
1010	REAL(kind=wp), DIMENSION(2) :: vl
1011	INTEGER, DIMENSION(1) :: cc,ds,dp
1012	ret = .false.
1013	CALL clear_dset(set)
1014	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1015	ALLOCATE(set%x(ds(1)), &
1016	set%data(ds(1)))
1017	cc(:) = 1
1018	DO i=1,ds(1)
1019	READ(666,*) vl(:)
1020	set%x(cc(1)) = vl(1)
1021	set%data(cc(1)) = vl(2)
1022	cc(1) = cc(1) + 1
1023	ENDDO
1024	READ(666,*,iostat=e) vl(1)
1025	IF (e == 0) THEN
1026	IF (debug) WRITE(,) 'ERROR: Extra value found...'
1027	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1028	ENDIF
1029	ret = .true.
1030	CLOSE(666)
1031	RETURN
1032	END FUNCTION ascii_rd_1d
1033
1034	FUNCTION ascii_rd_2d(path,set) RESULT(ret)
1035	!! Read a 2D data set from a ASCII file
1036	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1037	TYPE(dset2d), INTENT(out) :: set !! output 2D dataset
1038	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1039	INTEGER :: i,e
1040	REAL(kind=wp), DIMENSION(3) :: vl
1041	INTEGER, DIMENSION(2) :: cc,ds,dp
1042	ret = .false.
1043	CALL clear_dset(set)
1044	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1045	ALLOCATE(set%x(ds(1)),set%y(ds(2)), &
1046	set%data(ds(1),ds(2)))
1047	cc(:) = 1
1048	DO i=1,PRODUCT(ds)
1049	READ(666,*,iostat=e) vl ! Reads the line
1050	IF (e /= 0) THEN
1051	IF (debug) WRITE(*,'(a)') 'ERROR: Malformed file, line: ',i+2
1052	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1053	ENDIF
1054	! --- X coordinate
1055	set%x(cc(1)) = vl(1)
1056	! --- Y coordinate
1057	set%y(cc(2)) = vl(2)
1058	! --- Data
1059	set%data(cc(1),cc(2)) = vl(3)
1060	! - Update counters
1061	IF (mod(i,dp(1))==0) cc(1) = cc(1)+1 ; IF (cc(1) > ds(1)) cc(1)=1
1062	IF (mod(i,dp(2))==0) cc(2) = cc(2)+1 ; IF (cc(2) > ds(2)) cc(2)=1
1063	ENDDO
1064	READ(666,*,iostat=e) vl(1)
1065	IF (e == 0) THEN
1066	IF (debug) WRITE(*,'(a)') 'ERROR: Extra value found'
1067	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1068	ENDIF
1069	ret = .true.
1070	CLOSE(666)
1071	RETURN
1072	END FUNCTION ascii_rd_2d
1073
1074	FUNCTION ascii_rd_3d(path,set) RESULT(ret)
1075	!! Read a 3D data set from a ASCII file
1076	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1077	TYPE(dset3d), INTENT(out) :: set !! output 3D dataset
1078	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1079	INTEGER :: i,e
1080	REAL(kind=wp), DIMENSION(4) :: vl
1081	INTEGER, DIMENSION(3) :: cc,ds,dp
1082	ret = .false.
1083	CALL clear_dset(set)
1084	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1085	ALLOCATE(set%x(ds(1)),set%y(ds(2)),set%z(ds(3)), &
1086	set%data(ds(1),ds(2),ds(3)))
1087	cc(:) = 1
1088	DO i=1,PRODUCT(ds)
1089	READ(666,*,iostat=e) vl ! Reads the line
1090	IF (e /= 0) THEN
1091	IF (debug) WRITE(*,'(a)') 'ERROR: Malformed file, line: ',i+2
1092	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1093	ENDIF
1094	! --- X coordinate
1095	set%x(cc(1)) = vl(1)
1096	! --- Y coordinate
1097	set%y(cc(2)) = vl(2)
1098	! --- Z coordinate
1099	set%z(cc(3)) = vl(3)
1100	! --- Data
1101	set%data(cc(1),cc(2),cc(3)) = vl(4)
1102	! - Update counters
1103	IF (mod(i,dp(1))==0) cc(1) = cc(1)+1 ; IF (cc(1) > ds(1)) cc(1)=1
1104	IF (mod(i,dp(2))==0) cc(2) = cc(2)+1 ; IF (cc(2) > ds(2)) cc(2)=1
1105	IF (mod(i,dp(3))==0) cc(3) = cc(3)+1 ; IF (cc(3) > ds(3)) cc(3)=1
1106	ENDDO
1107	READ(666,*,iostat=e) vl(1)
1108	IF (e == 0) THEN
1109	IF (debug) WRITE(*,'(a)') 'ERROR: Extra value found'
1110	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1111	ENDIF
1112	ret = .true.
1113	CLOSE(666)
1114	RETURN
1115	END FUNCTION ascii_rd_3d
1116
1117	FUNCTION ascii_rd_4d(path,set) RESULT(ret)
1118	!! Read a 4D data set from a ASCII file
1119	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1120	TYPE(dset4d), INTENT(out) :: set !! output 4D dataset
1121	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1122	INTEGER :: i,e
1123	REAL(kind=wp), DIMENSION(5) :: vl
1124	INTEGER, DIMENSION(4) :: cc,ds,dp
1125	ret = .false.
1126	CALL clear_dset(set)
1127	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1128	ALLOCATE(set%x(ds(1)),set%y(ds(2)),set%z(ds(3)),set%t(ds(4)), &
1129	set%data(ds(1),ds(2),ds(3),ds(4)))
1130	cc(:) = 1
1131	DO i=1,PRODUCT(ds)
1132	READ(666,*,iostat=e) vl ! Reads the line
1133	IF (e /= 0) THEN
1134	IF (debug) WRITE(*,'(a)') 'ERROR: Malformed file, line: ',i+2
1135	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1136	ENDIF
1137	! --- X coordinate
1138	set%x(cc(1)) = vl(1)
1139	! --- Y coordinate
1140	set%y(cc(2)) = vl(2)
1141	! --- Z coordinate
1142	set%z(cc(3)) = vl(3)
1143	! --- T coordinate
1144	set%t(cc(4)) = vl(4)
1145	! --- Data
1146	set%data(cc(1),cc(2),cc(3),cc(4)) = vl(5)
1147	! - Update counters
1148	IF (mod(i,dp(1))==0) cc(1) = cc(1)+1 ; IF (cc(1) > ds(1)) cc(1)=1
1149	IF (mod(i,dp(2))==0) cc(2) = cc(2)+1 ; IF (cc(2) > ds(2)) cc(2)=1
1150	IF (mod(i,dp(3))==0) cc(3) = cc(3)+1 ; IF (cc(3) > ds(3)) cc(3)=1
1151	IF (mod(i,dp(4))==0) cc(4) = cc(4)+1 ; IF (cc(4) > ds(4)) cc(4)=1
1152	ENDDO
1153	READ(666,*,iostat=e) vl(1)
1154	IF (e == 0) THEN
1155	IF (debug) WRITE(*,'(a)') 'ERROR: Extra value found'
1156	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1157	ENDIF
1158	CLOSE(666)
1159	ret = .true.
1160	RETURN
1161	END FUNCTION ascii_rd_4d
1162
1163	FUNCTION ascii_rd_5d(path,set) RESULT(ret)
1164	!! Read a 5D data set from a ASCII file
1165	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1166	TYPE(dset5d), INTENT(out) :: set !! output 5D dataset
1167	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1168	INTEGER :: i,e
1169	REAL(kind=wp), DIMENSION(6) :: vl
1170	INTEGER, DIMENSION(5) :: cc,ds,dp
1171	ret = .false.
1172	CALL clear_dset(set)
1173	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1174	ALLOCATE(set%x(ds(1)),set%y(ds(2)),set%z(ds(3)),set%t(ds(4)),set%w(ds(5)), &
1175	set%data(ds(1),ds(2),ds(3),ds(4),ds(5)))
1176	cc(:) = 1
1177	DO i=1,PRODUCT(ds)
1178	READ(666,*,iostat=e) vl ! Reads the line
1179	IF (e /= 0) THEN
1180	IF (debug) WRITE(*,'(a)') 'ERROR: Malformed file, line: ',i+2
1181	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1182	ENDIF
1183	! --- X coordinate
1184	set%x(cc(1)) = vl(1)
1185	! --- Y coordinate
1186	set%y(cc(2)) = vl(2)
1187	! --- Z coordinate
1188	set%z(cc(3)) = vl(3)
1189	! --- T coordinate
1190	set%t(cc(4)) = vl(4)
1191	! --- W coordinate
1192	set%w(cc(5)) = vl(5)
1193	! --- Data
1194	set%data(cc(1),cc(2),cc(3),cc(4),cc(5)) = vl(6)
1195	! - Update counters
1196	IF (mod(i,dp(1)) == 0) cc(1)=cc(1)+1 ; IF (cc(1) > ds(1)) cc(1)=1
1197	IF (mod(i,dp(2)) == 0) cc(2)=cc(2)+1 ; IF (cc(2) > ds(2)) cc(2)=1
1198	IF (mod(i,dp(3)) == 0) cc(3)=cc(3)+1 ; IF (cc(3) > ds(3)) cc(3)=1
1199	IF (mod(i,dp(4)) == 0) cc(4)=cc(4)+1 ; IF (cc(4) > ds(4)) cc(4)=1
1200	IF (mod(i,dp(5)) == 0) cc(5)=cc(5)+1 ; IF (cc(5) > ds(5)) cc(5)=1
1201	ENDDO
1202	READ(666,*,iostat=e) vl(1)
1203	IF (e == 0) THEN
1204	IF (debug) WRITE(*,'(a)') 'ERROR: Extra value found'
1205	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1206	ENDIF
1207	CLOSE(666)
1208	ret = .true.
1209	RETURN
1210	END FUNCTION ascii_rd_5d
1211
1212	SUBROUTINE clr_1d_set(set)
1213	!! Clear the given 1D data set
1214	TYPE(DSET1D), INTENT(inout) :: set !! dataset object to clear
1215	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1216	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1217	RETURN
1218	END SUBROUTINE clr_1d_set
1219
1220	SUBROUTINE clr_2d_set(set)
1221	!! Clear the given 2D data set
1222	TYPE(DSET2D), INTENT(inout) :: set !! dataset object to clear
1223	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1224	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1225	IF (ALLOCATED(set%y)) DEALLOCATE(set%y)
1226	RETURN
1227	END SUBROUTINE clr_2d_set
1228
1229	SUBROUTINE clr_3d_set(set)
1230	!! Clear the given 3D data set
1231	TYPE(DSET3D), INTENT(inout) :: set !! dataset object to clear
1232	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1233	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1234	IF (ALLOCATED(set%y)) DEALLOCATE(set%y)
1235	IF (ALLOCATED(set%z)) DEALLOCATE(set%z)
1236	RETURN
1237	END SUBROUTINE clr_3d_set
1238
1239	SUBROUTINE clr_4d_set(set)
1240	!! Clear the given 4D data set
1241	TYPE(DSET4D), INTENT(inout) :: set !! dataset object to clear
1242	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1243	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1244	IF (ALLOCATED(set%y)) DEALLOCATE(set%y)
1245	IF (ALLOCATED(set%z)) DEALLOCATE(set%z)
1246	IF (ALLOCATED(set%t)) DEALLOCATE(set%t)
1247	RETURN
1248	END SUBROUTINE clr_4d_set
1249
1250	SUBROUTINE clr_5d_set(set)
1251	!! Clear the given 5D data set
1252	TYPE(DSET5D), INTENT(inout) :: set !! dataset object to clear
1253	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1254	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1255	IF (ALLOCATED(set%y)) DEALLOCATE(set%y)
1256	IF (ALLOCATED(set%z)) DEALLOCATE(set%z)
1257	IF (ALLOCATED(set%t)) DEALLOCATE(set%t)
1258	IF (ALLOCATED(set%w)) DEALLOCATE(set%w)
1259	RETURN
1260	END SUBROUTINE clr_5d_set
1261
1262	END MODULE DATASETS

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

Download in other formats:

Original Format