Context Navigation

datasets.F90 @ 3094

Last change on this file since 3094 was 1897, checked in by jvatant, 7 years ago

Making Titan's hazy again - part II
+ Major updates of J.Burgalat YAMMS library and optical coupling, including :
++ Added the routines for haze optics inside YAMMS
++ Calling rad. transf. with interactive haze is plugged
in but should stay unactive as long as the microphysics is
in test phase : cf "uncoupl_optic_haze" flag : true for now !
++ Also some sanity checks for negative tendencies and
some others upkeep of YAMMS model
+ Also added a temporary CPP key USE_QTEST in physiq_mod
that enables to have microphysical tendencies separated
from dynamics for debugging and test phases
-- JVO and JB

File size: 68.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,2017
39
40	MODULE DATASETS
41	!! dataset definitions module
42	!!
43	!! This module defines simple derived types that encapsulate data set variables. For a N-dimension
44	!! dataset, N vectors of \(n_{i}\) elements and a single N-dimensional array of
45	!! \(\prod_{i}^{N} n_{i}\) elements are defined.
46	!!
47	!! If the data set is to be used within [[lintdset(module)]] module then each coordinate values of the
48	!! dataset must be sorted either in ascending or descending order with no duplicates. The module does
49	!! not provide functionnality to sort and to check such requirements.
50	!!
51	!! The module also provides two interfaces to initialize data sets from input data file which can be
52	!! a NetCDF file (if NetCDF support is enabled at compilation time) or an ASCII file. In the latter
53	!! case, the file must contain a header which is formatted as follows:
54	!!
55	!! - The first line must contain only one value which is the number of coordinates (__N__)
56	!! - The second line must contain all the dimension sizes (that is __N__ values)
57	!! - Each other lines should contain __N+1__ columns with, for the __N__ first columns, the
58	!! coordinates values and finally the point value in the last column.
59	!!
60	!! @note
61	!! Note that for ASCII file, data must be ordered so first dimensions vary first. Same requirement
62	!! is needed for NetCDF file but in most cases, it is implicitly done (if dimensions are ordered).
63	USE NETCDF
64	IMPLICIT NONE
65
66	PRIVATE
67
68	PUBLIC :: read_dset, write_dset, clear_dset, is_in, has_data, debug
69
70	LOGICAL :: debug = .false. !! A control flag to enable verbose mode
71
72	!> Initialize a data set from either an ASCII or a NetCDF file
73	!!
74	!! Whatever the kind of file, this interface assumes that you know the dimensionnality of the
75	!! extracted variable. If file content is not compatible with the kind of data set given as
76	!! output argument, or if some I/O error occured, the dataset is cleared.
77	!! @note
78	!! Netcdf reader interface is available only if the library has been compiled with NetCDF support.
79	INTERFACE read_dset
80	MODULE PROCEDURE ncdf_rd_1d,ncdf_rd_2d,ncdf_rd_3d,ncdf_rd_4d,ncdf_rd_5d
81	#if HAVE_NC4_FTN
82	MODULE PROCEDURE ncdf4_rd_1d,ncdf4_rd_2d,ncdf4_rd_3d,ncdf4_rd_4d,ncdf4_rd_5d
83	#endif
84	MODULE PROCEDURE ascii_rd_1d,ascii_rd_2d,ascii_rd_3d,ascii_rd_4d,ascii_rd_5d
85	END INTERFACE
86
87	!> Write a dataset in a netcdf (classic) file.
88	INTERFACE write_dset
89	MODULE PROCEDURE nc_wr_1d,nc_wr_2d,nc_wr_3d,nc_wr_4d,nc_wr_5d
90	END INTERFACE
91
92	!> Clear the given data set
93	INTERFACE clear_dset
94	MODULE PROCEDURE clr_1d_set,clr_2d_set,clr_3d_set,clr_4d_set,clr_5d_set
95	END INTERFACE
96
97	!> Check if given point is within the data set
98	INTERFACE is_in
99	MODULE PROCEDURE is_in_1d,is_in_2d,is_in_3d,is_in_4d,is_in_5d
100	END INTERFACE
101
102	!> Private interface to netcdf informations getters
103	INTERFACE get_nc_info
104	MODULE PROCEDURE get_nc3_info
105	#if HAVE_NC4_FTN
106	MODULE PROCEDURE get_nc4_info
107	#endif
108	END INTERFACE
109
110	INTERFACE has_data
111	MODULE PROCEDURE has_d_1d,has_d_2d,has_d_3d,has_d_4d,has_d_5d
112	END INTERFACE
113
114
115	TYPE, PUBLIC :: DSET1D
116	!! A 1D data set
117	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
118	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
119	CHARACTER(len=NF90_MAX_NAME) :: xname = "X" !! Name of the X coordinate
120	CHARACTER(len=NF90_MAX_NAME) :: dname = "data" !! Name of the data block
121	END TYPE DSET1D
122
123	TYPE, PUBLIC :: DSET2D
124	!! A 2D data set
125	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
126	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: y !! Y coordinate tabulated values
127	REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
128	CHARACTER(len=NF90_MAX_NAME) :: xname = "X" !! Name of the X coordinate
129	CHARACTER(len=NF90_MAX_NAME) :: yname = "Y" !! Name of the Y coordinate
130	CHARACTER(len=NF90_MAX_NAME) :: dname = "data" !! Name of the data block
131	END TYPE DSET2D
132
133	TYPE, PUBLIC :: DSET3D
134	!! A 3D data set
135	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
136	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: y !! Y coordinate tabulated values
137	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: z !! Z coordinate tabulated values
138	REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
139	CHARACTER(len=NF90_MAX_NAME) :: xname = "X" !! Name of the X coordinate
140	CHARACTER(len=NF90_MAX_NAME) :: yname = "Y" !! Name of the Y coordinate
141	CHARACTER(len=NF90_MAX_NAME) :: zname = "Z" !! Name of the Z coordinate
142	CHARACTER(len=NF90_MAX_NAME) :: dname = "data" !! Name of the data block
143	END TYPE DSET3D
144
145	TYPE, PUBLIC :: DSET4D
146	!! A 4D data set
147	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
148	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: y !! Y coordinate tabulated values
149	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: z !! Z coordinate tabulated values
150	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: t !! T coordinate tabulated values
151	REAL(kind=8), DIMENSION(:,:,:,:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
152	CHARACTER(len=NF90_MAX_NAME) :: xname = "X" !! Name of the X coordinate
153	CHARACTER(len=NF90_MAX_NAME) :: yname = "Y" !! Name of the Y coordinate
154	CHARACTER(len=NF90_MAX_NAME) :: zname = "Z" !! Name of the Z coordinate
155	CHARACTER(len=NF90_MAX_NAME) :: tname = "T" !! Name of the T coordinate
156	CHARACTER(len=NF90_MAX_NAME) :: dname = "data" !! Name of the data block
157	END TYPE DSET4D
158
159	TYPE, PUBLIC :: DSET5D
160	!! A 5D data set
161	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: x !! X coordinate tabulated values
162	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: y !! Y coordinate tabulated values
163	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: z !! Z coordinate tabulated values
164	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: t !! T coordinate tabulated values
165	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: w !! W coordinate tabulated values
166	REAL(kind=8), DIMENSION(:,:,:,:,:), ALLOCATABLE :: data !! Tabulated function's value at each coordinate
167	CHARACTER(len=NF90_MAX_NAME) :: xname = "X" !! Name of the X coordinate
168	CHARACTER(len=NF90_MAX_NAME) :: yname = "Y" !! Name of the Y coordinate
169	CHARACTER(len=NF90_MAX_NAME) :: zname = "Z" !! Name of the Z coordinate
170	CHARACTER(len=NF90_MAX_NAME) :: tname = "T" !! Name of the T coordinate
171	CHARACTER(len=NF90_MAX_NAME) :: wname = "W" !! Name of the W coordinate
172	CHARACTER(len=NF90_MAX_NAME) :: dname = "data" !! Name of the data block
173	END TYPE DSET5D
174
175	CONTAINS
176
177
178	FUNCTION is_in_1d(set,x) RESULT(ret)
179	!! Check if point is in the 1D data set
180	TYPE(DSET1D), INTENT(in) :: set !! Dataset object to search in
181	REAL(kind=8), INTENT(in) :: x !! coordinate of the point to check
182	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
183	REAL(kind=8) :: l,u
184	ret=.true.
185	l = set%x(1) ; u= set%x(size(set%x))
186	IF ((x>=l.EQV.x<=u).OR.(x<=l.EQV.x>=u)) ret=.false.
187	RETURN
188	END FUNCTION is_in_1d
189
190	FUNCTION is_in_2d(set,x,y) RESULT(ret)
191	!! Check if point is in the 2D data set
192	TYPE(DSET2D), INTENT(in) :: set !! Dataset object to search in
193	REAL(kind=8), INTENT(in) :: x !! X coordinate of the point to check
194	REAL(kind=8), INTENT(in) :: y !! Y coordinate of the point to check
195	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
196	REAL(kind=8) :: l,u
197	ret=.false.
198	l = set%x(1) ; u= set%x(size(set%x))
199	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
200	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
201	l = set%y(1) ; u= set%y(size(set%y))
202	IF ((y>l.AND.y>u).OR.(y<l.AND.y<u)) RETURN
203	ret=.true.
204	RETURN
205	END FUNCTION is_in_2d
206
207	FUNCTION is_in_3d(set,x,y,z) RESULT(ret)
208	!! Check if point is in the 3D data set
209	TYPE(DSET3D), INTENT(in) :: set !! Dataset object to search in
210	REAL(kind=8), INTENT(in) :: x !! X coordinate of the point to check
211	REAL(kind=8), INTENT(in) :: y !! Y coordinate of the point to check
212	REAL(kind=8), INTENT(in) :: z !! Z coordinate of the point to check
213	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
214	REAL(kind=8) :: l,u
215	ret=.false.
216	l = set%x(1) ; u= set%x(size(set%x))
217	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
218	l = set%y(1) ; u= set%y(size(set%y))
219	IF ((y>l.AND.y>u).OR.(y<l.AND.y<u)) RETURN
220	l = set%z(1) ; u= set%z(size(set%z))
221	IF ((z>l.AND.z>u).OR.(z<l.AND.z<u)) RETURN
222	ret=.true.
223	RETURN
224	END FUNCTION is_in_3d
225
226	FUNCTION is_in_4d(set,x,y,z,t) RESULT(ret)
227	!! Check if point is in the 4D data set
228	TYPE(DSET4D), INTENT(in) :: set !! Dataset object to search in
229	REAL(kind=8), INTENT(in) :: x !! X coordinate of the point to check
230	REAL(kind=8), INTENT(in) :: y !! Y coordinate of the point to check
231	REAL(kind=8), INTENT(in) :: z !! Z coordinate of the point to check
232	REAL(kind=8), INTENT(in) :: t !! T coordinate of the point to check
233	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
234	REAL(kind=8) :: l,u
235	ret=.false.
236	l = set%x(1) ; u= set%x(size(set%x))
237	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
238	l = set%y(1) ; u= set%y(size(set%y))
239	IF ((y>l.AND.y>u).OR.(y<l.AND.y<u)) RETURN
240	l = set%z(1) ; u= set%z(size(set%z))
241	IF ((z>l.AND.z>u).OR.(z<l.AND.z<u)) RETURN
242	l = set%t(1) ; u= set%t(size(set%t))
243	IF ((t>l.AND.t>u).OR.(t<l.AND.t<u)) RETURN
244	ret=.true.
245	RETURN
246	END FUNCTION is_in_4d
247
248	FUNCTION is_in_5d(set,x,y,z,t,w) RESULT(ret)
249	!! Check if point is in the 4D data set
250	TYPE(DSET5D), INTENT(in) :: set !! Dataset object to search in
251	REAL(kind=8), INTENT(in) :: x !! X coordinate of the point to check
252	REAL(kind=8), INTENT(in) :: y !! Y coordinate of the point to check
253	REAL(kind=8), INTENT(in) :: z !! Z coordinate of the point to check
254	REAL(kind=8), INTENT(in) :: t !! T coordinate of the point to check
255	REAL(kind=8), INTENT(in) :: w !! W coordinate of the point to check
256	LOGICAL :: ret !! .true. if the point is in the data set, .false. otherwise
257	REAL(kind=8) :: l,u
258	ret=.false.
259	l = set%x(1) ; u= set%x(size(set%x))
260	IF ((x>l.AND.x>u).OR.(x<l.AND.x<u)) RETURN
261	l = set%y(1) ; u= set%y(size(set%y))
262	IF ((y>l.AND.y>u).OR.(y<l.AND.y<u)) RETURN
263	l = set%z(1) ; u= set%z(size(set%z))
264	IF ((z>l.AND.z>u).OR.(z<l.AND.z<u)) RETURN
265	l = set%t(1) ; u= set%t(size(set%t))
266	IF ((t>l.AND.t>u).OR.(t<l.AND.t<u)) RETURN
267	l = set%w(1) ; u= set%w(size(set%w))
268	IF ((w>l.AND.w>u).OR.(w<l.AND.w<u)) RETURN
269	ret=.true.
270	RETURN
271	END FUNCTION is_in_5d
272
273	FUNCTION ascii_header(path,ds,dp) RESULT(ret)
274	!! Read ASCII file header
275	!!
276	!! The method assumes the header is on two lines :
277	!!
278	!! - the first line must contain a single value which is the number of
279	!! dimensions.
280	!! - the second must contain N values with the size of each dimensions (N
281	!! referring to the first line number).
282	!!
283	!! The file remains open in the logical unit 666 unless some error occured.
284	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
285	INTEGER, INTENT(out), DIMENSION(:) :: ds !! Size of each dimensions
286	INTEGER, INTENT(out), DIMENSION(:) :: dp !! Product of each lower dimension size (1 for 1st)
287	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
288	INTEGER :: nd,i,e
289	CHARACTER(len=15) :: i2s
290	INQUIRE(666,OPENED=ret)
291	IF (ret) THEN
292	WRITE(,) 'ERROR: LUN 666 already used...'
293	ret = .false. ; RETURN
294	ENDIF
295	ret = .false.
296	OPEN(666,file=TRIM(path))
297	READ(666,*) nd
298	IF (nd /= SIZE(ds)) THEN
299	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
300	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be "//TRIM(i2s)//"D"
301	RETURN
302	ENDIF
303	READ(666,*,iostat=e) ds
304	IF (e /= 0) THEN
305	WRITE(*,'(a)') 'ERROR: Cannot get dimensions size'
306	CLOSE(666) ; RETURN
307	ENDIF
308	dp(1) = 1
309	DO i=2,nd
310	dp(i)=PRODUCT(ds(:i-1))
311	ENDDO
312	ret = .true.
313	END FUNCTION ascii_header
314
315	FUNCTION nc_get_dim_by_name(fid,dn,values,ds,did) RESULT(ret)
316	!! Get informations about a dimension.
317	!!
318	!! The method gets the values and size of a given dimension.
319	!!
320	!! Errors occur if:
321	!!
322	!! - the given name is not a dimension name.
323	!! - the size of the dimension is less than 2 (the method assumes dimension has an afferent variable with values).
324	!! - the method can not retrieve the values of the afferent variable.
325	INTEGER, INTENT(in) :: fid
326	!! Id of the NetCDF file
327	CHARACTER(len=NF90_MAX_NAME), INTENT(in) :: dn
328	!! Name of the dimension.
329	REAL(kind=8), INTENT(out), DIMENSION(:), ALLOCATABLE :: values
330	!! Values of the dimension
331	INTEGER, INTENT(out) :: ds
332	!! Size of __values__
333	INTEGER, INTENT(out) :: did
334	!! Id of the NetCDF dimension
335	LOGICAL :: ret
336	!! .true. if no error(s) occured, .false. otherwise
337	INTEGER :: vid,err
338	CHARACTER(len=15) :: i2s
339	ret = .false.
340	! --- Get dimension informations
341	err = NF90_INQ_DIMID(fid,dn,did)
342	IF (err /= NF90_NOERR) RETURN
343	err = NF90_INQUIRE_DIMENSION(fid,did,len=ds)
344	IF (err /= NF90_NOERR) RETURN
345	IF (ds < 2) RETURN
346	err = NF90_INQ_VARID(fid,TRIM(dn),vid)
347	IF (err /= NF90_NOERR) RETURN
348	ALLOCATE(values(ds))
349	err = NF90_GET_VAR(fid,vid,values)
350	IF (err /= NF90_NOERR) RETURN
351	ret = .true.
352	END FUNCTION nc_get_dim_by_name
353
354	FUNCTION nc_get_dim_by_id(fid,did,values,ds,dn) RESULT(ret)
355	!! Get informations about a dimension.
356	!!
357	!! The method gets the values and size of a given dimension.
358	!!
359	!! Errors occur if:
360	!!
361	!! - the given id is not a dimension identifier.
362	!! - the size of the dimension is less than 2 (the method assumes dimension has an afferent variable with values).
363	!! - the method can not retrieve the values of the afferent variable.
364	INTEGER, INTENT(in) :: fid
365	!! Id of the NetCDF file
366	INTEGER, INTENT(in) :: did
367	!! Id of the NetCDF dimension
368	REAL(kind=8), INTENT(out), DIMENSION(:), ALLOCATABLE :: values
369	!! Values of the dimension
370	INTEGER, INTENT(out) :: ds
371	!! Size of __values__
372	CHARACTER(len=NF90_MAX_NAME), INTENT(out) :: dn
373	!! Name of the dimension.
374	LOGICAL :: ret
375	!! .true. if no error(s) occured, .false. otherwise
376	INTEGER :: vid,err
377	CHARACTER(len=15) :: i2s
378	ret = .false.
379	! --- Get dimension informations
380	IF (NF90_INQUIRE_DIMENSION(fid,did,dn,ds) /= NF90_NOERR) RETURN
381	IF (ds < 2) RETURN
382	IF (NF90_INQ_VARID(fid,TRIM(dn),vid) /= NF90_NOERR) RETURN
383	ALLOCATE(values(ds))
384	IF (NF90_GET_VAR(fid,vid,values) /= NF90_NOERR) RETURN
385	ret = .true.
386	END FUNCTION nc_get_dim_by_id
387
388	FUNCTION get_nc3_info(path,variable,ncid,vid,dimids,verbose) RESULT(ret)
389	!! Get variable informations from NetCDF file
390	!!
391	!! The method attempts to read the given NetCDF file header and retrieves
392	!! variable's informations from it:
393	!!
394	!! - parent file ID
395	!! - variable ID
396	!! - variable's dimensions ID
397	!! - variable's dimensions sizes
398	!!
399	!! The method always opens and closes the file. If the file cannot be opened,
400	!! __ncid__ is set to -1.
401	CHARACTER(len=*), INTENT(in) :: path
402	!! Path of the NetCDF file
403	CHARACTER(len=*), INTENT(in) :: variable
404	!! Name of the NetCDF variable to extract
405	INTEGER, INTENT(out) :: ncid
406	!! NetCDF file ID
407	INTEGER, INTENT(out) :: vid
408	!! NetCDF Variable ID
409	INTEGER, INTENT(out), DIMENSION(:), ALLOCATABLE :: dimids
410	!! Id of the variable's dimensions. __dimids__ is not allocated on error.
411	LOGICAL, INTENT(in), OPTIONAL :: verbose
412	!! True to print out message on error (default to False).
413	LOGICAL :: ret
414	!! .true. if no errors occured, .false. otherwise
415	INTEGER :: fid,ty,nc,err
416	LOGICAL :: zlog
417	zlog = .false. ; IF (PRESENT(verbose)) zlog = verbose
418	ret = .false.
419	ncid = -1
420	! Opens file
421	IF (NF90_OPEN(TRIM(path),NF90_NOWRITE,fid) /= NF90_NOERR) THEN
422	IF (zlog) WRITE(*,'(a)') 'ERROR: Cannot open '//trim(path)
423	RETURN
424	ENDIF
425	ncid = fid
426	! Searches for variable
427	IF (NF90_INQ_VARID(ncid,TRIM(variable),vid) /= NF90_NOERR) THEN
428	IF (zlog) WRITE(*,'(a)') 'Cannot find '//TRIM(variable)//' in '//trim(path)
429	nc = NF90_CLOSE(fid)
430	RETURN
431	ENDIF
432	! Get variable infos
433	! 1st call to get type and number of dimensions)
434	IF (NF90_INQUIRE_VARIABLE(ncid,vid,xtype=ty,ndims=nc) /= NF90_NOERR) THEN
435	IF (zlog) WRITE(*,'(a)') 'Cannot access to '//TRIM(variable)//' informations'
436	nc = NF90_CLOSE(fid)
437	RETURN
438	ELSE
439	! Checks type
440	IF (ty == NF90_CHAR) THEN
441	IF (zlog) WRITE(*,'(a)') 'Inconsistent variable type (should be numeric)'
442	nc = NF90_CLOSE(fid)
443	RETURN
444	ENDIF
445	ALLOCATE(dimids(nc))
446	ENDIF
447	! Gets variable's dimensions informations
448	! first get dimensions id
449	IF (NF90_INQUIRE_VARIABLE(ncid,vid,dimids=dimids) /= NF90_NOERR) THEN
450	IF (zlog) WRITE(*,'(a)') 'Cannot access to '//TRIM(variable)//' informations'
451	nc = NF90_CLOSE(fid)
452	DEALLOCATE(dimids)
453	RETURN
454	ENDIF
455	nc = NF90_CLOSE(fid)
456	ret = .true.
457	END FUNCTION get_nc3_info
458
459	#if HAVE_NC4_FTN
460	FUNCTION get_nc4_info(path,variable,group,ncid,vid,dimids,verbose) RESULT(ret)
461	!! Get variable informations from NetCDF4 file
462	!!
463	!! The method attempts to read the given NetCDF file header and retrieves
464	!! variable's informations from it :
465	!!
466	!! - parent group/file ID
467	!! - variable ID
468	!! - variable's dimensions ID
469	!! - variable's dimensions sizes
470	!!
471	!! The method always opens and closes the file. If the file cannot be opened,
472	!! __ncid__ is set to -1.
473	CHARACTER(len=*), INTENT(in) :: path
474	!! Path of the NetCDF file
475	CHARACTER(len=*), INTENT(in) :: variable
476	!! Name of the NetCDF variable to extract
477	CHARACTER(len=*), INTENT(in) :: group
478	!! Fullname of the variable's parent group (should be set to empty string for root group)
479	INTEGER, INTENT(out) :: ncid
480	!! NetCDF group ID
481	INTEGER, INTENT(out) :: vid
482	!! NetCDF Variable ID
483	INTEGER, INTENT(out), DIMENSION(:), ALLOCATABLE :: dimids
484	!! Id of the variable's dimensions. On error, __dimids__ is not allocated.
485	LOGICAL, INTENT(in), OPTIONAL :: verbose
486	!! True to print out message on error (default to False).
487	LOGICAL :: ret
488	!! .true. if no errors occured, .false. otherwise
489	INTEGER :: fid,ty,nc,err
490	LOGICAL :: zlog
491	zlog = .false. ; IF (PRESENT(verbose)) zlog = verbose
492	ret = .false.
493	ncid = -1
494	! Opens file
495	IF (NF90_OPEN(TRIM(path),NF90_NOWRITE,fid) /= NF90_NOERR) THEN
496	IF (zlog) WRITE(*,'(a)') 'ERROR: Cannot open '//trim(path)
497	RETURN
498	ENDIF
499	! Searches for group based on its fullname
500	IF (LEN_TRIM(group) == 0) THEN
501	nc = NF90_NOERR ; ncid = fid
502	ELSE
503	nc = NF90_INQ_GRP_FULL_NCID(fid,TRIM(group),ncid)
504	ENDIF
505	SELECT CASE(nc)
506	! NF90_ENOGRP is missing from Netcdf-Fortran4.2 : its value=-125
507	CASE(-125)
508	IF (zlog) WRITE(*,'(a)') TRIM(group)//' does not exist in '//TRIM(path)
509	nc = NF90_CLOSE(fid) ; RETURN
510	CASE(NF90_ENOTNC4,NF90_ESTRICTNC3)
511	IF (zlog) WRITE(*,'(a)') TRIM(path)//' is not a NetCDF-4 file with "group" feature'
512	nc = NF90_CLOSE(fid) ; RETURN
513	CASE(NF90_EHDFERR)
514	IF (zlog) WRITE(*,'(a)') "Too bad, an HDF5 error has been reported..."
515	nc = NF90_CLOSE(fid) ; RETURN
516	END SELECT
517	! Searches for variable
518	IF (NF90_INQ_VARID(ncid,TRIM(variable),vid) /= NF90_NOERR) THEN
519	IF (zlog) WRITE(*,'(a)') 'Cannot find '//TRIM(variable)//' in '//trim(path)
520	nc = NF90_CLOSE(fid)
521	RETURN
522	ENDIF
523	! Get variable infos
524	! 1st call to get type and number of dimensions)
525	IF (NF90_INQUIRE_VARIABLE(ncid,vid,xtype=ty,ndims=nc) /= NF90_NOERR) THEN
526	IF (zlog) WRITE(*,'(a)') 'Cannot access to '//TRIM(variable)//' informations'
527	nc = NF90_CLOSE(fid)
528	RETURN
529	ELSE
530	! Checks type
531	IF (ty == NF90_CHAR) THEN
532	IF (zlog) WRITE(*,'(a)') 'Inconsistent variable type (should be numeric)'
533	nc = NF90_CLOSE(fid)
534	RETURN
535	ENDIF
536	ALLOCATE(dimids(nc))
537	ENDIF
538	! Gets variable's dimensions informations
539	! first get dimensions id
540	IF (NF90_INQUIRE_VARIABLE(ncid,vid,dimids=dimids) /= NF90_NOERR) THEN
541	IF (zlog) WRITE(*,'(a)') 'Cannot access to '//TRIM(variable)//' informations'
542	nc = NF90_CLOSE(fid)
543	DEALLOCATE(dimids)
544	RETURN
545	ENDIF
546	nc = NF90_CLOSE(fid)
547	ret = .true.
548	END FUNCTION get_nc4_info
549	#endif
550
551	!-------------------------
552	! NetCDF data file readers
553	!-------------------------
554
555	FUNCTION ncdf_rd_1d(path,variable,set) RESULT(ret)
556	!! Read a 1D data set from a NetCDF file
557	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
558	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
559	TYPE(DSET1D), INTENT(out) :: set !! Output dataset object
560	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
561	INTEGER :: fi,vi,nd,iret
562	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
563	CHARACTER(len=NF90_MAX_NAME) :: dn
564	CHARACTER(len=15) :: i2s
565	ret = .false.
566	! --- Reset the data set
567	CALL clear_dset(set)
568	! --- Check NetCDF file info
569	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
570	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
571	! --- Check dimension size
572	nd = SIZE(di)
573	IF (nd /= 1) THEN
574	IF (debug) THEN
575	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
576	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
577	ENDIF
578	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
579	ENDIF
580	! --- Get coordinates values
581	ALLOCATE(ds(nd))
582	! ------ X coordinate
583	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
584	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
585	ENDIF
586	! --- Read data
587	ALLOCATE(set%data(ds(1)))
588	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
589	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
590	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
591	ENDIF
592	nd = NF90_CLOSE(fi)
593	set%dname = variable
594	ret = .true.
595	RETURN
596	END FUNCTION ncdf_rd_1d
597
598	FUNCTION ncdf_rd_2d(path,variable,set) RESULT(ret)
599	!! Read a 2D data set from a NetCDF file
600	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
601	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
602	TYPE(DSET2D), INTENT(out) :: set !! Output dataset object
603	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
604	INTEGER :: fi,vi,nd,iret
605	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
606	CHARACTER(len=15) :: i2s
607	ret = .false.
608	! --- Reset the data set
609	CALL clear_dset(set)
610	! --- Check NetCDF file info
611	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
612	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
613	! --- Check dimension size
614	nd = SIZE(di)
615	IF (nd /= 2) THEN
616	IF (debug) THEN
617	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
618	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
619	ENDIF
620	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
621	ENDIF
622	! --- Get dimensions informations
623	ALLOCATE(ds(nd))
624	! ------ X coordinate
625	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
626	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
627	ENDIF
628	! ------ Y coordinate
629	IF (.NOT.nc_get_dim_by_id(fi,di(2),set%y,ds(2),set%yname)) THEN
630	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
631	ENDIF
632	! --- Read data
633	ALLOCATE(set%data(ds(1),ds(2)))
634	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
635	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
636	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
637	ENDIF
638	nd = NF90_CLOSE(fi)
639	set%dname = variable
640	ret = .true.
641	RETURN
642	END FUNCTION ncdf_rd_2d
643
644	FUNCTION ncdf_rd_3d(path,variable,set) RESULT(ret)
645	!! Read a 3D data set from a NetCDF file
646	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
647	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
648	TYPE(DSET3D), INTENT(out) :: set !! Output dataset object
649	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
650	INTEGER :: fi,vi,nd,iret
651	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
652	CHARACTER(len=15) :: i2s
653	ret = .false.
654	! --- Reset the data set
655	CALL clear_dset(set)
656	! --- Check NetCDF file info
657	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
658	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
659	! --- Check dimension size
660	nd = SIZE(di)
661	IF (nd /= 3) THEN
662	IF (debug) THEN
663	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
664	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
665	ENDIF
666	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
667	ENDIF
668	! --- Get dimensions informations
669	ALLOCATE(ds(nd))
670	! ------ X coordinate
671	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
672	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
673	ENDIF
674	! ------ Y coordinate
675	IF (.NOT.nc_get_dim_by_id(fi,di(2),set%y,ds(2),set%yname)) THEN
676	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
677	ENDIF
678	! ------ Z coordinate
679	IF (.NOT.nc_get_dim_by_id(fi,di(3),set%z,ds(3),set%zname)) THEN
680	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
681	ENDIF
682	! --- Read data
683	ALLOCATE(set%data(ds(1),ds(2),ds(3)))
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	set%dname = variable
690	ret = .true.
691	RETURN
692	END FUNCTION ncdf_rd_3d
693
694	FUNCTION ncdf_rd_4d(path,variable,set) RESULT(ret)
695	!! Read a 4D data set from a NetCDF file
696	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
697	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
698	TYPE(DSET4D), INTENT(out) :: set !! Output dataset object
699	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
700	INTEGER :: fi,vi,nd,iret
701	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
702	CHARACTER(len=15) :: i2s
703	ret = .false.
704	! --- Reset the data set
705	CALL clear_dset(set)
706	! --- Check NetCDF file info
707	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
708	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
709	! --- Check dimension size
710	nd = SIZE(di)
711	IF (nd /= 4) THEN
712	IF (debug) THEN
713	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
714	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
715	ENDIF
716	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
717	ENDIF
718	! --- Get dimensions informations
719	ALLOCATE(ds(nd))
720	! ------ X coordinate
721	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
722	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
723	ENDIF
724	! ------ Y coordinate
725	IF (.NOT.nc_get_dim_by_id(fi,di(2),set%y,ds(2),set%yname)) THEN
726	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
727	ENDIF
728	! ------ Z coordinate
729	IF (.NOT.nc_get_dim_by_id(fi,di(3),set%z,ds(3),set%zname)) THEN
730	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
731	ENDIF
732	! ------ T coordinate
733	IF (.NOT.nc_get_dim_by_id(fi,di(4),set%t,ds(4),set%tname)) THEN
734	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
735	ENDIF
736	! --- Read data
737	ALLOCATE(set%data(ds(1),ds(2),ds(3),ds(4)))
738	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
739	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
740	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
741	ENDIF
742	nd = NF90_CLOSE(fi)
743	set%dname = variable
744	ret = .true.
745	RETURN
746	END FUNCTION ncdf_rd_4d
747
748	FUNCTION ncdf_rd_5d(path,variable,set) RESULT(ret)
749	!! Read a 5D data set from a NetCDF file
750	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
751	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
752	TYPE(DSET5D), INTENT(out) :: set !! Output dataset object
753	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
754	INTEGER :: fi,vi,nd,iret
755	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
756	CHARACTER(len=15) :: i2s
757	ret = .false.
758	! --- Reset the data set
759	CALL clear_dset(set)
760	! --- Check NetCDF file info
761	IF (.NOT.get_nc_info(path,variable,fi,vi,di)) RETURN
762	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
763	! --- Check dimension size
764	nd = SIZE(di)
765	IF (nd /= 5) THEN
766	IF (debug) THEN
767	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
768	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
769	ENDIF
770	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
771	ENDIF
772	! --- Get dimensions informations
773	ALLOCATE(ds(nd))
774	! ------ X coordinate
775	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
776	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
777	ENDIF
778	! ------ Y coordinate
779	IF (.NOT.nc_get_dim_by_id(fi,di(2),set%y,ds(2),set%yname)) THEN
780	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
781	ENDIF
782	! ------ Z coordinate
783	IF (.NOT.nc_get_dim_by_id(fi,di(3),set%z,ds(3),set%zname)) THEN
784	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
785	ENDIF
786	! ------ T coordinate
787	IF (.NOT.nc_get_dim_by_id(fi,di(4),set%t,ds(4),set%tname)) THEN
788	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
789	ENDIF
790	! ------ W coordinate
791	IF (.NOT.nc_get_dim_by_id(fi,di(5),set%w,ds(5),set%wname)) THEN
792	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
793	ENDIF
794	! --- Read data
795	ALLOCATE(set%data(ds(1),ds(2),ds(3),ds(4),ds(5)))
796	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
797	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
798	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
799	ENDIF
800	nd = NF90_CLOSE(fi)
801	set%dname = variable
802	ret = .true.
803	RETURN
804	END FUNCTION ncdf_rd_5d
805
806	! NetCDF4 with groups !
807	#if HAVE_NC4_FTN
808
809	FUNCTION ncdf4_rd_1d(path,group,variable,set) RESULT(ret)
810	!! Read a 1D data set from a NetCDF4 file
811	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
812	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
813	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
814	TYPE(DSET1D), INTENT(out) :: set !! Output dataset
815	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
816	INTEGER :: fi,vi,nd,iret
817	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
818	CHARACTER(len=NF90_MAX_NAME) :: dn
819	CHARACTER(len=15) :: i2s
820	ret = .false.
821	! --- Reset the data set
822	CALL clear_dset(set)
823	! --- Check NetCDF file info
824	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) RETURN
825	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
826	! --- Check dimension size
827	nd = SIZE(di)
828	IF (nd /= 1) THEN
829	IF (debug) THEN
830	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
831	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
832	ENDIF
833	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
834	ENDIF
835	! --- Get coordinates values
836	ALLOCATE(ds(nd))
837	! ------ X coordinate
838	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
839	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
840	ENDIF
841	! --- Read data
842	ALLOCATE(set%data(ds(1)))
843	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
844	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
845	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
846	ENDIF
847	nd = NF90_CLOSE(fi)
848	set%dname = variable
849	ret = .true.
850	RETURN
851	END FUNCTION ncdf4_rd_1d
852
853	FUNCTION ncdf4_rd_2d(path,group,variable,set) RESULT(ret)
854	!! Read a 2D data set from a NetCDF4 file
855	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
856	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
857	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
858	TYPE(DSET2D), INTENT(out) :: set !! Output dataset
859	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
860	INTEGER :: fi,vi,nd,iret
861	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
862	CHARACTER(len=15) :: i2s
863	ret = .false.
864	! --- Reset the data set
865	CALL clear_dset(set)
866	! --- Check NetCDF file info
867	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) rETURN
868	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
869	! --- Check dimension size
870	nd = SIZE(di)
871	IF (nd /= 2) THEN
872	IF (debug) THEN
873	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
874	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
875	ENDIF
876	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
877	ENDIF
878	! --- Get dimensions informations
879	ALLOCATE(ds(nd))
880	! ------ X coordinate
881	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
882	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
883	ENDIF
884	! ------ Y coordinate
885	IF (.NOT.nc_get_dim_by_id(fi,di(2),set%y,ds(2),set%yname)) THEN
886	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
887	ENDIF
888	! --- Read data
889	ALLOCATE(set%data(ds(1),ds(2)))
890	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
891	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
892	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
893	ENDIF
894	nd = NF90_CLOSE(fi)
895	set%dname = variable
896	ret = .true.
897	RETURN
898	END FUNCTION ncdf4_rd_2d
899
900	FUNCTION ncdf4_rd_3d(path,group,variable,set) RESULT(ret)
901	!! Read a 3D data set from a NetCDF4 file
902	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
903	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
904	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
905	TYPE(DSET3D), INTENT(out) :: set !! Output dataset
906	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
907	INTEGER :: fi,vi,nd,iret
908	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
909	CHARACTER(len=15) :: i2s
910	ret = .false.
911	! --- Reset the data set
912	CALL clear_dset(set)
913	! --- Check NetCDF file info
914	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) RETURN
915	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
916	! --- Check dimension size
917	nd = SIZE(di)
918	IF (nd /= 3) THEN
919	IF (debug) THEN
920	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
921	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
922	ENDIF
923	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
924	ENDIF
925	! --- Get dimensions informations
926	ALLOCATE(ds(nd))
927	! ------ X coordinate
928	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
929	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
930	ENDIF
931	! ------ Y coordinate
932	IF (.NOT.nc_get_dim_by_id(fi,di(2),set%y,ds(2),set%yname)) THEN
933	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
934	ENDIF
935	! ------ Z coordinate
936	IF (.NOT.nc_get_dim_by_id(fi,di(3),set%z,ds(3),set%zname)) THEN
937	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
938	ENDIF
939	! --- Read data
940	ALLOCATE(set%data(ds(1),ds(2),ds(3)))
941	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
942	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
943	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
944	ENDIF
945	nd = NF90_CLOSE(fi)
946	set%dname = variable
947	ret = .true.
948	RETURN
949	END FUNCTION ncdf4_rd_3d
950
951	FUNCTION ncdf4_rd_4d(path,group,variable,set) RESULT(ret)
952	!! Read a 4D data set from a NetCDF4 file
953	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
954	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
955	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
956	TYPE(DSET4D), INTENT(out) :: set !! Output dataset
957	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
958	INTEGER :: fi,vi,nd,iret
959	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
960	CHARACTER(len=15) :: i2s
961	ret = .false.
962	! --- Reset the data set
963	CALL clear_dset(set)
964	! --- Check NetCDF file info
965	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) RETURN
966	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
967	! --- Check dimension size
968	nd = SIZE(di)
969	IF (nd /= 4) THEN
970	IF (debug) THEN
971	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
972	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
973	ENDIF
974	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
975	ENDIF
976	! --- Get dimensions informations
977	ALLOCATE(ds(nd))
978	! ------ X coordinate
979	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
980	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
981	ENDIF
982	! ------ Y coordinate
983	IF (.NOT.nc_get_dim_by_id(fi,di(2),set%y,ds(2),set%yname)) THEN
984	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
985	ENDIF
986	! ------ Z coordinate
987	IF (.NOT.nc_get_dim_by_id(fi,di(3),set%z,ds(3),set%zname)) THEN
988	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
989	ENDIF
990	! ------ T coordinate
991	IF (.NOT.nc_get_dim_by_id(fi,di(4),set%t,ds(4),set%tname)) THEN
992	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
993	ENDIF
994	! --- Read data
995	ALLOCATE(set%data(ds(1),ds(2),ds(3),ds(4)))
996	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
997	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
998	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
999	ENDIF
1000	nd = NF90_CLOSE(fi)
1001	set%dname = variable
1002	ret = .true.
1003	RETURN
1004	END FUNCTION ncdf4_rd_4d
1005
1006	FUNCTION ncdf4_rd_5d(path,group,variable,set) RESULT(ret)
1007	!! Read a 5D data set from a NetCDF4 file
1008	CHARACTER(len=*), INTENT(in) :: path !! Path of the NetCDF file
1009	CHARACTER(len=*), INTENT(in) :: group !! Parent group of the variable (empty for root)
1010	CHARACTER(len=*), INTENT(in) :: variable !! Name of the NetCDF variable to extract
1011	TYPE(DSET5D), INTENT(out) :: set !! Output dataset
1012	LOGICAL :: ret !! .true. if no errors occured, .false. otherwise
1013	INTEGER :: fi,vi,nd,iret
1014	INTEGER, DIMENSION(:), ALLOCATABLE :: di,ds
1015	CHARACTER(len=15) :: i2s
1016	ret = .false.
1017	! --- Reset the data set
1018	CALL clear_dset(set)
1019	! --- Check NetCDF file info
1020	IF (.NOT.get_nc_info(path,group,variable,fi,vi,di)) RETURN
1021	iret = NF90_OPEN(path,NF90_NOWRITE,fi)
1022	! --- Check dimension size
1023	nd = SIZE(di)
1024	IF (nd /= 5) THEN
1025	IF (debug) THEN
1026	WRITE(i2s,*) nd ; i2s=ADJUSTL(i2s)
1027	WRITE(*,'(a)') "ERROR: Incompatible size, DSET should be"//TRIM(i2s)//"D"
1028	ENDIF
1029	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
1030	ENDIF
1031	! --- Get dimensions informations
1032	ALLOCATE(ds(nd))
1033	! ------ X coordinate
1034	IF (.NOT.nc_get_dim_by_id(fi,di(1),set%x,ds(1),set%xname)) THEN
1035	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
1036	ENDIF
1037	! ------ Y coordinate
1038	IF (.NOT.nc_get_dim_by_id(fi,di(2),set%y,ds(2),set%yname)) THEN
1039	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
1040	ENDIF
1041	! ------ Z coordinate
1042	IF (.NOT.nc_get_dim_by_id(fi,di(3),set%z,ds(3),set%zname)) THEN
1043	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
1044	ENDIF
1045	! ------ T coordinate
1046	IF (.NOT.nc_get_dim_by_id(fi,di(4),set%t,ds(4),set%tname)) THEN
1047	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
1048	ENDIF
1049	! ------ W coordinate
1050	IF (.NOT.nc_get_dim_by_id(fi,di(5),set%w,ds(5),set%wname)) THEN
1051	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
1052	CALL clear_dset(set) ; RETURN
1053	ENDIF
1054	! --- Read data
1055	ALLOCATE(set%data(ds(1),ds(2),ds(3),ds(4),ds(5)))
1056	IF (NF90_GET_VAR(fi,vi,set%data) /= NF90_NOERR) THEN
1057	IF (debug) WRITE(*,'(a)') "ERROR:"//TRIM(variable)//": Cannot get values"
1058	nd = NF90_CLOSE(fi) ; CALL clear_dset(set) ; RETURN
1059	ENDIF
1060	nd = NF90_CLOSE(fi)
1061	set%dname = variable
1062	ret = .true.
1063	RETURN
1064	END FUNCTION ncdf4_rd_5d
1065	#endif
1066
1067	!------------------------
1068	! ASCII data file readers
1069	!------------------------
1070
1071	FUNCTION ascii_rd_1d(path,set) RESULT(ret)
1072	!! Read a 1D data set from a ASCII file
1073	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1074	TYPE(dset1d), INTENT(out) :: set !! output 1D dataset
1075	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1076	INTEGER :: i,e
1077	REAL(kind=8), DIMENSION(2) :: vl
1078	INTEGER, DIMENSION(1) :: cc,ds,dp
1079	ret = .false.
1080	CALL clear_dset(set)
1081	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1082	ALLOCATE(set%x(ds(1)), &
1083	set%data(ds(1)))
1084	cc(:) = 1
1085	DO i=1,ds(1)
1086	READ(666,*) vl(:)
1087	set%x(cc(1)) = vl(1)
1088	set%data(cc(1)) = vl(2)
1089	cc(1) = cc(1) + 1
1090	ENDDO
1091	READ(666,*,iostat=e) vl(1)
1092	IF (e == 0) THEN
1093	IF (debug) WRITE(,) 'ERROR: Extra value found...'
1094	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1095	ENDIF
1096	ret = .true.
1097	CLOSE(666)
1098	RETURN
1099	END FUNCTION ascii_rd_1d
1100
1101	FUNCTION ascii_rd_2d(path,set) RESULT(ret)
1102	!! Read a 2D data set from a ASCII file
1103	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1104	TYPE(dset2d), INTENT(out) :: set !! output 2D dataset
1105	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1106	INTEGER :: i,e
1107	REAL(kind=8), DIMENSION(3) :: vl
1108	INTEGER, DIMENSION(2) :: cc,ds,dp
1109	ret = .false.
1110	CALL clear_dset(set)
1111	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1112	ALLOCATE(set%x(ds(1)),set%y(ds(2)), &
1113	set%data(ds(1),ds(2)))
1114	cc(:) = 1
1115	DO i=1,PRODUCT(ds)
1116	READ(666,*,iostat=e) vl ! Reads the line
1117	IF (e /= 0) THEN
1118	IF (debug) WRITE(*,'(a)') 'ERROR: Malformed file, line: ',i+2
1119	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1120	ENDIF
1121	! --- X coordinate
1122	set%x(cc(1)) = vl(1)
1123	! --- Y coordinate
1124	set%y(cc(2)) = vl(2)
1125	! --- Data
1126	set%data(cc(1),cc(2)) = vl(3)
1127	! - Update counters
1128	IF (mod(i,dp(1))==0) cc(1) = cc(1)+1 ; IF (cc(1) > ds(1)) cc(1)=1
1129	IF (mod(i,dp(2))==0) cc(2) = cc(2)+1 ; IF (cc(2) > ds(2)) cc(2)=1
1130	ENDDO
1131	READ(666,*,iostat=e) vl(1)
1132	IF (e == 0) THEN
1133	IF (debug) WRITE(*,'(a)') 'ERROR: Extra value found'
1134	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1135	ENDIF
1136	ret = .true.
1137	CLOSE(666)
1138	RETURN
1139	END FUNCTION ascii_rd_2d
1140
1141	FUNCTION ascii_rd_3d(path,set) RESULT(ret)
1142	!! Read a 3D data set from a ASCII file
1143	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1144	TYPE(dset3d), INTENT(out) :: set !! output 3D dataset
1145	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1146	INTEGER :: i,e
1147	REAL(kind=8), DIMENSION(4) :: vl
1148	INTEGER, DIMENSION(3) :: cc,ds,dp
1149	ret = .false.
1150	CALL clear_dset(set)
1151	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1152	ALLOCATE(set%x(ds(1)),set%y(ds(2)),set%z(ds(3)), &
1153	set%data(ds(1),ds(2),ds(3)))
1154	cc(:) = 1
1155	DO i=1,PRODUCT(ds)
1156	READ(666,*,iostat=e) vl ! Reads the line
1157	IF (e /= 0) THEN
1158	IF (debug) WRITE(*,'(a)') 'ERROR: Malformed file, line: ',i+2
1159	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1160	ENDIF
1161	! --- X coordinate
1162	set%x(cc(1)) = vl(1)
1163	! --- Y coordinate
1164	set%y(cc(2)) = vl(2)
1165	! --- Z coordinate
1166	set%z(cc(3)) = vl(3)
1167	! --- Data
1168	set%data(cc(1),cc(2),cc(3)) = vl(4)
1169	! - Update counters
1170	IF (mod(i,dp(1))==0) cc(1) = cc(1)+1 ; IF (cc(1) > ds(1)) cc(1)=1
1171	IF (mod(i,dp(2))==0) cc(2) = cc(2)+1 ; IF (cc(2) > ds(2)) cc(2)=1
1172	IF (mod(i,dp(3))==0) cc(3) = cc(3)+1 ; IF (cc(3) > ds(3)) cc(3)=1
1173	ENDDO
1174	READ(666,*,iostat=e) vl(1)
1175	IF (e == 0) THEN
1176	IF (debug) WRITE(*,'(a)') 'ERROR: Extra value found'
1177	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1178	ENDIF
1179	ret = .true.
1180	CLOSE(666)
1181	RETURN
1182	END FUNCTION ascii_rd_3d
1183
1184	FUNCTION ascii_rd_4d(path,set) RESULT(ret)
1185	!! Read a 4D data set from a ASCII file
1186	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1187	TYPE(dset4d), INTENT(out) :: set !! output 4D dataset
1188	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1189	INTEGER :: i,e
1190	REAL(kind=8), DIMENSION(5) :: vl
1191	INTEGER, DIMENSION(4) :: cc,ds,dp
1192	ret = .false.
1193	CALL clear_dset(set)
1194	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1195	ALLOCATE(set%x(ds(1)),set%y(ds(2)),set%z(ds(3)),set%t(ds(4)), &
1196	set%data(ds(1),ds(2),ds(3),ds(4)))
1197	cc(:) = 1
1198	DO i=1,PRODUCT(ds)
1199	READ(666,*,iostat=e) vl ! Reads the line
1200	IF (e /= 0) THEN
1201	IF (debug) WRITE(*,'(a)') 'ERROR: Malformed file, line: ',i+2
1202	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1203	ENDIF
1204	! --- X coordinate
1205	set%x(cc(1)) = vl(1)
1206	! --- Y coordinate
1207	set%y(cc(2)) = vl(2)
1208	! --- Z coordinate
1209	set%z(cc(3)) = vl(3)
1210	! --- T coordinate
1211	set%t(cc(4)) = vl(4)
1212	! --- Data
1213	set%data(cc(1),cc(2),cc(3),cc(4)) = vl(5)
1214	! - Update counters
1215	IF (mod(i,dp(1))==0) cc(1) = cc(1)+1 ; IF (cc(1) > ds(1)) cc(1)=1
1216	IF (mod(i,dp(2))==0) cc(2) = cc(2)+1 ; IF (cc(2) > ds(2)) cc(2)=1
1217	IF (mod(i,dp(3))==0) cc(3) = cc(3)+1 ; IF (cc(3) > ds(3)) cc(3)=1
1218	IF (mod(i,dp(4))==0) cc(4) = cc(4)+1 ; IF (cc(4) > ds(4)) cc(4)=1
1219	ENDDO
1220	READ(666,*,iostat=e) vl(1)
1221	IF (e == 0) THEN
1222	IF (debug) WRITE(*,'(a)') 'ERROR: Extra value found'
1223	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1224	ENDIF
1225	CLOSE(666)
1226	ret = .true.
1227	RETURN
1228	END FUNCTION ascii_rd_4d
1229
1230	FUNCTION ascii_rd_5d(path,set) RESULT(ret)
1231	!! Read a 5D data set from a ASCII file
1232	CHARACTER(len=*), INTENT(in) :: path !! Path of the ASCII file to read
1233	TYPE(dset5d), INTENT(out) :: set !! output 5D dataset
1234	LOGICAL :: ret !! .true. if no error occured, .false. otherwise
1235	INTEGER :: i,e
1236	REAL(kind=8), DIMENSION(6) :: vl
1237	INTEGER, DIMENSION(5) :: cc,ds,dp
1238	ret = .false.
1239	CALL clear_dset(set)
1240	IF (.NOT.ascii_header(path,ds,dp)) RETURN
1241	ALLOCATE(set%x(ds(1)),set%y(ds(2)),set%z(ds(3)),set%t(ds(4)),set%w(ds(5)), &
1242	set%data(ds(1),ds(2),ds(3),ds(4),ds(5)))
1243	cc(:) = 1
1244	DO i=1,PRODUCT(ds)
1245	READ(666,*,iostat=e) vl ! Reads the line
1246	IF (e /= 0) THEN
1247	IF (debug) WRITE(*,'(a)') 'ERROR: Malformed file, line: ',i+2
1248	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1249	ENDIF
1250	! --- X coordinate
1251	set%x(cc(1)) = vl(1)
1252	! --- Y coordinate
1253	set%y(cc(2)) = vl(2)
1254	! --- Z coordinate
1255	set%z(cc(3)) = vl(3)
1256	! --- T coordinate
1257	set%t(cc(4)) = vl(4)
1258	! --- W coordinate
1259	set%w(cc(5)) = vl(5)
1260	! --- Data
1261	set%data(cc(1),cc(2),cc(3),cc(4),cc(5)) = vl(6)
1262	! - Update counters
1263	IF (mod(i,dp(1)) == 0) cc(1)=cc(1)+1 ; IF (cc(1) > ds(1)) cc(1)=1
1264	IF (mod(i,dp(2)) == 0) cc(2)=cc(2)+1 ; IF (cc(2) > ds(2)) cc(2)=1
1265	IF (mod(i,dp(3)) == 0) cc(3)=cc(3)+1 ; IF (cc(3) > ds(3)) cc(3)=1
1266	IF (mod(i,dp(4)) == 0) cc(4)=cc(4)+1 ; IF (cc(4) > ds(4)) cc(4)=1
1267	IF (mod(i,dp(5)) == 0) cc(5)=cc(5)+1 ; IF (cc(5) > ds(5)) cc(5)=1
1268	ENDDO
1269	READ(666,*,iostat=e) vl(1)
1270	IF (e == 0) THEN
1271	IF (debug) WRITE(*,'(a)') 'ERROR: Extra value found'
1272	ret = .false. ; call clear_dset(set) ; CLOSE(666) ; RETURN
1273	ENDIF
1274	CLOSE(666)
1275	ret = .true.
1276	RETURN
1277	END FUNCTION ascii_rd_5d
1278
1279	SUBROUTINE clr_1d_set(set)
1280	!! Clear the given 1D data set
1281	TYPE(DSET1D), INTENT(inout) :: set !! dataset object to clear
1282	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1283	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1284	RETURN
1285	END SUBROUTINE clr_1d_set
1286
1287	SUBROUTINE clr_2d_set(set)
1288	!! Clear the given 2D data set
1289	TYPE(DSET2D), INTENT(inout) :: set !! dataset object to clear
1290	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1291	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1292	IF (ALLOCATED(set%y)) DEALLOCATE(set%y)
1293	RETURN
1294	END SUBROUTINE clr_2d_set
1295
1296	SUBROUTINE clr_3d_set(set)
1297	!! Clear the given 3D data set
1298	TYPE(DSET3D), INTENT(inout) :: set !! dataset object to clear
1299	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1300	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1301	IF (ALLOCATED(set%y)) DEALLOCATE(set%y)
1302	IF (ALLOCATED(set%z)) DEALLOCATE(set%z)
1303	RETURN
1304	END SUBROUTINE clr_3d_set
1305
1306	SUBROUTINE clr_4d_set(set)
1307	!! Clear the given 4D data set
1308	TYPE(DSET4D), INTENT(inout) :: set !! dataset object to clear
1309	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1310	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1311	IF (ALLOCATED(set%y)) DEALLOCATE(set%y)
1312	IF (ALLOCATED(set%z)) DEALLOCATE(set%z)
1313	IF (ALLOCATED(set%t)) DEALLOCATE(set%t)
1314	RETURN
1315	END SUBROUTINE clr_4d_set
1316
1317	SUBROUTINE clr_5d_set(set)
1318	!! Clear the given 5D data set
1319	TYPE(DSET5D), INTENT(inout) :: set !! dataset object to clear
1320	IF (ALLOCATED(set%data)) DEALLOCATE(set%data)
1321	IF (ALLOCATED(set%x)) DEALLOCATE(set%x)
1322	IF (ALLOCATED(set%y)) DEALLOCATE(set%y)
1323	IF (ALLOCATED(set%z)) DEALLOCATE(set%z)
1324	IF (ALLOCATED(set%t)) DEALLOCATE(set%t)
1325	IF (ALLOCATED(set%w)) DEALLOCATE(set%w)
1326	RETURN
1327	END SUBROUTINE clr_5d_set
1328
1329	FUNCTION has_d_1d(dset) RESULT(yes)
1330	!! Check whether or not the dataset has data.
1331	TYPE(DSET1D), INTENT(in) :: dset !! Dataset to check
1332	LOGICAL :: yes !! return status
1333	yes = (ALLOCATED(dset%data).AND. &
1334	ALLOCATED(dset%x))
1335	END FUNCTION has_d_1d
1336
1337	FUNCTION has_d_2d(dset) RESULT(yes)
1338	!! Check whether or not the dataset has data.
1339	TYPE(DSET2D), INTENT(in) :: dset !! Dataset to check
1340	LOGICAL :: yes !! return status
1341	yes = (ALLOCATED(dset%data).AND. &
1342	ALLOCATED(dset%x) .AND. &
1343	ALLOCATED(dset%y))
1344	END FUNCTION has_d_2d
1345
1346	FUNCTION has_d_3d(dset) RESULT(yes)
1347	!! Check whether or not the dataset has data.
1348	TYPE(DSET3D), INTENT(in) :: dset !! Dataset to check
1349	LOGICAL :: yes !! return status
1350	yes = (ALLOCATED(dset%data).AND. &
1351	ALLOCATED(dset%x) .AND. &
1352	ALLOCATED(dset%y) .AND. &
1353	ALLOCATED(dset%z))
1354	END FUNCTION has_d_3d
1355
1356	FUNCTION has_d_4d(dset) RESULT(yes)
1357	!! Check whether or not the dataset has data.
1358	TYPE(DSET4D), INTENT(in) :: dset !! Dataset to check
1359	LOGICAL :: yes !! return status
1360	yes = (ALLOCATED(dset%data).AND. &
1361	ALLOCATED(dset%x) .AND. &
1362	ALLOCATED(dset%y) .AND. &
1363	ALLOCATED(dset%z) .AND. &
1364	ALLOCATED(dset%t))
1365	END FUNCTION has_d_4d
1366
1367	FUNCTION has_d_5d(dset) RESULt(yes)
1368	!! Check whether or not the dataset has data.
1369	TYPE(DSET5D), INTENT(in) :: dset !! Dataset to check
1370	LOGICAL :: yes !! return status
1371	yes = (ALLOCATED(dset%data).AND. &
1372	ALLOCATED(dset%x) .AND. &
1373	ALLOCATED(dset%y) .AND. &
1374	ALLOCATED(dset%z) .AND. &
1375	ALLOCATED(dset%t) .AND. &
1376	ALLOCATED(dset%w))
1377	END FUNCTION has_d_5d
1378
1379	FUNCTION nc_wr_1d(path,dset) RESULT(ret)
1380	CHARACTER(len=*), INTENT(in) :: path !! Path of the netcdf file.
1381	TYPE(DSET1D), INTENT(in) :: dset !! Dataset to write
1382	LOGICAL :: ret !! Return status
1383	INTEGER :: fi,vi,nd,ds,iret
1384	INTEGER, DIMENSION(:), ALLOCATABLE :: di
1385	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: tv
1386	CHARACTER(len=15) :: i2s
1387	LOGICAL :: hxd,ok
1388	INTEGER :: xid,vid
1389	ret = has_data(dset)
1390	IF (.NOT.ret) RETURN
1391	ret = .false.
1392	INQUIRE(FILE=TRIM(path),EXIST=ok)
1393	IF (ok) THEN
1394	IF (NF90_OPEN(TRIM(path), NF90_WRITE, fi) /= NF90_NOERR) RETURN
1395	ELSE
1396	IF (NF90_CREATE(TRIM(path), NF90_NOCLOBBER, fi) /= NF90_NOERR) RETURN
1397	ENDIF
1398	iret = NF90_CLOSE(fi)
1399	! if variable already exist: get out of here !
1400	IF (get_nc_info(path,dset%dname,fi,vi,di)) RETURN
1401
1402	iret = NF90_OPEN(path,NF90_WRITE,fi)
1403
1404	hxd = nc_get_dim_by_name(fi,dset%xname,tv,ds,xid)
1405	IF (hxd) THEN
1406	IF(.NOT.compare(dset%x,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1407	ENDIF
1408
1409	IF (.NOT.hxd) THEN
1410	ret = nc_set_dim(fi,dset%xname,dset%x,xid)
1411	IF (.NOT.ret) THEN ; iret= NF90_CLOSE(fi) ; RETURN ; ENDIF
1412	ENDIF
1413	ret = .false.
1414	iret = nf90_redef(fi)
1415	iret = nf90_def_var(fi, TRIM(dset%dname), nf90_double,(/xid/),vid)
1416	IF (iret /= 0) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1417	iret = nf90_enddef(fi)
1418	iret = nf90_put_var(fi,vid,dset%data)
1419	ret = (iret == 0)
1420	iret=NF90_CLOSE(fi)
1421	ret = .true.
1422	RETURN
1423	END FUNCTION nc_wr_1d
1424
1425	FUNCTION nc_wr_2d(path,dset) RESULT(ret)
1426	CHARACTER(len=*), INTENT(in) :: path !! Path of the netcdf file.
1427	TYPE(DSET2D), INTENT(in) :: dset !! Dataset to write
1428	LOGICAL :: ret !! Return status
1429	INTEGER :: fi,vi,nd,ds,iret,nferr
1430	INTEGER, DIMENSION(:), ALLOCATABLE :: di
1431	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: tv
1432	CHARACTER(len=15) :: i2s
1433	LOGICAL :: hxd,hyd,ok
1434	INTEGER :: xid,yid,vid
1435	ret = has_data(dset)
1436	IF (.NOT.ret) RETURN
1437	ret = .false.
1438	INQUIRE(FILE=TRIM(path),EXIST=ok)
1439	IF (ok) THEN
1440	IF (NF90_OPEN(TRIM(path), NF90_WRITE, fi) /= NF90_NOERR) RETURN
1441	ELSE
1442	IF (NF90_CREATE(TRIM(path), NF90_NOCLOBBER, fi) /= NF90_NOERR) RETURN
1443	ENDIF
1444	iret = NF90_CLOSE(fi)
1445
1446	! if variable already exist: get out of here !
1447	IF (get_nc_info(path,dset%dname,fi,vi,di)) RETURN
1448
1449	iret = NF90_OPEN(path,NF90_WRITE,fi)
1450
1451	hxd = nc_get_dim_by_name(fi,dset%xname,tv,ds,xid)
1452	IF (hxd) THEN
1453	IF(.NOT.compare(dset%x,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1454	ENDIF
1455	hyd = nc_get_dim_by_name(fi,dset%yname,tv,ds,yid)
1456	IF (hyd) THEN
1457	IF(.NOT.compare(dset%y,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1458	ENDIF
1459
1460	IF (.NOT.hxd) THEN
1461	ret = nc_set_dim(fi,dset%xname,dset%x,xid)
1462	IF (.NOT.ret) THEN ; iret= NF90_CLOSE(fi) ; RETURN ; ENDIF
1463	ENDIF
1464	IF (.NOT.hyd) THEN
1465	ret = nc_set_dim(fi,dset%yname,dset%y,yid)
1466	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1467	ENDIF
1468	ret = .false.
1469	iret = nf90_redef(fi)
1470	iret = nf90_def_var(fi, TRIM(dset%dname), nf90_double,(/xid,yid/),vid)
1471	IF (iret /= 0) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1472	iret = nf90_enddef(fi)
1473	iret = nf90_put_var(fi,vid,dset%data)
1474	ret = (iret == 0)
1475	iret=NF90_CLOSE(fi)
1476	ret = .true.
1477	RETURN
1478	END FUNCTION nc_wr_2d
1479
1480	FUNCTION nc_wr_3d(path,dset) RESULT(ret)
1481	CHARACTER(len=*), INTENT(in) :: path !! Path of the netcdf file.
1482	TYPE(DSET3D), INTENT(in) :: dset !! Dataset to write
1483	LOGICAL :: ret !! Return status
1484	INTEGER :: fi,vi,nd,ds,iret
1485	INTEGER, DIMENSION(:), ALLOCATABLE :: di
1486	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: tv
1487	CHARACTER(len=15) :: i2s
1488	LOGICAL :: hxd,hyd,hzd,ok
1489	INTEGER :: xid,yid,zid,vid
1490	ret = has_data(dset)
1491	IF (.NOT.ret) RETURN
1492	ret = .false.
1493	INQUIRE(FILE=TRIM(path),EXIST=ok)
1494	IF (ok) THEN
1495	IF (NF90_OPEN(TRIM(path), NF90_WRITE, fi) /= NF90_NOERR) RETURN
1496	ELSE
1497	IF (NF90_CREATE(TRIM(path), NF90_NOCLOBBER, fi) /= NF90_NOERR) RETURN
1498	ENDIF
1499	iret = NF90_CLOSE(fi)
1500
1501	! if variable already exist: get out of here !
1502	IF (get_nc_info(path,dset%dname,fi,vi,di)) RETURN
1503
1504	iret = NF90_OPEN(path,NF90_WRITE,fi)
1505
1506	hxd = nc_get_dim_by_name(fi,dset%xname,tv,ds,xid)
1507	IF (hxd) THEN
1508	IF(.NOT.compare(dset%x,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1509	ENDIF
1510	hyd = nc_get_dim_by_name(fi,dset%yname,tv,ds,yid)
1511	IF (hyd) THEN
1512	IF(.NOT.compare(dset%y,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1513	ENDIF
1514	hzd = nc_get_dim_by_name(fi,dset%zname,tv,ds,zid)
1515	IF (hzd) THEN
1516	IF(.NOT.compare(dset%z,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1517	ENDIF
1518
1519	IF (.NOT.hxd) THEN
1520	ret = nc_set_dim(fi,dset%xname,dset%x,xid)
1521	IF (.NOT.ret) THEN ; iret= NF90_CLOSE(fi) ; RETURN ; ENDIF
1522	ENDIF
1523	IF (.NOT.hyd) THEN
1524	ret = nc_set_dim(fi,dset%yname,dset%y,yid)
1525	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1526	ENDIF
1527	IF (.NOT.hzd) THEN
1528	ret = nc_set_dim(fi,dset%zname,dset%z,zid)
1529	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1530	ENDIF
1531	ret = .false.
1532	iret = nf90_redef(fi)
1533	iret = nf90_def_var(fi, TRIM(dset%dname), nf90_double,(/xid,yid,zid/),vid)
1534	IF (iret /= 0) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1535	iret = nf90_enddef(fi)
1536	iret = nf90_put_var(fi,vid,dset%data)
1537	ret = (iret == 0)
1538	iret=NF90_CLOSE(fi)
1539	ret = .true.
1540	RETURN
1541	END FUNCTION nc_wr_3d
1542
1543	FUNCTION nc_wr_4d(path,dset) RESULT(ret)
1544	CHARACTER(len=*), INTENT(in) :: path !! Path of the netcdf file.
1545	TYPE(DSET4D), INTENT(in) :: dset !! Dataset to write
1546	LOGICAL :: ret !! Return status
1547	INTEGER :: fi,vi,nd,ds,iret
1548	INTEGER, DIMENSION(:), ALLOCATABLE :: di
1549	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: tv
1550	CHARACTER(len=15) :: i2s
1551	LOGICAL :: hxd,hyd,hzd,htd,ok
1552	INTEGER :: xid,yid,zid,tid,vid
1553	ret = has_data(dset)
1554	IF (.NOT.ret) RETURN
1555	ret = .false.
1556	INQUIRE(FILE=TRIM(path),EXIST=ok)
1557	IF (ok) THEN
1558	IF (NF90_OPEN(TRIM(path), NF90_WRITE, fi) /= NF90_NOERR) RETURN
1559	ELSE
1560	IF (NF90_CREATE(TRIM(path), NF90_NOCLOBBER, fi) /= NF90_NOERR) RETURN
1561	ENDIF
1562	iret = NF90_CLOSE(fi)
1563
1564	! if variable already exist: get out of here !
1565	IF (get_nc_info(path,dset%dname,fi,vi,di)) RETURN
1566
1567	iret = NF90_OPEN(path,NF90_WRITE,fi)
1568
1569	hxd = nc_get_dim_by_name(fi,dset%xname,tv,ds,xid)
1570	IF (hxd) THEN
1571	IF(.NOT.compare(dset%x,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1572	ENDIF
1573	hyd = nc_get_dim_by_name(fi,dset%yname,tv,ds,yid)
1574	IF (hyd) THEN
1575	IF(.NOT.compare(dset%y,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1576	ENDIF
1577	hzd = nc_get_dim_by_name(fi,dset%zname,tv,ds,zid)
1578	IF (hzd) THEN
1579	IF(.NOT.compare(dset%z,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1580	ENDIF
1581	htd = nc_get_dim_by_name(fi,dset%tname,tv,ds,tid)
1582	IF (htd) THEN
1583	IF(.NOT.compare(dset%t,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1584	ENDIF
1585
1586	IF (.NOT.hxd) THEN
1587	ret = nc_set_dim(fi,dset%xname,dset%x,xid)
1588	IF (.NOT.ret) THEN ; iret= NF90_CLOSE(fi) ; RETURN ; ENDIF
1589	ENDIF
1590	IF (.NOT.hyd) THEN
1591	ret = nc_set_dim(fi,dset%yname,dset%y,yid)
1592	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1593	ENDIF
1594	IF (.NOT.hzd) THEN
1595	ret = nc_set_dim(fi,dset%zname,dset%z,zid)
1596	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1597	ENDIF
1598	IF (.NOT.htd) THEN
1599	ret = nc_set_dim(fi,dset%tname,dset%t,tid)
1600	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1601	ENDIF
1602	ret = .false.
1603	iret = nf90_redef(fi)
1604	iret = nf90_def_var(fi, TRIM(dset%dname), nf90_double,(/xid,yid,zid,tid/),vid)
1605	IF (iret /= 0) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1606	iret = nf90_enddef(fi)
1607	iret = nf90_put_var(fi,vid,dset%data)
1608	ret = (iret == 0)
1609	iret=NF90_CLOSE(fi)
1610	ret = .true.
1611	RETURN
1612	END FUNCTION nc_wr_4d
1613
1614	FUNCTION nc_wr_5d(path,dset) RESULT(ret)
1615	CHARACTER(len=*), INTENT(in) :: path !! Path of the netcdf file.
1616	TYPE(DSET5D), INTENT(in) :: dset !! Dataset to write
1617	LOGICAL :: ret !! Return status
1618	INTEGER :: fi,vi,nd,ds,iret
1619	INTEGER, DIMENSION(:), ALLOCATABLE :: di
1620	REAL(kind=8), DIMENSION(:), ALLOCATABLE :: tv
1621	CHARACTER(len=15) :: i2s
1622	LOGICAL :: hxd,hyd,hzd,htd,hwd,ok
1623	INTEGER :: xid,yid,zid,tid,wid,vid
1624	ret = has_data(dset)
1625	IF (.NOT.ret) RETURN
1626	ret = .false.
1627	INQUIRE(FILE=TRIM(path),EXIST=ok)
1628	IF (ok) THEN
1629	IF (NF90_OPEN(TRIM(path), NF90_WRITE, fi) /= NF90_NOERR) RETURN
1630	ELSE
1631	IF (NF90_CREATE(TRIM(path), NF90_NOCLOBBER, fi) /= NF90_NOERR) RETURN
1632	ENDIF
1633	iret = NF90_CLOSE(fi)
1634
1635	! if variable already exist: get out of here !
1636	IF (get_nc_info(path,dset%dname,fi,vi,di)) RETURN
1637
1638	iret = NF90_OPEN(path,NF90_WRITE,fi)
1639
1640	hxd = nc_get_dim_by_name(fi,dset%xname,tv,ds,xid)
1641	IF (hxd) THEN
1642	IF(.NOT.compare(dset%x,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1643	ENDIF
1644	hyd = nc_get_dim_by_name(fi,dset%yname,tv,ds,yid)
1645	IF (hyd) THEN
1646	IF(.NOT.compare(dset%y,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1647	ENDIF
1648	hzd = nc_get_dim_by_name(fi,dset%zname,tv,ds,zid)
1649	IF (hzd) THEN
1650	IF(.NOT.compare(dset%z,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1651	ENDIF
1652	htd = nc_get_dim_by_name(fi,dset%tname,tv,ds,tid)
1653	IF (htd) THEN
1654	IF(.NOT.compare(dset%t,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1655	ENDIF
1656	hwd = nc_get_dim_by_name(fi,dset%wname,tv,ds,wid)
1657	IF (hwd) THEN
1658	IF (.NOT.compare(dset%w,tv)) THEN ; ret = .false. ; RETURN ; ENDIF
1659	ENDIF
1660
1661	IF (.NOT.hxd) THEN
1662	ret = nc_set_dim(fi,dset%xname,dset%x,xid)
1663	IF (.NOT.ret) THEN ; iret= NF90_CLOSE(fi) ; RETURN ; ENDIF
1664	ENDIF
1665	IF (.NOT.hyd) THEN
1666	ret = nc_set_dim(fi,dset%yname,dset%y,yid)
1667	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1668	ENDIF
1669	IF (.NOT.hzd) THEN
1670	ret = nc_set_dim(fi,dset%zname,dset%z,zid)
1671	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1672	ENDIF
1673	IF (.NOT.htd) THEN
1674	ret = nc_set_dim(fi,dset%tname,dset%t,tid)
1675	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1676	ENDIF
1677	IF (.NOT.hwd) THEN
1678	ret = nc_set_dim(fi,dset%wname,dset%w,wid)
1679	IF (.NOT.ret) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1680	ENDIF
1681	ret = .false.
1682	iret = nf90_redef(fi)
1683	iret = nf90_def_var(fi, TRIM(dset%dname), nf90_double,(/xid,yid,zid,tid,wid/),vid)
1684	IF (iret /= 0) THEN ; iret = NF90_CLOSE(fi) ; RETURN ; ENDIF
1685	iret = nf90_enddef(fi)
1686	iret = nf90_put_var(fi,vid,dset%data)
1687	ret = (iret == 0)
1688	iret=NF90_CLOSE(fi)
1689	ret = .true.
1690	RETURN
1691	END FUNCTION nc_wr_5d
1692
1693	FUNCTION nc_set_dim(fid,name,values,dimid) RESULT(ret)
1694	!! Set a new dimension (and its associated values) in a NetCDF file.
1695	INTEGER, INTENT(in) :: fid !! Netcdf file id.
1696	CHARACTER(len=*), INTENT(in) :: name !! Name of the dimension to set.
1697	REAL(kind=8), DIMENSION(:), INTENT(in) :: values !! Associated values.
1698	INTEGER, INTENT(out) :: dimid !! Output dimension id.
1699	LOGICAL :: ret !! Return status
1700	INTEGER :: iret,vid
1701	ret = .false.
1702	iret = nf90_redef(fid)
1703	iret = nf90_def_dim(fid, TRIM(name), size(values), dimid)
1704	WRITE(*,'(a)') "nc_set_dim: "//TRIM(name)//" --> "//TRIM(nf90_strerror(iret))
1705	IF (iret /=0) RETURN
1706	iret = nf90_def_var(fid, TRIM(name), nf90_double,(/dimid/),vid)
1707	IF (iret /=0) RETURN
1708	iret = nf90_enddef(fid)
1709	iret = nf90_put_var(fid, vid,values)
1710	ret = .true.
1711	END FUNCTION nc_set_dim
1712
1713	FUNCTION compare(vec1,vec2,tol) RESULT(ret)
1714	!! Compare two vector of double.
1715	!!
1716	!! The test is checked against the difference (element wise) of the two vector compared to
1717	!! to a given tolerance.
1718	REAL(kind=8), DIMENSION(:), INTENT(in) :: vec1 !! First vector to compare
1719	REAL(kind=8), DIMENSION(:), INTENT(in) :: vec2 !! Second vector to compare
1720	REAL(kind=8), INTENT(in), OPTIONAL :: tol !! Tolerance to apply for the comparison
1721	LOGICAL :: ret !! .true. if both vectors are equal, .false. otherwise.
1722	REAL(kind=8) :: ztol
1723	INTEGER :: i
1724	ztol = 1d-10 ; IF (PRESENT(tol)) ztol = abs(tol)
1725	ret = .false.
1726	IF (size(vec1) /= size(vec2)) RETURN
1727
1728	DO i=1,size(vec1) ; IF (abs(vec1(i)-vec2(i)) > ztol) RETURN ; ENDDO
1729	ret = .true.
1730	END FUNCTION compare
1731
1732
1733
1734	END MODULE DATASETS

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

Download in other formats:

Original Format