mo_netcdf.fypp

!> \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

#:include "common.fypp"

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

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/

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

  use netcdf, only : &
          nf90_open, nf90_close, nf90_strerror, nf90_def_dim, nf90_def_var, &
          nf90_put_var, nf90_get_var, nf90_put_att, nf90_get_att, nf90_inq_attname, &
          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, &
          NF90_BYTE, NF90_SHORT, NF90_INT, NF90_INT64, NF90_FLOAT, NF90_DOUBLE, NF90_CHAR, &
          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

  implicit none

  ! --------------------------------------------------------------------------------------
  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'
  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

  type, abstract :: NcBase

    integer(i4) :: id

  contains

    procedure(getNameInterface), deferred :: getName
    procedure(getParentInterface), deferred :: getParent

  end type NcBase

  type, abstract, extends(NcBase) :: NcAttributable

  contains

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

#:for kind, type in REAL_KINDS_TYPES + INT_KINDS_TYPES
  #:for rank in [0, 1]
    procedure, private :: setAttribute_${rank}$d_${kind}$
    generic, public :: setAttribute => setAttribute_${rank}$d_${kind}$
    procedure, private :: getAttribute_${rank}$d_${kind}$
    generic, public :: getAttribute => getAttribute_${rank}$d_${kind}$
  #:endfor
#:endfor
#:for kind, type in CHAR_KINDS_TYPES
  #:for rank in [0]
    procedure, private :: setAttribute_${rank}$d_${kind}$
    generic, public :: setAttribute => setAttribute_${rank}$d_${kind}$
    procedure, private :: getAttribute_${rank}$d_${kind}$
    generic, public :: getAttribute => getAttribute_${rank}$d_${kind}$
  #:endfor
#:endfor

  end type NcAttributable

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

  type, extends(NcAttributable) :: NcGroup

  contains

    ! 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
    procedure, public  :: setGroup
    procedure, public  :: setDimension
    procedure, public  :: setCoordinate
    procedure, private :: set_scrip_dimension
    procedure, private :: set_1d_coordinate_variable
    procedure, private :: setVariableWithTypes
    procedure, private :: setVariableWithNames
    procedure, private :: setVariableWithIds

    generic, public :: setVariable => &
            setVariableWithNames, &
            setVariableWithTypes, &
            setVariableWithIds

  end type NcGroup

  interface NcGroup
    procedure newNcGroup
  end interface NcGroup

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

  type, extends(NcGroup) :: NcDataset

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

  contains

    procedure, public :: sync
    procedure, public :: close

  end type NcDataset

  interface NcDataset
    procedure newNcDataset
  end interface NcDataset

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

  type, extends(NcBase) :: NcDimension

    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
  end type NcDimension

  interface NcDimension
    procedure newNcDimension
  end interface NcDimension
  ! --------------------------------------------------------------------------------------

  type, extends(NcAttributable) :: NcVariable
    type(NcGroup) :: parent   !> The variables's parent

  contains

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

#: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
    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}$
#:endfor

    procedure, public :: getNoDimensions

    procedure, public :: getDimensions => getVariableDimensions

    procedure, public :: getRank => getVariableRank

    procedure, public :: getShape => getVariableShape

    procedure, public :: getDtype => getVariableDtype

    procedure, public :: isUnlimited => isUnlimitedVariable


  end type NcVariable

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

  ! abstract interfaces
  interface
    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
  end interface

  interface operator (==)
    procedure equalNcBases
  end interface operator (==)

contains

  function newNcDataset(fname, fmode, cmode) result(out)
    character(*), intent(in) :: fname
    character(1), intent(in) :: fmode
    character(*), intent(inout), optional :: cmode
    integer(i4) :: status
    type(NcDataset) :: out

    select case(fmode)
    case("w")
      status = nf90_create(trim(fname), getCreationMode(cmode), out%id)
    case("r")
      status = nf90_open(trim(fname), NF90_NOWRITE, out%id)
    case("a")
      status = nf90_open(trim(fname), NF90_WRITE, out%id)
    case default
      write(*, *) "Mode argument must be in 'w','r','a' ! "
      stop 1
    end select
    call check(status, "Failed to open file: " // fname)

    out%fname = fname
    out%mode = fmode
  end function newNcDataset

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

    out%id = id
    out%parent = parent
  end function newNcVariable

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

    out%id = id
    out%parent = parent
  end function newNcDimension

  function newNcGroup(id) result(out)
    integer(i4), intent(in) :: id
    type(NcGroup) :: out

    out%id = id
  end function newNcGroup

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

    call check(nf90_sync(self%id), "Failed to sync file: " // self%fname)
  end subroutine sync

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

    call check(nf90_close(self%id), "Failed to close file: " // self%fname)
  end subroutine close

  function setGroup(self, name)
    class(NcGroup), intent(inout) :: self
    character(*), intent(in) :: name
    integer(i4) :: id
    type(NcGroup) :: setGroup

    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
    integer(i4) :: id
    type(NcGroup) :: getGroupParent

    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
    character(256) :: getGroupName

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

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

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

  function getDimensionParent(self)
    class(NcDimension), intent(in) :: self
    type(NcGroup) :: getDimensionParent

    getDimensionParent = self%parent
  end function getDimensionParent

  function getVariableParent(self)
    class(NcVariable), intent(in) :: self
    type(NcGroup) :: getVariableParent

    getVariableParent = self%parent
  end function getVariableParent

  function getVariableIds(self)
    class(NcGroup), intent(in) :: self
    integer(i4), dimension(:), allocatable :: getVariableIds
    integer(i4) :: tmp

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

  function getVariables(self)
    class(NcGroup), intent(in) :: self
    type(NcVariable), dimension(:), allocatable :: getVariables
    integer(i4), dimension(:), allocatable :: varids
    integer(i4) :: i, nvars

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

    varids = self%getVariableIds()
    do i = 1, size(varids)
      getVariables(i) = NcVariable(varids(i), self)
    end do

  end function getVariables

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

    call check(nf90_inquire_dimension(self%parent%id, self%id, name = getDimensionName), &
            "Failed to inquire dimension name")
  end function getDimensionName

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

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

  function isDatasetUnlimited(self)
    class(NcGroup), intent(in) :: self
    logical :: isDatasetUnlimited
    integer(i4) :: dimid

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

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

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

    if (dimid == -1) then
      write(*, *) "Dataset has no unlimited dimension"
      stop 1
    end if

    getUnlimitedDimension = self%getDimension(dimid)
  end function getUnlimitedDimension

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

    out = (left%id == right%id)
  end function equalNcBases

  function isUnlimitedDimension(self)
    class(NcDimension), intent(in) :: self
    logical :: isUnlimitedDimension

    isUnlimitedDimension = .false.
    if (self%parent%isUnlimited()) then
      isUnlimitedDimension = (self == self%parent%getUnlimitedDimension())
    end if
  end function isUnlimitedDimension

  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)
    class(NcGroup), intent(in) :: self
    character(*)  , intent(in) :: name
    integer(i4), intent(in), optional :: length

    type(NcDimension) :: ncDim
    integer(i4) :: dimLength
    integer(i4) :: id

    dimLength = NF90_UNLIMITED
    if (present(length)) then
      if (length > 0_i4) then
        dimLength = length
      end if
    end if

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

    ncDim = NcDimension(id, self)

  end function setDimension

  function setCoordinate(self, name, length, bounds, reference, attribute_names, attribute_values, &
                        centersDim1, centersDim2, cornersDim1, cornersDim2, subDimSizes, units) result(ncDim)
    class(NcGroup), intent(in) :: self
    character(*)  , intent(in) :: name
    integer(i4), intent(in), optional :: length
    real(dp)      , intent(in), optional, dimension(:) :: bounds
    integer(i4)   , intent(in), optional :: reference
    character(*) , intent(in), optional, dimension(:) :: attribute_names
    character(*) , intent(in), optional, dimension(:) :: attribute_values
    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
    integer(i4)   , intent(in), optional, dimension(:) :: subDimSizes
    character(*), intent(in), optional :: units

    type(NcDimension) :: ncDim
    type(NcVariable) :: ncVar
    integer(i4) :: iAtt

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

      if (present(bounds)) then
        call self%set_1D_coordinate_variable(name, ncDim, bounds, reference, ncVar)
        ! set attributes
        ! already set attributes
        if (present(attribute_names) .and. present(attribute_values)) then
          do iAtt = 1, size(attribute_names)
            call ncVar%setAttribute(trim(attribute_names(iAtt)), &
                    trim(attribute_values(iAtt)))
          end do
        end if
      end if
    end if

  end function setCoordinate

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

    hasVariable = (nf90_inq_varid(self%id, name, tmpid) == NF90_NOERR)
  end function hasVariable

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

    hasDimension = (nf90_inq_dimid(self%id, name, tmpid) == NF90_NOERR)
  end function hasDimension

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

    hasGroup = (nf90_inq_ncid(self%id, name, tmpid) == NF90_NOERR)
  end function hasGroup

  function setVariableWithIds(self, name, dtype, dimensions, contiguous, &
          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
    integer(i4), intent(in) :: dimensions(:)
    logical, intent(in), optional :: contiguous, shuffle, fletcher32
    integer(i4), intent(in), optional :: endianness, deflate_level, cache_size, &
            cache_nelems, cache_preemption, chunksizes(:)
    type(NcVariable) :: setVariableWithIds
    integer(i4) :: varid, status

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

  function setVariableWithNames(self, name, dtype, dimensions, contiguous, &
          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
    integer(i4), intent(in), optional :: endianness, deflate_level, cache_size, &
            cache_nelems, cache_preemption, chunksizes(:)
    type(NcVariable) :: setVariableWithNames
    type(NcDimension) :: dim
    integer(i4) :: i, dimids(size(dimensions))

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

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

  function setVariableWithTypes(self, name, dtype, dimensions, contiguous, &
          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
    integer(i4), intent(in), optional :: endianness, deflate_level, cache_size, &
            cache_nelems, cache_preemption, chunksizes(:)
    type(NcVariable) :: setVariableWithTypes
    type(NcDimension) :: dim
    integer(i4) :: i, dimids(size(dimensions))

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

    setVariableWithTypes = setVariableWithIds(self, name, dtype, dimids, contiguous, &
            chunksizes, deflate_level, shuffle, fletcher32, endianness, &
            cache_size, cache_nelems, cache_preemption)
  end function setVariableWithTypes

  function getDimensionById(self, id)
    class(NcGroup), intent(in) :: self
    integer(i4) :: id
    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)
  end function getDimensionById

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

    call check(nf90_inq_dimid(self%id, name, id), &
            "Could not inquire dimension: " // name)
    getDimensionByName = self%getDimensionById(id)
  end function getDimensionByName

  function getGroupByName(self, name)
    class(NcGroup), intent(in) :: self
    character(*), intent(in) :: name
    type(NcGroup) :: getGroupByName
    integer(i4) :: id

    call check(nf90_inq_ncid(self%id, name, id), &
            "Could not inquire variable: " // name)
    getGroupByName = NcGroup(id)
  end function getGroupByName

  function getVariableByName(self, name)
    class(NcGroup), intent(in) :: self
    character(*), intent(in) :: name
    type(NcVariable) :: getVariableByName
    integer(i4) :: id

    call check(nf90_inq_varid(self%id, name, id), &
            "Could not inquire variable: " // name)
    getVariableByName = NcVariable(id, self)

  end function getVariableByName

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

    call check(nf90_inquire_variable(self%parent%id, self%id, name = getVariableName), &
            "Could not inquire variable name")
  end function getVariableName

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

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

  function getVariableDimensions(self)
    class(NcVariable), intent(in) :: self
    type(NcDimension), allocatable :: getVariableDimensions(:)
    integer(i4), allocatable :: dimids(:)
    integer(i4) :: ii, ndims

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

    do ii = 1, ndims
      getVariableDimensions (ii) = self%parent%getDimension(dimids(ii))
    end do
  end function getVariableDimensions

  function getVariableShape(self)
    class(NcVariable), intent(in) :: self
    integer(i4), allocatable :: getVariableShape(:)
    type(NcDimension), allocatable :: dims(:)
    integer(i4) :: ii, ndims

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

    dims = self%getDimensions()
    do ii = 1, size(dims)
      getVariableShape(ii) = dims(ii)%getLength()
    end do
  end function getVariableShape

  function getVariableRank(self)
    class(NcVariable), intent(in) :: self
    integer(i4) :: getVariableRank

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

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

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

  function isUnlimitedVariable(self)
    class(NcVariable), intent(in) :: self
    logical :: isUnlimitedVariable
    type(NcDimension), allocatable :: dims(:)
    type(NcDimension) :: dim
    integer(i4) :: ii

    allocate(dims(self%getNoDimensions()))

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

    do ii = 1, size(dims)
      dim = dims(ii)
      if (dim%isUnlimited()) then
        isUnlimitedVariable = .true.
      end if
    end do
  end function isUnlimitedVariable

  logical function hasAttribute(self, name, xtype, len, attnum)
    class(NcAttributable), intent(in) :: self
    character(*), intent(in) :: name
    integer(i4), intent(out), optional :: xtype
    integer(i4), intent(out), optional :: len
    integer(i4), intent(out), optional :: attnum
    integer(i4) :: status

    select type (self)
    class is (NcGroup)
      print*, 'hasAttribute group'
      status = nf90_inquire_attribute(self%id, NF90_GLOBAL, name, xtype=xtype, len=len, attnum=attnum)
    class is (NcVariable)
      print*, 'hasAttribute var'
      status = nf90_inquire_attribute(self%parent%id, self%id, name, xtype=xtype, len=len, attnum=attnum)
    end select

    hasAttribute = (status == NF90_NOERR)
  end function hasAttribute

  function getAttributeNames(self) result(attributeNames)
    class(NcAttributable), intent(in) :: self
    character(256), dimension(:), allocatable :: attributeNames

    integer(i4), parameter :: maxNames = 100_i4
    integer(i4) :: nAtts
    integer(i4) :: status

    ! assume a maximum number of 100 attributes that are checked
    allocate(attributeNames(maxNames))
    nAtts = 0_i4
    do while (nAtts < maxNames)
      select type (self)
      class is (NcGroup)
        print*, 'is group', nAtts + 1_i4
        status = nf90_inq_attname(self%id, NF90_GLOBAL, nAtts + 1_i4, attributeNames(nAtts + 1_i4))
      class is (NcVariable)
        print*, 'is variable', nAtts + 1_i4
        status = nf90_inq_attname(self%parent%id, self%id, nAtts + 1_i4, attributeNames(nAtts + 1_i4))
      end select
      ! if the status is negative, exit loop, else increase counter
      if (status /= NF90_NOERR) then
        exit
      else
        nAtts = nAtts + 1_i4
      end if
    end do
    ! select only valid names
    attributeNames = attributeNames(1:nAtts)

  end function getAttributeNames

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

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

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

#:enddef setAttribute_template

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

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

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

#:enddef getAttribute_template

#: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

  function getAttributableIds(self)
    class(NcAttributable), intent(in) :: self
    integer(i4) :: getAttributableIds(2)
    select type(self)
    class is (NcGroup)
      getAttributableIds(1) = self%id
      getAttributableIds(2) = NF90_GLOBAL
    class is (NcVariable)
      getAttributableIds(1) = self%parent%id
      getAttributableIds(2) = self%id
    end select
  end function getAttributableIds

  subroutine renameAttribute(self, oldname, newname)
    class(NcAttributable), intent(inout) :: self
    character(len = *), intent(in) :: oldname, newname
    integer(i4) :: ids(2)
    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
    character(len = *), intent(in) :: name
    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
    character(len = *), intent(in) :: name
    call check(nf90_rename_dim(self%parent%id, self%id, name), "Failed to rename dimension: " // self%getName())
  end subroutine renameDimension

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

    if (.not. self%hasAttribute(CF_FILL_VALUE)) then
      call self%setAttribute(CF_FILL_VALUE, fvalue)
    end if

  end subroutine setVariableFillValue_${kind}$

#:enddef setFillValue_template

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

    if (self%hasAttribute(CF_FILL_VALUE)) then
      call self%getAttribute(CF_FILL_VALUE, fvalue)
    else
      $: '      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]
    end if

  end subroutine getVariableFillValue_${kind}$

#:enddef getFillValue_template

#:def getCFAttributes_template(kind, type)