Context Navigation

easy_netcdf.F90 @ 4484

Last change on this file since 4484 was 3908, checked in by idelkadi, 3 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

Inclusion of the ecrad directory containing the sources of the ECRAD code
- interface routine : radiation_scheme.F90
Adaptation of compilation scripts :
- compilation under CPP key CPP_ECRAD
- compilation with option "-rad ecard" or "-ecard true"
- The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
Runing LMDZ simulations with ecrad, you need :
- logical key iflag_rrtm = 2 in physiq.def
- namelist_ecrad (DefLists?)
- the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...

File size: 83.9 KB

Line
1	! easy_netcdf.F90 - Module providing convenient NetCDF read/write capability
2	!
3	! (C) Copyright 2014- ECMWF.
4	!
5	! This software is licensed under the terms of the Apache Licence Version 2.0
6	! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7	!
8	! In applying this licence, ECMWF does not waive the privileges and immunities
9	! granted to it by virtue of its status as an intergovernmental organisation
10	! nor does it submit to any jurisdiction.
11	!
12	! Author: Robin Hogan
13	! Email: r.j.hogan@ecmwf.int
14	!
15	! Modifications
16	! 2017-04-28 R. Hogan Fix "reshape" when writing 3D array
17	! 2017-10-23 A. Bozzo Reading 4-D array
18	! 2018-03-14 R. Hogan Fix "reshape" properly this time
19	! 2019-01-04 R. Hogan Allow reading and writing a slice of a larger array
20	! 2019-01-07 R. Hogan HDF5 writing support, allowing larger files, provided NC_NETCDF4 defined
21	! 2019-01-16 R. Hogan Revised interpretation of "iverbose"
22	! 2019-06-17 R. Hogan Pass through deflate_level and shuffle to variable definition
23
24	module easy_netcdf
25
26	use netcdf
27	use parkind1, only : jprb, jpib, jprm, jprd
28	use radiation_io, only : nulout, nulerr, my_abort => radiation_abort
29
30	implicit none
31	public
32
33	!---------------------------------------------------------------------
34	! An object of this type provides convenient read or write access to
35	! a NetCDF file
36	type netcdf_file
37	integer :: ncid ! NetCDF file ID
38	integer :: iverbose ! Verbosity: 0 = report only fatal errors,
39	! 1 = ...and warnings,
40	! 2 = ...and when opening files,
41	! 3 = ...and when reading/writing variables,
42	! 4 = ...and variable attributes and when writing dimensions,
43	! 5 = ...and debugging information
44	logical :: do_transpose_2d = .false. ! Transpose 2D arrays on read/write?
45	logical :: is_write_mode = .false. ! .false. for read, .true. for write
46	logical :: is_define_mode = .true. ! .true. if in NetCDF define mode
47	logical :: is_double_precision = .false. ! Write reals in double precision?
48	logical :: do_permute_3d = .false. ! Permute 3D arrays on write?
49	logical :: do_permute_4d = .false. ! Permute 3D arrays on write?
50	integer :: i_permute_3d(3) = (/1,2,3/) ! New order of dimensions
51	integer :: i_permute_4d(4) = (/1,2,3,4/) ! New order of dimensions
52	character(len=511) :: file_name
53	contains
54	procedure :: open => open_netcdf_file
55	procedure :: create => create_netcdf_file
56	procedure :: close => close_netcdf_file
57	procedure :: get_real_scalar
58	procedure :: get_real_vector
59	procedure :: get_integer_vector
60	procedure :: get_real_matrix
61	procedure :: get_real_array3
62	procedure :: get_real_scalar_indexed
63	procedure :: get_real_vector_indexed
64	procedure :: get_real_matrix_indexed
65	procedure :: get_real_array3_indexed
66	procedure :: get_real_array4
67	generic :: get => get_real_scalar, get_real_vector, &
68	& get_real_matrix, get_real_array3, &
69	& get_real_array4, get_integer_vector, &
70	& get_real_scalar_indexed, get_real_vector_indexed, &
71	& get_real_matrix_indexed, get_real_array3_indexed
72	procedure :: get_real_scalar_attribute
73	procedure :: get_string_attribute
74	generic :: get_attribute => get_real_scalar_attribute, &
75	& get_string_attribute
76	procedure :: get_global_attribute
77
78	procedure :: define_dimension
79	procedure :: define_variable
80	procedure :: put_attribute
81	procedure :: put_global_attributes
82	procedure :: put_global_attribute
83	procedure :: put_real_scalar
84	procedure :: put_real_vector
85	procedure :: put_real_matrix
86	procedure :: put_real_array3
87	procedure :: put_real_scalar_indexed
88	procedure :: put_real_vector_indexed
89	procedure :: put_real_matrix_indexed
90	generic :: put => put_real_scalar, put_real_vector, &
91	& put_real_matrix, put_real_array3, &
92	& put_real_scalar_indexed, put_real_vector_indexed, &
93	& put_real_matrix_indexed
94	procedure :: set_verbose
95	procedure :: transpose_matrices
96	procedure :: double_precision
97	procedure :: permute_3d_arrays
98	procedure :: get_rank
99	procedure :: exists
100	procedure :: get_outer_dimension
101	procedure :: attribute_exists
102	procedure :: global_attribute_exists
103	procedure, private :: get_array_dimensions
104	procedure, private :: get_variable_id
105	procedure, private :: end_define_mode
106	procedure, private :: print_variable_attributes
107	end type netcdf_file
108
109	contains
110
111	! --- GENERIC SUBROUTINES ---
112
113	!---------------------------------------------------------------------
114	! Open a NetCDF file with name "file_name", optionally specifying the
115	! verbosity level (0-5) and if the file is for writing (the default
116	! is read-only)
117	subroutine open_netcdf_file(this, file_name, iverbose, is_write_mode, is_hdf5_file)
118	class(netcdf_file) :: this
119	character(len=*), intent(in) :: file_name
120	integer, intent(in), optional :: iverbose
121	logical, intent(in), optional :: is_write_mode
122	logical, intent(in), optional :: is_hdf5_file ! Only for write mode
123
124	integer :: istatus
125	integer :: i_write_mode
126
127
128	! Store verbosity level in object
129	if (present(iverbose)) then
130	this%iverbose = iverbose
131	else
132	! By default announce files being opened and closed, but not
133	! variables read/written
134	this%iverbose = 2
135	end if
136
137	! Store read/write mode in object
138	if (present(is_write_mode)) then
139	this%is_write_mode = is_write_mode
140	else
141	this%is_write_mode = .false.
142	end if
143
144	! By default we don't transpose 2D arrays on read/write
145	this%do_transpose_2d = .false.
146
147	! Store filename
148	this%file_name = file_name
149
150	! Open file according to write mode
151	if (.not. this%is_write_mode) then
152	istatus = nf90_open(file_name, NF90_NOWRITE, this%ncid)
153	if (this%iverbose >= 2) then
154	write(nulout,'(a,a)') 'Reading NetCDF file ', file_name
155	!write(nulout,'(a,a,a,i0,a)') 'Reading NetCDF file ', file_name, ' (ID=', this%ncid, ')'
156	end if
157	this%is_define_mode = .false.
158	else
159	i_write_mode = NF90_CLOBBER
160	! Check if HDF5 file is to be written (which can be larger)
161	if (present(is_hdf5_file)) then
162	if (is_hdf5_file) then
163	#ifdef NC_NETCDF4
164	i_write_mode = ior(i_write_mode, NF90_HDF5)
165	#else
166	if (this%iverbose >= 1) then
167	write(nulout,'(a,a)') 'Warning: cannot use HDF5 format for writing ', file_name, &
168	& ' unless compiled with NC_NETCDF4 defined'
169	end if
170	#endif
171	end if
172	end if
173
174	istatus = nf90_create(file_name, i_write_mode, this%ncid)
175	if (this%iverbose >= 2) then
176	write(nulout,'(a,a)') 'Writing NetCDF file ', file_name
177	end if
178	this%is_define_mode = .true.
179	end if
180
181	! Check the file opened correctly
182	if (istatus /= NF90_NOERR) then
183	write(nulerr,'(a,a,a,a)') '*** Error opening NetCDF file ', &
184	& file_name, ': ', trim(nf90_strerror(istatus))
185	call my_abort('Error opening NetCDF file')
186	end if
187
188	end subroutine open_netcdf_file
189
190
191	!---------------------------------------------------------------------
192	! Open a NetCDF file for writing
193	subroutine create_netcdf_file(this, file_name, iverbose, is_hdf5_file)
194	class(netcdf_file) :: this
195	character(len=*), intent(in) :: file_name
196	integer, intent(in), optional :: iverbose
197	logical, intent(in), optional :: is_hdf5_file
198
199	integer :: istatus
200	integer :: i_write_mode
201
202	if (present(iverbose)) then
203	this%iverbose = iverbose
204	else
205	this%iverbose = 2
206	end if
207
208	this%do_transpose_2d = .false.
209
210	i_write_mode = NF90_CLOBBER
211	! Check if HDF5 file is to be written (which can be large)
212	if (present(is_hdf5_file)) then
213	if (is_hdf5_file) then
214	#ifdef NC_NETCDF4
215	i_write_mode = ior(i_write_mode, NF90_HDF5)
216	#else
217	if (this%iverbose >= 1) then
218	write(nulout,'(a,a)') 'Warning: cannot use HDF5 format for writing ', file_name, &
219	& ' unless compiled with NC_NETCDF4 defined'
220	end if
221	#endif
222	end if
223	end if
224
225	istatus = nf90_create(file_name, i_write_mode, this%ncid)
226	if (this%iverbose >= 2) then
227	write(nulout,'(a,a)') 'Writing NetCDF file ', file_name
228	!write(nulout,'(a,a,a,i0,a)') 'Writing NetCDF file ', file_name, ' (ID=', this%ncid, ')'
229	end if
230	this%is_define_mode = .true.
231
232	if (istatus /= NF90_NOERR) then
233	write(nulerr,'(a,a,a)') '*** Error opening NetCDF file ', file_name, &
234	& ': ', trim(nf90_strerror(istatus))
235	stop
236	end if
237	this%file_name = file_name
238
239	end subroutine create_netcdf_file
240
241
242	!---------------------------------------------------------------------
243	! Close the NetCDF file
244	subroutine close_netcdf_file(this)
245	class(netcdf_file) :: this
246	integer :: istatus
247
248	if (this%iverbose >= 3) then
249	write(nulout,'(a,a)') 'Closing NetCDF file ', trim(this%file_name)
250	end if
251
252	istatus = nf90_close(this%ncid)
253	if (istatus /= NF90_NOERR) then
254	write(nulerr,'(a,a,a,a)') '*** Error closing NetCDF file ', &
255	& trim(this%file_name), ': ', trim(nf90_strerror(istatus))
256	stop
257	end if
258
259	end subroutine close_netcdf_file
260
261
262	!---------------------------------------------------------------------
263	! Set the verbosity level from 0 to 5, where the codes have the
264	! following meaning: 0=errors only, 1=warning, 2=info, 3=progress,
265	! 4=detailed, 5=debug
266	subroutine set_verbose(this, ival)
267	class(netcdf_file) :: this
268	integer, optional :: ival
269
270	if (present(ival)) then
271	this%iverbose = ival
272	else
273	this%iverbose = 2
274	end if
275
276	end subroutine set_verbose
277
278
279	!---------------------------------------------------------------------
280	! Specify whether floating-point arrays should be written in double precision
281	subroutine double_precision(this, is_double)
282	class(netcdf_file) :: this
283	logical, optional :: is_double
284
285	if (present(is_double)) then
286	this%is_double_precision = is_double
287	else
288	this%is_double_precision = .true.
289	end if
290
291	end subroutine double_precision
292
293
294	!---------------------------------------------------------------------
295	! Specify whether 2D arrays should be transposed on read/write
296	subroutine transpose_matrices(this, do_transpose)
297	class(netcdf_file) :: this
298	logical, optional :: do_transpose
299
300	if (present(do_transpose)) then
301	this%do_transpose_2d = do_transpose
302	else
303	this%do_transpose_2d = .true.
304	end if
305
306	end subroutine transpose_matrices
307
308
309	!---------------------------------------------------------------------
310	! Specify that 3D arrays should be permuted on write, with the new
311	! dimension order in the input argument "ipermute" (e.g. 3,2,1)
312	subroutine permute_3d_arrays(this, ipermute)
313	class(netcdf_file) :: this
314	integer, intent(in) :: ipermute(3)
315
316	this%do_permute_3d = .true.
317	this%i_permute_3d = ipermute
318
319	end subroutine permute_3d_arrays
320
321
322	!---------------------------------------------------------------------
323	! Specify that 4D arrays should be permuted on write, with the new
324	! dimension order in the input argument "ipermute" (e.g. 4,3,2,1)
325	subroutine permute_4d_arrays(this, ipermute)
326	class(netcdf_file) :: this
327	integer, intent(in) :: ipermute(4)
328
329	this%do_permute_4d = .true.
330	this%i_permute_4d = ipermute
331
332	end subroutine permute_4d_arrays
333
334
335	! --- PRIVATE SUBROUTINES ---
336
337	!---------------------------------------------------------------------
338	! Return the NetCDF variable ID for variable "var_name", or abort if
339	! not present
340	subroutine get_variable_id(this, var_name, ivarid)
341	class(netcdf_file) :: this
342	character(len=*), intent(in) :: var_name
343	integer, intent(out) :: ivarid
344
345	integer :: istatus
346
347	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
348	if (istatus /= NF90_NOERR) then
349	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
350	& var_name, ': ', trim(nf90_strerror(istatus))
351	call my_abort('Error reading NetCDF file')
352	end if
353
354	end subroutine get_variable_id
355
356
357	!---------------------------------------------------------------------
358	! Return the array dimensions of variable with specified ID, along
359	! with the number of dimensions and optionally the total number of
360	! elements, or abort if variable not present
361	subroutine get_array_dimensions(this, ivarid, ndims, ndimlens, ntotal)
362	class(netcdf_file) :: this
363	integer, intent(in) :: ivarid
364	integer, intent(out) :: ndims
365	integer, intent(out) :: ndimlens(NF90_MAX_VAR_DIMS)
366	integer(kind=jpib), intent(out), optional :: ntotal
367
368	integer :: j, istatus
369	integer :: dimids(NF90_MAX_VAR_DIMS)
370
371	istatus = nf90_inquire_variable(this%ncid, ivarid, &
372	& ndims=ndims, dimids=dimids)
373	if (istatus /= NF90_NOERR) then
374	write(nulerr,'(a,i0,a,a)') '*** Error inquiring about NetCDF variable with id ', &
375	& ivarid, ': ', trim(nf90_strerror(istatus))
376	call my_abort('Error reading NetCDF file')
377	end if
378
379	ndimlens(:) = 0
380	do j = 1,ndims
381	istatus = nf90_inquire_dimension(this%ncid, dimids(j), len=ndimlens(j))
382	if (istatus /= NF90_NOERR) then
383	write(nulerr,'(a,i0,a,i0,a,a)') '*** Error reading length of dimension ', &
384	& j, ' of NetCDF variable with id ', ivarid, ': ', &
385	& trim(nf90_strerror(istatus))
386	call my_abort('Error reading NetCDF file')
387	end if
388	end do
389
390	if (present(ntotal)) then
391	ntotal = 1
392	do j = 1, ndims
393	ntotal = ntotal * ndimlens(j)
394	end do
395	end if
396
397	end subroutine get_array_dimensions
398
399
400	!---------------------------------------------------------------------
401	! End define mode, if in define mode, and check the error code
402	! (errors are possible if variables are too large for the format,
403	! for example)
404	subroutine end_define_mode(this)
405	class(netcdf_file) :: this
406	integer :: istatus
407	if (this%is_define_mode) then
408	istatus = nf90_enddef(this%ncid)
409	if (istatus /= NF90_NOERR) then
410	write(nulerr,'(a,a,a,a)') '*** Error ending define mode when writing ', &
411	& trim(this%file_name), ': ', trim(nf90_strerror(istatus))
412	call my_abort('Error writing NetCDF file')
413	end if
414	this%is_define_mode = .false.
415	end if
416	end subroutine end_define_mode
417
418
419	! --- READING SUBROUTINES ---
420
421	!---------------------------------------------------------------------
422	! Return the number of dimensions of variable with name var_name, or
423	! -1 if the variable is not found
424	function get_rank(this, var_name) result(ndims)
425	class(netcdf_file) :: this
426	character(len=*), intent(in) :: var_name
427
428	integer :: ndims
429	integer :: ivarid
430	integer :: ndimlens(NF90_MAX_VAR_DIMS)
431	integer :: istatus
432
433	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
434	if (istatus /= NF90_NOERR) then
435	if (istatus == NF90_ENOTVAR) then
436	if (this%iverbose >= 1) then
437	write(nulout,'(a,a,a)') ' Warning: variable ', var_name, ' not found'
438	end if
439	else
440	write(nulerr,'(a,a,a)') '*** Error inquiring about variable ', &
441	& var_name, ': ', trim(nf90_strerror(istatus))
442	call my_abort('Error reading NetCDF file')
443	end if
444	ndims = -1
445	else
446	call this%get_array_dimensions(ivarid, ndims, ndimlens)
447	end if
448
449	end function get_rank
450
451
452	!---------------------------------------------------------------------
453	! Return the length of the slowest-varying dimension of variable
454	! with name var_name, or -1 if the variable is not found
455	function get_outer_dimension(this, var_name) result(n)
456	class(netcdf_file) :: this
457	character(len=*), intent(in) :: var_name
458
459	integer :: n
460	integer :: ndims
461	integer :: ivarid
462	integer :: ndimlens(NF90_MAX_VAR_DIMS)
463	integer :: istatus
464
465	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
466	if (istatus /= NF90_NOERR) then
467	if (istatus == NF90_ENOTVAR) then
468	if (this%iverbose >= 1) then
469	write(nulout,'(a,a,a)') ' Warning: variable ', var_name, ' not found'
470	end if
471	else
472	write(nulerr,'(a,a,a)') '*** Error inquiring about variable ', &
473	& var_name, ': ', trim(nf90_strerror(istatus))
474	call my_abort('Error reading NetCDF file')
475	end if
476	n = -1
477	else
478	call this%get_array_dimensions(ivarid, ndims, ndimlens)
479	n = ndimlens(ndims)
480	end if
481
482	end function get_outer_dimension
483
484
485	!---------------------------------------------------------------------
486	! Return true if the variable is present, false otherwise
487	function exists(this, var_name) result(is_present)
488	class(netcdf_file) :: this
489	character(len=*), intent(in) :: var_name
490
491	logical :: is_present
492
493	integer :: ivarid
494	integer :: istatus
495
496	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
497	if (istatus /= NF90_NOERR) then
498	is_present = .false.
499	else
500	is_present = .true.
501	end if
502
503	end function exists
504
505
506	!---------------------------------------------------------------------
507	! Return true if the attribute is present, false otherwise. If
508	! argument "len" is provided then return false if len is smaller
509	! than the length of the attribute. This is useful if you have a
510	! fixed array size and want to check whether the attribute will fit
511	! into it.
512	function attribute_exists(this, var_name, attr_name, len) result(is_present)
513	class(netcdf_file) :: this
514	character(len=*), intent(in) :: var_name, attr_name
515	integer, optional, intent(in) :: len
516
517	logical :: is_present
518	integer :: i_attr_len, ivarid
519	integer :: istatus
520
521	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
522	if (istatus /= NF90_NOERR) then
523	is_present = .false.
524	else
525	istatus = nf90_inquire_attribute(this%ncid, ivarid, attr_name, &
526	& len=i_attr_len)
527	if (istatus /= NF90_NOERR) then
528	is_present = .false.
529	else
530	is_present = .true.
531	if (present(len)) then
532	if (i_attr_len > len) then
533	is_present = .false.
534	end if
535	end if
536	end if
537	end if
538
539	end function attribute_exists
540
541
542	!---------------------------------------------------------------------
543	! Return true if the global attribute is present, false otherwise.
544	! If argument "len" is provided then return false if len is smaller
545	! than the length of the attribute. This is useful if you have a
546	! fixed array size and want to check whether the attribute will fit
547	! into it.
548	function global_attribute_exists(this, attr_name, len) result(is_present)
549	class(netcdf_file) :: this
550	character(len=*), intent(in) :: attr_name
551	integer, optional, intent(in) :: len
552
553	logical :: is_present
554	integer :: i_attr_len
555	integer :: istatus
556
557	istatus = nf90_inquire_attribute(this%ncid, NF90_GLOBAL, attr_name, &
558	& len=i_attr_len)
559	if (istatus /= NF90_NOERR) then
560	is_present = .false.
561	else
562	is_present = .true.
563	if (present(len)) then
564	if (i_attr_len > len) then
565	is_present = .false.
566	end if
567	end if
568	end if
569
570	end function global_attribute_exists
571
572
573	!---------------------------------------------------------------------
574	! The method "get" will read either a scalar, vector or matrix
575	! depending on the rank of the output argument. This version reads a
576	! scalar.
577	subroutine get_real_scalar(this, var_name, scalar)
578	class(netcdf_file) :: this
579	character(len=*), intent(in) :: var_name
580	real(jprb), intent(out) :: scalar
581
582	integer :: istatus
583	integer :: ivarid, ndims
584	integer :: ndimlens(NF90_MAX_VAR_DIMS)
585	integer :: j, ntotal
586
587	! Inquire the ID, shape & size of the variable
588	call this%get_variable_id(var_name, ivarid)
589	call this%get_array_dimensions(ivarid, ndims, ndimlens)
590
591	! Compute number of elements of the variable in the file
592	ntotal = 1
593	do j = 1, ndims
594	ntotal = ntotal * ndimlens(j)
595	end do
596
597	if (this%iverbose >= 3) then
598	write(nulout,'(a,a)',advance='no') ' Reading ', var_name
599	call this%print_variable_attributes(ivarid,nulout)
600	end if
601
602	! Abort if the number of elements is anything other than 1
603	if (ntotal /= 1) then
604	write(nulerr,'(a,a,a,i0,a)') '*** Error reading NetCDF variable ', &
605	& var_name, ' with total length ', ntotal, ' as a scalar'
606	call my_abort('Error reading NetCDF file')
607	end if
608
609	! Read variable
610	istatus = nf90_get_var(this%ncid, ivarid, scalar)
611	if (istatus /= NF90_NOERR) then
612	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
613	& var_name, ' as a scalar: ', trim(nf90_strerror(istatus))
614	call my_abort('Error reading NetCDF file')
615	end if
616
617	end subroutine get_real_scalar
618
619
620	!---------------------------------------------------------------------
621	! Read a scalar from a larger array, where "index" indexes the most
622	! slowly varying dimension
623	subroutine get_real_scalar_indexed(this, var_name, index, scalar)
624	class(netcdf_file) :: this
625	character(len=*), intent(in) :: var_name
626	integer, intent(in) :: index
627	real(jprb), intent(out) :: scalar
628
629	integer :: istatus
630	integer :: ivarid, ndims
631	integer :: ndimlens(NF90_MAX_VAR_DIMS)
632	integer :: vstart(NF90_MAX_VAR_DIMS)
633	integer :: j, ntotal
634
635	! Inquire the ID, shape & size of the variable
636	call this%get_variable_id(var_name, ivarid)
637	call this%get_array_dimensions(ivarid, ndims, ndimlens)
638
639	! Compute number of elements of the slice in the file,
640	! i.e. excluding the slowest varying dimension, indexed by "index"
641	ntotal = 1
642	do j = 1, ndims-1
643	ntotal = ntotal * ndimlens(j)
644	end do
645
646	if (this%iverbose >= 3) then
647	write(nulout,'(a,a,i0,a,a)') ' Reading element ', index, ' of ', var_name
648	end if
649
650	! Abort if the number of elements is anything other than 1
651	if (ntotal /= 1) then
652	write(nulerr,'(a,a,a,i0,a)') '*** Error reading NetCDF variable ', &
653	& var_name, ' with slice length ', ntotal, ' as a scalar'
654	call my_abort('Error reading NetCDF file')
655	end if
656
657	if (index < 1 .or. index > ndimlens(ndims)) then
658	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error reading element ', index, &
659	& ' of NetCDF variable ', &
660	& var_name, ' with outer dimension ', ndimlens(ndims)
661	call my_abort('Error reading NetCDF file')
662	end if
663
664	! Read variable
665	vstart(1:ndims-1) = 1
666	vstart(ndims) = index
667	istatus = nf90_get_var(this%ncid, ivarid, scalar, start=vstart)
668	if (istatus /= NF90_NOERR) then
669	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
670	& var_name, ' as a scalar: ', trim(nf90_strerror(istatus))
671	call my_abort('Error reading NetCDF file')
672	end if
673
674	end subroutine get_real_scalar_indexed
675
676
677	!---------------------------------------------------------------------
678	! Read a 1D array into "vector", which must be allocatable and will
679	! be reallocated if necessary
680	subroutine get_real_vector(this, var_name, vector)
681	class(netcdf_file) :: this
682	character(len=*), intent(in) :: var_name
683	real(jprb), allocatable, intent(out) :: vector(:)
684
685	integer :: n ! Length of vector
686	integer :: istatus
687	integer :: ivarid, ndims
688	integer :: ndimlens(NF90_MAX_VAR_DIMS)
689	integer :: j
690
691	call this%get_variable_id(var_name, ivarid)
692	call this%get_array_dimensions(ivarid, ndims, ndimlens)
693
694	! Ensure variable has only one dimension in the file
695	n = 1
696	do j = 1, ndims
697	n = n * ndimlens(j)
698	if (j > 1 .and. ndimlens(j) > 1) then
699	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
700	& var_name, &
701	& ' as a vector: all dimensions above the first must be singletons'
702	call my_abort('Error reading NetCDF file')
703	end if
704	end do
705
706	! Reallocate if necessary
707	if (allocated(vector)) then
708	if (size(vector) /= n) then
709	if (this%iverbose >= 1) then
710	write(nulout,'(a,a)') ' Warning: resizing vector to read ', var_name
711	end if
712	deallocate(vector)
713	allocate(vector(n))
714	end if
715	else
716	allocate(vector(n))
717	end if
718
719	if (this%iverbose >= 3) then
720	write(nulout,'(a,a,a,i0,a)',advance='no') ' Reading ', var_name, '(', n, ')'
721	call this%print_variable_attributes(ivarid,nulout)
722	end if
723
724	! Read variable
725	istatus = nf90_get_var(this%ncid, ivarid, vector)
726	if (istatus /= NF90_NOERR) then
727	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
728	& var_name, ' as a vector: ', trim(nf90_strerror(istatus))
729	call my_abort('Error reading NetCDF file')
730	end if
731
732	end subroutine get_real_vector
733
734
735	!---------------------------------------------------------------------
736	! Read a 1D integer array into "vector", which must be allocatable
737	! and will be reallocated if necessary
738	subroutine get_integer_vector(this, var_name, vector)
739	class(netcdf_file) :: this
740	character(len=*), intent(in) :: var_name
741	integer, allocatable, intent(out) :: vector(:)
742
743	integer :: n ! Length of vector
744	integer :: istatus
745	integer :: ivarid, ndims
746	integer :: ndimlens(NF90_MAX_VAR_DIMS)
747	integer :: j
748
749	call this%get_variable_id(var_name, ivarid)
750	call this%get_array_dimensions(ivarid, ndims, ndimlens)
751
752	! Ensure variable has only one dimension in the file
753	n = 1
754	do j = 1, ndims
755	n = n * ndimlens(j)
756	if (j > 1 .and. ndimlens(j) > 1) then
757	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
758	& var_name, &
759	& ' as a vector: all dimensions above the first must be singletons'
760	call my_abort('Error reading NetCDF file')
761	end if
762	end do
763
764	! Reallocate if necessary
765	if (allocated(vector)) then
766	if (size(vector) /= n) then
767	if (this%iverbose >= 1) then
768	write(nulout,'(a,a)') ' Warning: resizing vector to read ', var_name
769	end if
770	deallocate(vector)
771	allocate(vector(n))
772	end if
773	else
774	allocate(vector(n))
775	end if
776
777	if (this%iverbose >= 3) then
778	write(nulout,'(a,a,a,i0,a)',advance='no') ' Reading ', var_name, '(', n, ')'
779	call this%print_variable_attributes(ivarid,nulout)
780	end if
781
782	! Read variable
783	istatus = nf90_get_var(this%ncid, ivarid, vector)
784	if (istatus /= NF90_NOERR) then
785	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
786	& var_name, ' as an integer vector: ', trim(nf90_strerror(istatus))
787	call my_abort('Error reading NetCDF file')
788	end if
789
790	end subroutine get_integer_vector
791
792
793	!---------------------------------------------------------------------
794	! Read a vector of data from a larger array; the vector must be
795	! allocatable and will be reallocated if necessary
796	subroutine get_real_vector_indexed(this, var_name, index, vector)
797	class(netcdf_file) :: this
798	character(len=*), intent(in) :: var_name
799	integer, intent(in) :: index
800	real(jprb), allocatable, intent(out) :: vector(:)
801
802	integer :: n ! Length of vector
803	integer :: istatus
804	integer :: ivarid, ndims
805	integer :: ndimlens(NF90_MAX_VAR_DIMS)
806	integer :: vstart(NF90_MAX_VAR_DIMS)
807	integer :: vcount(NF90_MAX_VAR_DIMS)
808	integer :: j
809
810	call this%get_variable_id(var_name, ivarid)
811	call this%get_array_dimensions(ivarid, ndims, ndimlens)
812
813	! Ensure variable has only one dimension aside from the last one
814	n = 1
815	do j = 1, ndims-1
816	n = n * ndimlens(j)
817	if (j > 1 .and. ndimlens(j) > 1) then
818	write(nulerr,'(a,a,a)') '*** Error reading 1D slice from NetCDF variable ', &
819	& var_name, &
820	& ': all dimensions except the first and last must be singletons'
821	call my_abort('Error reading NetCDF file')
822	end if
823	end do
824
825	! Reallocate if necessary
826	if (allocated(vector)) then
827	if (size(vector) /= n) then
828	if (this%iverbose >= 1) then
829	write(nulout,'(a,i0,a,a)') ' Warning: resizing vector to length ', n, &
830	& ' to read slice of ', var_name
831	end if
832	deallocate(vector)
833	allocate(vector(n))
834	end if
835	else
836	allocate(vector(n))
837	end if
838
839	if (this%iverbose >= 3) then
840	write(nulout,'(a,i0,a,a)') ' Reading column ', index, ' of ', var_name
841	end if
842
843	if (index < 1 .or. index > ndimlens(ndims)) then
844	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error reading element ', index, &
845	& ' of NetCDF variable ', &
846	& var_name, ' with outer dimension ', ndimlens(ndims)
847	call my_abort('Error reading NetCDF file')
848	end if
849
850	! Read variable
851	vstart(1:ndims-1) = 1
852	vstart(ndims) = index
853	vcount(1:ndims-1) = ndimlens(1:ndims-1)
854	vcount(ndims) = 1
855	istatus = nf90_get_var(this%ncid, ivarid, vector, &
856	& start=vstart, count=vcount)
857	if (istatus /= NF90_NOERR) then
858	write(nulerr,'(a,a,a,a)') '*** Error reading 1D slice of NetCDF variable ', &
859	& var_name, ': ', trim(nf90_strerror(istatus))
860	call my_abort('Error reading NetCDF file')
861	end if
862
863	end subroutine get_real_vector_indexed
864
865
866	!---------------------------------------------------------------------
867	! Read 2D array into "matrix", which must be allocatable and will be
868	! reallocated if necessary. Whether to transpose is specifed by the
869	! final optional argument, but can also be specified by the
870	! do_transpose_2d class data member.
871	subroutine get_real_matrix(this, var_name, matrix, do_transp)
872	class(netcdf_file) :: this
873	character(len=*), intent(in) :: var_name
874	real(jprb), allocatable, intent(out) :: matrix(:,:)
875	logical, optional, intent(in):: do_transp ! Transpose data?
876
877	real(jprb), allocatable :: tmp_matrix(:,:)
878	integer :: ndimlen1, ndimlen2
879	integer :: istatus
880	integer :: ivarid, ndims
881	integer :: ndimlens(NF90_MAX_VAR_DIMS)
882	integer :: j, ntotal
883	logical :: do_transpose
884
885	! Decide whether to transpose the array
886	if (present(do_transp)) then
887	do_transpose = do_transp
888	else
889	do_transpose = this%do_transpose_2d
890	end if
891
892	call this%get_variable_id(var_name, ivarid)
893	call this%get_array_dimensions(ivarid, ndims, ndimlens)
894
895	! Ensure the variable has no more than two non-singleton
896	! dimensions
897	ntotal = 1
898	do j = 1, ndims
899	ntotal = ntotal * ndimlens(j)
900	if (j > 2 .and. ndimlens(j) > 1) then
901	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
902	& var_name, &
903	& ' as a matrix: all dimensions above the second must be singletons'
904	call my_abort('Error reading NetCDF file')
905	end if
906	end do
907
908	! Work out dimension lengths
909	if (ndims >= 2) then
910	ndimlen1 = ndimlens(1)
911	ndimlen2 = ntotal/ndimlen1
912	else
913	ndimlen1 = ntotal
914	ndimlen2 = 1
915	end if
916
917	if (do_transpose) then
918	! Read and transpose
919	allocate(tmp_matrix(ndimlen1, ndimlen2))
920
921	! Reallocate if necessary
922	if (allocated(matrix)) then
923	if (size(matrix,1) /= ndimlen2 .or. size(matrix,2) /= ndimlen1) then
924	if (this%iverbose >= 1) then
925	write(nulout,'(a,a)') ' Warning: resizing matrix to read ', var_name
926	end if
927	deallocate(matrix)
928	allocate(matrix(ndimlen2, ndimlen1))
929	end if
930	else
931	allocate(matrix(ndimlen2, ndimlen1))
932	end if
933
934	if (this%iverbose >= 3) then
935	write(nulout,'(a,a,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, '(', &
936	& ndimlen2, ',', ndimlen1, ')'
937	call this%print_variable_attributes(ivarid,nulout)
938	end if
939
940	istatus = nf90_get_var(this%ncid, ivarid, tmp_matrix)
941	matrix = transpose(tmp_matrix)
942	deallocate(tmp_matrix)
943	else
944	! Read data without transposition
945
946	! Reallocate if necessary
947	if (allocated(matrix)) then
948	if (size(matrix,1) /= ndimlen1 .or. size(matrix,2) /= ndimlen2) then
949	if (this%iverbose >= 1) then
950	write(nulout,'(a,a)') ' Warning: resizing matrix to read ', var_name
951	end if
952	allocate(matrix(ndimlen1, ndimlen2))
953	end if
954	else
955	allocate(matrix(ndimlen1, ndimlen2))
956	end if
957
958	if (this%iverbose >= 3) then
959	write(nulout,'(a,a,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, '(', &
960	& ndimlen1, ',', ndimlen2, ')'
961	call this%print_variable_attributes(ivarid,nulout)
962	end if
963
964	istatus = nf90_get_var(this%ncid, ivarid, matrix)
965	end if
966
967	if (istatus /= NF90_NOERR) then
968	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
969	& var_name, ' as a matrix: ', trim(nf90_strerror(istatus))
970	call my_abort('Error reading NetCDF file')
971	end if
972
973	end subroutine get_real_matrix
974
975
976	!---------------------------------------------------------------------
977	! Read matrix of data from a larger array, which must be allocatable
978	! and will be reallocated if necessary. Whether to transpose is
979	! specifed by the final optional argument, but can also be specified
980	! by the do_transpose_2d class data member.
981	subroutine get_real_matrix_indexed(this, var_name, index, matrix, do_transp)
982	class(netcdf_file) :: this
983	character(len=*), intent(in) :: var_name
984	integer, intent(in) :: index
985	real(jprb), allocatable, intent(out) :: matrix(:,:)
986	logical, optional, intent(in):: do_transp ! Transpose data?
987
988	real(jprb), allocatable :: tmp_matrix(:,:)
989	integer :: ndimlen1, ndimlen2
990	integer :: istatus
991	integer :: ivarid, ndims
992	integer :: ndimlens(NF90_MAX_VAR_DIMS)
993	integer :: vstart(NF90_MAX_VAR_DIMS)
994	integer :: vcount(NF90_MAX_VAR_DIMS)
995	integer :: j, ntotal
996	logical :: do_transpose
997
998	! Decide whether to transpose the array
999	if (present(do_transp)) then
1000	do_transpose = do_transp
1001	else
1002	do_transpose = this%do_transpose_2d
1003	end if
1004
1005	call this%get_variable_id(var_name, ivarid)
1006	call this%get_array_dimensions(ivarid, ndims, ndimlens)
1007
1008	! Ensure the variable has no more than two non-singleton
1009	! dimensions aside from the last one
1010	ntotal = 1
1011	do j = 1, ndims-1
1012	ntotal = ntotal * ndimlens(j)
1013	if (j > 2 .and. ndimlens(j) > 1) then
1014	write(nulerr,'(a,a,a)') '*** Error reading 2D slice from NetCDF variable ', &
1015	& var_name, &
1016	& ': all dimensions except the first, second and last must be singletons'
1017	call my_abort('Error reading NetCDF file')
1018	end if
1019	end do
1020
1021	if (index < 1 .or. index > ndimlens(ndims)) then
1022	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error reading element ', index, &
1023	& ' of NetCDF variable ', &
1024	& var_name, ' with outer dimension ', ndimlens(ndims)
1025	call my_abort('Error reading NetCDF file')
1026	end if
1027
1028	! Work out dimension lengths
1029	if (ndims >= 2) then
1030	ndimlen1 = ndimlens(1)
1031	ndimlen2 = ntotal/ndimlen1
1032	else
1033	ndimlen1 = ntotal
1034	ndimlen2 = 1
1035	end if
1036
1037	vstart(1:ndims-1) = 1
1038	vstart(ndims) = index
1039	vcount(1:ndims-1) = ndimlens(1:ndims-1)
1040	vcount(ndims) = 1
1041
1042	if (do_transpose) then
1043	! Read and transpose
1044	allocate(tmp_matrix(ndimlen1, ndimlen2))
1045
1046	! Reallocate if necessary
1047	if (allocated(matrix)) then
1048	if (size(matrix,1) /= ndimlen2 .or. size(matrix,2) /= ndimlen1) then
1049	if (this%iverbose >= 1) then
1050	write(nulout,'(a,a)') ' Warning: resizing matrix to read ', var_name
1051	end if
1052	allocate(matrix(ndimlen2, ndimlen1))
1053	end if
1054	else
1055	allocate(matrix(ndimlen2, ndimlen1))
1056	end if
1057
1058	if (this%iverbose >= 3) then
1059	write(nulout,'(a,i0,a,a,a,i0,a,i0,a)') ' Reading slice ', index, &
1060	& ' of ', var_name, ' as ', ndimlen2, 'x', ndimlen1, ' array'
1061	end if
1062
1063	istatus = nf90_get_var(this%ncid, ivarid, tmp_matrix, &
1064	& start=vstart, count=vcount)
1065	matrix = transpose(tmp_matrix)
1066	deallocate(tmp_matrix)
1067	else
1068	! Read data without transposition
1069
1070	! Reallocate if necessary
1071	if (allocated(matrix)) then
1072	if (size(matrix,1) /= ndimlen1 .or. size(matrix,2) /= ndimlen2) then
1073	if (this%iverbose >= 1) then
1074	write(nulout,'(a,a)') ' Warning: resizing matrix to read ', var_name
1075	end if
1076	allocate(matrix(ndimlen1, ndimlen2))
1077	end if
1078	else
1079	allocate(matrix(ndimlen1, ndimlen2))
1080	end if
1081
1082	if (this%iverbose >= 3) then
1083	write(nulout,'(a,i0,a,a,a,i0,a,i0,a)') ' Reading slice ', index, &
1084	& ' of ', var_name, ' as ', ndimlen1, 'x', ndimlen2, ' array'
1085	end if
1086
1087	istatus = nf90_get_var(this%ncid, ivarid, matrix, &
1088	& start=vstart, count=vcount)
1089	end if
1090
1091	if (istatus /= NF90_NOERR) then
1092	write(nulerr,'(a,a,a,a)') '*** Error reading 2D slice of NetCDF variable ', &
1093	& var_name, ': ', trim(nf90_strerror(istatus))
1094	call my_abort('Error reading NetCDF file')
1095	end if
1096
1097	end subroutine get_real_matrix_indexed
1098
1099
1100	!---------------------------------------------------------------------
1101	! Read 3D array into "var", which must be allocatable and will be
1102	! reallocated if necessary. Whether to pemute is specifed by the
1103	! final optional argument
1104	subroutine get_real_array3(this, var_name, var, ipermute)
1105	class(netcdf_file) :: this
1106	character(len=*), intent(in) :: var_name
1107	real(jprb), allocatable, intent(out) :: var(:,:,:)
1108	integer, optional, intent(in) :: ipermute(3)
1109
1110	real(jprb), allocatable :: var_permute(:,:,:)
1111	integer :: ndimlen1, ndimlen2, ndimlen3
1112	integer :: istatus
1113	integer :: ivarid, ndims
1114	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1115	integer :: j, ntotal
1116	integer :: n_dimlens_permuted(3)
1117	integer :: i_permute_3d(3)
1118	logical :: do_permute
1119
1120	! Decide whether to permute
1121	if (present(ipermute)) then
1122	do_permute = .true.
1123	i_permute_3d = ipermute
1124	else
1125	do_permute = this%do_permute_3d
1126	i_permute_3d = this%i_permute_3d
1127	end if
1128
1129	call this%get_variable_id(var_name, ivarid)
1130	call this%get_array_dimensions(ivarid, ndims, ndimlens)
1131
1132	! Ensure the variable has no more than three non-singleton
1133	! dimensions
1134	ntotal = 1
1135	do j = 1, ndims
1136	ntotal = ntotal * ndimlens(j)
1137	if (j > 3 .and. ndimlens(j) > 1) then
1138	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
1139	& var_name, &
1140	& ' as a 3D array: all dimensions above the third must be singletons'
1141	call my_abort('Error reading NetCDF file')
1142	end if
1143	end do
1144
1145	! Work out dimension lengths
1146	if (ndims >= 3) then
1147	ndimlen1 = ndimlens(1)
1148	ndimlen2 = ndimlens(2)
1149	ndimlen3 = ntotal/(ndimlen1*ndimlen2)
1150	else if (ndims >= 2) then
1151	ndimlen1 = ndimlens(1)
1152	ndimlen2 = ntotal/ndimlen1
1153	ndimlen3 = 1
1154	else
1155	ndimlen1 = ntotal
1156	ndimlen2 = 1
1157	ndimlen3 = 1
1158	end if
1159
1160	! Deallocate if necessary
1161	if (allocated(var)) then
1162	deallocate(var)
1163	end if
1164
1165	if (do_permute) then
1166	! Read and permute
1167	allocate(var_permute(ndimlen1, ndimlen2, ndimlen3))
1168	n_dimlens_permuted(i_permute_3d) = ndimlens(1:3)
1169
1170	! Reallocate if necessary
1171	if (allocated(var)) then
1172	if (size(var,1) /= n_dimlens_permuted(1) &
1173	& .or. size(var,2) /= n_dimlens_permuted(2) &
1174	& .or. size(var,3) /= n_dimlens_permuted(3)) then
1175	if (this%iverbose >= 1) then
1176	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1177	end if
1178	deallocate(var)
1179	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1180	& n_dimlens_permuted(3)))
1181	end if
1182	else
1183	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1184	& n_dimlens_permuted(3)))
1185	end if
1186
1187	if (this%iverbose >= 3) then
1188	write(nulout,'(a,a,a,i0,i0,i0,a)',advance='no') ' Reading ', var_name, &
1189	& ' (permuted dimensions ', i_permute_3d, ')'
1190	call this%print_variable_attributes(ivarid,nulout)
1191	end if
1192
1193	istatus = nf90_get_var(this%ncid, ivarid, var_permute)
1194	var = reshape(var_permute, n_dimlens_permuted, order=i_permute_3d)
1195	deallocate(var_permute)
1196
1197	else
1198	! Read data without permutation
1199
1200	! Reallocate if necessary
1201	if (allocated(var)) then
1202	if (size(var,1) /= ndimlen1 &
1203	& .or. size(var,2) /= ndimlen2 &
1204	& .or. size(var,3) /= ndimlen3) then
1205	if (this%iverbose >= 1) then
1206	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1207	end if
1208	deallocate(var)
1209	allocate(var(ndimlen1, ndimlen2, ndimlen3))
1210	end if
1211	else
1212	allocate(var(ndimlen1, ndimlen2, ndimlen3))
1213	end if
1214
1215	if (this%iverbose >= 3) then
1216	write(nulout,'(a,a,a,i0,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, &
1217	& '(', ndimlen1, ',', ndimlen2, ',', ndimlen3, ')'
1218	call this%print_variable_attributes(ivarid,nulout)
1219	end if
1220
1221	istatus = nf90_get_var(this%ncid, ivarid, var)
1222	end if
1223
1224	if (istatus /= NF90_NOERR) then
1225	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
1226	& var_name, ' as a 3D array: ', trim(nf90_strerror(istatus))
1227	call my_abort('Error reading NetCDF file')
1228	end if
1229
1230	end subroutine get_real_array3
1231
1232
1233	!---------------------------------------------------------------------
1234	! Read 3D array of data from a larger-dimensional array, which must
1235	! be allocatable and will be reallocated if necessary. Whether to
1236	! pemute is specifed by the final optional argument
1237	subroutine get_real_array3_indexed(this, var_name, index, var, ipermute)
1238	class(netcdf_file) :: this
1239	character(len=*), intent(in) :: var_name
1240	integer, intent(in) :: index
1241	real(jprb), allocatable, intent(out) :: var(:,:,:)
1242	integer, optional, intent(in) :: ipermute(3)
1243
1244	real(jprb), allocatable :: var_permute(:,:,:)
1245	integer :: ndimlen1, ndimlen2, ndimlen3
1246	integer :: istatus
1247	integer :: ivarid, ndims
1248	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1249	integer :: vstart(NF90_MAX_VAR_DIMS)
1250	integer :: vcount(NF90_MAX_VAR_DIMS)
1251	integer :: j, ntotal
1252	integer :: n_dimlens_permuted(3)
1253	integer :: i_permute_3d(3)
1254	logical :: do_permute
1255
1256	! Decide whether to permute
1257	if (present(ipermute)) then
1258	do_permute = .true.
1259	i_permute_3d = ipermute
1260	else
1261	do_permute = this%do_permute_3d
1262	i_permute_3d = this%i_permute_3d
1263	end if
1264
1265	call this%get_variable_id(var_name, ivarid)
1266	call this%get_array_dimensions(ivarid, ndims, ndimlens)
1267
1268	! Ensure the variable has no more than three non-singleton
1269	! dimensions aside from the last one
1270	ntotal = 1
1271	do j = 1, ndims-1
1272	ntotal = ntotal * ndimlens(j)
1273	if (j > 3 .and. ndimlens(j) > 1) then
1274	write(nulerr,'(a,a,a)') '*** Error reading 3D slice from NetCDF variable ', &
1275	& var_name, &
1276	& ': all dimensions above the third must be singletons'
1277	call my_abort('Error reading NetCDF file')
1278	end if
1279	end do
1280
1281	if (index < 1 .or. index > ndimlens(ndims)) then
1282	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error reading element ', index, &
1283	& ' of NetCDF variable ', &
1284	& var_name, ' with outer dimension ', ndimlens(ndims)
1285	call my_abort('Error reading NetCDF file')
1286	end if
1287
1288	! Work out dimension lengths
1289	if (ndims >= 3) then
1290	ndimlen1 = ndimlens(1)
1291	ndimlen2 = ndimlens(2)
1292	ndimlen3 = ntotal/(ndimlen1*ndimlen2)
1293	else if (ndims >= 2) then
1294	ndimlen1 = ndimlens(1)
1295	ndimlen2 = ntotal/ndimlen1
1296	ndimlen3 = 1
1297	else
1298	ndimlen1 = ntotal
1299	ndimlen2 = 1
1300	ndimlen3 = 1
1301	end if
1302
1303	vstart(1:ndims-1) = 1
1304	vstart(ndims) = index
1305	vcount(1:ndims-1) = ndimlens(1:ndims-1)
1306	vcount(ndims) = 1
1307
1308	if (do_permute) then
1309	! Read and permute
1310	allocate(var_permute(ndimlen1, ndimlen2, ndimlen3))
1311	n_dimlens_permuted(i_permute_3d) = ndimlens(1:3)
1312
1313	! Reallocate if necessary
1314	if (allocated(var)) then
1315	if (size(var,1) /= n_dimlens_permuted(1) &
1316	& .or. size(var,2) /= n_dimlens_permuted(2) &
1317	& .or. size(var,3) /= n_dimlens_permuted(3)) then
1318	if (this%iverbose >= 1) then
1319	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1320	end if
1321	deallocate(var)
1322	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1323	& n_dimlens_permuted(3)))
1324	end if
1325	else
1326	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1327	& n_dimlens_permuted(3)))
1328	end if
1329
1330	if (this%iverbose >= 3) then
1331	write(nulout,'(a,i0,a,a,a,i0,i0,i0,a)') ' Reading slice ', index, &
1332	& ' of ', var_name, &
1333	& ' (permuted dimensions ', i_permute_3d, ')'
1334	end if
1335
1336	istatus = nf90_get_var(this%ncid, ivarid, var_permute, &
1337	& start=vstart, count=vcount)
1338	var = reshape(var_permute, n_dimlens_permuted, order=i_permute_3d)
1339	deallocate(var_permute)
1340
1341	else
1342	! Read data without permutation
1343
1344	! Reallocate if necessary
1345	if (allocated(var)) then
1346	if (size(var,1) /= ndimlen1 &
1347	& .or. size(var,2) /= ndimlen2 &
1348	& .or. size(var,3) /= ndimlen3) then
1349	if (this%iverbose >= 1) then
1350	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1351	end if
1352	deallocate(var)
1353	allocate(var(ndimlen1, ndimlen2, ndimlen3))
1354	end if
1355	else
1356	allocate(var(ndimlen1, ndimlen2, ndimlen3))
1357	end if
1358
1359	if (this%iverbose >= 3) then
1360	write(nulout,'(a,i0,a,a,a,i0,a,i0,a,i0,a)') ' Reading slice ', index, &
1361	& ' of ', var_name, ' as ', ndimlen1, 'x', ndimlen2, 'x', &
1362	& ndimlen3, 'array'
1363	end if
1364
1365	istatus = nf90_get_var(this%ncid, ivarid, var, &
1366	& start=vstart, count=vcount)
1367	end if
1368
1369	if (istatus /= NF90_NOERR) then
1370	write(nulerr,'(a,a,a,a)') '*** Error reading 3D slice of NetCDF variable ', &
1371	& var_name, ': ', trim(nf90_strerror(istatus))
1372	call my_abort('Error reading NetCDF file')
1373	end if
1374
1375	end subroutine get_real_array3_indexed
1376
1377
1378	!---------------------------------------------------------------------
1379	! Read 4D array into "var", which must be allocatable and will be
1380	! reallocated if necessary. Whether to pemute is specifed by the
1381	! final optional argument
1382	subroutine get_real_array4(this, var_name, var, ipermute)
1383	class(netcdf_file) :: this
1384	character(len=*), intent(in) :: var_name
1385	real(jprb), allocatable, intent(out) :: var(:,:,:,:)
1386	integer, optional, intent(in) :: ipermute(4)
1387
1388	real(jprb), allocatable :: var_permute(:,:,:,:)
1389	integer :: ndimlen1, ndimlen2, ndimlen3, ndimlen4
1390	integer :: istatus
1391	integer :: ivarid, ndims
1392	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1393	integer :: j, ntotal
1394	integer :: n_dimlens_permuted(4)
1395	integer :: i_permute_4d(4)
1396	logical :: do_permute
1397
1398	! Decide whether to permute
1399	if (present(ipermute)) then
1400	do_permute = .true.
1401	i_permute_4d = ipermute
1402	else
1403	do_permute = this%do_permute_4d
1404	i_permute_4d = this%i_permute_4d
1405	end if
1406
1407	call this%get_variable_id(var_name, ivarid)
1408	call this%get_array_dimensions(ivarid, ndims, ndimlens)
1409
1410	! Ensure the variable has no more than three non-singleton
1411	! dimensions
1412	ntotal = 1
1413	do j = 1, ndims
1414	ntotal = ntotal * ndimlens(j)
1415	if (j > 4 .and. ndimlens(j) > 1) then
1416	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
1417	& var_name, &
1418	& ' as a 4D array: all dimensions above the third must be singletons'
1419	call my_abort('Error reading NetCDF file')
1420	end if
1421	end do
1422
1423	! Work out dimension lengths
1424	if (ndims >= 4) then
1425	ndimlen1 = ndimlens(1)
1426	ndimlen2 = ndimlens(2)
1427	ndimlen3 = ndimlens(3)
1428	ndimlen4 = ntotal/(ndimlen1ndimlen2ndimlen3)
1429	else if (ndims >= 3) then
1430	ndimlen1 = ndimlens(1)
1431	ndimlen2 = ndimlens(2)
1432	ndimlen3 = ntotal/(ndimlen1*ndimlen2)
1433	else if (ndims >= 2) then
1434	ndimlen1 = ndimlens(1)
1435	ndimlen2 = ntotal/ndimlen1
1436	ndimlen3 = 1
1437	else
1438	ndimlen1 = ntotal
1439	ndimlen2 = 1
1440	ndimlen3 = 1
1441	end if
1442
1443	! Deallocate if necessary
1444	if (allocated(var)) then
1445	deallocate(var)
1446	end if
1447
1448	if (do_permute) then
1449	! Read and permute - not tested
1450	allocate(var_permute(ndimlen1, ndimlen2, ndimlen3, ndimlen4))
1451	n_dimlens_permuted(i_permute_4d) = ndimlens(1:4)
1452
1453	! Reallocate if necessary
1454	if (allocated(var)) then
1455	if (size(var,1) /= n_dimlens_permuted(1) &
1456	& .or. size(var,2) /= n_dimlens_permuted(2) &
1457	& .or. size(var,3) /= n_dimlens_permuted(3) &
1458	& .or. size(var,4) /= n_dimlens_permuted(4)) then
1459	if (this%iverbose >= 1) then
1460	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1461	end if
1462	deallocate(var)
1463	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1464	& n_dimlens_permuted(3), n_dimlens_permuted(4)))
1465	end if
1466	else
1467	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1468	& n_dimlens_permuted(3), n_dimlens_permuted(4)))
1469	end if
1470
1471	if (this%iverbose >= 3) then
1472	write(nulout,'(a,a,a,i0,i0,i0,a)',advance='no') ' Reading ', var_name, &
1473	& ' (permuted dimensions ', i_permute_4d, ')'
1474	call this%print_variable_attributes(ivarid,nulout)
1475	end if
1476
1477	istatus = nf90_get_var(this%ncid, ivarid, var_permute)
1478	var = reshape(var_permute, n_dimlens_permuted, order=i_permute_4d)
1479	deallocate(var_permute)
1480
1481	else
1482	! Read data without permutation
1483
1484	! Reallocate if necessary
1485	if (allocated(var)) then
1486	if (size(var,1) /= ndimlen1 &
1487	& .or. size(var,2) /= ndimlen2 &
1488	& .or. size(var,3) /= ndimlen3 &
1489	& .or. size(var,4) /= ndimlen4) then
1490	if (this%iverbose >= 1) then
1491	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1492	end if
1493	deallocate(var)
1494	allocate(var(ndimlen1, ndimlen2, ndimlen3, ndimlen4))
1495	end if
1496	else
1497	allocate(var(ndimlen1, ndimlen2, ndimlen3, ndimlen4))
1498	end if
1499
1500	if (this%iverbose >= 3) then
1501	write(nulout,'(a,a,a,i0,a,i0,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, &
1502	& '(', ndimlen1, ',', ndimlen2, ',', ndimlen3,',', ndimlen4, ')'
1503	call this%print_variable_attributes(ivarid,nulout)
1504	end if
1505
1506	istatus = nf90_get_var(this%ncid, ivarid, var)
1507	end if
1508
1509	if (istatus /= NF90_NOERR) then
1510	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
1511	& var_name, ' as a 4D array: ', trim(nf90_strerror(istatus))
1512	call my_abort('Error reading NetCDF file')
1513	end if
1514
1515	end subroutine get_real_array4
1516
1517
1518	!---------------------------------------------------------------------
1519	! Get attribute as a character string
1520	subroutine get_string_attribute(this, var_name, attr_name, attr_str)
1521	class(netcdf_file) :: this
1522
1523	character(len=*), intent(in) :: var_name, attr_name
1524	character(len=*), intent(inout) :: attr_str
1525
1526	integer :: i_attr_len, ivarid
1527	integer :: istatus
1528	integer :: j
1529
1530	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
1531	if (istatus /= NF90_NOERR) then
1532	write(nulerr,'(a,a,a,a)') '*** Error inquiring about variable ', var_name, &
1533	& ': ', trim(nf90_strerror(istatus))
1534	call my_abort('Error reading NetCDF file')
1535	end if
1536	istatus = nf90_inquire_attribute(this%ncid, ivarid, attr_name, &
1537	& len = i_attr_len)
1538	if (istatus /= NF90_NOERR) then
1539	write(nulerr,'(a,a,a,a)') '*** Error reading size of attribute ', attr_name, &
1540	& ': ', trim(nf90_strerror(istatus))
1541	call my_abort('Error reading NetCDF file')
1542	end if
1543
1544	! Allocatable character strings not supported on enough compilers
1545	! yet
1546	! if (allocated(attr_str)) then
1547	! deallocate(attr_str)
1548	! end if
1549	! allocate(character(len=i_attr_len) :: attr_str)
1550	if (len(attr_str) < i_attr_len) then
1551	write(nulerr,'(a,a)') '*** Error: not enough space to read attribute ', attr_name
1552	call my_abort('Error reading NetCDF file')
1553	end if
1554
1555	istatus = nf90_get_att(this%ncid, ivarid, attr_name, attr_str)
1556	if (istatus /= NF90_NOERR) then
1557	write(nulerr,'(a,a,a,a)') '*** Error reading attribute ', attr_name, &
1558	& ': ', trim(nf90_strerror(istatus))
1559	call my_abort('Error reading NetCDF file')
1560	end if
1561
1562	! Pad with blanks since nf90_get_att does not do this
1563	do j = i_attr_len+1,len(attr_str)
1564	attr_str(j:j) = ' '
1565	end do
1566
1567	end subroutine get_string_attribute
1568
1569
1570
1571	!---------------------------------------------------------------------
1572	! Get attribute as a real scalar
1573	subroutine get_real_scalar_attribute(this, var_name, attr_name, attr)
1574	class(netcdf_file) :: this
1575
1576	character(len=*), intent(in) :: var_name, attr_name
1577	real(jprb), intent(out) :: attr
1578
1579	integer :: i_attr_len, ivarid
1580	integer :: istatus
1581	integer :: j
1582
1583	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
1584	if (istatus /= NF90_NOERR) then
1585	write(nulerr,'(a,a,a,a)') '*** Error inquiring about variable ', var_name, &
1586	& ': ', trim(nf90_strerror(istatus))
1587	call my_abort('Error reading NetCDF file')
1588	end if
1589	istatus = nf90_get_att(this%ncid, ivarid, attr_name, attr)
1590	if (istatus /= NF90_NOERR) then
1591	write(nulerr,'(a,a,a,a)') '*** Error reading attribute ', attr_name, &
1592	& ': ', trim(nf90_strerror(istatus))
1593	call my_abort('Error reading NetCDF file')
1594	end if
1595
1596	end subroutine get_real_scalar_attribute
1597
1598
1599	!---------------------------------------------------------------------
1600	! Get a global attribute as a character string
1601	subroutine get_global_attribute(this, attr_name, attr_str)
1602	class(netcdf_file) :: this
1603
1604	character(len=*), intent(in) :: attr_name
1605	character(len=*), intent(inout) :: attr_str
1606
1607	integer :: i_attr_len
1608	integer :: istatus
1609	integer :: j
1610
1611	istatus = nf90_inquire_attribute(this%ncid, NF90_GLOBAL, attr_name, &
1612	& len = i_attr_len)
1613	if (istatus /= NF90_NOERR) then
1614	write(nulerr,'(a,a,a,a)') '*** Error reading size of global attribute ', attr_name, &
1615	& ': ', trim(nf90_strerror(istatus))
1616	call my_abort('Error reading NetCDF file')
1617	end if
1618
1619	! Allocatable character strings not supported one enough compilers
1620	! yet
1621	! if (allocated(attr_str)) then
1622	! deallocate(attr_str)
1623	! end if
1624	! allocate(character(len=i_attr_len) :: attr_str)
1625	if (len(attr_str) < i_attr_len) then
1626	write(nulerr,'(a,a)') '*** Error: not enough space to read global attribute ', attr_name
1627	call my_abort('Error reading NetCDF file')
1628	end if
1629
1630	istatus = nf90_get_att(this%ncid, NF90_GLOBAL, attr_name, attr_str)
1631	if (istatus /= NF90_NOERR) then
1632	write(nulerr,'(a,a,a,a)') '*** Error reading global attribute ', attr_name, &
1633	& ': ', trim(nf90_strerror(istatus))
1634	call my_abort('Error reading NetCDF file')
1635	end if
1636
1637	! Pad with blanks since nf90_get_att does not do this
1638	do j = i_attr_len+1,len(attr_str)
1639	attr_str(j:j) = ' '
1640	end do
1641
1642	end subroutine get_global_attribute
1643
1644
1645	!---------------------------------------------------------------------
1646	! Print a variable's long_name, units and comment, according to
1647	! verbosity level
1648	subroutine print_variable_attributes(this, ivarid, iunit)
1649	class(netcdf_file) :: this
1650	integer, intent(in) :: ivarid ! NetCDF ID of variable
1651	integer, intent(in) :: iunit ! Unit to print information to
1652
1653	character(len=4000) :: attr_str
1654	integer :: i_attr_len
1655	integer :: istatus
1656	integer :: j
1657
1658	if (this%iverbose >= 4) then
1659	istatus = nf90_get_att(this%ncid, ivarid, 'long_name', attr_str)
1660	if (istatus == NF90_NOERR) then
1661	write(iunit, '(a)') ':'
1662	write(iunit, '(a,a)', advance='no') ' ', trim(attr_str)
1663	istatus = nf90_get_att(this%ncid, ivarid, 'units', attr_str)
1664	if (istatus == NF90_NOERR) then
1665	if (trim(attr_str) == '1') then
1666	write(iunit, '(a)') ' (dimensionless)'
1667	else
1668	write(iunit, '(a,a,a)') ' (', trim(attr_str), ')'
1669	end if
1670	else
1671	! No units present
1672	write(iunit, '(1x)')
1673	end if
1674	if (this%iverbose >= 5) then
1675	istatus = nf90_get_att(this%ncid, ivarid, 'comment', attr_str)
1676	if (istatus == NF90_NOERR) then
1677	write(iunit, '(a,a,a)') 'comment="', trim(attr_str), '"'
1678	end if
1679	end if
1680	end if
1681	else
1682	! No long_name present
1683	write(iunit, '(1x)')
1684	end if
1685
1686	end subroutine print_variable_attributes
1687
1688
1689	! --- WRITING SUBROUTINES ---
1690
1691	!---------------------------------------------------------------------
1692	! Define a dimension with name dim_name of length n (or 0 to
1693	! indicate the unlimited dimension)
1694	subroutine define_dimension(this, dim_name, n)
1695	class(netcdf_file) :: this
1696	character(len=*), intent(in) :: dim_name
1697	integer, intent(in) :: n
1698	integer :: idimid, istatus
1699
1700	istatus = nf90_def_dim(this%ncid, dim_name, n, idimid)
1701	if (istatus /= NF90_NOERR) then
1702	write(nulerr,'(a,a,a,a)') '*** Error defining ', dim_name, &
1703	& ' as a dimension: ', trim(nf90_strerror(istatus))
1704	call my_abort('Error writing NetCDF file')
1705	end if
1706
1707	if (this%iverbose >= 4) then
1708	write(nulout,'(a,a,a,i0)') ' Defining dimension ',trim(dim_name), &
1709	& ' of length ', n
1710	end if
1711
1712	end subroutine define_dimension
1713
1714
1715	!---------------------------------------------------------------------
1716	! Define a variable with name var_name, and specify its shape via
1717	! the dim1_name, dim2_name and dim3_name optional arguments, which
1718	! are strings referring to existing dimensions. If none are
1719	! specified, the variable will be a scalar. The optional arguments
1720	! long_name, units and comment write string attributes with these
1721	! names.
1722	subroutine define_variable(this, var_name, dim1_name, dim2_name, dim3_name, &
1723	& long_name, units_str, comment_str, standard_name, is_double, &
1724	& data_type_name, fill_value, deflate_level, shuffle, chunksizes)
1725	class(netcdf_file) :: this
1726	character(len=*), intent(in) :: var_name
1727	character(len=*), intent(in), optional :: long_name, units_str, comment_str, standard_name
1728	character(len=*), intent(in), optional :: dim1_name, dim2_name, dim3_name
1729	logical, intent(in), optional :: is_double
1730	character(len=*), intent(in), optional :: data_type_name
1731	real(jprb), intent(in), optional :: fill_value
1732	integer, intent(in), optional :: deflate_level ! Compression: 0 (none) to 9 (most)
1733	logical, intent(in), optional :: shuffle ! Shuffle bytes before compression
1734	integer, dimension(:), intent(in), optional :: chunksizes
1735
1736	integer :: istatus, ndims, ivarid
1737	integer, dimension(NF90_MAX_VAR_DIMS) :: idimids
1738	integer :: data_type
1739
1740	if (present(dim1_name)) then
1741	! Variable is at least one dimensional
1742	ndims = 1
1743	istatus = nf90_inq_dimid(this%ncid, dim1_name, idimids(1))
1744	if (istatus /= NF90_NOERR) then
1745	write(nulerr,'(a,a,a,a)') '*** Error inquiring ID of dimension ', &
1746	& dim1_name, ': ', trim(nf90_strerror(istatus))
1747	call my_abort('Error writing NetCDF file')
1748	end if
1749	if (present(dim2_name)) then
1750	! Variable is at least two dimensional
1751	ndims = 2
1752	istatus = nf90_inq_dimid(this%ncid, dim2_name, idimids(2))
1753	if (istatus /= NF90_NOERR) then
1754	write(nulerr,'(a,a,a)') '*** Error inquiring ID of dimension ', &
1755	& dim2_name, ': ', trim(nf90_strerror(istatus))
1756	call my_abort('Error writing NetCDF file')
1757	end if
1758	if (present(dim3_name)) then
1759	! Variable is at least three dimensional
1760	ndims = 3
1761	istatus = nf90_inq_dimid(this%ncid, dim3_name, idimids(3))
1762	if (istatus /= NF90_NOERR) then
1763	write(nulerr,'(a,a,a,a)') '*** Error inquiring ID of dimension ', &
1764	& dim3_name, ': ', trim(nf90_strerror(istatus))
1765	call my_abort('Error writing NetCDF file')
1766	end if
1767	end if
1768	end if
1769	else
1770	! Variable is a scalar
1771	ndims = 0
1772	end if
1773
1774	! Read output precision from optional argument "is_double" if
1775	! present, otherwise from default output precision for this file
1776	data_type = NF90_FLOAT ! Default
1777	if (present(data_type_name)) then
1778	if (data_type_name == 'double') then
1779	data_type = NF90_DOUBLE
1780	else if (data_type_name == 'byte') then
1781	data_type = NF90_BYTE
1782	else if (data_type_name == 'short') then
1783	data_type = NF90_SHORT
1784	else if (data_type_name == 'int') then
1785	data_type = NF90_INT
1786	else if (data_type_name == 'float') then
1787	data_type = NF90_FLOAT
1788	else
1789	write(nulerr,'(a,a,a)') '*** Error: netCDF data type "', data_type_name, '" not supported'
1790	call my_abort('Error writing NetCDF file')
1791	end if
1792	else if (present(is_double)) then
1793	data_type = NF90_DOUBLE
1794	end if
1795
1796	! Define variable
1797	#ifdef NC_NETCDF4
1798	istatus = nf90_def_var(this%ncid, var_name, data_type, idimids(1:ndims), &
1799	& ivarid, deflate_level=deflate_level, shuffle=shuffle, chunksizes=chunksizes)
1800	#else
1801	istatus = nf90_def_var(this%ncid, var_name, data_type, idimids(1:ndims), ivarid)
1802	#endif
1803	if (istatus /= NF90_NOERR) then
1804	write(nulerr,'(a,a,a,a)') '*** Error defining variable ', var_name, &
1805	& ': ', trim(nf90_strerror(istatus))
1806	call my_abort('Error writing NetCDF file')
1807	end if
1808
1809	! Add optional attributes
1810	if (present(long_name)) then
1811	istatus = nf90_put_att(this%ncid, ivarid, "long_name", long_name)
1812	if (this%iverbose >= 4) then
1813	write(nulout,'(a,a,a,a)') ' Defining ',trim(var_name),': ',long_name
1814	end if
1815	else
1816	if (this%iverbose >= 4) then
1817	write(nulout,'(a,a)') ' Defining ',trim(var_name)
1818	end if
1819	end if
1820	if (present(units_str)) then
1821	istatus = nf90_put_att(this%ncid, ivarid, "units", units_str)
1822	end if
1823	if (present(standard_name)) then
1824	istatus = nf90_put_att(this%ncid, ivarid, "standard_name", standard_name)
1825	end if
1826	if (present(comment_str)) then
1827	istatus = nf90_put_att(this%ncid, ivarid, "comment", comment_str)
1828	end if
1829
1830	if (present(fill_value)) then
1831	#ifdef NC_NETCDF4
1832	if (data_type == NF90_DOUBLE) then
1833	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, real(fill_value,jprd))
1834	else if (data_type == NF90_FLOAT) then
1835	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, real(fill_value,jprm))
1836	else if (data_type == NF90_INT) then
1837	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, int(fill_value,4))
1838	else if (data_type == NF90_SHORT) then
1839	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, int(fill_value,2))
1840	else if (data_type == NF90_BYTE) then
1841	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, int(fill_value,1))
1842	end if
1843	#else
1844	if (data_type == NF90_DOUBLE) then
1845	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", real(fill_value,jprd))
1846	else if (data_type == NF90_FLOAT) then
1847	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", real(fill_value,jprm))
1848	else if (data_type == NF90_INT) then
1849	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", int(fill_value,4))
1850	else if (data_type == NF90_SHORT) then
1851	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", int(fill_value,2))
1852	else if (data_type == NF90_BYTE) then
1853	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", int(fill_value,1))
1854	end if
1855	#endif
1856	end if
1857
1858	end subroutine define_variable
1859
1860
1861	!---------------------------------------------------------------------
1862	! Put CF-compliant global attributes into the file
1863	subroutine put_global_attributes(this, title_str, inst_str, source_str, &
1864	& comment_str, references_str, creator_name, creator_email_str, &
1865	& contributor_name, project_str, conventions_str)
1866	class(netcdf_file) :: this
1867
1868	character(len=*), intent(in), optional :: title_str
1869	character(len=*), intent(in), optional :: inst_str
1870	character(len=*), intent(in), optional :: source_str
1871	character(len=*), intent(in), optional :: creator_name, creator_email_str
1872	character(len=*), intent(in), optional :: contributor_name, project_str
1873	character(len=*), intent(in), optional :: comment_str, conventions_str
1874	character(len=*), intent(in), optional :: references_str
1875
1876	character(len=32) :: date_time_str
1877	character(len=4000) :: command_line_str
1878	character(len=4032) :: history_str
1879
1880	integer :: time_vals(8)
1881	integer :: i ! status
1882
1883	call date_and_time(values=time_vals)
1884	call get_command(command_line_str)
1885
1886	write(date_time_str,"(i0.4,'-',i0.2,'-',i0.2,' ',i0.2,':',i0.2,':',i0.2)") &
1887	& time_vals(1), time_vals(2), time_vals(3), time_vals(5), time_vals(6), time_vals(7)
1888
1889	history_str = trim(date_time_str) // ': ' // trim(command_line_str)
1890
1891	if (present(title_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "title", title_str)
1892	if (present(inst_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "institution", inst_str)
1893	if (present(source_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "source", source_str)
1894	if (present(creator_name))i=nf90_put_att(this%ncid, NF90_GLOBAL, "creator_name", creator_name)
1895	if (present(creator_email_str))i=nf90_put_att(this%ncid, NF90_GLOBAL, "creator_email", creator_email_str)
1896	if (present(contributor_name))i=nf90_put_att(this%ncid, NF90_GLOBAL, "contributor_name", contributor_name)
1897
1898	i = nf90_put_att(this%ncid, NF90_GLOBAL, "history", history_str)
1899
1900	if (present(project_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "project", project_str)
1901	if (present(comment_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "comment", comment_str)
1902	if (present(references_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, &
1903	& "references", references_str)
1904	if (present(conventions_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, &
1905	& "conventions", conventions_str)
1906
1907	end subroutine put_global_attributes
1908
1909
1910	!---------------------------------------------------------------------
1911	! Put a non-standard global attribute into the file
1912	subroutine put_global_attribute(this, attr_name, attr_str)
1913	class(netcdf_file) :: this
1914
1915	character(len=*), intent(in) :: attr_name, attr_str
1916
1917	integer :: istatus
1918
1919	istatus = nf90_put_att(this%ncid, NF90_GLOBAL, trim(attr_name), trim(attr_str))
1920
1921	if (istatus /= NF90_NOERR) then
1922	write(nulerr,'(a,a,a,a)') '*** Error writing global attribute ', attr_name, &
1923	& ': ', trim(nf90_strerror(istatus))
1924	call my_abort('Error writing NetCDF file')
1925	end if
1926
1927	end subroutine put_global_attribute
1928
1929
1930	!---------------------------------------------------------------------
1931	! Put a non-standard variable attribute into the file
1932	subroutine put_attribute(this, var_name, attr_name, attr_str)
1933	class(netcdf_file) :: this
1934
1935	character(len=*), intent(in) :: var_name, attr_name, attr_str
1936
1937	integer :: istatus, ivarid
1938
1939	call this%get_variable_id(var_name, ivarid)
1940
1941	istatus = nf90_put_att(this%ncid, ivarid, trim(attr_name), trim(attr_str))
1942
1943	if (istatus /= NF90_NOERR) then
1944	write(nulerr,'(a,a,a,a)') '*** Error writing attribute ', attr_name, &
1945	& ': ', trim(nf90_strerror(istatus))
1946	call my_abort('Error writing NetCDF file')
1947	end if
1948
1949	end subroutine put_attribute
1950
1951
1952	!---------------------------------------------------------------------
1953	! The "put" method saves a scalar, vector or matrix into the
1954	! variable with name var_name, according to the rank of the var
1955	! argument. This version saves a scalar.
1956	subroutine put_real_scalar(this, var_name, var)
1957	class(netcdf_file) :: this
1958	character(len=*), intent(in) :: var_name
1959	real(jprb), intent(in) :: var
1960
1961	integer :: ivarid, ndims, istatus
1962	integer(kind=jpib) :: ntotal
1963	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1964
1965	! If we are in define mode, exit define mode
1966	call this%end_define_mode()
1967
1968	! Check the variable is a scalar
1969	call this%get_variable_id(var_name, ivarid)
1970	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
1971	if (ntotal /= 1) then
1972	write(nulerr,'(a,a,a,i0)') '*** Error: attempt to write scalar to ', &
1973	& var_name, ' which has total length ', ntotal
1974	call my_abort('Error writing NetCDF file')
1975	end if
1976
1977	! Save the scalar
1978	istatus = nf90_put_var(this%ncid, ivarid, var)
1979	if (istatus /= NF90_NOERR) then
1980	write(nulerr,'(a,a,a,a)') '*** Error writing scalar ', var_name, ': ', &
1981	& trim(nf90_strerror(istatus))
1982	call my_abort('Error writing NetCDF file')
1983	end if
1984
1985	end subroutine put_real_scalar
1986
1987
1988	!---------------------------------------------------------------------
1989	! Save a scalar.
1990	subroutine put_real_scalar_indexed(this, var_name, index, var)
1991	class(netcdf_file) :: this
1992	character(len=*), intent(in) :: var_name
1993	real(jprb), intent(in) :: var
1994	integer, intent(in) :: index
1995
1996	integer :: ivarid, ndims, istatus
1997	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1998	integer :: vstart(NF90_MAX_VAR_DIMS)
1999
2000	! If we are in define mode, exit define mode
2001	call this%end_define_mode()
2002
2003	! Check the variable is a scalar
2004	call this%get_variable_id(var_name, ivarid)
2005	call this%get_array_dimensions(ivarid, ndims, ndimlens)
2006	if (index < 1 .or. index > ndimlens(ndims)) then
2007	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write scalar to element ', &
2008	& index, ' of ', var_name, ' which has outer dimension ', ndimlens(ndims)
2009	call my_abort('Error writing NetCDF file')
2010	end if
2011
2012	! Save the scalar
2013	vstart(1:ndims-1) = 1
2014	vstart(ndims) = index
2015	istatus = nf90_put_var(this%ncid, ivarid, var, start=vstart)
2016	if (istatus /= NF90_NOERR) then
2017	write(nulerr,'(a,a,a,a)') '*** Error writing scalar to ', var_name, ': ', &
2018	& trim(nf90_strerror(istatus))
2019	call my_abort('Error writing NetCDF file')
2020	end if
2021
2022	end subroutine put_real_scalar_indexed
2023
2024
2025	!---------------------------------------------------------------------
2026	! Save a vector with name var_name in the file
2027	subroutine put_real_vector(this, var_name, var)
2028	class(netcdf_file) :: this
2029	character(len=*), intent(in) :: var_name
2030	real(jprb), intent(in) :: var(:)
2031
2032	integer :: ivarid, ndims, istatus
2033	integer(kind=jpib) :: ntotal
2034	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2035
2036	call this%end_define_mode()
2037
2038	! Check the vector is of the right length
2039	call this%get_variable_id(var_name, ivarid)
2040	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2041	if (ntotal /= size(var,kind=jpib)) then
2042	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write vector of length ', &
2043	& size(var), ' to ', var_name, ' which has total length ', ntotal
2044	call my_abort('Error writing NetCDF file')
2045	end if
2046
2047	! Save the vector
2048	istatus = nf90_put_var(this%ncid, ivarid, var)
2049	if (istatus /= NF90_NOERR) then
2050	write(nulerr,'(a,a,a,a)') '*** Error writing vector ', var_name, ': ', &
2051	& trim(nf90_strerror(istatus))
2052	call my_abort('Error writing NetCDF file')
2053	end if
2054
2055	end subroutine put_real_vector
2056
2057
2058	!---------------------------------------------------------------------
2059	! Save a vector slice with name var_name in the file
2060	subroutine put_real_vector_indexed(this, var_name, index, var)
2061	class(netcdf_file) :: this
2062	character(len=*), intent(in) :: var_name
2063	real(jprb), intent(in) :: var(:)
2064	integer, intent(in) :: index
2065
2066	integer :: ivarid, ndims, istatus
2067	integer(kind=jpib) :: ntotal
2068	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2069	integer :: vstart(NF90_MAX_VAR_DIMS)
2070	integer :: vcount(NF90_MAX_VAR_DIMS)
2071
2072	call this%end_define_mode()
2073
2074	! Check the vector is of the right length
2075	call this%get_variable_id(var_name, ivarid)
2076	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2077	ntotal = ntotal / ndimlens(ndims)
2078	if (ntotal /= size(var,kind=jpib)) then
2079	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write vector of length ', &
2080	& size(var), ' to slice of ', var_name, ' which has length ', ntotal
2081	call my_abort('Error writing NetCDF file')
2082	end if
2083	if (index < 1 .or. index > ndimlens(ndims)) then
2084	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write vector to slice ', &
2085	& index, ' of ', var_name, ' which has outer dimension ', ndimlens(ndims)
2086	call my_abort('Error writing NetCDF file')
2087	end if
2088
2089	! Save the vector
2090	vstart(1:ndims-1) = 1
2091	vstart(ndims) = index
2092	vcount(1:ndims-1) = ndimlens(1:ndims-1)
2093	vcount(ndims) = 1
2094	istatus = nf90_put_var(this%ncid, ivarid, var, start=vstart, count=vcount)
2095	if (istatus /= NF90_NOERR) then
2096	write(nulerr,'(a,a,a,a)') '*** Error writing vector to ', var_name, ': ', &
2097	& trim(nf90_strerror(istatus))
2098	call my_abort('Error writing NetCDF file')
2099	end if
2100
2101	end subroutine put_real_vector_indexed
2102
2103
2104	!---------------------------------------------------------------------
2105	! Save a matrix with name var_name in the file, transposing its
2106	! dimensions if either optional argument transp is .true., or the
2107	! transpose_matrices method has already been called.
2108	subroutine put_real_matrix(this, var_name, var, do_transp)
2109	class(netcdf_file) :: this
2110	character(len=*), intent(in) :: var_name
2111	real(jprb), intent(in) :: var(:,:)
2112	real(jprb), allocatable :: var_transpose(:,:)
2113	logical, optional, intent(in):: do_transp
2114
2115	integer :: ivarid, ndims, nvarlen, istatus
2116	integer(kind=jpib) :: ntotal
2117	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2118
2119	logical :: do_transpose
2120
2121	if (present(do_transp)) then
2122	do_transpose = do_transp
2123	else
2124	do_transpose = this%do_transpose_2d
2125	end if
2126
2127	call this%end_define_mode()
2128
2129	call this%get_variable_id(var_name, ivarid)
2130	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2131
2132	nvarlen = size(var,1)*size(var,2)
2133
2134	! Check the total size of the variable to be stored (but receiving
2135	! ntotal is zero then there must be an unlimited dimension)
2136	if (ntotal /= size(var,kind=jpib) .and. ntotal /= 0) then
2137	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write matrix of total size ', &
2138	& nvarlen, ' to ', var_name, ' which has total size ', ntotal
2139	call my_abort('Error writing NetCDF file')
2140	end if
2141
2142	if (do_transpose) then
2143	! Save the matrix with transposition
2144	if (this%iverbose >= 3) then
2145	write(nulout,'(a,a,a)') ' Writing ', var_name, &
2146	& ' (transposing dimensions)'
2147	end if
2148	allocate(var_transpose(size(var,2), size(var,1)))
2149	var_transpose = transpose(var)
2150	istatus = nf90_put_var(this%ncid, ivarid, var_transpose)
2151	deallocate(var_transpose)
2152	else
2153	! Save the matrix without transposition
2154	if (this%iverbose >= 3) then
2155	write(nulout,'(a,a)') ' Writing ', var_name
2156	end if
2157	istatus = nf90_put_var(this%ncid, ivarid, var)
2158	end if
2159
2160	if (istatus /= NF90_NOERR) then
2161	write(nulerr,'(a,a,a,a)') '*** Error writing matrix ', var_name, &
2162	& ': ', trim(nf90_strerror(istatus))
2163	call my_abort('Error writing NetCDF file')
2164	end if
2165
2166	end subroutine put_real_matrix
2167
2168
2169	!---------------------------------------------------------------------
2170	! Save a matrix slice with name var_name in the file, transposing its
2171	! dimensions if either optional argument transp is .true., or the
2172	! transpose_matrices method has already been called.
2173	subroutine put_real_matrix_indexed(this, var_name, index, var, do_transp)
2174	class(netcdf_file) :: this
2175	character(len=*), intent(in) :: var_name
2176	real(jprb), intent(in) :: var(:,:)
2177	integer, intent(in) :: index
2178
2179	real(jprb), allocatable :: var_transpose(:,:)
2180	logical, optional, intent(in):: do_transp
2181
2182	integer :: ivarid, ndims, nvarlen, istatus
2183	integer(kind=jpib) :: ntotal
2184	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2185	integer :: vstart(NF90_MAX_VAR_DIMS)
2186	integer :: vcount(NF90_MAX_VAR_DIMS)
2187
2188	logical :: do_transpose
2189
2190	if (present(do_transp)) then
2191	do_transpose = do_transp
2192	else
2193	do_transpose = this%do_transpose_2d
2194	end if
2195
2196	call this%end_define_mode()
2197
2198	call this%get_variable_id(var_name, ivarid)
2199	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2200
2201	nvarlen = size(var,1)*size(var,2)
2202
2203	! Check the total size of the variable to be stored (but receiving
2204	! ntotal is zero then there must be an unlimited dimension)
2205	ntotal = ntotal / ndimlens(ndims)
2206	if (ntotal /= size(var,kind=jpib) .and. ntotal /= 0) then
2207	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write matrix of total size ', &
2208	& nvarlen, ' to ', var_name, ' which has total size ', ntotal
2209	call my_abort('Error writing NetCDF file')
2210	end if
2211
2212	vstart(1:ndims-1) = 1
2213	vstart(ndims) = index
2214	vcount(1:ndims-1) = ndimlens(1:ndims-1)
2215	vcount(ndims) = 1
2216
2217	if (do_transpose) then
2218	! Save the matrix with transposition
2219	if (this%iverbose >= 3) then
2220	write(nulout,'(a,i0,a,a,a)') ' Writing slice ', index, ' of ', var_name, &
2221	& ' (transposing dimensions)'
2222	end if
2223	allocate(var_transpose(size(var,2), size(var,1)))
2224	var_transpose = transpose(var)
2225	istatus = nf90_put_var(this%ncid, ivarid, var_transpose, start=vstart, count=vcount)
2226	deallocate(var_transpose)
2227	else
2228	! Save the matrix without transposition
2229	if (this%iverbose >= 3) then
2230	write(nulout,'(a,i0,a,a)') ' Writing slice ', index, ' of ', var_name
2231	end if
2232	istatus = nf90_put_var(this%ncid, ivarid, var, start=vstart, count=vcount)
2233	end if
2234
2235	if (istatus /= NF90_NOERR) then
2236	write(nulerr,'(a,a,a,a)') '*** Error writing matrix ', var_name, &
2237	& ': ', trim(nf90_strerror(istatus))
2238	call my_abort('Error writing NetCDF file')
2239	end if
2240
2241	end subroutine put_real_matrix_indexed
2242
2243
2244	!---------------------------------------------------------------------
2245	! Save a 3D array with name var_name in the file. The optional
2246	! argument permute specifies that the dimensions should first be
2247	! permuted according to the three integers therein (or if
2248	! permute_3d_arrays has already been called). ipermute is
2249	! interpretted such that if OLD and NEW are 3-element vectors
2250	! containing the size of each dimension in memory and in the written
2251	! file, respectively, then NEW=OLD(ipermute).
2252	subroutine put_real_array3(this, var_name, var, ipermute)
2253	class(netcdf_file) :: this
2254	character(len=*), intent(in) :: var_name
2255	real(jprb), intent(in) :: var(:,:,:)
2256	real(jprb), allocatable :: var_permute(:,:,:)
2257	integer, optional, intent(in) :: ipermute(3)
2258
2259	integer :: ivarid, ndims, nvarlen, istatus
2260	integer(kind=jpib) :: ntotal
2261	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2262
2263	logical :: do_permute ! Do we permute?
2264	integer :: i_permute_3d(3)
2265	integer :: n_dimlens_permuted(3)
2266	integer :: i_order(3)
2267
2268	! Decide whether to permute
2269	if (present(ipermute)) then
2270	do_permute = .true.
2271	i_permute_3d = ipermute
2272	else
2273	do_permute = this%do_permute_3d
2274	i_permute_3d = this%i_permute_3d
2275	end if
2276
2277	call this%end_define_mode()
2278
2279	! Check total size
2280	call this%get_variable_id(var_name, ivarid)
2281	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2282	nvarlen = size(var,1)size(var,2)size(var,3)
2283	if (ntotal /= size(var,kind=jpib)) then
2284	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write array of total size ', &
2285	& nvarlen, ' to ', var_name, ' which has total size ', ntotal
2286	call my_abort('Error writing NetCDF file')
2287	end if
2288
2289	if (do_permute) then
2290	! Save array after permuting dimensions
2291	if (this%iverbose >= 3) then
2292	write(nulout,'(a,a,a,i0,i0,i0,a)') ' Writing ', var_name, &
2293	& ' (permuted dimensions: ', i_permute_3d, ')'
2294	end if
2295	n_dimlens_permuted = (/ size(var,i_permute_3d(1)), &
2296	& size(var,i_permute_3d(2)), &
2297	& size(var,i_permute_3d(3)) /)
2298	if (this%iverbose >= 4) then
2299	write(nulout,'(a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a)') ' (', &
2300	& n_dimlens_permuted(1), ',', n_dimlens_permuted(2), &
2301	& ',', n_dimlens_permuted(3), ') -> (', ndimlens(1), &
2302	& ',', ndimlens(2), ',', ndimlens(3), ')'
2303	end if
2304	allocate(var_permute(n_dimlens_permuted(1), &
2305	& n_dimlens_permuted(2), n_dimlens_permuted(3)))
2306	! Due to the odd way that ORDER works in Fortran RESHAPE, we
2307	! need to do this:
2308	i_order(i_permute_3d) = (/ 1, 2, 3 /)
2309	var_permute = reshape(var, n_dimlens_permuted, order=i_order)
2310	istatus = nf90_put_var(this%ncid, ivarid, var_permute)
2311	deallocate(var_permute)
2312	else
2313	! Save array without permuting dimensions
2314	if (this%iverbose >= 3) then
2315	write(nulout,'(a,a)') ' Writing ', var_name
2316	end if
2317	istatus = nf90_put_var(this%ncid, ivarid, var)
2318	end if
2319
2320	if (istatus /= NF90_NOERR) then
2321	write(nulerr,'(a,a,a,a)') '*** Error writing array ', var_name, &
2322	& ': ', trim(nf90_strerror(istatus))
2323	call my_abort('Error writing NetCDF file')
2324	end if
2325
2326	end subroutine put_real_array3
2327
2328	end module easy_netcdf

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

Download in other formats:

Original Format