mo_netcdf.fypp 41.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
!> \file mo_netcdf.f90

!> \brief NetCDF Fortran 90 interface wrapper

!> \details A wrapper around the NetCDF Fortran 90 interface.
!
!> \authors David Schaefer
!> \date Jun 2015

10
11
12
#:include "common.fypp"

#:set KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES + CHAR_KINDS_TYPES
13
#:set KINDS = REAL_KINDS + INT_KINDS
14
15
#:set RANKS = range(0, MAXRANK)

16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
module mo_netcdf

  ! This module provides a thin wrapper around the NetCDF Fortran 90 interface,
  ! following a object-oriented approach.

  ! Written  David Schaefer, Jun 2015
  ! Modified Matthias Cuntz, Jan 2016 - compiled with PGI Fortran rev 15.9 - no automatic allocation of left-hand-side
  ! Modified Ricardo Torres, Feb 2017 - add derived type NcGroup and NcAttributable. NcAttributable is the base derived type,
  !                                     NcGroup and NcVariable are extended from it. NcDataset is extended from NcGroup. No more
  !                                     duplicated routines to set attributes.

  ! License
  ! -------
  ! GNU Lesser General Public License http://www.gnu.org/licenses/

31
32
  use mo_kind, only: ${', '.join(KINDS)}$
  use mo_utils, only: ne
33
34
  use ieee_arithmetic, only : ieee_is_nan

35
36
  use netcdf, only : &
          nf90_open, nf90_close, nf90_strerror, nf90_def_dim, nf90_def_var, &
37
          nf90_put_var, nf90_get_var, nf90_put_att, nf90_get_att, nf90_inq_attname, &
38
39
40
41
42
          nf90_inquire, nf90_inq_dimid, nf90_inquire_dimension, &
          nf90_inq_varid, nf90_inq_varids, nf90_inquire_variable, nf90_inquire_attribute, &
          nf90_inq_ncid, nf90_inq_grp_parent, nf90_inq_grpname, nf90_def_grp, &
          nf90_rename_dim, nf90_rename_var, nf90_rename_att, nf90_sync, &
          NF90_OPEN, NF90_NETCDF4, NF90_CREATE, NF90_WRITE, NF90_NOWRITE, &
Robert Schweppe's avatar
Robert Schweppe committed
43
          NF90_BYTE, NF90_SHORT, NF90_INT, NF90_INT64, NF90_FLOAT, NF90_DOUBLE, NF90_CHAR, &
44
45
46
          NF90_FILL_BYTE, NF90_FILL_SHORT, NF90_FILL_INT, NF90_FILL_FLOAT, NF90_FILL_DOUBLE, &
          NF90_NOERR, NF90_UNLIMITED, NF90_GLOBAL, NF90_SHARE, NF90_HDF5, &
          NF90_64BIT_OFFSET, NF90_CLASSIC_MODEL
47
48
49
50

  implicit none

  ! --------------------------------------------------------------------------------------
51
52
53
54
  character(10), parameter :: CF_FILL_VALUE = '_FillValue'
  character(11), parameter :: CF_VALID_RANGE = 'valid_range'
  character(9), parameter :: CF_VALID_MIN = 'valid_min'
  character(9), parameter :: CF_VALID_MAX = 'valid_max'
55
56
57
58
59
  integer(i4), parameter :: CF_USE_FILL_VALUE = 1_i4
  integer(i4), parameter :: CF_USE_VALID_MIN = 2_i4
  integer(i4), parameter :: CF_USE_VALID_MAX = 3_i4
  integer(i4), parameter :: CF_USE_VALID_RANGE = 4_i4
  integer(i4), parameter :: CF_USE_NAN = 5_i4
60
61
62

  type, abstract :: NcBase

63
    integer(i4) :: id
64

65
  contains
66

67
68
    procedure(getNameInterface), deferred :: getName
    procedure(getParentInterface), deferred :: getParent
69
70
71
72
73

  end type NcBase

  type, abstract, extends(NcBase) :: NcAttributable

74
75
76
77
78
  contains

    procedure, public :: hasAttribute
    procedure, public :: renameAttribute
    procedure, private :: getAttributableIds
79
    procedure, public :: getAttributeNames
80

81
82
#:for kind, type in REAL_KINDS_TYPES + INT_KINDS_TYPES
  #:for rank in [0, 1]
83
84
    procedure, private :: setAttribute_${rank}$d_${kind}$
    generic, public :: setAttribute => setAttribute_${rank}$d_${kind}$
85
86
    procedure, private :: getAttribute_${rank}$d_${kind}$
    generic, public :: getAttribute => getAttribute_${rank}$d_${kind}$
87
88
89
90
  #:endfor
#:endfor
#:for kind, type in CHAR_KINDS_TYPES
  #:for rank in [0]
91
92
    procedure, private :: setAttribute_${rank}$d_${kind}$
    generic, public :: setAttribute => setAttribute_${rank}$d_${kind}$
93
94
    procedure, private :: getAttribute_${rank}$d_${kind}$
    generic, public :: getAttribute => getAttribute_${rank}$d_${kind}$
95
96
  #:endfor
#:endfor
97

98
99
100
101
102
103
104
105
  end type NcAttributable

  ! --------------------------------------------------------------------------------------

  type, extends(NcAttributable) :: NcGroup

  contains

106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
    ! getter
    procedure, private :: getVariableIds
    procedure, public :: getVariables
    procedure, public :: getUnlimitedDimension
    procedure, public :: getNoVariables

    procedure, private :: getDimensionByName
    procedure, private :: getDimensionById

    procedure, public :: getParent => getGroupParent
    procedure, public :: getName => getGroupName
    procedure, public :: getGroup => getGroupByName
    procedure, public :: getVariable => getVariableByName
    generic, public :: getDimension => &
            getDimensionById, &
            getDimensionByName

    ! checker
    procedure, public :: hasVariable
    procedure, public :: hasDimension
    procedure, public :: hasGroup
    procedure, public :: isUnlimited => isDatasetUnlimited

    ! setter
130
    procedure, public  :: setGroup
131
132
133
134
    procedure, public  :: setDimension
    procedure, public  :: setCoordinate
    procedure, private :: set_scrip_dimension
    procedure, private :: set_1d_coordinate_variable
135
136
137
138
139
140
141
142
    procedure, private :: setVariableWithTypes
    procedure, private :: setVariableWithNames
    procedure, private :: setVariableWithIds

    generic, public :: setVariable => &
            setVariableWithNames, &
            setVariableWithTypes, &
            setVariableWithIds
143
144
145
146

  end type NcGroup

  interface NcGroup
147
    procedure newNcGroup
148
149
150
151
152
153
  end interface NcGroup

  ! --------------------------------------------------------------------------------------

  type, extends(NcGroup) :: NcDataset

154
155
    character(256) :: fname !> Filename of the opened dataset
    character(1) :: mode  !> File open mode
156

157
  contains
158

159
160
    procedure, public :: sync
    procedure, public :: close
161
162
163
164

  end type NcDataset

  interface NcDataset
165
    procedure newNcDataset
166
167
  end interface NcDataset

168
  ! --------------------------------------------------------------------------------------
169
170
171

  type, extends(NcBase) :: NcDimension

172
173
174
175
176
177
178
    type(NcGroup) :: parent  !> The dimension's parent
  contains
    procedure, public :: renameDimension
    procedure, public :: getParent => getDimensionParent
    procedure, public :: getName => getDimensionName
    procedure, public :: getLength => getDimensionLength
    procedure, public :: isUnlimited => isUnlimitedDimension
179
  end type NcDimension
180

181
  interface NcDimension
182
    procedure newNcDimension
183
  end interface NcDimension
184
  ! --------------------------------------------------------------------------------------
185

186
  type, extends(NcAttributable) :: NcVariable
187
188
189
190
191
192
193
194
195
    type(NcGroup) :: parent   !> The variables's parent

  contains

    procedure, public :: renameVariable
    procedure, public :: getParent => getVariableParent
    procedure, public :: getName => getVariableName
    procedure, private :: getSlicingShape

196
197
198
199
200
201
202
#:for kind, type in REAL_KINDS_TYPES + INT_KINDS_TYPES
  #:for rank in RANKS
    procedure, private :: setData_${rank}$d_${kind}$
    generic, public :: setData => setData_${rank}$d_${kind}$
    procedure, private :: getData_${rank}$d_${kind}$
    generic, public :: getData => getData_${rank}$d_${kind}$
  #:endfor
203
204
205
206
207
208
    procedure, private :: getCFAttributes_${kind}$
    generic, public :: getCFAttributes => getCFAttributes_${kind}$
    procedure, private :: setVariableFillValue_${kind}$
    generic, public :: setFillValue => setVariableFillValue_${kind}$
    procedure, private :: getVariableFillValue_${kind}$
    generic, public :: getFillValue => getVariableFillValue_${kind}$
209
#:endfor
210
211
212
213
214
215
216
217
218
219
220
221
222
223

    procedure, public :: getNoDimensions

    procedure, public :: getDimensions => getVariableDimensions

    procedure, public :: getRank => getVariableRank

    procedure, public :: getShape => getVariableShape

    procedure, public :: getDtype => getVariableDtype

    procedure, public :: isUnlimited => isUnlimitedVariable


224
225
226
227
228
229
230
231
232
  end type NcVariable

  interface NcVariable
    procedure newNcVariable
  end interface NcVariable
  ! --------------------------------------------------------------------------------------

  ! abstract interfaces
  interface
233
234
235
236
237
238
239
240
241
242
243
    function getNameInterface(self)
      import NcBase
      class(NcBase), intent(in) :: self
      character(len = 256) :: getNameInterface
    end function getNameInterface

    function getParentInterface(self)
      import NcBase, NcGroup
      class(NcBase), intent(in) :: self
      type(NcGroup) :: getParentInterface
    end function getParentInterface
244
245
246
  end interface

  interface operator (==)
247
    procedure equalNcBases
248
249
250
251
252
  end interface operator (==)

contains

  function newNcDataset(fname, fmode, cmode) result(out)
253
254
    character(*), intent(in) :: fname
    character(1), intent(in) :: fmode
255
    character(*), intent(inout), optional :: cmode
256
    integer(i4) :: status
257
    type(NcDataset) :: out
258
259
260

    select case(fmode)
    case("w")
261
      status = nf90_create(trim(fname), getCreationMode(cmode), out%id)
262
    case("r")
263
      status = nf90_open(trim(fname), NF90_NOWRITE, out%id)
264
    case("a")
265
      status = nf90_open(trim(fname), NF90_WRITE, out%id)
266
    case default
267
268
      write(*, *) "Mode argument must be in 'w','r','a' ! "
      stop 1
269
    end select
270
    call check(status, "Failed to open file: " // fname)
271
272

    out%fname = fname
273
    out%mode = fmode
274
275
276
  end function newNcDataset

  function newNcVariable(id, parent) result(out)
277
    integer(i4), intent(in) :: id
278
    type(NcGroup), intent(in) :: parent
279
    type(NcVariable) :: out
280

281
    out%id = id
282
283
284
285
    out%parent = parent
  end function newNcVariable

  function newNcDimension(id, parent) result(out)
286
    integer(i4), intent(in) :: id
287
    type(NcGroup), intent(in) :: parent
288
    type(NcDimension) :: out
289

290
    out%id = id
291
292
293
294
    out%parent = parent
  end function newNcDimension

  function newNcGroup(id) result(out)
295
    integer(i4), intent(in) :: id
296
    type(NcGroup) :: out
297
298
299
300
301
302
303

    out%id = id
  end function newNcGroup

  subroutine sync(self)
    class(NcDataset) :: self

304
    call check(nf90_sync(self%id), "Failed to sync file: " // self%fname)
305
306
307
308
309
  end subroutine sync

  subroutine close(self)
    class(NcDataset) :: self

310
    call check(nf90_close(self%id), "Failed to close file: " // self%fname)
311
312
313
314
  end subroutine close

  function setGroup(self, name)
    class(NcGroup), intent(inout) :: self
315
    character(*), intent(in) :: name
316
    integer(i4) :: id
317
    type(NcGroup) :: setGroup
318
319
320
321
322
323
324

    call check(nf90_def_grp(self%id, name, id), "Failed to create new group: " // name)
    setGroup = NcGroup(id)
  end function setGroup

  function getGroupParent(self)
    class(NcGroup), intent(in) :: self
325
    integer(i4) :: id
326
    type(NcGroup) :: getGroupParent
327
328
329
330
331
332
333

    call check(nf90_inq_grp_parent(self%id, id), "Failed to get parent group of: " // self%getName())
    getGroupParent = NcGroup(id)
  end function getGroupParent

  function getGroupName(self)
    class(NcGroup), intent(in) :: self
334
    character(256) :: getGroupName
335
336
337
338
339
340

    call check(nf90_inq_grpname(self%id, getGroupName), "Failed to inquire group name")
  end function getGroupName

  function getNoVariables(self)
    class(NcGroup), intent(in) :: self
341
    integer(i4) :: getNoVariables
342

343
    call check(nf90_inquire(self%id, nvariables = getNoVariables), "Failed inquire number of variables")
344
345
346
347
  end function getNoVariables

  function getDimensionParent(self)
    class(NcDimension), intent(in) :: self
348
    type(NcGroup) :: getDimensionParent
349
350
351
352
353
354

    getDimensionParent = self%parent
  end function getDimensionParent

  function getVariableParent(self)
    class(NcVariable), intent(in) :: self
355
    type(NcGroup) :: getVariableParent
356
357
358
359
360

    getVariableParent = self%parent
  end function getVariableParent

  function getVariableIds(self)
361
    class(NcGroup), intent(in) :: self
362
363
    integer(i4), dimension(:), allocatable :: getVariableIds
    integer(i4) :: tmp
364
365
366
367
368
369

    allocate(getVariableIds(self%getNoVariables()))
    call check(nf90_inq_varids(self%id, tmp, getVariableIds), "Failed to inquire variable ids")
  end function getVariableIds

  function getVariables(self)
370
    class(NcGroup), intent(in) :: self
371
    type(NcVariable), dimension(:), allocatable :: getVariables
372
373
    integer(i4), dimension(:), allocatable :: varids
    integer(i4) :: i, nvars
374
375
376
377
378

    nvars = self%getNoVariables()
    allocate(getVariables(nvars), varids(nvars))

    varids = self%getVariableIds()
379
380
    do i = 1, size(varids)
      getVariables(i) = NcVariable(varids(i), self)
381
382
383
384
385
386
    end do

  end function getVariables

  function getDimensionName(self)
    class(NcDimension), intent(in) :: self
387
    character(len = 256) :: getDimensionName
388

389
390
    call check(nf90_inquire_dimension(self%parent%id, self%id, name = getDimensionName), &
            "Failed to inquire dimension name")
391
392
393
394
  end function getDimensionName

  function getDimensionLength(self)
    class(NcDimension), intent(in) :: self
395
    integer(i4) :: getDimensionLength
396

397
398
    call check(nf90_inquire_dimension(self%parent%id, self%id, len = getDimensionLength), &
            "Failed to inquire dimension: " // self%getName())
399
400
401
  end function getDimensionLength

  function isDatasetUnlimited(self)
402
403
    class(NcGroup), intent(in) :: self
    logical :: isDatasetUnlimited
404
    integer(i4) :: dimid
405

406
407
    call check(nf90_inquire(self%id, unlimitedDimId = dimid), &
            "Failed to inquire group " // self%getName())
408
    isDatasetUnlimited = (dimid /= -1)
409
410
411
  end function isDatasetUnlimited

  function getUnlimitedDimension(self)
412
413
    class(NcGroup), intent(in) :: self
    type(NcDimension) :: getUnlimitedDimension
414
    integer(i4) :: dimid
415

416
417
    call check(nf90_inquire(self%id, unlimitedDimId = dimid), &
            "Failed to inquire group " // self%getName())
418

419
    if (dimid == -1) then
420
421
      write(*, *) "Dataset has no unlimited dimension"
      stop 1
422
423
424
425
426
427
428
    end if

    getUnlimitedDimension = self%getDimension(dimid)
  end function getUnlimitedDimension

  function equalNcBases(left, right) result(out)
    class(NcBase), intent(in) :: left, right
429
    logical :: out
430

431
    out = (left%id == right%id)
432
433
434
435
  end function equalNcBases

  function isUnlimitedDimension(self)
    class(NcDimension), intent(in) :: self
436
    logical :: isUnlimitedDimension
437
438
439

    isUnlimitedDimension = .false.
    if (self%parent%isUnlimited()) then
440
      isUnlimitedDimension = (self == self%parent%getUnlimitedDimension())
441
442
443
    end if
  end function isUnlimitedDimension

444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
  function set_scrip_dimension(self, centersDim1, centersDim2, cornersDim1, cornersDim2, subDimSizes, units) &
          result(ncDim)
    class(NcGroup), intent(in) :: self
    real(dp)      , intent(in), dimension(:) :: centersDim1
    real(dp)      , intent(in), dimension(:) :: centersDim2
    real(dp)      , intent(in), dimension(:,:) :: cornersDim1
    real(dp)      , intent(in), dimension(:,:) :: cornersDim2
    integer(i4)   , intent(in), dimension(:) :: subDimSizes
    character(256), intent(in) :: units
    type(NcDimension) :: ncDim

    type(NcDimension)          :: cornerDim, rankDim
    type(NcVariable)           :: ncVar
    integer(i4), allocatable, dimension(:) :: imask_data

    ! set the new ncDimension (integer values and name)
    ncDim = self%setDimension('grid_size', size(centersDim1))
    cornerDim = self%setDimension('grid_corners', size(cornersDim1, 1))
    rankDim = self%setDimension('grid_rank', size(subDimSizes))
    ! here we set the reference to ncDimension for labelled ncDimension which in fact is a variable
    ncVar = self%setVariable('grid_center_lon', "f64", [ncDim])
    call ncVar%setData(centersDim1)
    call ncVar%setAttribute('units', trim(units))
    ncVar = self%setVariable('grid_center_lat', "f64", [ncDim])
    call ncVar%setData(centersDim2)
    call ncVar%setAttribute('units', trim(units))
    ncVar = self%setVariable('grid_corner_lon', "f64", [cornerDim, ncDim])
    call ncVar%setData(cornersDim1)
    call ncVar%setAttribute('units', trim(units))
    ncVar = self%setVariable('grid_corner_lat', "f64", [cornerDim, ncDim])
    call ncVar%setData(cornersDim2)
    call ncVar%setAttribute('units', trim(units))
    ncVar = self%setVariable('grid_dims', "i32", [rankDim])
    call ncVar%setData(subDimSizes)
    ! set all values to 1 (True) for mask
    ncVar = self%setVariable('grid_imask', "i32", [ncDim])
    allocate(imask_data(size(centersDim1)))
    imask_data = 1_i4
    call ncVar%setData(imask_data)
    deallocate(imask_data)
    call ncVar%setAttribute('units', 'unitless')

  end function set_scrip_dimension

  subroutine set_1D_coordinate_variable(self, name, ncDim, bounds, referenceArg, ncVar)
    class(NcGroup), intent(in) :: self
    character(*)  , intent(in) :: name
    type(NcDimension), intent(in) :: ncDim
    real(dp)      , intent(in), dimension(:) :: bounds
    integer(i4)   , intent(in), optional :: referenceArg
    type(NcVariable):: ncVar

    type(NcDimension) :: bndsDim
    integer(i8) :: dimLength, iBound
    integer(i4) :: reference, iAtt
    character(256) :: dimBoundName
    real(dp), allocatable, dimension(:, :) :: boundData

    ! init
    dimLength = size(bounds, kind=i8)
    reference = 1_i4
    if (present(referenceArg)) then
      reference = referenceArg
    end if
    ! here we set the reference to ncDimension for labelled ncDimension which in fact is a variable
    ncVar = self%setVariable(name, "f64", [ncDim])
    ! write the data based on the type of reference
    select case(reference)
    case(0_i4)
      ! set the start values
      call ncVar%setData(bounds(1_i8:dimLength - 1_i8))
    case(1_i4)
      ! set the center values
      call ncVar%setData((bounds(2_i8:dimLength) + bounds(1_i8:dimLength-1_i8)) / 2.0_dp)
    case(2_4)
      ! set the end values
      call ncVar%setData(bounds(2_i8:dimLength))
    case default
      write(*,*) "reference id for set_Dimension is ", reference, ", must be 0, 1 or 2."
      stop 1
    end select

    ! --- bounds ---
    ! create dimension name for bounds
    dimBoundName = trim(name) // "_bnds"
    ! set the attribute
    call ncVar%setAttribute('bounds', trim(dimBoundName))
    ! set the dimensions used for the bounds array
    if (self%hasDimension("bnds")) then
      ! add it to our bounds of ncDimensions for the current array
      bndsDim = self%getDimension("bnds")
    else
      bndsDim = self%setDimension("bnds", 2_i4)
    end if
    ncVar = self%setVariable(dimBoundName, "f64", [bndsDim, ncDim])

    ! allocate array for data
    allocate(boundData(2_i8, dimLength-1_i8))
    do iBound = 1_i8, dimLength-1_i8
      boundData(1_i8, iBound) = bounds(iBound)
      boundData(2_i8, iBound) = bounds(iBound + 1_i8)
    end do
    call ncVar%setData(boundData)
    deallocate(boundData)

  end subroutine set_1D_coordinate_variable

  function setDimension(self, name, length) result(ncDim)
552
    class(NcGroup), intent(in) :: self
553
    character(*)  , intent(in) :: name
554
    integer(i4), intent(in), optional :: length
555

556
557
558
    type(NcDimension) :: ncDim
    integer(i4) :: dimLength
    integer(i4) :: id
559

560
    dimLength = NF90_UNLIMITED
561
    if (present(length)) then
562
563
      if (length > 0_i4) then
        dimLength = length
564
565
566
      end if
    end if

567
    call check(nf90_def_dim(self%id, name, dimLength, id), &
568
569
         "Failed to create dimension: " // name)

570
    ncDim = NcDimension(id, self)
571

572
  end function setDimension
573

574
575
  function setCoordinate(self, name, length, bounds, reference, attribute_names, attribute_values, &
                        centersDim1, centersDim2, cornersDim1, cornersDim2, subDimSizes, units) result(ncDim)
576
    class(NcGroup), intent(in) :: self
577
    character(*)  , intent(in) :: name
578
    integer(i4), intent(in), optional :: length
579
    real(dp)      , intent(in), optional, dimension(:) :: bounds
580
    integer(i4)   , intent(in), optional :: reference
581
582
    character(*) , intent(in), optional, dimension(:) :: attribute_names
    character(*) , intent(in), optional, dimension(:) :: attribute_values
583
584
585
586
    real(dp)      , intent(in), optional, dimension(:) :: centersDim1
    real(dp)      , intent(in), optional, dimension(:) :: centersDim2
    real(dp)      , intent(in), optional, dimension(:,:) :: cornersDim1
    real(dp)      , intent(in), optional, dimension(:,:) :: cornersDim2
587
    integer(i4)   , intent(in), optional, dimension(:) :: subDimSizes
588
    character(*), intent(in), optional :: units
589

590
591
592
    type(NcDimension) :: ncDim
    type(NcVariable) :: ncVar
    integer(i4) :: iAtt
593
594
595

    if (present(centersDim1) .and. present(centersDim2) .and. present(cornersDim1) .and. present(cornersDim2) &
             .and. present(subDimSizes) .and. present(units)) then
596
      ncDim = self%set_scrip_dimension(centersDim1, centersDim2, cornersDim1, cornersDim2, subDimSizes, units)
597
    else
598
      ! set the new ncDimension (integer values and name)
599
      ncDim = self%setDimension(name, length)
600
601

      if (present(bounds)) then
602
        call self%set_1D_coordinate_variable(name, ncDim, bounds, reference, ncVar)
603
604
605
606
        ! set attributes
        ! already set attributes
        if (present(attribute_names) .and. present(attribute_values)) then
          do iAtt = 1, size(attribute_names)
607
            call ncVar%setAttribute(trim(attribute_names(iAtt)), &
608
609
610
611
612
                    trim(attribute_values(iAtt)))
          end do
        end if
      end if
    end if
613

614
  end function setCoordinate
615
616
617

  function hasVariable(self, name)
    class(NcGroup), intent(in) :: self
618
619
    character(*), intent(in) :: name
    logical :: hasVariable
620
    integer(i4) :: tmpid
621

622
    hasVariable = (nf90_inq_varid(self%id, name, tmpid) == NF90_NOERR)
623
624
625
626
  end function hasVariable

  function hasDimension(self, name)
    class(NcGroup), intent(in) :: self
627
628
    character(*), intent(in) :: name
    logical :: hasDimension
629
    integer(i4) :: tmpid
630

631
    hasDimension = (nf90_inq_dimid(self%id, name, tmpid) == NF90_NOERR)
632
633
634
635
  end function hasDimension

  function hasGroup(self, name)
    class(NcGroup), intent(in) :: self
636
637
    character(*), intent(in) :: name
    logical :: hasGroup
638
    integer(i4) :: tmpid
639

640
    hasGroup = (nf90_inq_ncid(self%id, name, tmpid) == NF90_NOERR)
641
642
643
  end function hasGroup

  function setVariableWithIds(self, name, dtype, dimensions, contiguous, &
644
645
646
647
648
          chunksizes, deflate_level, shuffle, fletcher32, endianness, &
          cache_size, cache_nelems, cache_preemption)
    class(NcGroup), intent(in) :: self
    character(*), intent(in) :: name
    character(*), intent(in) :: dtype
649
    integer(i4), intent(in) :: dimensions(:)
650
    logical, intent(in), optional :: contiguous, shuffle, fletcher32
651
    integer(i4), intent(in), optional :: endianness, deflate_level, cache_size, &
652
653
            cache_nelems, cache_preemption, chunksizes(:)
    type(NcVariable) :: setVariableWithIds
654
    integer(i4) :: varid, status
655
656

    status = nf90_def_var(self%id, name, getDtypeFromString(dtype), dimensions, varid, contiguous, &
657
658
            chunksizes, deflate_level, shuffle, fletcher32, endianness, &
            cache_size, cache_nelems, cache_preemption)
659
660
661
662
663
    call check(status, "Failed to create variable: " // name)
    setVariableWithIds = NcVariable(varid, self)
  end function setVariableWithIds

  function setVariableWithNames(self, name, dtype, dimensions, contiguous, &
664
665
666
667
668
669
670
671
          chunksizes, deflate_level, shuffle, fletcher32, endianness, &
          cache_size, cache_nelems, cache_preemption)

    class(NcGroup), intent(in) :: self
    character(*), intent(in) :: name
    character(*), intent(in) :: dtype
    character(*), intent(in) :: dimensions(:)
    logical, intent(in), optional :: contiguous, shuffle, fletcher32
672
    integer(i4), intent(in), optional :: endianness, deflate_level, cache_size, &
673
674
675
            cache_nelems, cache_preemption, chunksizes(:)
    type(NcVariable) :: setVariableWithNames
    type(NcDimension) :: dim
676
    integer(i4) :: i, dimids(size(dimensions))
677
678
679
680

    do i = 1, size(dimensions)
      dim = self%getDimension(dimensions(i))
      dimids(i) = dim%id
681
682
683
    end do

    setVariableWithNames = setVariableWithIds(self, name, dtype, dimids, contiguous, &
684
685
            chunksizes, deflate_level, shuffle, fletcher32, endianness, &
            cache_size, cache_nelems, cache_preemption)
686
687
688
  end function setVariableWithNames

  function setVariableWithTypes(self, name, dtype, dimensions, contiguous, &
689
690
691
692
693
694
695
          chunksizes, deflate_level, shuffle, fletcher32, endianness, &
          cache_size, cache_nelems, cache_preemption)
    class(NcGroup), intent(in) :: self
    character(*), intent(in) :: name
    character(*), intent(in) :: dtype
    type(NcDimension), intent(in) :: dimensions(:)
    logical, intent(in), optional :: contiguous, shuffle, fletcher32
696
    integer(i4), intent(in), optional :: endianness, deflate_level, cache_size, &
697
698
699
            cache_nelems, cache_preemption, chunksizes(:)
    type(NcVariable) :: setVariableWithTypes
    type(NcDimension) :: dim
700
    integer(i4) :: i, dimids(size(dimensions))
701
702
703
704

    do i = 1, size(dimensions)
      dim = dimensions(i)
      dimids(i) = dim%id
705
706
707
    end do

    setVariableWithTypes = setVariableWithIds(self, name, dtype, dimids, contiguous, &
708
709
            chunksizes, deflate_level, shuffle, fletcher32, endianness, &
            cache_size, cache_nelems, cache_preemption)
710
711
712
713
  end function setVariableWithTypes

  function getDimensionById(self, id)
    class(NcGroup), intent(in) :: self
714
    integer(i4) :: id
715
716
717
718
719
720
721
    type(NcDimension) :: getDimensionById
    character(32) :: msg, name

    write(msg, *) id
    call check(nf90_inquire_dimension(self%id, id, name), &
            "Could not inquire dimension: " // msg)
    getDimensionById = NcDimension(id, self)
722
723
724
725
  end function getDimensionById

  function getDimensionByName(self, name)
    class(NcGroup), intent(in) :: self
726
727
    character(*) :: name
    type(NcDimension) :: getDimensionByName
728
    integer(i4) :: id
729

730
731
    call check(nf90_inq_dimid(self%id, name, id), &
            "Could not inquire dimension: " // name)
732
733
734
735
736
    getDimensionByName = self%getDimensionById(id)
  end function getDimensionByName

  function getGroupByName(self, name)
    class(NcGroup), intent(in) :: self
737
738
    character(*), intent(in) :: name
    type(NcGroup) :: getGroupByName
739
    integer(i4) :: id
740
741

    call check(nf90_inq_ncid(self%id, name, id), &
742
            "Could not inquire variable: " // name)
743
744
745
746
747
    getGroupByName = NcGroup(id)
  end function getGroupByName

  function getVariableByName(self, name)
    class(NcGroup), intent(in) :: self
748
749
    character(*), intent(in) :: name
    type(NcVariable) :: getVariableByName
750
    integer(i4) :: id
751
752

    call check(nf90_inq_varid(self%id, name, id), &
753
            "Could not inquire variable: " // name)
754
755
756
757
758
759
    getVariableByName = NcVariable(id, self)

  end function getVariableByName

  function getVariableName(self)
    class(NcVariable), intent(in) :: self
760
    character(len = 256) :: getVariableName
761

762
763
    call check(nf90_inquire_variable(self%parent%id, self%id, name = getVariableName), &
            "Could not inquire variable name")
764
765
766
767
  end function getVariableName

  function getNoDimensions(self)
    class(NcVariable), intent(in) :: self
768
    integer(i4) :: getNoDimensions
769

770
771
    call check(nf90_inquire_variable(self%parent%id, self%id, ndims = getNoDimensions), &
            "Could not inquire variable: " // self%getName())
772
773
774
  end function getNoDimensions

  function getVariableDimensions(self)
775
    class(NcVariable), intent(in) :: self
776
    type(NcDimension), allocatable :: getVariableDimensions(:)
777
778
    integer(i4), allocatable :: dimids(:)
    integer(i4) :: ii, ndims
779
780
781

    ndims = self%getNoDimensions()
    allocate(dimids(ndims), getVariableDimensions(ndims))
782
783
    call check(nf90_inquire_variable(self%parent%id, self%id, dimids = dimids), &
            "Could not inquire variable: " // self%getName())
784

785
786
    do ii = 1, ndims
      getVariableDimensions (ii) = self%parent%getDimension(dimids(ii))
787
788
789
790
    end do
  end function getVariableDimensions

  function getVariableShape(self)
791
    class(NcVariable), intent(in) :: self
792
    integer(i4), allocatable :: getVariableShape(:)
793
    type(NcDimension), allocatable :: dims(:)
794
    integer(i4) :: ii, ndims
795
796
797
798
799

    ndims = self%getNoDimensions()
    allocate(getVariableShape(ndims), dims(ndims))

    dims = self%getDimensions()
800
801
    do ii = 1, size(dims)
      getVariableShape(ii) = dims(ii)%getLength()
802
803
804
805
    end do
  end function getVariableShape

  function getVariableRank(self)
806
    class(NcVariable), intent(in) :: self
807
    integer(i4) :: getVariableRank
808
809
810
811
812
813

    getVariableRank = size(self%getDimensions())
  end function getVariableRank

  function getVariableDtype(self)
    class(NcVariable), intent(in) :: self
814
    integer(i4) :: dtype
815
    character(3) :: getVariableDtype
816

817
818
    call check(nf90_inquire_variable(self%parent%id, self%id, xtype = dtype), &
            "Could not inquire variable: " // self%getName())
819
820
821
822
    getVariableDtype = getDtypeFromInteger(dtype)
  end function getVariableDtype

  function isUnlimitedVariable(self)
823
824
    class(NcVariable), intent(in) :: self
    logical :: isUnlimitedVariable
825
    type(NcDimension), allocatable :: dims(:)
826
    type(NcDimension) :: dim
827
    integer(i4) :: ii
828
829
830
831
832
833

    allocate(dims(self%getNoDimensions()))

    isUnlimitedVariable = .false.
    dims = self%getDimensions()

834
835
836
837
838
    do ii = 1, size(dims)
      dim = dims(ii)
      if (dim%isUnlimited()) then
        isUnlimitedVariable = .true.
      end if
839
840
841
    end do
  end function isUnlimitedVariable

Robert Schweppe's avatar
Robert Schweppe committed
842
  logical function hasAttribute(self, name, xtype, len, attnum)
843
    class(NcAttributable), intent(in) :: self
844
    character(*), intent(in) :: name
Robert Schweppe's avatar
Robert Schweppe committed
845
846
847
    integer(i4), intent(out), optional :: xtype
    integer(i4), intent(out), optional :: len
    integer(i4), intent(out), optional :: attnum
848
    integer(i4) :: status
849
850
851

    select type (self)
    class is (NcGroup)
852
      print*, 'hasAttribute group'
Robert Schweppe's avatar
Robert Schweppe committed
853
      status = nf90_inquire_attribute(self%id, NF90_GLOBAL, name, xtype=xtype, len=len, attnum=attnum)
854
    class is (NcVariable)
855
      print*, 'hasAttribute var'
Robert Schweppe's avatar
Robert Schweppe committed
856
      status = nf90_inquire_attribute(self%parent%id, self%id, name, xtype=xtype, len=len, attnum=attnum)
857
858
    end select

859
    hasAttribute = (status == NF90_NOERR)
860
861
  end function hasAttribute

862
863
864
865
866
  function getAttributeNames(self) result(attributeNames)
    class(NcAttributable), intent(in) :: self
    character(256), dimension(:), allocatable :: attributeNames

    integer(i4), parameter :: maxNames = 100_i4
867
    integer(i4) :: nAtts
Robert Schweppe's avatar
Robert Schweppe committed
868
    integer(i4) :: status
869
870

    ! assume a maximum number of 100 attributes that are checked
Robert Schweppe's avatar
Robert Schweppe committed
871
    allocate(attributeNames(maxNames))
872
    nAtts = 0_i4
873
874
875
    do while (nAtts < maxNames)
      select type (self)
      class is (NcGroup)
876
        print*, 'is group', nAtts + 1_i4
Robert Schweppe's avatar
Robert Schweppe committed
877
        status = nf90_inq_attname(self%id, NF90_GLOBAL, nAtts + 1_i4, attributeNames(nAtts + 1_i4))
878
      class is (NcVariable)
879
        print*, 'is variable', nAtts + 1_i4
Robert Schweppe's avatar
Robert Schweppe committed
880
        status = nf90_inq_attname(self%parent%id, self%id, nAtts + 1_i4, attributeNames(nAtts + 1_i4))
881
      end select
Robert Schweppe's avatar
Robert Schweppe committed
882
883
      ! if the status is negative, exit loop, else increase counter
      if (status /= NF90_NOERR) then
884
885
886
887
888
889
890
891
892
893
        exit
      else
        nAtts = nAtts + 1_i4
      end if
    end do
    ! select only valid names
    attributeNames = attributeNames(1:nAtts)

  end function getAttributeNames

894
895
#:def setAttribute_template(kind, type, rank)
  subroutine setAttribute_${rank}$d_${kind}$(self, name, data)
896
    class(NcAttributable), intent(in) :: self
897
898
    character(len=*), intent(in) :: name
    ${type}$, intent(in) :: data${ranksuffix(rank)}$
899
    integer(i4) :: ids(2)
900

901
    ids = self%getAttributableIds()
902
    call check(nf90_put_att(ids(1), ids(2), name, data), &
903
            "Failed to write attribute: " // name)
904

905
  end subroutine setAttribute_${rank}$d_${kind}$
906

907
#:enddef setAttribute_template
908

909
910
#:def getAttribute_template(kind, type, rank)
  subroutine getAttribute_${rank}$d_${kind}$(self, name, avalue)
911
    class(NcAttributable), intent(in) :: self
912
913
    character(len=*), intent(in) :: name
    ${type}$, intent(out) :: avalue${ranksuffix(rank)}$
914
    integer(i4) :: length, ids(2)
915
916

    ids = self%getAttributableIds()
917
918
    call check(nf90_inquire_attribute(ids(1), ids(2), name, len = length), &
            "Could not inquire attribute " // name)
919
    call check(nf90_get_att(ids(1), ids(2), name, avalue), &
920
            "Could not read attribute " // name)
921

922
  end subroutine getAttribute_${rank}$d_${kind}$
923

924
#:enddef getAttribute_template
925

926
927
928
929
930
931
932
933
934
935
936
937
#:for kind, type in REAL_KINDS_TYPES + INT_KINDS_TYPES
  #:for rank in [0, 1]
    $:setAttribute_template(kind, type, rank)
    $:getAttribute_template(kind, type, rank)
  #:endfor
#:endfor
#:for kind, type in CHAR_KINDS_TYPES
  #:for rank in [0]
    $:setAttribute_template(kind, type, rank)
    $:getAttribute_template(kind, type, rank)
  #:endfor
#:endfor
938
939
940

  function getAttributableIds(self)
    class(NcAttributable), intent(in) :: self
941
    integer(i4) :: getAttributableIds(2)
942
943
    select type(self)
    class is (NcGroup)
944
945
      getAttributableIds(1) = self%id
      getAttributableIds(2) = NF90_GLOBAL
946
    class is (NcVariable)
947
948
      getAttributableIds(1) = self%parent%id
      getAttributableIds(2) = self%id
949
950
951
952
953
    end select
  end function getAttributableIds

  subroutine renameAttribute(self, oldname, newname)
    class(NcAttributable), intent(inout) :: self
954
    character(len = *), intent(in) :: oldname, newname
955
    integer(i4) :: ids(2)
956
957
958
959
960
961
    ids = self%getAttributableIds()
    call check(nf90_rename_att(ids(1), ids(2), oldname, newname), "Failed to rename attribute: " // oldname)
  end subroutine renameAttribute

  subroutine renameVariable(self, name)
    class(NcVariable), intent(inout) :: self
962
    character(len = *), intent(in) :: name
963
964
965
966
967
    call check(nf90_rename_var(self%parent%id, self%id, name), "Failed to rename variable: " // self%getName())
  end subroutine renameVariable

  subroutine renameDimension(self, name)
    class(NcDimension), intent(inout) :: self
968
    character(len = *), intent(in) :: name
969
970
971
    call check(nf90_rename_dim(self%parent%id, self%id, name), "Failed to rename dimension: " // self%getName())
  end subroutine renameDimension

972
973
#:def setFillValue_template(kind, type)
  subroutine setVariableFillValue_${kind}$(self, fvalue)
974
    class(NcVariable), intent(inout) :: self
975
    ${type}$, intent(in) :: fvalue
976

977
978
    if (.not. self%hasAttribute(CF_FILL_VALUE)) then
      call self%setAttribute(CF_FILL_VALUE, fvalue)
979
980
    end if

981
  end subroutine setVariableFillValue_${kind}$
982

983
#:enddef setFillValue_template
984

985
986
987
988
#:def getFillValue_template(kind, type)
  subroutine getVariableFillValue_${kind}$(self, fvalue)
    class(NcVariable), intent(in) :: self
    ${type}$, intent(out) :: fvalue
989

990
991
    if (self%hasAttribute(CF_FILL_VALUE)) then
      call self%getAttribute(CF_FILL_VALUE, fvalue)
992
    else
993
      $: '      fvalue = ' + {'i1': 'NF90_FILL_BYTE', 'i2': 'NF90_FILL_SHORT', 'i4': 'NF90_FILL_INT', 'i8': 'NF90_FILL_INT', 'sp': 'NF90_FILL_FLOAT', 'dp': 'NF90_FILL_DOUBLE', }[kind]
994
995
    end if

996
  end subroutine getVariableFillValue_${kind}$
997

998
#:enddef getFillValue_template
999

1000
#:def getCFAttributes_template(kind, type)
For faster browsing, not all history is shown. View entire blame