Context Navigation

easy_netcdf.F90 @ 5159

Last change on this file since 5159 was 5159, checked in by abarral, 5 months ago
Put dimensions.h and paramet.h into modules
File size: 104.1 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 = -1! 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 :: is_open
58	procedure :: get_real_scalar
59	procedure :: get_int_scalar
60	procedure :: get_real_vector
61	procedure :: get_int_vector
62	procedure :: get_real_matrix
63	procedure :: get_real_array3
64	procedure :: get_real_scalar_indexed
65	procedure :: get_real_vector_indexed
66	procedure :: get_real_matrix_indexed
67	procedure :: get_real_array3_indexed
68	procedure :: get_real_array4
69	procedure :: get_char_vector
70	procedure :: get_char_matrix
71	generic :: get => get_real_scalar, get_int_scalar, &
72	& get_real_vector, get_int_vector, &
73	& get_real_matrix, get_real_array3, &
74	& get_real_array4, &
75	& get_real_scalar_indexed, get_real_vector_indexed, &
76	& get_real_matrix_indexed, get_real_array3_indexed, &
77	& get_char_vector, get_char_matrix
78	procedure :: get_real_scalar_attribute
79	procedure :: get_string_attribute
80	generic :: get_attribute => get_real_scalar_attribute, &
81	& get_string_attribute
82	procedure :: get_global_attribute
83
84	procedure :: define_dimension
85	procedure :: define_variable
86	procedure :: put_attribute
87	procedure :: put_global_attributes
88	procedure :: put_global_attribute
89	procedure :: put_real_scalar
90	procedure :: put_real_vector
91	procedure :: put_int_vector
92	procedure :: put_real_matrix
93	procedure :: put_real_array3
94	procedure :: put_real_scalar_indexed
95	procedure :: put_real_vector_indexed
96	procedure :: put_real_matrix_indexed
97	generic :: put => put_real_scalar, put_real_vector, &
98	& put_real_matrix, put_real_array3, &
99	& put_real_scalar_indexed, put_real_vector_indexed, &
100	& put_real_matrix_indexed, put_int_vector
101	procedure :: set_verbose
102	procedure :: transpose_matrices
103	procedure :: double_precision
104	procedure :: permute_3d_arrays
105	procedure :: get_rank
106	procedure :: exists
107	procedure :: get_outer_dimension
108	procedure :: attribute_exists
109	procedure :: global_attribute_exists
110	#ifdef NC_NETCDF4
111	procedure :: copy_dimensions
112	#endif
113	procedure :: copy_variable_definition
114	procedure :: copy_variable
115	procedure, private :: get_array_dimensions
116	procedure, private :: get_variable_id
117	procedure, private :: end_define_mode
118	procedure, private :: print_variable_attributes
119	end type netcdf_file
120
121	contains
122
123	! --- GENERIC SUBROUTINES ---
124
125	!---------------------------------------------------------------------
126	! Open a NetCDF file with name "file_name", optionally specifying the
127	! verbosity level (0-5) and if the file is for writing (the default
128	! is read-only)
129	subroutine open_netcdf_file(this, file_name, iverbose, is_write_mode, is_hdf5_file)
130	class(netcdf_file) :: this
131	character(len=*), intent(in) :: file_name
132	integer, intent(in), optional :: iverbose
133	logical, intent(in), optional :: is_write_mode
134	logical, intent(in), optional :: is_hdf5_file ! Only for write mode
135
136	integer :: istatus
137	integer :: i_write_mode
138
139
140	! Store verbosity level in object
141	if (present(iverbose)) then
142	this%iverbose = iverbose
143	else
144	! By default announce files being opened and closed, but not
145	! variables read/written
146	this%iverbose = 2
147	end if
148
149	! Store read/write mode in object
150	if (present(is_write_mode)) then
151	this%is_write_mode = is_write_mode
152	else
153	this%is_write_mode = .false.
154	end if
155
156	! By default we don't transpose 2D arrays on read/write
157	this%do_transpose_2d = .false.
158
159	! Store filename
160	this%file_name = file_name
161
162	! Open file according to write mode
163	if (.not. this%is_write_mode) then
164	istatus = nf90_open(file_name, NF90_NOWRITE, this%ncid)
165	if (this%iverbose >= 2) then
166	write(nulout,'(a,a)') 'Reading NetCDF file ', file_name
167	!write(nulout,'(a,a,a,i0,a)') 'Reading NetCDF file ', file_name, ' (ID=', this%ncid, ')'
168	end if
169	this%is_define_mode = .false.
170	else
171	i_write_mode = NF90_CLOBBER
172	! Check if HDF5 file is to be written (which can be larger)
173	if (present(is_hdf5_file)) then
174	if (is_hdf5_file) then
175	#ifdef NC_NETCDF4
176	i_write_mode = ior(i_write_mode, NF90_HDF5)
177	#else
178	if (this%iverbose >= 1) then
179	write(nulout,'(a,a)') 'Warning: cannot use HDF5 format for writing ', file_name, &
180	& ' unless compiled with NC_NETCDF4 defined'
181	end if
182	#endif
183	end if
184	end if
185
186	istatus = nf90_create(file_name, i_write_mode, this%ncid)
187	if (this%iverbose >= 2) then
188	write(nulout,'(a,a)') 'Writing NetCDF file ', file_name
189	end if
190	this%is_define_mode = .true.
191	end if
192
193	! Check the file opened correctly
194	if (istatus /= NF90_NOERR) then
195	write(nulerr,'(a,a,a,a)') '*** Error opening NetCDF file ', &
196	& file_name, ': ', trim(nf90_strerror(istatus))
197	call my_abort('Error opening NetCDF file')
198	end if
199
200	end subroutine open_netcdf_file
201
202
203	!---------------------------------------------------------------------
204	! Open a NetCDF file for writing
205	subroutine create_netcdf_file(this, file_name, iverbose, is_hdf5_file)
206	class(netcdf_file) :: this
207	character(len=*), intent(in) :: file_name
208	integer, intent(in), optional :: iverbose
209	logical, intent(in), optional :: is_hdf5_file
210
211	integer :: istatus
212	integer :: i_write_mode
213
214	if (present(iverbose)) then
215	this%iverbose = iverbose
216	else
217	this%iverbose = 2
218	end if
219
220	this%do_transpose_2d = .false.
221
222	i_write_mode = NF90_CLOBBER
223	! Check if HDF5 file is to be written (which can be large)
224	if (present(is_hdf5_file)) then
225	if (is_hdf5_file) then
226	#ifdef NC_NETCDF4
227	i_write_mode = ior(i_write_mode, NF90_HDF5)
228	#else
229	if (this%iverbose >= 1) then
230	write(nulout,'(a,a)') 'Warning: cannot use HDF5 format for writing ', file_name, &
231	& ' unless compiled with NC_NETCDF4 defined'
232	end if
233	#endif
234	end if
235	end if
236
237	istatus = nf90_create(file_name, i_write_mode, this%ncid)
238	if (this%iverbose >= 2) then
239	write(nulout,'(a,a)') 'Writing NetCDF file ', file_name
240	!write(nulout,'(a,a,a,i0,a)') 'Writing NetCDF file ', file_name, ' (ID=', this%ncid, ')'
241	end if
242	this%is_define_mode = .true.
243
244	if (istatus /= NF90_NOERR) then
245	write(nulerr,'(a,a,a)') '*** Error opening NetCDF file ', file_name, &
246	& ': ', trim(nf90_strerror(istatus))
247	stop
248	end if
249	this%file_name = file_name
250
251	end subroutine create_netcdf_file
252
253
254	!---------------------------------------------------------------------
255	! Close the NetCDF file
256	subroutine close_netcdf_file(this)
257	class(netcdf_file) :: this
258	integer :: istatus
259
260	if (this%iverbose >= 3) then
261	write(nulout,'(a,a)') 'Closing NetCDF file ', trim(this%file_name)
262	end if
263
264	istatus = nf90_close(this%ncid)
265	if (istatus /= NF90_NOERR) then
266	write(nulerr,'(a,a,a,a)') '*** Error closing NetCDF file ', &
267	& trim(this%file_name), ': ', trim(nf90_strerror(istatus))
268	stop
269	end if
270
271	this%ncid = -1
272
273	end subroutine close_netcdf_file
274
275
276	!---------------------------------------------------------------------
277	! Set the verbosity level from 0 to 5, where the codes have the
278	! following meaning: 0=errors only, 1=warning, 2=info, 3=progress,
279	! 4=detailed, 5=debug
280	subroutine set_verbose(this, ival)
281	class(netcdf_file) :: this
282	integer, optional :: ival
283
284	if (present(ival)) then
285	this%iverbose = ival
286	else
287	this%iverbose = 2
288	end if
289
290	end subroutine set_verbose
291
292
293	!---------------------------------------------------------------------
294	! Specify whether floating-point arrays should be written in double precision
295	subroutine double_precision(this, is_double)
296	class(netcdf_file) :: this
297	logical, optional :: is_double
298
299	if (present(is_double)) then
300	this%is_double_precision = is_double
301	else
302	this%is_double_precision = .true.
303	end if
304
305	end subroutine double_precision
306
307
308	!---------------------------------------------------------------------
309	! Specify whether 2D arrays should be transposed on read/write
310	subroutine transpose_matrices(this, do_transpose)
311	class(netcdf_file) :: this
312	logical, optional :: do_transpose
313
314	if (present(do_transpose)) then
315	this%do_transpose_2d = do_transpose
316	else
317	this%do_transpose_2d = .true.
318	end if
319
320	end subroutine transpose_matrices
321
322
323	!---------------------------------------------------------------------
324	! Specify that 3D arrays should be permuted on write, with the new
325	! dimension order in the input argument "ipermute" (e.g. 3,2,1)
326	subroutine permute_3d_arrays(this, ipermute)
327	class(netcdf_file) :: this
328	integer, intent(in) :: ipermute(3)
329
330	this%do_permute_3d = .true.
331	this%i_permute_3d = ipermute
332
333	end subroutine permute_3d_arrays
334
335
336	!---------------------------------------------------------------------
337	! Specify that 4D arrays should be permuted on write, with the new
338	! dimension order in the input argument "ipermute" (e.g. 4,3,2,1)
339	subroutine permute_4d_arrays(this, ipermute)
340	class(netcdf_file) :: this
341	integer, intent(in) :: ipermute(4)
342
343	this%do_permute_4d = .true.
344	this%i_permute_4d = ipermute
345
346	end subroutine permute_4d_arrays
347
348
349	! --- PRIVATE SUBROUTINES ---
350
351	!---------------------------------------------------------------------
352	! Return the NetCDF variable ID for variable "var_name", or abort if
353	! not present
354	subroutine get_variable_id(this, var_name, ivarid)
355	class(netcdf_file) :: this
356	character(len=*), intent(in) :: var_name
357	integer, intent(out) :: ivarid
358
359	integer :: istatus
360
361	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
362	if (istatus /= NF90_NOERR) then
363	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
364	& var_name, ': ', trim(nf90_strerror(istatus))
365	call my_abort('Error reading NetCDF file')
366	end if
367
368	end subroutine get_variable_id
369
370
371	!---------------------------------------------------------------------
372	! Return the array dimensions of variable with specified ID, along
373	! with the number of dimensions and optionally the total number of
374	! elements, or abort if variable not present
375	subroutine get_array_dimensions(this, ivarid, ndims, ndimlens, ntotal)
376	class(netcdf_file) :: this
377	integer, intent(in) :: ivarid
378	integer, intent(out) :: ndims
379	integer, intent(out) :: ndimlens(NF90_MAX_VAR_DIMS)
380	integer(kind=jpib), intent(out), optional :: ntotal
381
382	integer :: j, istatus
383	integer :: idimids(NF90_MAX_VAR_DIMS)
384
385	istatus = nf90_inquire_variable(this%ncid, ivarid, &
386	& ndims=ndims, dimids=idimids)
387	if (istatus /= NF90_NOERR) then
388	write(nulerr,'(a,i0,a,a)') '*** Error inquiring about NetCDF variable with id ', &
389	& ivarid, ': ', trim(nf90_strerror(istatus))
390	call my_abort('Error reading NetCDF file')
391	end if
392
393	ndimlens(:) = 0
394	DO j = 1,ndims
395	istatus = nf90_inquire_dimension(this%ncid, idimids(j), len=ndimlens(j))
396	if (istatus /= NF90_NOERR) then
397	write(nulerr,'(a,i0,a,i0,a,a)') '*** Error reading length of dimension ', &
398	& j, ' of NetCDF variable with id ', ivarid, ': ', &
399	& trim(nf90_strerror(istatus))
400	call my_abort('Error reading NetCDF file')
401	end if
402	end do
403
404	if (present(ntotal)) then
405	ntotal = 1
406	DO j = 1, ndims
407	ntotal = ntotal * ndimlens(j)
408	end do
409	end if
410
411	end subroutine get_array_dimensions
412
413
414	!---------------------------------------------------------------------
415	! End define mode, if in define mode, and check the error code
416	! (errors are possible if variables are too large for the format,
417	! for example)
418	subroutine end_define_mode(this)
419	class(netcdf_file) :: this
420	integer :: istatus
421	if (this%is_define_mode) then
422	istatus = nf90_enddef(this%ncid)
423	if (istatus /= NF90_NOERR) then
424	write(nulerr,'(a,a,a,a)') '*** Error ending define mode when writing ', &
425	& trim(this%file_name), ': ', trim(nf90_strerror(istatus))
426	call my_abort('Error writing NetCDF file')
427	end if
428	this%is_define_mode = .false.
429	end if
430	end subroutine end_define_mode
431
432
433	! --- READING SUBROUTINES ---
434
435	!---------------------------------------------------------------------
436	! Return true if file is open, false otherwise
437	function is_open(this)
438	class(netcdf_file) :: this
439	logical :: is_open
440	is_open = (this%ncid >= 0)
441	end function is_open
442
443	!---------------------------------------------------------------------
444	! Return the number of dimensions of variable with name var_name, or
445	! -1 if the variable is not found
446	function get_rank(this, var_name) result(ndims)
447	class(netcdf_file) :: this
448	character(len=*), intent(in) :: var_name
449
450	integer :: ndims
451	integer :: ivarid
452	integer :: ndimlens(NF90_MAX_VAR_DIMS)
453	integer :: istatus
454
455	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
456	if (istatus /= NF90_NOERR) then
457	if (istatus == NF90_ENOTVAR) then
458	if (this%iverbose >= 1) then
459	write(nulout,'(a,a,a)') ' Warning: variable ', var_name, ' not found'
460	end if
461	else
462	write(nulerr,'(a,a,a)') '*** Error inquiring about variable ', &
463	& var_name, ': ', trim(nf90_strerror(istatus))
464	call my_abort('Error reading NetCDF file')
465	end if
466	ndims = -1
467	else
468	call this%get_array_dimensions(ivarid, ndims, ndimlens)
469	end if
470
471	end function get_rank
472
473
474	!---------------------------------------------------------------------
475	! Return the length of the slowest-varying dimension of variable
476	! with name var_name, or -1 if the variable is not found
477	function get_outer_dimension(this, var_name) result(n)
478	class(netcdf_file) :: this
479	character(len=*), intent(in) :: var_name
480
481	integer :: n
482	integer :: ndims
483	integer :: ivarid
484	integer :: ndimlens(NF90_MAX_VAR_DIMS)
485	integer :: istatus
486
487	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
488	if (istatus /= NF90_NOERR) then
489	if (istatus == NF90_ENOTVAR) then
490	if (this%iverbose >= 1) then
491	write(nulout,'(a,a,a)') ' Warning: variable ', var_name, ' not found'
492	end if
493	else
494	write(nulerr,'(a,a,a)') '*** Error inquiring about variable ', &
495	& var_name, ': ', trim(nf90_strerror(istatus))
496	call my_abort('Error reading NetCDF file')
497	end if
498	n = -1
499	else
500	call this%get_array_dimensions(ivarid, ndims, ndimlens)
501	n = ndimlens(ndims)
502	end if
503
504	end function get_outer_dimension
505
506
507	!---------------------------------------------------------------------
508	! Return true if the variable is present, false otherwise
509	function exists(this, var_name) result(is_present)
510	class(netcdf_file) :: this
511	character(len=*), intent(in) :: var_name
512
513	logical :: is_present
514
515	integer :: ivarid
516	integer :: istatus
517
518	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
519	if (istatus /= NF90_NOERR) then
520	is_present = .false.
521	else
522	is_present = .true.
523	end if
524
525	end function exists
526
527
528	!---------------------------------------------------------------------
529	! Return true if the attribute is present, false otherwise. If
530	! argument "len" is provided then return false if len is smaller
531	! than the length of the attribute. This is useful if you have a
532	! fixed array size and want to check whether the attribute will fit
533	! into it.
534	function attribute_exists(this, var_name, attr_name, len) result(is_present)
535	class(netcdf_file) :: this
536	character(len=*), intent(in) :: var_name, attr_name
537	integer, optional, intent(in) :: len
538
539	logical :: is_present
540	integer :: i_attr_len, ivarid
541	integer :: istatus
542
543	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
544	if (istatus /= NF90_NOERR) then
545	is_present = .false.
546	else
547	istatus = nf90_inquire_attribute(this%ncid, ivarid, attr_name, &
548	& len=i_attr_len)
549	if (istatus /= NF90_NOERR) then
550	is_present = .false.
551	else
552	is_present = .true.
553	if (present(len)) then
554	if (i_attr_len > len) then
555	is_present = .false.
556	end if
557	end if
558	end if
559	end if
560
561	end function attribute_exists
562
563
564	!---------------------------------------------------------------------
565	! Return true if the global attribute is present, false otherwise.
566	! If argument "len" is provided then return false if len is smaller
567	! than the length of the attribute. This is useful if you have a
568	! fixed array size and want to check whether the attribute will fit
569	! into it.
570	function global_attribute_exists(this, attr_name, len) result(is_present)
571	class(netcdf_file) :: this
572	character(len=*), intent(in) :: attr_name
573	integer, optional, intent(in) :: len
574
575	logical :: is_present
576	integer :: i_attr_len
577	integer :: istatus
578
579	istatus = nf90_inquire_attribute(this%ncid, NF90_GLOBAL, attr_name, &
580	& len=i_attr_len)
581	if (istatus /= NF90_NOERR) then
582	is_present = .false.
583	else
584	is_present = .true.
585	if (present(len)) then
586	if (i_attr_len > len) then
587	is_present = .false.
588	end if
589	end if
590	end if
591
592	end function global_attribute_exists
593
594
595	!---------------------------------------------------------------------
596	! The method "get" will read either a scalar, vector or matrix
597	! depending on the rank of the output argument. This version reads a
598	! scalar.
599	subroutine get_real_scalar(this, var_name, scalar)
600	class(netcdf_file) :: this
601	character(len=*), intent(in) :: var_name
602	real(jprb), intent(out) :: scalar
603
604	integer :: istatus
605	integer :: ivarid, ndims
606	integer :: ndimlens(NF90_MAX_VAR_DIMS)
607	integer :: j, ntotal
608
609	! Inquire the ID, shape & size of the variable
610	call this%get_variable_id(var_name, ivarid)
611	call this%get_array_dimensions(ivarid, ndims, ndimlens)
612
613	! Compute number of elements of the variable in the file
614	ntotal = 1
615	DO j = 1, ndims
616	ntotal = ntotal * ndimlens(j)
617	end do
618
619	if (this%iverbose >= 3) then
620	write(nulout,'(a,a)',advance='no') ' Reading ', var_name
621	call this%print_variable_attributes(ivarid,nulout)
622	end if
623
624	! Abort if the number of elements is anything other than 1
625	if (ntotal /= 1) then
626	write(nulerr,'(a,a,a,i0,a)') '*** Error reading NetCDF variable ', &
627	& var_name, ' with total length ', ntotal, ' as a scalar'
628	call my_abort('Error reading NetCDF file')
629	end if
630
631	! Read variable
632	istatus = nf90_get_var(this%ncid, ivarid, scalar)
633	if (istatus /= NF90_NOERR) then
634	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
635	& var_name, ' as a scalar: ', trim(nf90_strerror(istatus))
636	call my_abort('Error reading NetCDF file')
637	end if
638
639	end subroutine get_real_scalar
640
641
642	!---------------------------------------------------------------------
643	! Read an integer scalar
644	subroutine get_int_scalar(this, var_name, scalar)
645	class(netcdf_file) :: this
646	character(len=*), intent(in) :: var_name
647	integer, intent(out):: scalar
648
649	integer :: istatus
650	integer :: ivarid, ndims
651	integer :: ndimlens(NF90_MAX_VAR_DIMS)
652	integer :: j, ntotal
653
654	! Inquire the ID, shape & size of the variable
655	call this%get_variable_id(var_name, ivarid)
656	call this%get_array_dimensions(ivarid, ndims, ndimlens)
657
658	! Compute number of elements of the variable in the file
659	ntotal = 1
660	DO j = 1, ndims
661	ntotal = ntotal * ndimlens(j)
662	end do
663
664	if (this%iverbose >= 3) then
665	write(nulout,'(a,a)',advance='no') ' Reading ', var_name
666	call this%print_variable_attributes(ivarid,nulout)
667	end if
668
669	! Abort if the number of elements is anything other than 1
670	if (ntotal /= 1) then
671	write(nulerr,'(a,a,a,i0,a)') '*** Error reading NetCDF variable ', &
672	& var_name, ' with total length ', ntotal, ' as a scalar'
673	call my_abort('Error reading NetCDF file')
674	end if
675
676	! Read variable
677	istatus = nf90_get_var(this%ncid, ivarid, scalar)
678	if (istatus /= NF90_NOERR) then
679	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
680	& var_name, ' as a scalar: ', trim(nf90_strerror(istatus))
681	call my_abort('Error reading NetCDF file')
682	end if
683
684	end subroutine get_int_scalar
685
686
687	!---------------------------------------------------------------------
688	! Read a scalar from a larger array, where "index" indexes the most
689	! slowly varying dimension
690	subroutine get_real_scalar_indexed(this, var_name, scalar, index)
691	class(netcdf_file) :: this
692	character(len=*), intent(in) :: var_name
693	integer, intent(in) :: index
694	real(jprb), intent(out) :: scalar
695
696	integer :: istatus
697	integer :: ivarid, ndims
698	integer :: ndimlens(NF90_MAX_VAR_DIMS)
699	integer :: vstart(NF90_MAX_VAR_DIMS)
700	integer :: j, ntotal
701
702	! Inquire the ID, shape & size of the variable
703	call this%get_variable_id(var_name, ivarid)
704	call this%get_array_dimensions(ivarid, ndims, ndimlens)
705
706	! Compute number of elements of the slice in the file,
707	! i.e. excluding the slowest varying dimension, indexed by "index"
708	ntotal = 1
709	DO j = 1, ndims-1
710	ntotal = ntotal * ndimlens(j)
711	end do
712
713	if (this%iverbose >= 3) then
714	write(nulout,'(a,a,i0,a,a)') ' Reading element ', index, ' of ', var_name
715	end if
716
717	! Abort if the number of elements is anything other than 1
718	if (ntotal /= 1) then
719	write(nulerr,'(a,a,a,i0,a)') '*** Error reading NetCDF variable ', &
720	& var_name, ' with slice length ', ntotal, ' as a scalar'
721	call my_abort('Error reading NetCDF file')
722	end if
723
724	if (index < 1 .or. index > ndimlens(ndims)) then
725	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error reading element ', index, &
726	& ' of NetCDF variable ', &
727	& var_name, ' with outer dimension ', ndimlens(ndims)
728	call my_abort('Error reading NetCDF file')
729	end if
730
731	! Read variable
732	vstart(1:ndims-1) = 1
733	vstart(ndims) = index
734	istatus = nf90_get_var(this%ncid, ivarid, scalar, start=vstart)
735	if (istatus /= NF90_NOERR) then
736	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
737	& var_name, ' as a scalar: ', trim(nf90_strerror(istatus))
738	call my_abort('Error reading NetCDF file')
739	end if
740
741	end subroutine get_real_scalar_indexed
742
743
744	!---------------------------------------------------------------------
745	! Read a 1D real array into "vector", which must be allocatable and
746	! will be reallocated if necessary
747	subroutine get_real_vector(this, var_name, vector)
748	class(netcdf_file) :: this
749	character(len=*), intent(in) :: var_name
750	real(jprb), allocatable, intent(out) :: vector(:)
751
752	integer :: n ! Length of vector
753	integer :: istatus
754	integer :: ivarid, ndims
755	integer :: ndimlens(NF90_MAX_VAR_DIMS)
756	integer :: j
757
758	call this%get_variable_id(var_name, ivarid)
759	call this%get_array_dimensions(ivarid, ndims, ndimlens)
760
761	! Ensure variable has only one dimension in the file
762	n = 1
763	DO j = 1, ndims
764	n = n * ndimlens(j)
765	if (j > 1 .and. ndimlens(j) > 1) then
766	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
767	& var_name, &
768	& ' as a vector: all dimensions above the first must be singletons'
769	call my_abort('Error reading NetCDF file')
770	end if
771	end do
772
773	! Reallocate if necessary
774	if (allocated(vector)) then
775	if (size(vector) /= n) then
776	if (this%iverbose >= 1) then
777	write(nulout,'(a,a)') ' Warning: resizing vector to read ', var_name
778	end if
779	deallocate(vector)
780	allocate(vector(n))
781	end if
782	else
783	allocate(vector(n))
784	end if
785
786	if (this%iverbose >= 3) then
787	write(nulout,'(a,a,a,i0,a)',advance='no') ' Reading ', var_name, '(', n, ')'
788	call this%print_variable_attributes(ivarid,nulout)
789	end if
790
791	! Read variable
792	istatus = nf90_get_var(this%ncid, ivarid, vector)
793	if (istatus /= NF90_NOERR) then
794	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
795	& var_name, ' as a vector: ', trim(nf90_strerror(istatus))
796	call my_abort('Error reading NetCDF file')
797	end if
798
799	end subroutine get_real_vector
800
801
802	!---------------------------------------------------------------------
803	! Read a 1D character array into "vector", which must be allocatable
804	! and will be reallocated if necessary
805	subroutine get_char_vector(this, var_name, vector)
806	class(netcdf_file) :: this
807	character(len=*), intent(in) :: var_name
808	character(len=1), allocatable, intent(out) :: vector(:)
809
810	integer :: n ! Length of vector
811	integer :: istatus
812	integer :: ivarid, ndims
813	integer :: ndimlens(NF90_MAX_VAR_DIMS)
814	integer :: j
815
816	call this%get_variable_id(var_name, ivarid)
817	call this%get_array_dimensions(ivarid, ndims, ndimlens)
818
819	! Ensure variable has only one dimension in the file
820	n = 1
821	DO j = 1, ndims
822	n = n * ndimlens(j)
823	if (j > 1 .and. ndimlens(j) > 1) then
824	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
825	& var_name, &
826	& ' as a vector: all dimensions above the first must be singletons'
827	call my_abort('Error reading NetCDF file')
828	end if
829	end do
830
831	! Reallocate if necessary
832	if (allocated(vector)) then
833	if (size(vector) /= n) then
834	if (this%iverbose >= 1) then
835	write(nulout,'(a,a)') ' Warning: resizing vector to read ', var_name
836	end if
837	deallocate(vector)
838	allocate(vector(n))
839	end if
840	else
841	allocate(vector(n))
842	end if
843
844	if (this%iverbose >= 3) then
845	write(nulout,'(a,a,a,i0,a)',advance='no') ' Reading ', var_name, '(', n, ')'
846	call this%print_variable_attributes(ivarid,nulout)
847	end if
848
849	! Read variable
850	istatus = nf90_get_var(this%ncid, ivarid, vector)
851	if (istatus /= NF90_NOERR) then
852	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
853	& var_name, ' as a vector of chars: ', trim(nf90_strerror(istatus))
854	call my_abort('Error reading NetCDF file')
855	end if
856
857	end subroutine get_char_vector
858
859
860	!---------------------------------------------------------------------
861	! Read a 1D integer array into "vector", which must be allocatable
862	! and will be reallocated if necessary
863	subroutine get_int_vector(this, var_name, vector)
864
865	class(netcdf_file) :: this
866	character(len=*), intent(in) :: var_name
867	integer, allocatable, intent(out) :: vector(:)
868
869	integer :: n ! Length of vector
870	integer :: istatus
871	integer :: ivarid, ndims
872	integer :: ndimlens(NF90_MAX_VAR_DIMS)
873	integer :: j
874
875	call this%get_variable_id(var_name, ivarid)
876	call this%get_array_dimensions(ivarid, ndims, ndimlens)
877
878	! Ensure variable has only one dimension in the file
879	n = 1
880	DO j = 1, ndims
881	n = n * ndimlens(j)
882	if (j > 1 .and. ndimlens(j) > 1) then
883	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
884	& var_name, &
885	& ' as a vector: all dimensions above the first must be singletons'
886	call my_abort('Error reading NetCDF file')
887	end if
888	end do
889
890	! Reallocate if necessary
891	if (allocated(vector)) then
892	if (size(vector) /= n) then
893	if (this%iverbose >= 1) then
894	write(nulout,'(a,a)') ' Warning: resizing vector to read ', var_name
895	end if
896	deallocate(vector)
897	allocate(vector(n))
898	end if
899	else
900	allocate(vector(n))
901	end if
902
903	if (this%iverbose >= 3) then
904	write(nulout,'(a,a,a,i0,a)',advance='no') ' Reading ', var_name, '(', n, ')'
905	call this%print_variable_attributes(ivarid,nulout)
906	end if
907
908	! Read variable
909	istatus = nf90_get_var(this%ncid, ivarid, vector)
910	if (istatus /= NF90_NOERR) then
911	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
912	& var_name, ' as an integer vector: ', trim(nf90_strerror(istatus))
913	call my_abort('Error reading NetCDF file')
914	end if
915
916	end subroutine get_int_vector
917
918	!---------------------------------------------------------------------
919	! Read a vector of data from a larger array; the vector must be
920	! allocatable and will be reallocated if necessary
921	subroutine get_real_vector_indexed(this, var_name, vector, index)
922	class(netcdf_file) :: this
923	character(len=*), intent(in) :: var_name
924	integer, intent(in) :: index
925	real(jprb), allocatable, intent(out) :: vector(:)
926
927	integer :: n ! Length of vector
928	integer :: istatus
929	integer :: ivarid, ndims
930	integer :: ndimlens(NF90_MAX_VAR_DIMS)
931	integer :: vstart(NF90_MAX_VAR_DIMS)
932	integer :: vcount(NF90_MAX_VAR_DIMS)
933	integer :: j
934
935	call this%get_variable_id(var_name, ivarid)
936	call this%get_array_dimensions(ivarid, ndims, ndimlens)
937
938	! Ensure variable has only one dimension aside from the last one
939	n = 1
940	DO j = 1, ndims-1
941	n = n * ndimlens(j)
942	if (j > 1 .and. ndimlens(j) > 1) then
943	write(nulerr,'(a,a,a)') '*** Error reading 1D slice from NetCDF variable ', &
944	& var_name, &
945	& ': all dimensions except the first and last must be singletons'
946	call my_abort('Error reading NetCDF file')
947	end if
948	end do
949
950	! Reallocate if necessary
951	if (allocated(vector)) then
952	if (size(vector) /= n) then
953	if (this%iverbose >= 1) then
954	write(nulout,'(a,i0,a,a)') ' Warning: resizing vector to length ', n, &
955	& ' to read slice of ', var_name
956	end if
957	deallocate(vector)
958	allocate(vector(n))
959	end if
960	else
961	allocate(vector(n))
962	end if
963
964	if (this%iverbose >= 3) then
965	write(nulout,'(a,i0,a,a)') ' Reading column ', index, ' of ', var_name
966	end if
967
968	if (index < 1 .or. index > ndimlens(ndims)) then
969	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error reading element ', index, &
970	& ' of NetCDF variable ', &
971	& var_name, ' with outer dimension ', ndimlens(ndims)
972	call my_abort('Error reading NetCDF file')
973	end if
974
975	! Read variable
976	vstart(1:ndims-1) = 1
977	vstart(ndims) = index
978	vcount(1:ndims-1) = ndimlens(1:ndims-1)
979	vcount(ndims) = 1
980	istatus = nf90_get_var(this%ncid, ivarid, vector, &
981	& start=vstart, count=vcount)
982	if (istatus /= NF90_NOERR) then
983	write(nulerr,'(a,a,a,a)') '*** Error reading 1D slice of NetCDF variable ', &
984	& var_name, ': ', trim(nf90_strerror(istatus))
985	call my_abort('Error reading NetCDF file')
986	end if
987
988	end subroutine get_real_vector_indexed
989
990
991	!---------------------------------------------------------------------
992	! Read 2D array into "matrix", which must be allocatable and will be
993	! reallocated if necessary. Whether to transpose is specifed by the
994	! final optional argument, but can also be specified by the
995	! do_transpose_2d class data member.
996	subroutine get_real_matrix(this, var_name, matrix, do_transp)
997	class(netcdf_file) :: this
998	character(len=*), intent(in) :: var_name
999	real(jprb), allocatable, intent(out) :: matrix(:,:)
1000	logical, optional, intent(in):: do_transp ! Transpose data?
1001
1002	real(jprb), allocatable :: tmp_matrix(:,:)
1003	integer :: ndimlen1, ndimlen2
1004	integer :: istatus
1005	integer :: ivarid, ndims
1006	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1007	integer :: j, ntotal
1008	logical :: do_transpose
1009
1010	! Decide whether to transpose the array
1011	if (present(do_transp)) then
1012	do_transpose = do_transp
1013	else
1014	do_transpose = this%do_transpose_2d
1015	end if
1016
1017	call this%get_variable_id(var_name, ivarid)
1018	call this%get_array_dimensions(ivarid, ndims, ndimlens)
1019
1020	! Ensure the variable has no more than two non-singleton
1021	! dimensions
1022	ntotal = 1
1023	DO j = 1, ndims
1024	ntotal = ntotal * ndimlens(j)
1025	if (j > 2 .and. ndimlens(j) > 1) then
1026	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
1027	& var_name, &
1028	& ' as a matrix: all dimensions above the second must be singletons'
1029	call my_abort('Error reading NetCDF file')
1030	end if
1031	end do
1032
1033	! Work out dimension lengths
1034	if (ndims >= 2) then
1035	ndimlen1 = ndimlens(1)
1036	ndimlen2 = ntotal/ndimlen1
1037	else
1038	ndimlen1 = ntotal
1039	ndimlen2 = 1
1040	end if
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	deallocate(matrix)
1053	allocate(matrix(ndimlen2, ndimlen1))
1054	end if
1055	else
1056	allocate(matrix(ndimlen2, ndimlen1))
1057	end if
1058
1059	if (this%iverbose >= 3) then
1060	write(nulout,'(a,a,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, '(', &
1061	& ndimlen2, ',', ndimlen1, ')'
1062	call this%print_variable_attributes(ivarid,nulout)
1063	end if
1064
1065	istatus = nf90_get_var(this%ncid, ivarid, tmp_matrix)
1066	matrix = transpose(tmp_matrix)
1067	deallocate(tmp_matrix)
1068	else
1069	! Read data without transposition
1070
1071	! Reallocate if necessary
1072	if (allocated(matrix)) then
1073	if (size(matrix,1) /= ndimlen1 .or. size(matrix,2) /= ndimlen2) then
1074	if (this%iverbose >= 1) then
1075	write(nulout,'(a,a)') ' Warning: resizing matrix to read ', var_name
1076	end if
1077	allocate(matrix(ndimlen1, ndimlen2))
1078	end if
1079	else
1080	allocate(matrix(ndimlen1, ndimlen2))
1081	end if
1082
1083	if (this%iverbose >= 3) then
1084	write(nulout,'(a,a,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, '(', &
1085	& ndimlen1, ',', ndimlen2, ')'
1086	call this%print_variable_attributes(ivarid,nulout)
1087	end if
1088
1089	istatus = nf90_get_var(this%ncid, ivarid, matrix)
1090	end if
1091
1092	if (istatus /= NF90_NOERR) then
1093	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
1094	& var_name, ' as a matrix: ', trim(nf90_strerror(istatus))
1095	call my_abort('Error reading NetCDF file')
1096	end if
1097
1098	end subroutine get_real_matrix
1099
1100
1101	!---------------------------------------------------------------------
1102	! Read 2D array of characters into "matrix", which must be
1103	! allocatable and will be reallocated if necessary. Whether to
1104	! transpose is specifed by the final optional argument, but can also
1105	! be specified by the do_transpose_2d class data member.
1106	subroutine get_char_matrix(this, var_name, matrix, do_transp)
1107	class(netcdf_file) :: this
1108	character(len=*), intent(in) :: var_name
1109	character(len=1), allocatable, intent(inout) :: matrix(:,:)
1110	logical, optional, intent(in):: do_transp ! Transpose data?
1111
1112	character(len=1), allocatable:: tmp_matrix(:,:)
1113	integer :: ndimlen1, ndimlen2
1114	integer :: istatus
1115	integer :: ivarid, ndims
1116	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1117	integer :: vstart(NF90_MAX_VAR_DIMS)
1118	integer :: vcount(NF90_MAX_VAR_DIMS)
1119	integer :: j, ntotal
1120	logical :: do_transpose
1121
1122	! Decide whether to transpose the array
1123	if (present(do_transp)) then
1124	do_transpose = do_transp
1125	else
1126	do_transpose = this%do_transpose_2d
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 two non-singleton
1133	! dimensions
1134	ntotal = 1
1135	DO j = 1, ndims
1136	ntotal = ntotal * ndimlens(j)
1137	if (j > 2 .and. ndimlens(j) > 1) then
1138	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
1139	& var_name, &
1140	& ' as a matrix: all dimensions above the second 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 >= 2) then
1147	ndimlen1 = ndimlens(1)
1148	ndimlen2 = ntotal/ndimlen1
1149	else
1150	ndimlen1 = ntotal
1151	ndimlen2 = 1
1152	end if
1153
1154	if (do_transpose) then
1155	! Read and transpose
1156	allocate(tmp_matrix(ndimlen1, ndimlen2))
1157
1158	! Reallocate if necessary
1159	if (allocated(matrix)) then
1160	if (size(matrix,1) /= ndimlen2 .or. size(matrix,2) /= ndimlen1) then
1161	if (this%iverbose >= 1) then
1162	write(nulout,'(a,a)') ' Warning: resizing matrix to read ', var_name
1163	end if
1164	deallocate(matrix)
1165	allocate(matrix(ndimlen2, ndimlen1))
1166	end if
1167	else
1168	allocate(matrix(ndimlen2, ndimlen1))
1169	end if
1170
1171	if (this%iverbose >= 3) then
1172	write(nulout,'(a,a,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, '(', &
1173	& ndimlen2, ',', ndimlen1, ')'
1174	call this%print_variable_attributes(ivarid,nulout)
1175	end if
1176
1177	istatus = nf90_get_var(this%ncid, ivarid, tmp_matrix)
1178	matrix = transpose(tmp_matrix)
1179	deallocate(tmp_matrix)
1180	else
1181	! Read data without transposition
1182
1183	! Reallocate if necessary
1184	if (allocated(matrix)) then
1185	if (size(matrix,1) /= ndimlen1 .or. size(matrix,2) /= ndimlen2) then
1186	if (this%iverbose >= 1) then
1187	write(nulout,'(a,a)') ' Warning: resizing matrix to read ', var_name
1188	end if
1189	allocate(matrix(ndimlen1, ndimlen2))
1190	end if
1191	else
1192	allocate(matrix(ndimlen1, ndimlen2))
1193	end if
1194
1195	if (this%iverbose >= 3) then
1196	write(nulout,'(a,a,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, '(', &
1197	& ndimlen1, ',', ndimlen2, ')'
1198	call this%print_variable_attributes(ivarid,nulout)
1199	end if
1200
1201	vstart = 1
1202	vcount(1:2) = [ndimlen1,1]
1203
1204	DO j = 1,ndimlen2
1205	vstart(2) = j
1206	istatus = nf90_get_var(this%ncid, ivarid, matrix(:,j), start=vstart, count=vcount)
1207	end do
1208	end if
1209
1210	if (istatus /= NF90_NOERR) then
1211	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
1212	& var_name, ' as a matrix of characters: ', trim(nf90_strerror(istatus))
1213	call my_abort('Error reading NetCDF file')
1214	end if
1215
1216	end subroutine get_char_matrix
1217
1218
1219	!---------------------------------------------------------------------
1220	! Read matrix of data from a larger array, which must be allocatable
1221	! and will be reallocated if necessary. Whether to transpose is
1222	! specifed by the final optional argument, but can also be specified
1223	! by the do_transpose_2d class data member.
1224	subroutine get_real_matrix_indexed(this, var_name, matrix, index, do_transp)
1225	class(netcdf_file) :: this
1226	character(len=*), intent(in) :: var_name
1227	integer, intent(in) :: index
1228	real(jprb), allocatable, intent(out) :: matrix(:,:)
1229	logical, optional, intent(in):: do_transp ! Transpose data?
1230
1231	real(jprb), allocatable :: tmp_matrix(:,:)
1232	integer :: ndimlen1, ndimlen2
1233	integer :: istatus
1234	integer :: ivarid, ndims
1235	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1236	integer :: vstart(NF90_MAX_VAR_DIMS)
1237	integer :: vcount(NF90_MAX_VAR_DIMS)
1238	integer :: j, ntotal
1239	logical :: do_transpose
1240
1241	! Decide whether to transpose the array
1242	if (present(do_transp)) then
1243	do_transpose = do_transp
1244	else
1245	do_transpose = this%do_transpose_2d
1246	end if
1247
1248	call this%get_variable_id(var_name, ivarid)
1249	call this%get_array_dimensions(ivarid, ndims, ndimlens)
1250
1251	! Ensure the variable has no more than two non-singleton
1252	! dimensions aside from the last one
1253	ntotal = 1
1254	DO j = 1, ndims-1
1255	ntotal = ntotal * ndimlens(j)
1256	if (j > 2 .and. ndimlens(j) > 1) then
1257	write(nulerr,'(a,a,a)') '*** Error reading 2D slice from NetCDF variable ', &
1258	& var_name, &
1259	& ': all dimensions except the first, second and last must be singletons'
1260	call my_abort('Error reading NetCDF file')
1261	end if
1262	end do
1263
1264	if (index < 1 .or. index > ndimlens(ndims)) then
1265	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error reading element ', index, &
1266	& ' of NetCDF variable ', &
1267	& var_name, ' with outer dimension ', ndimlens(ndims)
1268	call my_abort('Error reading NetCDF file')
1269	end if
1270
1271	! Work out dimension lengths
1272	if (ndims >= 2) then
1273	ndimlen1 = ndimlens(1)
1274	ndimlen2 = ntotal/ndimlen1
1275	else
1276	ndimlen1 = ntotal
1277	ndimlen2 = 1
1278	end if
1279
1280	vstart(1:ndims-1) = 1
1281	vstart(ndims) = index
1282	vcount(1:ndims-1) = ndimlens(1:ndims-1)
1283	vcount(ndims) = 1
1284
1285	if (do_transpose) then
1286	! Read and transpose
1287	allocate(tmp_matrix(ndimlen1, ndimlen2))
1288
1289	! Reallocate if necessary
1290	if (allocated(matrix)) then
1291	if (size(matrix,1) /= ndimlen2 .or. size(matrix,2) /= ndimlen1) then
1292	if (this%iverbose >= 1) then
1293	write(nulout,'(a,a)') ' Warning: resizing matrix to read ', var_name
1294	end if
1295	allocate(matrix(ndimlen2, ndimlen1))
1296	end if
1297	else
1298	allocate(matrix(ndimlen2, ndimlen1))
1299	end if
1300
1301	if (this%iverbose >= 3) then
1302	write(nulout,'(a,i0,a,a,a,i0,a,i0,a)') ' Reading slice ', index, &
1303	& ' of ', var_name, ' as ', ndimlen2, 'x', ndimlen1, ' array'
1304	end if
1305
1306	istatus = nf90_get_var(this%ncid, ivarid, tmp_matrix, &
1307	& start=vstart, count=vcount)
1308	matrix = transpose(tmp_matrix)
1309	deallocate(tmp_matrix)
1310	else
1311	! Read data without transposition
1312
1313	! Reallocate if necessary
1314	if (allocated(matrix)) then
1315	if (size(matrix,1) /= ndimlen1 .or. size(matrix,2) /= ndimlen2) then
1316	if (this%iverbose >= 1) then
1317	write(nulout,'(a,a)') ' Warning: resizing matrix to read ', var_name
1318	end if
1319	allocate(matrix(ndimlen1, ndimlen2))
1320	end if
1321	else
1322	allocate(matrix(ndimlen1, ndimlen2))
1323	end if
1324
1325	if (this%iverbose >= 3) then
1326	write(nulout,'(a,i0,a,a,a,i0,a,i0,a)') ' Reading slice ', index, &
1327	& ' of ', var_name, ' as ', ndimlen1, 'x', ndimlen2, ' array'
1328	end if
1329
1330	istatus = nf90_get_var(this%ncid, ivarid, matrix, &
1331	& start=vstart, count=vcount)
1332	end if
1333
1334	if (istatus /= NF90_NOERR) then
1335	write(nulerr,'(a,a,a,a)') '*** Error reading 2D slice of NetCDF variable ', &
1336	& var_name, ': ', trim(nf90_strerror(istatus))
1337	call my_abort('Error reading NetCDF file')
1338	end if
1339
1340	end subroutine get_real_matrix_indexed
1341
1342
1343	!---------------------------------------------------------------------
1344	! Read 3D array into "var", which must be allocatable and will be
1345	! reallocated if necessary. Whether to pemute is specifed by the
1346	! final optional argument
1347	subroutine get_real_array3(this, var_name, var, ipermute)
1348	class(netcdf_file) :: this
1349	character(len=*), intent(in) :: var_name
1350	real(jprb), allocatable, intent(out) :: var(:,:,:)
1351	integer, optional, intent(in) :: ipermute(3)
1352
1353	real(jprb), allocatable :: var_permute(:,:,:)
1354	integer :: ndimlen1, ndimlen2, ndimlen3
1355	integer :: istatus
1356	integer :: ivarid, ndims
1357	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1358	integer :: j, ntotal
1359	integer :: n_dimlens_permuted(3)
1360	integer :: i_permute_3d(3)
1361	logical :: do_permute
1362
1363	! Decide whether to permute
1364	if (present(ipermute)) then
1365	do_permute = .true.
1366	i_permute_3d = ipermute
1367	else
1368	do_permute = this%do_permute_3d
1369	i_permute_3d = this%i_permute_3d
1370	end if
1371
1372	call this%get_variable_id(var_name, ivarid)
1373	call this%get_array_dimensions(ivarid, ndims, ndimlens)
1374
1375	! Ensure the variable has no more than three non-singleton
1376	! dimensions
1377	ntotal = 1
1378	DO j = 1, ndims
1379	ntotal = ntotal * ndimlens(j)
1380	if (j > 3 .and. ndimlens(j) > 1) then
1381	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
1382	& var_name, &
1383	& ' as a 3D array: all dimensions above the third must be singletons'
1384	call my_abort('Error reading NetCDF file')
1385	end if
1386	end do
1387
1388	! Work out dimension lengths
1389	if (ndims >= 3) then
1390	ndimlen1 = ndimlens(1)
1391	ndimlen2 = ndimlens(2)
1392	ndimlen3 = ntotal/(ndimlen1*ndimlen2)
1393	else if (ndims >= 2) then
1394	ndimlen1 = ndimlens(1)
1395	ndimlen2 = ntotal/ndimlen1
1396	ndimlen3 = 1
1397	else
1398	ndimlen1 = ntotal
1399	ndimlen2 = 1
1400	ndimlen3 = 1
1401	end if
1402
1403	! Deallocate if necessary
1404	if (allocated(var)) then
1405	deallocate(var)
1406	end if
1407
1408	if (do_permute) then
1409	! Read and permute
1410	allocate(var_permute(ndimlen1, ndimlen2, ndimlen3))
1411	n_dimlens_permuted(i_permute_3d) = ndimlens(1:3)
1412
1413	! Reallocate if necessary
1414	if (allocated(var)) then
1415	if (size(var,1) /= n_dimlens_permuted(1) &
1416	& .or. size(var,2) /= n_dimlens_permuted(2) &
1417	& .or. size(var,3) /= n_dimlens_permuted(3)) then
1418	if (this%iverbose >= 1) then
1419	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1420	end if
1421	deallocate(var)
1422	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1423	& n_dimlens_permuted(3)))
1424	end if
1425	else
1426	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1427	& n_dimlens_permuted(3)))
1428	end if
1429
1430	if (this%iverbose >= 3) then
1431	write(nulout,'(a,a,a,i0,i0,i0,a)',advance='no') ' Reading ', var_name, &
1432	& ' (permuting dimensions ', i_permute_3d, ')'
1433	call this%print_variable_attributes(ivarid,nulout)
1434	end if
1435
1436	istatus = nf90_get_var(this%ncid, ivarid, var_permute)
1437	var = reshape(var_permute, n_dimlens_permuted, order=i_permute_3d)
1438	deallocate(var_permute)
1439
1440	else
1441	! Read data without permutation
1442
1443	! Reallocate if necessary
1444	if (allocated(var)) then
1445	if (size(var,1) /= ndimlen1 &
1446	& .or. size(var,2) /= ndimlen2 &
1447	& .or. size(var,3) /= ndimlen3) then
1448	if (this%iverbose >= 1) then
1449	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1450	end if
1451	deallocate(var)
1452	allocate(var(ndimlen1, ndimlen2, ndimlen3))
1453	end if
1454	else
1455	allocate(var(ndimlen1, ndimlen2, ndimlen3))
1456	end if
1457
1458	if (this%iverbose >= 3) then
1459	write(nulout,'(a,a,a,i0,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, &
1460	& '(', ndimlen1, ',', ndimlen2, ',', ndimlen3, ')'
1461	call this%print_variable_attributes(ivarid,nulout)
1462	end if
1463
1464	istatus = nf90_get_var(this%ncid, ivarid, var)
1465	end if
1466
1467	if (istatus /= NF90_NOERR) then
1468	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
1469	& var_name, ' as a 3D array: ', trim(nf90_strerror(istatus))
1470	call my_abort('Error reading NetCDF file')
1471	end if
1472
1473	end subroutine get_real_array3
1474
1475
1476	!---------------------------------------------------------------------
1477	! Read 3D array of data from a larger-dimensional array, which must
1478	! be allocatable and will be reallocated if necessary. Whether to
1479	! pemute is specifed by the final optional argument
1480	subroutine get_real_array3_indexed(this, var_name, var, index, ipermute)
1481	class(netcdf_file) :: this
1482	character(len=*), intent(in) :: var_name
1483	integer, intent(in) :: index
1484	real(jprb), allocatable, intent(out) :: var(:,:,:)
1485	integer, optional, intent(in) :: ipermute(3)
1486
1487	real(jprb), allocatable :: var_permute(:,:,:)
1488	integer :: ndimlen1, ndimlen2, ndimlen3
1489	integer :: istatus
1490	integer :: ivarid, ndims
1491	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1492	integer :: vstart(NF90_MAX_VAR_DIMS)
1493	integer :: vcount(NF90_MAX_VAR_DIMS)
1494	integer :: j, ntotal
1495	integer :: n_dimlens_permuted(3)
1496	integer :: i_permute_3d(3)
1497	logical :: do_permute
1498
1499	! Decide whether to permute
1500	if (present(ipermute)) then
1501	do_permute = .true.
1502	i_permute_3d = ipermute
1503	else
1504	do_permute = this%do_permute_3d
1505	i_permute_3d = this%i_permute_3d
1506	end if
1507
1508	call this%get_variable_id(var_name, ivarid)
1509	call this%get_array_dimensions(ivarid, ndims, ndimlens)
1510
1511	! Ensure the variable has no more than three non-singleton
1512	! dimensions aside from the last one
1513	ntotal = 1
1514	DO j = 1, ndims-1
1515	ntotal = ntotal * ndimlens(j)
1516	if (j > 3 .and. ndimlens(j) > 1) then
1517	write(nulerr,'(a,a,a)') '*** Error reading 3D slice from NetCDF variable ', &
1518	& var_name, &
1519	& ': all dimensions above the third must be singletons'
1520	call my_abort('Error reading NetCDF file')
1521	end if
1522	end do
1523
1524	if (index < 1 .or. index > ndimlens(ndims)) then
1525	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error reading element ', index, &
1526	& ' of NetCDF variable ', &
1527	& var_name, ' with outer dimension ', ndimlens(ndims)
1528	call my_abort('Error reading NetCDF file')
1529	end if
1530
1531	! Work out dimension lengths
1532	if (ndims >= 3) then
1533	ndimlen1 = ndimlens(1)
1534	ndimlen2 = ndimlens(2)
1535	ndimlen3 = ntotal/(ndimlen1*ndimlen2)
1536	else if (ndims >= 2) then
1537	ndimlen1 = ndimlens(1)
1538	ndimlen2 = ntotal/ndimlen1
1539	ndimlen3 = 1
1540	else
1541	ndimlen1 = ntotal
1542	ndimlen2 = 1
1543	ndimlen3 = 1
1544	end if
1545
1546	vstart(1:ndims-1) = 1
1547	vstart(ndims) = index
1548	vcount(1:ndims-1) = ndimlens(1:ndims-1)
1549	vcount(ndims) = 1
1550
1551	if (do_permute) then
1552	! Read and permute
1553	allocate(var_permute(ndimlen1, ndimlen2, ndimlen3))
1554	n_dimlens_permuted(i_permute_3d) = ndimlens(1:3)
1555
1556	! Reallocate if necessary
1557	if (allocated(var)) then
1558	if (size(var,1) /= n_dimlens_permuted(1) &
1559	& .or. size(var,2) /= n_dimlens_permuted(2) &
1560	& .or. size(var,3) /= n_dimlens_permuted(3)) then
1561	if (this%iverbose >= 1) then
1562	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1563	end if
1564	deallocate(var)
1565	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1566	& n_dimlens_permuted(3)))
1567	end if
1568	else
1569	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1570	& n_dimlens_permuted(3)))
1571	end if
1572
1573	if (this%iverbose >= 3) then
1574	write(nulout,'(a,i0,a,a,a,i0,i0,i0,a)') ' Reading slice ', index, &
1575	& ' of ', var_name, &
1576	& ' (permuting dimensions ', i_permute_3d, ')'
1577	end if
1578
1579	istatus = nf90_get_var(this%ncid, ivarid, var_permute, &
1580	& start=vstart, count=vcount)
1581	var = reshape(var_permute, n_dimlens_permuted, order=i_permute_3d)
1582	deallocate(var_permute)
1583
1584	else
1585	! Read data without permutation
1586
1587	! Reallocate if necessary
1588	if (allocated(var)) then
1589	if (size(var,1) /= ndimlen1 &
1590	& .or. size(var,2) /= ndimlen2 &
1591	& .or. size(var,3) /= ndimlen3) then
1592	if (this%iverbose >= 1) then
1593	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1594	end if
1595	deallocate(var)
1596	allocate(var(ndimlen1, ndimlen2, ndimlen3))
1597	end if
1598	else
1599	allocate(var(ndimlen1, ndimlen2, ndimlen3))
1600	end if
1601
1602	if (this%iverbose >= 3) then
1603	write(nulout,'(a,i0,a,a,a,i0,a,i0,a,i0,a)') ' Reading slice ', index, &
1604	& ' of ', var_name, ' as ', ndimlen1, 'x', ndimlen2, 'x', &
1605	& ndimlen3, 'array'
1606	end if
1607
1608	istatus = nf90_get_var(this%ncid, ivarid, var, &
1609	& start=vstart, count=vcount)
1610	end if
1611
1612	if (istatus /= NF90_NOERR) then
1613	write(nulerr,'(a,a,a,a)') '*** Error reading 3D slice of NetCDF variable ', &
1614	& var_name, ': ', trim(nf90_strerror(istatus))
1615	call my_abort('Error reading NetCDF file')
1616	end if
1617
1618	end subroutine get_real_array3_indexed
1619
1620
1621	!---------------------------------------------------------------------
1622	! Read 4D array into "var", which must be allocatable and will be
1623	! reallocated if necessary. Whether to pemute is specifed by the
1624	! final optional argument
1625	subroutine get_real_array4(this, var_name, var, ipermute)
1626	class(netcdf_file) :: this
1627	character(len=*), intent(in) :: var_name
1628	real(jprb), allocatable, intent(out) :: var(:,:,:,:)
1629	integer, optional, intent(in) :: ipermute(4)
1630
1631	real(jprb), allocatable :: var_permute(:,:,:,:)
1632	integer :: ndimlen1, ndimlen2, ndimlen3, ndimlen4
1633	integer :: istatus
1634	integer :: ivarid, ndims
1635	integer :: ndimlens(NF90_MAX_VAR_DIMS)
1636	integer :: j, ntotal
1637	integer :: n_dimlens_permuted(4)
1638	integer :: i_permute_4d(4)
1639	logical :: do_permute
1640
1641	! Decide whether to permute
1642	if (present(ipermute)) then
1643	do_permute = .true.
1644	i_permute_4d = ipermute
1645	else
1646	do_permute = this%do_permute_4d
1647	i_permute_4d = this%i_permute_4d
1648	end if
1649
1650	call this%get_variable_id(var_name, ivarid)
1651	call this%get_array_dimensions(ivarid, ndims, ndimlens)
1652
1653	! Ensure the variable has no more than three non-singleton
1654	! dimensions
1655	ntotal = 1
1656	DO j = 1, ndims
1657	ntotal = ntotal * ndimlens(j)
1658	if (j > 4 .and. ndimlens(j) > 1) then
1659	write(nulerr,'(a,a,a)') '*** Error reading NetCDF variable ', &
1660	& var_name, &
1661	& ' as a 4D array: all dimensions above the third must be singletons'
1662	call my_abort('Error reading NetCDF file')
1663	end if
1664	end do
1665
1666	! Work out dimension lengths
1667	if (ndims >= 4) then
1668	ndimlen1 = ndimlens(1)
1669	ndimlen2 = ndimlens(2)
1670	ndimlen3 = ndimlens(3)
1671	ndimlen4 = ntotal/(ndimlen1ndimlen2ndimlen3)
1672	else if (ndims >= 3) then
1673	ndimlen1 = ndimlens(1)
1674	ndimlen2 = ndimlens(2)
1675	ndimlen3 = ntotal/(ndimlen1*ndimlen2)
1676	else if (ndims >= 2) then
1677	ndimlen1 = ndimlens(1)
1678	ndimlen2 = ntotal/ndimlen1
1679	ndimlen3 = 1
1680	else
1681	ndimlen1 = ntotal
1682	ndimlen2 = 1
1683	ndimlen3 = 1
1684	end if
1685
1686	! Deallocate if necessary
1687	if (allocated(var)) then
1688	deallocate(var)
1689	end if
1690
1691	if (do_permute) then
1692	! Read and permute - not tested
1693	allocate(var_permute(ndimlen1, ndimlen2, ndimlen3, ndimlen4))
1694	n_dimlens_permuted(i_permute_4d) = ndimlens(1:4)
1695
1696	! Reallocate if necessary
1697	if (allocated(var)) then
1698	if (size(var,1) /= n_dimlens_permuted(1) &
1699	& .or. size(var,2) /= n_dimlens_permuted(2) &
1700	& .or. size(var,3) /= n_dimlens_permuted(3) &
1701	& .or. size(var,4) /= n_dimlens_permuted(4)) then
1702	if (this%iverbose >= 1) then
1703	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1704	end if
1705	deallocate(var)
1706	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1707	& n_dimlens_permuted(3), n_dimlens_permuted(4)))
1708	end if
1709	else
1710	allocate(var(n_dimlens_permuted(1), n_dimlens_permuted(2), &
1711	& n_dimlens_permuted(3), n_dimlens_permuted(4)))
1712	end if
1713
1714	if (this%iverbose >= 3) then
1715	write(nulout,'(a,a,a,i0,i0,i0,a)',advance='no') ' Reading ', var_name, &
1716	& ' (permuting dimensions ', i_permute_4d, ')'
1717	call this%print_variable_attributes(ivarid,nulout)
1718	end if
1719
1720	istatus = nf90_get_var(this%ncid, ivarid, var_permute)
1721	var = reshape(var_permute, n_dimlens_permuted, order=i_permute_4d)
1722	deallocate(var_permute)
1723
1724	else
1725	! Read data without permutation
1726
1727	! Reallocate if necessary
1728	if (allocated(var)) then
1729	if (size(var,1) /= ndimlen1 &
1730	& .or. size(var,2) /= ndimlen2 &
1731	& .or. size(var,3) /= ndimlen3 &
1732	& .or. size(var,4) /= ndimlen4) then
1733	if (this%iverbose >= 1) then
1734	write(nulout,'(a,a)') ' Warning: resizing array to read ', var_name
1735	end if
1736	deallocate(var)
1737	allocate(var(ndimlen1, ndimlen2, ndimlen3, ndimlen4))
1738	end if
1739	else
1740	allocate(var(ndimlen1, ndimlen2, ndimlen3, ndimlen4))
1741	end if
1742
1743	if (this%iverbose >= 3) then
1744	write(nulout,'(a,a,a,i0,a,i0,a,i0,a,i0,a)',advance='no') ' Reading ', var_name, &
1745	& '(', ndimlen1, ',', ndimlen2, ',', ndimlen3,',', ndimlen4, ')'
1746	call this%print_variable_attributes(ivarid,nulout)
1747	end if
1748
1749	istatus = nf90_get_var(this%ncid, ivarid, var)
1750	end if
1751
1752	if (istatus /= NF90_NOERR) then
1753	write(nulerr,'(a,a,a,a)') '*** Error reading NetCDF variable ', &
1754	& var_name, ' as a 4D array: ', trim(nf90_strerror(istatus))
1755	call my_abort('Error reading NetCDF file')
1756	end if
1757
1758	end subroutine get_real_array4
1759
1760
1761	!---------------------------------------------------------------------
1762	! Get attribute as a character string
1763	subroutine get_string_attribute(this, var_name, attr_name, attr_str)
1764	class(netcdf_file) :: this
1765
1766	character(len=*), intent(in) :: var_name, attr_name
1767	character(len=*), intent(inout) :: attr_str
1768
1769	integer :: i_attr_len, ivarid
1770	integer :: istatus
1771	integer :: j
1772
1773	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
1774	if (istatus /= NF90_NOERR) then
1775	write(nulerr,'(a,a,a,a)') '*** Error inquiring about variable ', var_name, &
1776	& ': ', trim(nf90_strerror(istatus))
1777	call my_abort('Error reading NetCDF file')
1778	end if
1779	istatus = nf90_inquire_attribute(this%ncid, ivarid, attr_name, &
1780	& len = i_attr_len)
1781	if (istatus /= NF90_NOERR) then
1782	write(nulerr,'(a,a,a,a)') '*** Error reading size of attribute ', attr_name, &
1783	& ': ', trim(nf90_strerror(istatus))
1784	call my_abort('Error reading NetCDF file')
1785	end if
1786
1787	! Allocatable character strings not supported on enough compilers
1788	! yet
1789	! if (allocated(attr_str)) then
1790	! deallocate(attr_str)
1791	! end if
1792	! allocate(character(len=i_attr_len) :: attr_str)
1793	if (len(attr_str) < i_attr_len) then
1794	write(nulerr,'(a,a)') '*** Not enough space to read attribute ', attr_name
1795	call my_abort('Error reading NetCDF file')
1796	end if
1797
1798	istatus = nf90_get_att(this%ncid, ivarid, attr_name, attr_str)
1799	if (istatus /= NF90_NOERR) then
1800	write(nulerr,'(a,a,a,a)') '*** Error reading attribute ', attr_name, &
1801	& ': ', trim(nf90_strerror(istatus))
1802	call my_abort('Error reading NetCDF file')
1803	end if
1804
1805	! Pad with blanks since nf90_get_att does not do this
1806	DO j = i_attr_len+1,len(attr_str)
1807	attr_str(j:j) = ' '
1808	end do
1809
1810	end subroutine get_string_attribute
1811
1812
1813
1814	!---------------------------------------------------------------------
1815	! Get attribute as a real scalar
1816	subroutine get_real_scalar_attribute(this, var_name, attr_name, attr)
1817	class(netcdf_file) :: this
1818
1819	character(len=*), intent(in) :: var_name, attr_name
1820	real(jprb), intent(out) :: attr
1821
1822	integer :: ivarid
1823	integer :: istatus
1824
1825	istatus = nf90_inq_varid(this%ncid, var_name, ivarid)
1826	if (istatus /= NF90_NOERR) then
1827	write(nulerr,'(a,a,a,a)') '*** Error inquiring about variable ', var_name, &
1828	& ': ', trim(nf90_strerror(istatus))
1829	call my_abort('Error reading NetCDF file')
1830	end if
1831	istatus = nf90_get_att(this%ncid, ivarid, attr_name, attr)
1832	if (istatus /= NF90_NOERR) then
1833	write(nulerr,'(a,a,a,a)') '*** Error reading attribute ', attr_name, &
1834	& ': ', trim(nf90_strerror(istatus))
1835	call my_abort('Error reading NetCDF file')
1836	end if
1837
1838	end subroutine get_real_scalar_attribute
1839
1840
1841	!---------------------------------------------------------------------
1842	! Get a global attribute as a character string
1843	subroutine get_global_attribute(this, attr_name, attr_str)
1844	class(netcdf_file) :: this
1845
1846	character(len=*), intent(in) :: attr_name
1847	character(len=*), intent(inout) :: attr_str
1848
1849	integer :: i_attr_len
1850	integer :: istatus
1851	integer :: j
1852
1853	istatus = nf90_inquire_attribute(this%ncid, NF90_GLOBAL, attr_name, &
1854	& len = i_attr_len)
1855	if (istatus /= NF90_NOERR) then
1856	write(nulerr,'(a,a,a,a)') '*** Error reading size of global attribute ', attr_name, &
1857	& ': ', trim(nf90_strerror(istatus))
1858	call my_abort('Error reading NetCDF file')
1859	end if
1860
1861	! Allocatable character strings not supported one enough compilers
1862	! yet
1863	! if (allocated(attr_str)) then
1864	! deallocate(attr_str)
1865	! end if
1866	! allocate(character(len=i_attr_len) :: attr_str)
1867	if (len(attr_str) < i_attr_len) then
1868	write(nulerr,'(a,a)') '*** Not enough space to read global attribute ', attr_name
1869	call my_abort('Error reading NetCDF file')
1870	end if
1871
1872	istatus = nf90_get_att(this%ncid, NF90_GLOBAL, attr_name, attr_str)
1873	if (istatus /= NF90_NOERR) then
1874	write(nulerr,'(a,a,a,a)') '*** Error reading global attribute ', attr_name, &
1875	& ': ', trim(nf90_strerror(istatus))
1876	call my_abort('Error reading NetCDF file')
1877	end if
1878
1879	! Pad with blanks since nf90_get_att does not do this
1880	DO j = i_attr_len+1,len(attr_str)
1881	attr_str(j:j) = ' '
1882	end do
1883
1884	end subroutine get_global_attribute
1885
1886
1887	!---------------------------------------------------------------------
1888	! Print a variable's long_name, units and comment, according to
1889	! verbosity level
1890	subroutine print_variable_attributes(this, ivarid, iunit)
1891	class(netcdf_file) :: this
1892	integer, intent(in) :: ivarid ! NetCDF ID of variable
1893	integer, intent(in) :: iunit ! Unit to print information to
1894
1895	character(len=4000) :: attr_str
1896	integer :: istatus
1897
1898	if (this%iverbose >= 4) then
1899	istatus = nf90_get_att(this%ncid, ivarid, 'long_name', attr_str)
1900	if (istatus == NF90_NOERR) then
1901	write(iunit, '(a,a,a)', advance='no') ': "', trim(attr_str), '"'
1902	istatus = nf90_get_att(this%ncid, ivarid, 'units', attr_str)
1903	if (istatus == NF90_NOERR) then
1904	if (trim(attr_str) == '1') then
1905	write(iunit, '(a)') ' (dimensionless)'
1906	else
1907	write(iunit, '(a,a,a)') ' (', trim(attr_str), ')'
1908	end if
1909	else
1910	! No units present
1911	write(iunit, '(1x)')
1912	end if
1913	if (this%iverbose >= 5) then
1914	istatus = nf90_get_att(this%ncid, ivarid, 'comment', attr_str)
1915	if (istatus == NF90_NOERR) then
1916	write(iunit, '(a,a,a)') 'comment="', trim(attr_str), '"'
1917	end if
1918	end if
1919	end if
1920	else
1921	! No long_name present
1922	write(iunit, '(1x)')
1923	end if
1924
1925	end subroutine print_variable_attributes
1926
1927
1928	! --- WRITING SUBROUTINES ---
1929
1930	!---------------------------------------------------------------------
1931	! Define a dimension with name dim_name of length n (or 0 to
1932	! indicate the unlimited dimension)
1933	subroutine define_dimension(this, dim_name, n)
1934	class(netcdf_file) :: this
1935	character(len=*), intent(in) :: dim_name
1936	integer, intent(in) :: n
1937	integer :: idimid, istatus
1938
1939	istatus = nf90_def_dim(this%ncid, dim_name, n, idimid)
1940	if (istatus /= NF90_NOERR) then
1941	write(nulerr,'(a,a,a,a)') '*** Error defining ', dim_name, &
1942	& ' as a dimension: ', trim(nf90_strerror(istatus))
1943	call my_abort('Error writing NetCDF file')
1944	end if
1945
1946	if (this%iverbose >= 4) then
1947	write(nulout,'(a,a,a,i0)') ' Defining dimension ',trim(dim_name), &
1948	& ' of length ', n
1949	end if
1950
1951	end subroutine define_dimension
1952
1953
1954	!---------------------------------------------------------------------
1955	! Define a variable with name var_name, and specify its shape via
1956	! the dim1_name, dim2_name and dim3_name optional arguments, which
1957	! are strings referring to existing dimensions. If none are
1958	! specified, the variable will be a scalar. The optional arguments
1959	! long_name, units and comment write string attributes with these
1960	! names.
1961	subroutine define_variable(this, var_name, dim1_name, dim2_name, dim3_name, &
1962	& dim4_name, long_name, units_str, comment_str, &
1963	& standard_name, is_double, data_type_name, fill_value, &
1964	& deflate_level, shuffle, chunksizes, ndims)
1965	class(netcdf_file) :: this
1966	character(len=*), intent(in) :: var_name
1967	character(len=*), intent(in), optional :: long_name, units_str, comment_str, standard_name
1968	character(len=*), intent(in), optional :: dim1_name, dim2_name, dim3_name, dim4_name
1969	logical, intent(in), optional :: is_double
1970	character(len=*), intent(in), optional :: data_type_name
1971	real(jprb), intent(in), optional :: fill_value
1972	integer, intent(in), optional :: deflate_level ! Compression: 0 (none) to 9 (most)
1973	logical, intent(in), optional :: shuffle ! Shuffle bytes before compression
1974	integer, dimension(:), intent(in), optional :: chunksizes
1975	integer, intent(in), optional :: ndims
1976
1977	integer :: istatus, ndims_local, ndims_input, ivarid
1978	integer, dimension(NF90_MAX_VAR_DIMS) :: idimids
1979	integer :: data_type
1980
1981	! Sometimes a program may not know at compile time the exact
1982	! dimensions of a variable - if ndims is present then only up to
1983	! that many dimensions will be defined
1984	ndims_input = 4
1985	if (present(ndims)) then
1986	ndims_input = ndims
1987	end if
1988
1989	if (present(dim1_name) .and. ndims_input >= 1) then
1990	! Variable is at least one dimensional
1991	ndims_local = 1
1992	istatus = nf90_inq_dimid(this%ncid, dim1_name, idimids(1))
1993	if (istatus /= NF90_NOERR) then
1994	write(nulerr,'(a,a,a,a)') '*** Error inquiring ID of dimension ', &
1995	& dim1_name, ': ', trim(nf90_strerror(istatus))
1996	call my_abort('Error writing NetCDF file')
1997	end if
1998	if (present(dim2_name) .and. ndims_input >= 2) then
1999	! Variable is at least two dimensional
2000	ndims_local = 2
2001	istatus = nf90_inq_dimid(this%ncid, dim2_name, idimids(2))
2002	if (istatus /= NF90_NOERR) then
2003	write(nulerr,'(a,a,a)') '*** Error inquiring ID of dimension ', &
2004	& dim2_name, ': ', trim(nf90_strerror(istatus))
2005	call my_abort('Error writing NetCDF file')
2006	end if
2007	if (present(dim3_name) .and. ndims_input >= 3) then
2008	! Variable is at least three dimensional
2009	ndims_local = 3
2010	istatus = nf90_inq_dimid(this%ncid, dim3_name, idimids(3))
2011	if (istatus /= NF90_NOERR) then
2012	write(nulerr,'(a,a,a,a)') '*** Error inquiring ID of dimension ', &
2013	& dim3_name, ': ', trim(nf90_strerror(istatus))
2014	call my_abort('Error writing NetCDF file')
2015	end if
2016	if (present(dim4_name) .and. ndims_input >= 4) then
2017	! Variable is at least three dimensional
2018	ndims_local = 4
2019	istatus = nf90_inq_dimid(this%ncid, dim4_name, idimids(4))
2020	if (istatus /= NF90_NOERR) then
2021	write(nulerr,'(a,a,a,a)') '*** Error inquiring ID of dimension ', &
2022	& dim4_name, ': ', trim(nf90_strerror(istatus))
2023	call my_abort('Error writing NetCDF file')
2024	end if
2025	end if
2026	end if
2027	end if
2028	else
2029	! Variable is a scalar
2030	ndims_local = 0
2031	end if
2032
2033	! Read output precision from optional argument "is_double" if
2034	! present, otherwise from default output precision for this file
2035	if (present(data_type_name)) then
2036	if (data_type_name == 'double') then
2037	data_type = NF90_DOUBLE
2038	else if (data_type_name == 'byte') then
2039	data_type = NF90_BYTE
2040	else if (data_type_name == 'short') then
2041	data_type = NF90_SHORT
2042	else if (data_type_name == 'int') then
2043	data_type = NF90_INT
2044	else if (data_type_name == 'float') then
2045	data_type = NF90_FLOAT
2046	else
2047	write(nulerr,'(a,a,a)') '*** NetCDF data type "', data_type_name, '" not supported'
2048	call my_abort('Error writing NetCDF file')
2049	end if
2050	else if (present(is_double)) then
2051	if (is_double) then
2052	data_type = NF90_DOUBLE
2053	else
2054	data_type = NF90_FLOAT
2055	end if
2056	else if (this%is_double_precision) then
2057	data_type = NF90_DOUBLE
2058	else
2059	data_type = NF90_FLOAT
2060	end if
2061
2062	! Define variable
2063	#ifdef NC_NETCDF4
2064	istatus = nf90_def_var(this%ncid, var_name, data_type, idimids(1:ndims_local), &
2065	& ivarid, deflate_level=deflate_level, shuffle=shuffle, chunksizes=chunksizes)
2066	#else
2067	istatus = nf90_def_var(this%ncid, var_name, data_type, idimids(1:ndims_local), ivarid)
2068	#endif
2069	if (istatus /= NF90_NOERR) then
2070	write(nulerr,'(a,a,a,a)') '*** Error defining variable ', var_name, &
2071	& ': ', trim(nf90_strerror(istatus))
2072	call my_abort('Error writing NetCDF file')
2073	end if
2074
2075	! Add optional attributes
2076	if (present(long_name)) then
2077	istatus = nf90_put_att(this%ncid, ivarid, "long_name", long_name)
2078	if (this%iverbose >= 4) then
2079	write(nulout,'(a,a,a,a,a)', advance='no') ' Defining ',trim(var_name), &
2080	& ': "', long_name, '"'
2081	end if
2082	else
2083	if (this%iverbose >= 4) then
2084	write(nulout,'(a,a)', advance='no') ' Defining ',trim(var_name)
2085	end if
2086	end if
2087
2088	if (present(units_str)) then
2089	istatus = nf90_put_att(this%ncid, ivarid, "units", units_str)
2090	if (this%iverbose >= 4) then
2091	if (trim(units_str) == '1') then
2092	write(nulout, '(a)') ' (dimensionless)'
2093	else
2094	write(nulout, '(a,a,a)') ' (', trim(units_str), ')'
2095	end if
2096	end if
2097	else
2098	if (this%iverbose >= 4) then
2099	write(nulout, '(1x)')
2100	end if
2101	end if
2102
2103	if (present(standard_name)) then
2104	istatus = nf90_put_att(this%ncid, ivarid, "standard_name", standard_name)
2105	end if
2106	if (present(comment_str)) then
2107	istatus = nf90_put_att(this%ncid, ivarid, "comment", comment_str)
2108	end if
2109
2110	if (present(fill_value)) then
2111	#ifdef NC_NETCDF4
2112	if (data_type == NF90_DOUBLE) then
2113	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, real(fill_value,jprd))
2114	else if (data_type == NF90_FLOAT) then
2115	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, real(fill_value,jprm))
2116	else if (data_type == NF90_INT) then
2117	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, int(fill_value,4))
2118	else if (data_type == NF90_SHORT) then
2119	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, int(fill_value,2))
2120	else if (data_type == NF90_BYTE) then
2121	istatus = nf90_def_var_fill(this%ncid, ivarid, 0, int(fill_value,1))
2122	end if
2123	#else
2124	if (data_type == NF90_DOUBLE) then
2125	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", real(fill_value,jprd))
2126	else if (data_type == NF90_FLOAT) then
2127	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", real(fill_value,jprm))
2128	else if (data_type == NF90_INT) then
2129	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", int(fill_value,4))
2130	else if (data_type == NF90_SHORT) then
2131	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", int(fill_value,2))
2132	else if (data_type == NF90_BYTE) then
2133	istatus = nf90_put_att(this%ncid, ivarid, "_FillValue", int(fill_value,1))
2134	end if
2135	#endif
2136	end if
2137
2138	end subroutine define_variable
2139
2140
2141	!---------------------------------------------------------------------
2142	! Put CF-compliant global attributes into the file
2143	subroutine put_global_attributes(this, title_str, inst_str, source_str, &
2144	& comment_str, references_str, creator_name, creator_email_str, &
2145	& contributor_name, project_str, conventions_str, prior_history_str)
2146	class(netcdf_file) :: this
2147
2148	character(len=*), intent(in), optional :: title_str
2149	character(len=*), intent(in), optional :: inst_str
2150	character(len=*), intent(in), optional :: source_str
2151	character(len=*), intent(in), optional :: creator_name, creator_email_str
2152	character(len=*), intent(in), optional :: contributor_name, project_str
2153	character(len=*), intent(in), optional :: comment_str, conventions_str
2154	character(len=*), intent(in), optional :: references_str, prior_history_str
2155
2156	character(len=32) :: date_time_str
2157	character(len=4000) :: command_line_str
2158	character(len=4032) :: history_str
2159
2160	integer :: time_vals(8)
2161	integer :: i ! status
2162
2163	call date_and_time(values=time_vals)
2164	call get_command(command_line_str)
2165
2166	write(date_time_str,"(i0.4,'-',i0.2,'-',i0.2,' ',i0.2,':',i0.2,':',i0.2)") &
2167	& time_vals(1), time_vals(2), time_vals(3), time_vals(5), time_vals(6), time_vals(7)
2168
2169	if (present(prior_history_str)) then
2170	history_str = trim(prior_history_str) // new_line('a') &
2171	& // trim(date_time_str) // ': ' // trim(command_line_str)
2172	else
2173	history_str = trim(date_time_str) // ': ' // trim(command_line_str)
2174	end if
2175
2176	if (present(title_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "title", title_str)
2177	if (present(inst_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "institution", inst_str)
2178	if (present(source_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "source", source_str)
2179	if (present(creator_name))i=nf90_put_att(this%ncid, NF90_GLOBAL, "creator_name", creator_name)
2180	if (present(creator_email_str))i=nf90_put_att(this%ncid, NF90_GLOBAL, "creator_email", creator_email_str)
2181	if (present(contributor_name))i=nf90_put_att(this%ncid, NF90_GLOBAL, "contributor_name", contributor_name)
2182
2183	i = nf90_put_att(this%ncid, NF90_GLOBAL, "history", history_str)
2184
2185	if (present(project_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "project", project_str)
2186	if (present(comment_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, "comment", comment_str)
2187	if (present(references_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, &
2188	& "references", references_str)
2189	if (present(conventions_str)) i=nf90_put_att(this%ncid, NF90_GLOBAL, &
2190	& "conventions", conventions_str)
2191
2192	end subroutine put_global_attributes
2193
2194
2195	!---------------------------------------------------------------------
2196	! Put a non-standard global attribute into the file
2197	subroutine put_global_attribute(this, attr_name, attr_str)
2198	class(netcdf_file) :: this
2199
2200	character(len=*), intent(in) :: attr_name, attr_str
2201
2202	integer :: istatus
2203
2204	istatus = nf90_put_att(this%ncid, NF90_GLOBAL, trim(attr_name), trim(attr_str))
2205
2206	if (istatus /= NF90_NOERR) then
2207	write(nulerr,'(a,a,a,a)') '*** Error writing global attribute ', attr_name, &
2208	& ': ', trim(nf90_strerror(istatus))
2209	call my_abort('Error writing NetCDF file')
2210	end if
2211
2212	end subroutine put_global_attribute
2213
2214
2215	!---------------------------------------------------------------------
2216	! Put a non-standard variable attribute into the file
2217	subroutine put_attribute(this, var_name, attr_name, attr_str)
2218	class(netcdf_file) :: this
2219
2220	character(len=*), intent(in) :: var_name, attr_name, attr_str
2221
2222	integer :: istatus, ivarid
2223
2224	call this%get_variable_id(var_name, ivarid)
2225
2226	istatus = nf90_put_att(this%ncid, ivarid, trim(attr_name), trim(attr_str))
2227
2228	if (istatus /= NF90_NOERR) then
2229	write(nulerr,'(a,a,a,a)') '*** Error writing attribute ', attr_name, &
2230	& ': ', trim(nf90_strerror(istatus))
2231	call my_abort('Error writing NetCDF file')
2232	end if
2233
2234	end subroutine put_attribute
2235
2236
2237	!---------------------------------------------------------------------
2238	! The "put" method saves a scalar, vector or matrix into the
2239	! variable with name var_name, according to the rank of the var
2240	! argument. This version saves a scalar.
2241	subroutine put_real_scalar(this, var_name, var)
2242	class(netcdf_file) :: this
2243	character(len=*), intent(in) :: var_name
2244	real(jprb), intent(in) :: var
2245
2246	integer :: ivarid, ndims, istatus
2247	integer(kind=jpib) :: ntotal
2248	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2249
2250	! If we are in define mode, exit define mode
2251	call this%end_define_mode()
2252
2253	! Check the variable is a scalar
2254	call this%get_variable_id(var_name, ivarid)
2255	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2256	if (ntotal /= 1) then
2257	write(nulerr,'(a,a,a,i0)') '*** Error: attempt to write scalar to ', &
2258	& var_name, ' which has total length ', ntotal
2259	call my_abort('Error writing NetCDF file')
2260	end if
2261
2262	! Save the scalar
2263	istatus = nf90_put_var(this%ncid, ivarid, var)
2264	if (istatus /= NF90_NOERR) then
2265	write(nulerr,'(a,a,a,a)') '*** Error writing scalar ', var_name, ': ', &
2266	& trim(nf90_strerror(istatus))
2267	call my_abort('Error writing NetCDF file')
2268	end if
2269
2270	end subroutine put_real_scalar
2271
2272
2273	!---------------------------------------------------------------------
2274	! Save a scalar.
2275	subroutine put_real_scalar_indexed(this, var_name, index, var)
2276	class(netcdf_file) :: this
2277	character(len=*), intent(in) :: var_name
2278	real(jprb), intent(in) :: var
2279	integer, intent(in) :: index
2280
2281	integer :: ivarid, ndims, istatus
2282	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2283	integer :: vstart(NF90_MAX_VAR_DIMS)
2284
2285	! If we are in define mode, exit define mode
2286	call this%end_define_mode()
2287
2288	! Check the variable is a scalar
2289	call this%get_variable_id(var_name, ivarid)
2290	call this%get_array_dimensions(ivarid, ndims, ndimlens)
2291	if (index < 1 .or. index > ndimlens(ndims)) then
2292	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write scalar to element ', &
2293	& index, ' of ', var_name, ' which has outer dimension ', ndimlens(ndims)
2294	call my_abort('Error writing NetCDF file')
2295	end if
2296
2297	! Save the scalar
2298	vstart(1:ndims-1) = 1
2299	vstart(ndims) = index
2300	istatus = nf90_put_var(this%ncid, ivarid, var, start=vstart)
2301	if (istatus /= NF90_NOERR) then
2302	write(nulerr,'(a,a,a,a)') '*** Error writing scalar to ', var_name, ': ', &
2303	& trim(nf90_strerror(istatus))
2304	call my_abort('Error writing NetCDF file')
2305	end if
2306
2307	end subroutine put_real_scalar_indexed
2308
2309
2310	!---------------------------------------------------------------------
2311	! Save a vector with name var_name in the file
2312	subroutine put_real_vector(this, var_name, var)
2313	class(netcdf_file) :: this
2314	character(len=*), intent(in) :: var_name
2315	real(jprb), intent(in) :: var(:)
2316
2317	integer :: ivarid, ndims, istatus
2318	integer(kind=jpib) :: ntotal
2319	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2320
2321	call this%end_define_mode()
2322
2323	! Check the vector is of the right length
2324	call this%get_variable_id(var_name, ivarid)
2325	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2326	if (ntotal /= size(var,kind=jpib)) then
2327	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write vector of length ', &
2328	& size(var), ' to ', var_name, ' which has total length ', ntotal
2329	call my_abort('Error writing NetCDF file')
2330	end if
2331
2332	! Save the vector
2333	istatus = nf90_put_var(this%ncid, ivarid, var)
2334	if (istatus /= NF90_NOERR) then
2335	write(nulerr,'(a,a,a,a)') '*** Error writing vector ', var_name, ': ', &
2336	& trim(nf90_strerror(istatus))
2337	call my_abort('Error writing NetCDF file')
2338	end if
2339
2340	end subroutine put_real_vector
2341
2342
2343	!---------------------------------------------------------------------
2344	! Save an integer vector with name var_name in the file
2345	subroutine put_int_vector(this, var_name, var)
2346	class(netcdf_file) :: this
2347	character(len=*), intent(in) :: var_name
2348	integer, intent(in) :: var(:)
2349
2350	integer :: ivarid, ndims, istatus
2351	integer(kind=jpib) :: ntotal
2352	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2353
2354	call this%end_define_mode()
2355
2356	! Check the vector is of the right length
2357	call this%get_variable_id(var_name, ivarid)
2358	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2359	if (ntotal /= size(var,kind=jpib)) then
2360	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write vector of length ', &
2361	& size(var), ' to ', var_name, ' which has total length ', ntotal
2362	call my_abort('Error writing NetCDF file')
2363	end if
2364
2365	! Save the vector
2366	istatus = nf90_put_var(this%ncid, ivarid, var)
2367	if (istatus /= NF90_NOERR) then
2368	write(nulerr,'(a,a,a,a)') '*** Error writing vector ', var_name, ': ', &
2369	& trim(nf90_strerror(istatus))
2370	call my_abort('Error writing NetCDF file')
2371	end if
2372
2373	end subroutine put_int_vector
2374
2375
2376	!---------------------------------------------------------------------
2377	! Save a vector slice with name var_name in the file
2378	subroutine put_real_vector_indexed(this, var_name, var, index2, index3)
2379	class(netcdf_file) :: this
2380	character(len=*), intent(in) :: var_name
2381	real(jprb), intent(in) :: var(:)
2382	integer, intent(in) :: index2
2383	integer, intent(in), optional :: index3
2384
2385	integer :: ivarid, ndims, istatus
2386	integer(kind=jpib) :: ntotal
2387	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2388	integer :: vstart(NF90_MAX_VAR_DIMS)
2389	integer :: vcount(NF90_MAX_VAR_DIMS)
2390
2391	character(len=512) :: var_slice_name
2392	integer :: index_last
2393
2394	call this%end_define_mode()
2395
2396	! Check the vector is of the right length
2397	call this%get_variable_id(var_name, ivarid)
2398	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2399	ntotal = ntotal / ndimlens(ndims)
2400	if (present(index3)) then
2401	ntotal = ntotal / ndimlens(ndims-1)
2402	index_last = index3
2403	write(var_slice_name,'(a,a,i0,a,i0,a)') var_name, '(:,', index2, ',', index3, ')'
2404	else
2405	index_last = index2
2406	write(var_slice_name,'(a,a,i0,a)') var_name, '(:,', index2, ')'
2407	end if
2408
2409	if (ntotal /= size(var,kind=jpib)) then
2410	write(nulerr,'(a,i0,a,a,i0)') '*** Error: attempt to write vector of length ', &
2411	& size(var), ' to ', trim(var_slice_name), ' which has length ', ntotal
2412	call my_abort('Error writing NetCDF file')
2413	end if
2414	if (index_last < 1 .or. index_last > ndimlens(ndims)) then
2415	write(nulerr,'(a,a,a,i0)') '*** Error: attempt to write vector to ', &
2416	& trim(var_slice_name), ' which has outer dimension ', ndimlens(ndims)
2417	call my_abort('Error writing NetCDF file')
2418	end if
2419
2420	! Save the vector
2421	vstart(1:ndims-1) = 1
2422	vcount(1:ndims-1) = ndimlens(1:ndims-1)
2423	vcount(ndims) = 1
2424	if (present(index3)) then
2425	vstart(ndims) = index3
2426	vstart(ndims-1) = index2
2427	vcount(ndims-1) = 1
2428	else
2429	vstart(ndims) = index2
2430	end if
2431
2432	istatus = nf90_put_var(this%ncid, ivarid, var, start=vstart, count=vcount)
2433	if (istatus /= NF90_NOERR) then
2434	write(nulerr,'(a,a,a,a)') '*** Error writing vector to ', trim(var_slice_name), &
2435	& ': ', trim(nf90_strerror(istatus))
2436	call my_abort('Error writing NetCDF file')
2437	end if
2438
2439	end subroutine put_real_vector_indexed
2440
2441
2442	!---------------------------------------------------------------------
2443	! Save a matrix with name var_name in the file, transposing its
2444	! dimensions if either optional argument transp is .true., or the
2445	! transpose_matrices method has already been called.
2446	subroutine put_real_matrix(this, var_name, var, do_transp)
2447	class(netcdf_file) :: this
2448	character(len=*), intent(in) :: var_name
2449	real(jprb), intent(in) :: var(:,:)
2450	real(jprb), allocatable :: var_transpose(:,:)
2451	logical, optional, intent(in):: do_transp
2452
2453	integer :: ivarid, ndims, nvarlen, istatus
2454	integer(kind=jpib) :: ntotal
2455	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2456
2457	logical :: do_transpose
2458
2459	if (present(do_transp)) then
2460	do_transpose = do_transp
2461	else
2462	do_transpose = this%do_transpose_2d
2463	end if
2464
2465	call this%end_define_mode()
2466
2467	call this%get_variable_id(var_name, ivarid)
2468	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2469
2470	nvarlen = size(var,1)*size(var,2)
2471
2472	! Check the total size of the variable to be stored (but receiving
2473	! ntotal is zero then there must be an unlimited dimension)
2474	if (ntotal /= size(var,kind=jpib) .and. ntotal /= 0) then
2475	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write matrix of total size ', &
2476	& nvarlen, ' to ', var_name, ' which has total size ', ntotal
2477	call my_abort('Error writing NetCDF file')
2478	end if
2479
2480	if (do_transpose) then
2481	! Save the matrix with transposition
2482	if (this%iverbose >= 3) then
2483	write(nulout,'(a,a,a)') ' Writing ', var_name, &
2484	& ' (transposing dimensions)'
2485	end if
2486	allocate(var_transpose(size(var,2), size(var,1)))
2487	var_transpose = transpose(var)
2488	istatus = nf90_put_var(this%ncid, ivarid, var_transpose)
2489	deallocate(var_transpose)
2490	else
2491	! Save the matrix without transposition
2492	if (this%iverbose >= 3) then
2493	write(nulout,'(a,a)') ' Writing ', var_name
2494	end if
2495	istatus = nf90_put_var(this%ncid, ivarid, var)
2496	end if
2497
2498	if (istatus /= NF90_NOERR) then
2499	write(nulerr,'(a,a,a,a)') '*** Error writing matrix ', var_name, &
2500	& ': ', trim(nf90_strerror(istatus))
2501	call my_abort('Error writing NetCDF file')
2502	end if
2503
2504	end subroutine put_real_matrix
2505
2506
2507	!---------------------------------------------------------------------
2508	! Save a matrix slice with name var_name in the file, transposing its
2509	! dimensions if either optional argument transp is .true., or the
2510	! transpose_matrices method has already been called.
2511	subroutine put_real_matrix_indexed(this, var_name, var, index3, index4, do_transp)
2512	class(netcdf_file) :: this
2513	character(len=*), intent(in) :: var_name
2514	real(jprb), intent(in) :: var(:,:)
2515	integer, intent(in) :: index3
2516	integer, intent(in), optional :: index4
2517
2518	real(jprb), allocatable :: var_transpose(:,:)
2519	logical, optional, intent(in) :: do_transp
2520
2521	integer :: ivarid, ndims, nvarlen, istatus
2522	integer(kind=jpib) :: ntotal
2523	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2524	integer :: vstart(NF90_MAX_VAR_DIMS)
2525	integer :: vcount(NF90_MAX_VAR_DIMS)
2526
2527	character(len=512) :: var_slice_name
2528
2529	logical :: do_transpose
2530
2531	if (present(do_transp)) then
2532	do_transpose = do_transp
2533	else
2534	do_transpose = this%do_transpose_2d
2535	end if
2536
2537	call this%end_define_mode()
2538
2539	call this%get_variable_id(var_name, ivarid)
2540	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2541
2542	nvarlen = size(var,1)*size(var,2)
2543
2544	! Check the total size of the variable to be stored (but receiving
2545	! ntotal is zero then there must be an unlimited dimension)
2546	ntotal = ntotal / ndimlens(ndims)
2547	if (present(index4)) then
2548	ntotal = ntotal / ndimlens(ndims-1)
2549	write(var_slice_name,'(a,a,i0,a,i0,a)') var_name, '(:,:,', index3, ',', index4, ')'
2550	else
2551	write(var_slice_name,'(a,a,i0,a)') var_name, '(:,:,', index3, ')'
2552	end if
2553	if (ntotal /= size(var,kind=jpib) .and. ntotal /= 0) then
2554	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write matrix of total size ', &
2555	& nvarlen, ' to ', trim(var_slice_name), ' which has total size ', ntotal
2556	call my_abort('Error writing NetCDF file')
2557	end if
2558
2559	vstart(1:ndims-1) = 1
2560	vcount(1:ndims-1) = ndimlens(1:ndims-1)
2561	vcount(ndims) = 1
2562	if (present(index4)) then
2563	vstart(ndims) = index4
2564	vstart(ndims-1) = index3
2565	vcount(ndims-1) = 1
2566	else
2567	vstart(ndims) = index3
2568	end if
2569
2570	if (do_transpose) then
2571	! Save the matrix with transposition
2572	if (this%iverbose >= 3) then
2573	write(nulout,'(a,a,a)') ' Writing ', trim(var_slice_name), &
2574	& ' (transposing dimensions)'
2575	end if
2576	allocate(var_transpose(size(var,2), size(var,1)))
2577	var_transpose = transpose(var)
2578	istatus = nf90_put_var(this%ncid, ivarid, var_transpose, start=vstart, count=vcount)
2579	deallocate(var_transpose)
2580	else
2581	! Save the matrix without transposition
2582	if (this%iverbose >= 3) then
2583	write(nulout,'(a,a)') ' Writing ', trim(var_slice_name)
2584	end if
2585	istatus = nf90_put_var(this%ncid, ivarid, var, start=vstart, count=vcount)
2586	end if
2587
2588	if (istatus /= NF90_NOERR) then
2589	write(nulerr,'(a,a,a)') '*** Error writing ', trim(var_slice_name), &
2590	& ': ', trim(nf90_strerror(istatus))
2591	call my_abort('Error writing NetCDF file')
2592	end if
2593
2594	end subroutine put_real_matrix_indexed
2595
2596
2597	!---------------------------------------------------------------------
2598	! Save a 3D array with name var_name in the file. The optional
2599	! argument permute specifies that the dimensions should first be
2600	! permuted according to the three integers therein (or if
2601	! permute_3d_arrays has already been called). ipermute is
2602	! interpretted such that if OLD and NEW are 3-element vectors
2603	! containing the size of each dimension in memory and in the written
2604	! file, respectively, then NEW=OLD(ipermute).
2605	subroutine put_real_array3(this, var_name, var, ipermute)
2606	class(netcdf_file) :: this
2607	character(len=*), intent(in) :: var_name
2608	real(jprb), intent(in) :: var(:,:,:)
2609	real(jprb), allocatable :: var_permute(:,:,:)
2610	integer, optional, intent(in) :: ipermute(3)
2611
2612	integer :: ivarid, ndims, nvarlen, istatus
2613	integer(kind=jpib) :: ntotal
2614	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2615
2616	logical :: do_permute ! Do we permute?
2617	integer :: i_permute_3d(3)
2618	integer :: n_dimlens_permuted(3)
2619	integer :: i_order(3)
2620
2621	! Decide whether to permute
2622	if (present(ipermute)) then
2623	do_permute = .true.
2624	i_permute_3d = ipermute
2625	else
2626	do_permute = this%do_permute_3d
2627	i_permute_3d = this%i_permute_3d
2628	end if
2629
2630	call this%end_define_mode()
2631
2632	! Check total size
2633	call this%get_variable_id(var_name, ivarid)
2634	call this%get_array_dimensions(ivarid, ndims, ndimlens, ntotal)
2635	nvarlen = size(var,1)size(var,2)size(var,3)
2636	if (ntotal /= size(var,kind=jpib)) then
2637	write(nulerr,'(a,i0,a,a,a,i0)') '*** Error: attempt to write array of total size ', &
2638	& nvarlen, ' to ', var_name, ' which has total size ', ntotal
2639	call my_abort('Error writing NetCDF file')
2640	end if
2641
2642	if (do_permute) then
2643	! Save array after permuting dimensions
2644	if (this%iverbose >= 3) then
2645	write(nulout,'(a,a,a,i0,i0,i0,a)') ' Writing ', var_name, &
2646	& ' (permuting dimensions: ', i_permute_3d, ')'
2647	end if
2648	n_dimlens_permuted = (/ size(var,i_permute_3d(1)), &
2649	& size(var,i_permute_3d(2)), &
2650	& size(var,i_permute_3d(3)) /)
2651	!! FIX: This makes it look like the dimensions have stayed the same
2652	! if (this%iverbose >= 4) then
2653	! write(nulout,'(a,i0,a,i0,a,i0,a,i0,a,i0,a,i0,a)') ' (', &
2654	! & n_dimlens_permuted(1), ',', n_dimlens_permuted(2), &
2655	! & ',', n_dimlens_permuted(3), ') -> (', ndimlens(1), &
2656	! & ',', ndimlens(2), ',', ndimlens(3), ')'
2657	! end if
2658	allocate(var_permute(n_dimlens_permuted(1), &
2659	& n_dimlens_permuted(2), n_dimlens_permuted(3)))
2660	! Due to the odd way that ORDER works in Fortran RESHAPE, we
2661	! need to do this:
2662	i_order(i_permute_3d) = (/ 1, 2, 3 /)
2663	var_permute = reshape(var, n_dimlens_permuted, order=i_order)
2664	istatus = nf90_put_var(this%ncid, ivarid, var_permute)
2665	deallocate(var_permute)
2666	else
2667	! Save array without permuting dimensions
2668	if (this%iverbose >= 3) then
2669	write(nulout,'(a,a)') ' Writing ', var_name
2670	end if
2671	istatus = nf90_put_var(this%ncid, ivarid, var)
2672	end if
2673
2674	if (istatus /= NF90_NOERR) then
2675	write(nulerr,'(a,a,a,a)') '*** Error writing array ', var_name, &
2676	& ': ', trim(nf90_strerror(istatus))
2677	call my_abort('Error writing NetCDF file')
2678	end if
2679
2680	end subroutine put_real_array3
2681
2682
2683	#ifdef NC_NETCDF4
2684	!---------------------------------------------------------------------
2685	! Copy dimensions from "infile" to "this"
2686	subroutine copy_dimensions(this, infile)
2687	class(netcdf_file) :: this
2688	type(netcdf_file), intent(in) :: infile
2689
2690	integer :: jdim
2691	integer :: ndims
2692	integer :: idimids(1024)
2693	integer :: dimlen
2694	character(len=512) :: dimname
2695	integer :: istatus
2696	integer :: include_parents
2697
2698	include_parents = 0
2699
2700	istatus = nf90_inq_dimids(infile%ncid, ndims, idimids, include_parents)
2701	if (istatus /= NF90_NOERR) then
2702	write(nulerr,'(a,a)') '*** Error reading dimensions of NetCDF file: ', &
2703	trim(nf90_strerror(istatus))
2704	call my_abort('Error reading NetCDF file')
2705	end if
2706
2707	DO jdim = 1,ndims
2708	istatus = nf90_inquire_dimension(infile%ncid, idimids(jdim), &
2709	& name=dimname, len=dimlen)
2710	if (istatus /= NF90_NOERR) then
2711	write(nulerr,'(a,a)') '*** Error reading NetCDF dimension properties: ', &
2712	trim(nf90_strerror(istatus))
2713	call my_abort('Error reading NetCDF file')
2714	end if
2715	call this%define_dimension(trim(dimname), dimlen)
2716	end do
2717
2718	end subroutine copy_dimensions
2719	#endif
2720
2721	!---------------------------------------------------------------------
2722	! Copy variable definition and attributes from "infile" to "this"
2723	subroutine copy_variable_definition(this, infile, var_name)
2724	class(netcdf_file) :: this
2725	type(netcdf_file), intent(in) :: infile
2726	character(len=*), intent(in) :: var_name
2727
2728	#ifdef NC_NETCDF4
2729	integer :: deflate_level ! Compression: 0 (none) to 9 (most)
2730	logical :: shuffle ! Shuffle bytes before compression
2731	integer :: chunksizes(NF90_MAX_VAR_DIMS)
2732	#endif
2733	integer :: data_type
2734	integer :: ndims
2735	integer :: idimids_in(NF90_MAX_VAR_DIMS)
2736	integer :: idimids_out(NF90_MAX_VAR_DIMS)
2737	integer :: nattr
2738	character(len=512) :: attr_name
2739	character(len=512) :: dim_name
2740
2741	integer :: istatus
2742	integer :: ivarid_in, ivarid_out
2743	integer :: jattr, jdim
2744
2745	if (this%iverbose >= 4) then
2746	write(nulout,'(a,a)') ' Copying definition of ', trim(var_name)
2747	end if
2748
2749	! Get variable ID from name
2750	istatus = nf90_inq_varid(infile%ncid, var_name, ivarid_in)
2751	if (istatus /= NF90_NOERR) then
2752	write(nulerr,'(a,i0,a)') '*** Error inquiring about NetCDF variable "', &
2753	& var_name, '": ', trim(nf90_strerror(istatus))
2754	call my_abort('Error reading NetCDF file')
2755	end if
2756
2757	! Get variable properties
2758	#ifdef NC_NETCDF4
2759	istatus = nf90_inquire_variable(infile%ncid, ivarid_in, xtype=data_type, ndims=ndims, &
2760	& dimids=idimids_in, chunksizes=chunksizes, deflate_level=deflate_level, &
2761	& shuffle=shuffle, natts=nattr)
2762	#else
2763	istatus = nf90_inquire_variable(infile%ncid, ivarid_in, xtype=data_type, ndims=ndims, &
2764	& dimids=idimids_in, natts=nattr)
2765	#endif
2766	if (istatus /= NF90_NOERR) then
2767	write(nulerr,'(a,a)') '*** Error reading NetCDF variable properties: ', &
2768	trim(nf90_strerror(istatus))
2769	call my_abort('Error reading NetCDF file')
2770	end if
2771
2772	! Map dimension IDs
2773	DO jdim = 1,ndims
2774	istatus = nf90_inquire_dimension(infile%ncid, idimids_in(jdim), name=dim_name)
2775	if (istatus /= NF90_NOERR) then
2776	write(nulerr,'(a,a)') '*** Error reading NetCDF dimension name: ', &
2777	trim(nf90_strerror(istatus))
2778	call my_abort('Error reading NetCDF file')
2779	end if
2780
2781	istatus = nf90_inq_dimid(this%ncid, trim(dim_name), idimids_out(jdim))
2782	if (istatus /= NF90_NOERR) then
2783	write(nulerr,'(a,a)') '*** Error reading NetCDF dimension ID: ', &
2784	trim(nf90_strerror(istatus))
2785	call my_abort('Error reading NetCDF file')
2786	end if
2787	end do
2788
2789	! Create variable
2790	#ifdef NC_NETCDF4
2791	istatus = nf90_def_var(this%ncid, var_name, data_type, idimids_out(1:ndims), &
2792	& ivarid_out, deflate_level=deflate_level, shuffle=shuffle, chunksizes=chunksizes(1:ndims))
2793	#else
2794	istatus = nf90_def_var(this%ncid, var_name, data_type, idimids_out(1:ndims), ivarid_out)
2795	#endif
2796	if (istatus /= NF90_NOERR) then
2797	write(nulerr,'(a,a,a,a)') '*** Error defining variable "', var_name, &
2798	& '": ', trim(nf90_strerror(istatus))
2799	call my_abort('Error writing NetCDF file')
2800	end if
2801
2802	! Copy attributes
2803	DO jattr = 1,nattr
2804	istatus = nf90_inq_attname(infile%ncid, ivarid_in, jattr, attr_name)
2805	if (istatus /= NF90_NOERR) then
2806	write(nulerr,'(a,a)') '*** Error reading attribute: ', &
2807	& trim(nf90_strerror(istatus))
2808	call my_abort('Error reading NetCDF file')
2809	end if
2810	istatus = nf90_copy_att(infile%ncid, ivarid_in, trim(attr_name), &
2811	& this%ncid, ivarid_out)
2812
2813	end do
2814
2815	end subroutine copy_variable_definition
2816
2817
2818	!---------------------------------------------------------------------
2819	! Copy variable from "infile" to "this"
2820	subroutine copy_variable(this, infile, var_name)
2821	class(netcdf_file) :: this
2822	class(netcdf_file), intent(in) :: infile
2823	character(len=*), intent(in) :: var_name
2824
2825	integer :: ivarid_in, ivarid_out
2826	integer :: ndims
2827	integer :: ndimlens(NF90_MAX_VAR_DIMS)
2828	integer(kind=jpib) :: ntotal
2829	integer :: data_type
2830	integer :: istatus
2831
2832	! We use the Fortran-77 functions because they don't check that
2833	! the rank of the arguments is correct
2834	integer, external :: nf_get_var_double, nf_put_var_double
2835	integer, external :: nf_get_var_int, nf_put_var_int
2836
2837	real(kind=jprd), allocatable :: data_real(:)
2838	integer, allocatable :: data_int(:)
2839
2840	! If we are in define mode, exit define mode
2841	call this%end_define_mode()
2842
2843	if (this%iverbose >= 4) then
2844	write(nulout,'(a,a)') ' Copying ', trim(var_name)
2845	end if
2846
2847	call infile%get_variable_id(var_name, ivarid_in)
2848	call infile%get_array_dimensions(ivarid_in, ndims, ndimlens, ntotal)
2849	istatus = nf90_inquire_variable(infile%ncid, ivarid_in, xtype=data_type)
2850	if (istatus /= NF90_NOERR) then
2851	write(nulerr,'(a,a,a,a)') '*** Error reading variable "', var_name, '": ', &
2852	& trim(nf90_strerror(istatus))
2853	call my_abort('Error reading NetCDF file')
2854	end if
2855
2856	call infile%get_variable_id(var_name, ivarid_out)
2857	if (data_type == NF90_DOUBLE .or. data_type == NF90_FLOAT) then
2858	allocate(data_real(ntotal))
2859	!istatus = nf90_get_var(infile%ncid, ivarid_in, data_real(1))
2860	istatus = nf_get_var_double(infile%ncid, ivarid_in, data_real)
2861	if (istatus /= NF90_NOERR) then
2862	deallocate(data_real)
2863	write(nulerr,'(a,a,a,a)') '*** Error reading variable "', var_name, '": ', &
2864	& trim(nf90_strerror(istatus))
2865	call my_abort('Error reading NetCDF file')
2866	end if
2867
2868	!istatus = nf90_put_var(this%ncid, ivarid_out, data_real)
2869	istatus = nf_put_var_double(this%ncid, ivarid_out, data_real)
2870	deallocate(data_real)
2871	if (istatus /= NF90_NOERR) then
2872	write(nulerr,'(a,a,a,a)') '*** Error writing variable "', var_name, '": ', &
2873	& trim(nf90_strerror(istatus))
2874	call my_abort('Error writing NetCDF file')
2875	end if
2876
2877	else
2878	allocate(data_int(ntotal))
2879	!istatus = nf90_get_var(infile%ncid, ivarid_in, data_int)
2880	istatus = nf_get_var_int(infile%ncid, ivarid_in, data_int)
2881	if (istatus /= NF90_NOERR) then
2882	deallocate(data_int)
2883
2884	write(nulerr,'(a,a,a,a)') '*** Error reading variable "', var_name, '": ', &
2885	& trim(nf90_strerror(istatus))
2886	call my_abort('Error reading NetCDF file')
2887	end if
2888
2889	!istatus = nf90_put_var(this%ncid, ivarid_out, data_int)
2890	istatus = nf_put_var_int(this%ncid, ivarid_out, data_int)
2891	deallocate(data_int)
2892	if (istatus /= NF90_NOERR) then
2893	write(nulerr,'(a,a,a,a)') '*** Error writing variable "', var_name, '": ', &
2894	& trim(nf90_strerror(istatus))
2895	call my_abort('Error writing NetCDF file')
2896	end if
2897	end if
2898
2899	end subroutine copy_variable
2900
2901	end module easy_netcdf

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

Download in other formats:

Original Format