source: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad.v1.5.1/easy_netcdf.F90 @ 5080

Last change on this file since 5080 was 4489, checked in by lguez, 21 months ago

Merge LMDZ_ECRad branch back into trunk!

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
24module 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
121contains
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/(ndimlen1*ndimlen2*ndimlen3)
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
2901end module easy_netcdf
Note: See TracBrowser for help on using the repository browser.