netcdf_utils.fpp Source File

# include "define.inc"

!> Provides some convenience wrappers around netCDF procedures, mostly
!> for handling netCDF errors in a consistent manner, removing the
!> need to manually check the returned status code.
!>
!> @warning This module requires the preprocessor macro `NETCDF` to be
!>          defined. If it is not defined, you won't be able to `use`
!>          this module!
module netcdf_utils

#ifndef NETCDF
! Can't just do `#error` because `make depends` always preprocesses
! this file, so force a compile time error instead. `#warning` may get
! overlooked if it's much earlier in the build process
Error: You are attempting to compile netcdf_utils.fpp but are building without netcdf. Ensure you guard any 'use netcdf_utils' statements with '#ifdef NETCDF'

#else

  use netcdf, only: NF90_FLOAT, NF90_DOUBLE
  use netcdf, only: NF90_NOWRITE, NF90_CLOBBER, NF90_NOERR, NF90_MAX_NAME
  use netcdf, only: NF90_ENOTATT, NF90_ENOTVAR
  use netcdf, only: nf90_strerror
  use netcdf, only: nf90_inquire_variable
  use netcdf, only: nf90_inquire_dimension
  use netcdf, only: nf90_open, nf90_close
  use netcdf, only: nf90_inq_varid
  use netcdf, only: NF90_INT
  use netcdf, only: nf90_inq_varid, nf90_def_var
  use netcdf, only: nf90_inquire_attribute, nf90_put_att

  implicit none

  private

  public :: netcdf_error
  public :: get_netcdf_code_precision
  public :: check_netcdf_file_precision
  public :: netcdf_real, kind_nf, netcdf_int
  public :: ensure_netcdf_var_exists
  public :: ensure_netcdf_att_exists
  public :: ensure_netcdf_dim_exists
  public :: get_netcdf_dim_length
  public :: set_default_netcdf_compression
  public :: get_default_netcdf_compression
  public :: set_default_netcdf_compression_level
  public :: get_default_netcdf_compression_level

  !> Get the netCDF ID for a variable, creating it if it doesn't exist
  !> already
  interface ensure_netcdf_var_exists
    module procedure ensure_netcdf_var_exists_scalar
    module procedure ensure_netcdf_var_exists_onedim
    module procedure ensure_netcdf_var_exists_manydims
  end interface ensure_netcdf_var_exists

  !> Add an attribute to an existing netCDF variable
  interface ensure_netcdf_att_exists
    module procedure ensure_netcdf_att_exists_text
  end interface ensure_netcdf_att_exists

  !> Kind of netCDF types and constants
  integer, parameter :: kind_nf = kind (NF90_NOERR)
  !> netCDF type of default `real`
#ifdef SINGLE_PRECISION
  integer (kind_nf), parameter :: netcdf_real = NF90_FLOAT
#else
#ifdef QUAD_PRECISION
# error "Quad precision is incompatible with netCDF"
#else
  integer (kind_nf), parameter :: netcdf_real = NF90_DOUBLE
#endif
#endif
  !> netCDF type of default `integer`
  integer (kind_nf), parameter :: netcdf_int = NF90_INT

  ! Enables some internal print statements for debugging
  logical, parameter :: test = .false.

  ! The default use of compression
#ifdef GK_NETCDF_DEFAULT_COMPRESSION_ON
  logical :: use_compression_default = .true.
#else
  logical :: use_compression_default = .false.
#endif

  ! The default compression level. Can be between 1 (fast) and 9 (more)
  integer(kind_nf) :: deflate_level_default = 1

contains

  !> Returns the correct netCDF constant for the `real` kind in use
  function get_netcdf_code_precision () result (code_real)
    use constants, only: pi, kind_rs, kind_rd
    use file_utils, only: error_unit
    integer :: code_real

    ! second condition for Cray
    if ( (kind(pi)==kind_rs) .or. (kind_rs==kind_rd) ) then
       code_real = NF90_FLOAT
    else if (kind(pi)==kind_rd) then
       code_real = NF90_DOUBLE
    else
       write (error_unit(),*) &
            'ERROR: precision mismatch in get_netcdf_code_precision'
    end if
  end function get_netcdf_code_precision

  !> Check that an existing netCDF file uses the same `real` kind as
  !> the current simulation. Only prints a warning!
  subroutine check_netcdf_file_precision (ncid, filename)
    use file_utils, only: error_unit
    implicit none
    integer (kind_nf), intent (in), optional :: ncid
    character (*), intent (in), optional :: filename
    integer (kind_nf) :: file_real
    integer (kind_nf) :: ist, ncid_private, tid
    integer :: ierr

    ierr = error_unit()

    ist = NF90_NOERR
    file_real = -1

    if (present(ncid)) then
       if (present(filename)) then
          write (ierr,*) 'WARNING: in calling check_netcdf_file_precision'
          write (ierr,*) &
               'WARNING: both filename and ncid given -- filename ignored'
       end if
       ncid_private = ncid
    else
       if (present(filename)) then
          ist = nf90_open (filename, NF90_NOWRITE, ncid_private)
          if (test) write (error_unit(),*) &
               'opened netcdf file ', trim(filename), ' with ncid: ', &
               ncid_private, ' in check_netcdf_file_precision'
          if (ist /= NF90_NOERR) then
             call netcdf_error (ist, file=filename)
             return
          end if
       else
          write (ierr,*) 'ERROR: in calling check_netcdf_file_precision'
          write (ierr,*) 'ERROR: either filename or ncid should be given'
          return
       end if
    end if

    ist = nf90_inq_varid (ncid_private, 't0', tid)
    if (ist /= NF90_NOERR) call netcdf_error (ist, var='t0')

    ! get file_real
    if (ist == NF90_NOERR) then
       ist = nf90_inquire_variable (ncid_private, tid, xtype=file_real)
       if (ist /= NF90_NOERR) call netcdf_error (ist, ncid_private, tid)
    end if

    if (.not.present(ncid)) then
       ist = nf90_close (ncid_private)
       if (ist /= NF90_NOERR) call netcdf_error (ist, file=filename)
    end if

    ! check if file_real == code_real
    if (file_real /= netcdf_real) then
       write (ierr,*) 'WARNING: precision mismatch in input netcdf file and running code'
       if (file_real == NF90_FLOAT) then
          write (ierr,*) 'WARNING: file_real = NF90_FLOAT'
       else if (file_real == NF90_DOUBLE) then
          write (ierr,*) 'WARNING: file_real = NF90_DOUBLE'
       else
          write (ierr,*) 'WARNING: unknown file_real', file_real
       end if
       if (netcdf_real == NF90_FLOAT) then
          write (ierr,*) 'WARNING: code_real = NF90_FLOAT'
       else if (netcdf_real == NF90_DOUBLE) then
          write (ierr,*) 'WARNING: code_real = NF90_DOUBLE'
       else
          write (ierr,*) 'WARNING: unknown code_real'
       end if
    end if
  end subroutine check_netcdf_file_precision

  !> Convert a netCDF error code to a nice error message. Writes to
  !> [[file_utils:error_unit]]
  !>
  !> All the arguments except `istatus` are optional and are used to
  !> give more detailed information about the error
  !>
  !> The `abort` argument will abort the simulation
  subroutine netcdf_error &
       (istatus, ncid, varid, dimid, file, dim, var, att, message, abort)

    use mp, only: mp_abort
    use file_utils, only: error_unit
    use netcdf, only: NF90_GLOBAL
    implicit none
    !> The netCDF error code to convert to a message
    integer (kind_nf), intent (in) :: istatus
    !> netCDF ID of the file
    integer (kind_nf), intent (in), optional :: ncid
    !> netCDF ID of the variable
    integer (kind_nf), intent (in), optional :: varid
    !> netCDF ID of the dimension
    integer (kind_nf), intent (in), optional :: dimid
    !> Name of the file
    character (*), intent (in), optional :: file
    !> Name of the dimension
    character (*), intent (in), optional :: dim
    !> Name of the variable
    character (*), intent (in), optional :: var
    !> Name of the attribute
    character (*), intent (in), optional :: att
    !> Custom text to append to the error message
    character (*), intent (in), optional :: message
    !> Immediately abort the program if present and true
    logical, intent (in), optional :: abort
    integer (kind_nf) :: ist
    integer :: ierr
    character (20) :: varname, dimname

    ierr = error_unit()

    write (ierr, '(2a)', advance='no') 'ERROR: ', trim (nf90_strerror (istatus))

    if (present(file)) &
         write (ierr, '(2a)', advance='no') ' in file: ', trim (file)

    if (present(dim)) &
         write (ierr, '(2a)', advance='no') ' in dimension: ', trim (dim)

    if (present(var)) &
         write (ierr, '(2a)', advance='no') ' in variable: ', trim (var)

    if (present(varid)) then
       if (present(ncid)) then
          if (present(att) ) then
             if (varid == NF90_GLOBAL) then
                write (ierr, '(2a)') ' in global attribute: ', trim(att)
             else
                write (ierr, '(2a)') ' with the attribute: ', trim(att)
             end if
          else
             ist = nf90_inquire_variable (ncid, varid, varname)
             if (ist == NF90_NOERR) then
                write (ierr, '(a,i8,2a)', advance='no') ' in varid: ', varid, &
                     & ' variable name: ', trim (varname)
             else
                write (ierr, *) ''
                write (ierr, '(3a,i8,a,i8)', advance='no') 'ERROR in netcdf_error: ', &
                     trim (nf90_strerror(ist)), ' in varid: ', varid, &
                     ', ncid: ', ncid
             end if
          end if
       else
          write (ierr, *) ''
          write (ierr, '(2a)', advance='no') 'ERROR in netcdf_error: ', &
               & 'ncid missing while varid present in the argument'
       end if
    end if

    if (present(dimid)) then
       if (present(ncid)) then
          ist = nf90_inquire_dimension (ncid, dimid, dimname)
          if (ist == NF90_NOERR) then
             write (ierr, '(a,i8,2a)', advance='no') ' in dimid: ', dimid, &
                  & ' dimension name: ', trim (dimname)
          else
             write (ierr, *) ''
             write (ierr, '(3a,i8,a,i8)', advance='no') 'ERROR in netcdf_error: ', &
                  trim (nf90_strerror(ist)), ' in dimid: ', dimid, &
                  ', ncid: ', ncid
          end if
       else
          write (ierr, *) ''
          write (ierr, '(2a)', advance='no') 'ERROR in netcdf_error: ', &
               & 'ncid missing while dimid present in the argument'
       end if
    end if

    if (present(message)) write (ierr, '(a)', advance='no') trim(message)

    ! append line-break
    write(ierr,*)

    ! maybe this switching is not necessary.
    ! if error is detected, the program should abort immediately
    if(present(abort)) then
       if(abort) then
          call mp_abort('Aborted by netcdf_error')
       endif
    endif
  end subroutine netcdf_error

  !> Get the netCDF ID for a variable, creating it if it doesn't exist
  !> already. Aborts if any errors are detected.
  subroutine ensure_netcdf_var_exists_scalar(file_id, var_name, var_type, var_id)
    !> ID of the file or group to look for or create the variable under
    integer(kind_nf), intent(in) :: file_id
    !> Name of the variable
    character(len=*), intent(in) :: var_name
    !> The netCDF type of the variable
    integer(kind_nf), intent(in) :: var_type
    !> The netCDF ID of the variable under `file_id`
    integer(kind_nf), intent(out) :: var_id

    ! Error code of the netCDF calls
    integer(kind_nf) :: status

    status = nf90_inq_varid(file_id, var_name, var_id)
    ! Variable doesn't exist, so let's create it
    if (status == NF90_ENOTVAR) then
       status = nf90_def_var(file_id, var_name, var_type, var_id)
    end if
    ! Something went wrong with one of the previous two calls
    if (status /= NF90_NOERR) then
      call netcdf_error(status, var=var_name, varid=var_id, &
                        message="(ensure_netcdf_var_exists)", abort=.true.)
    end if
  end subroutine ensure_netcdf_var_exists_scalar

  !> Get the netCDF ID for a variable, creating it if it doesn't exist
  !> already. Aborts if any errors are detected.
  subroutine ensure_netcdf_var_exists_onedim(file_id, var_name, var_type, dim_id, var_id, &
       compress, deflate_level)
    !> ID of the file or group to look for or create the variable under
    integer(kind_nf), intent(in) :: file_id
    !> Name of the variable
    character(len=*), intent(in) :: var_name
    !> The netCDF type of the variable
    integer(kind_nf), intent(in) :: var_type
    !> Array of dimension IDs
    integer(kind_nf), intent(in) :: dim_id
    !> The netCDF ID of the variable under `file_id`
    integer(kind_nf), intent(out) :: var_id
    !> Do we want to use compression for this variable
    logical, intent(in), optional :: compress
    !> The compression level to use for this variable
    integer(kind_nf), intent(in), optional :: deflate_level

    call ensure_netcdf_var_exists_manydims(file_id, var_name, var_type, [dim_id], &
         var_id, compress, deflate_level)
  end subroutine ensure_netcdf_var_exists_onedim

  !> Get the netCDF ID for a variable, creating it if it doesn't exist
  !> already. Aborts if any errors are detected.
  subroutine ensure_netcdf_var_exists_manydims(file_id, var_name, var_type, dim_id, var_id, &
       compress, deflate_level)
    use optionals, only: get_option_with_default
    !> ID of the file or group to look for or create the variable under
    integer(kind_nf), intent(in) :: file_id
    !> Name of the variable
    character(len=*), intent(in) :: var_name
    !> The netCDF type of the variable
    integer(kind_nf), intent(in) :: var_type
    !> Array of dimension IDs
    integer(kind_nf), dimension(:), intent(in) :: dim_id
    !> The netCDF ID of the variable under `file_id`
    integer(kind_nf), intent(out) :: var_id
    !> Do we want to use compression for this variable
    logical, intent(in), optional :: compress
    !> The compression level to use for this variable
    integer(kind_nf), intent(in), optional :: deflate_level

    ! Error code of the netCDF calls
    integer(kind_nf) :: status
    ! Do we want to use compression?
    logical :: use_compression
    ! What compression level do we want to use?
    integer(kind_nf) :: compression_level

    status = nf90_inq_varid(file_id, var_name, var_id)
    ! Variable doesn't exist, so let's create it
    if (status == NF90_ENOTVAR) then
       use_compression = get_option_with_default(compress, use_compression_default)
       compression_level = 0
       if (use_compression) &
            compression_level = get_option_with_default(deflate_level, deflate_level_default)

       status = nf90_def_var(file_id, var_name, var_type, dim_id, var_id, &
            deflate_level = compression_level, &
            shuffle = use_compression &
            )
    end if
    ! Something went wrong with one of the previous two calls
    if (status /= NF90_NOERR) then
      call netcdf_error(status, var=var_name, varid=var_id, &
                        message="(ensure_netcdf_var_exists)", abort=.true.)
    end if
  end subroutine ensure_netcdf_var_exists_manydims

  !> Add an attribute to an existing netCDF variable. Aborts if any
  !> errors are detected.
  subroutine ensure_netcdf_att_exists_text(file_id, var_id, att_name, att_text)
    !> ID of the file or group to look for or create the variable under
    integer(kind_nf), intent(in) :: file_id
    !> The netCDF ID of the variable under `file_id`
    integer(kind_nf), intent(in) :: var_id
    !> Name of the attribute
    character(len=*), intent(in) :: att_name
    !> Attribute text to store
    character(len=*), intent(in) :: att_text

    ! Error code of the netCDF calls
    integer(kind_nf) :: status

    status = nf90_inquire_attribute(file_id, var_id, att_name)
    ! Attribute doesn't exist, so let's create it
    if (status == NF90_ENOTATT) then
      status = nf90_put_att(file_id, var_id, att_name, att_text)
    end if
    ! Something went wrong with one of the previous two calls
    if (status /= NF90_NOERR) then
      call netcdf_error(status, file_id, var_id, att=att_name, &
                        message="(ensure_netcdf_att_exists)", abort=.true.)
    end if
  end subroutine ensure_netcdf_att_exists_text

  !> Get the netCDF ID for a dimension, creating it if it doesn't
  !> exist already. Aborts if any errors are detected.
  subroutine ensure_netcdf_dim_exists(file_id, dim_name, dim_size, dim_id)
    use netcdf, only: nf90_inq_dimid, nf90_def_dim, NF90_NOERR, NF90_EBADDIM
    !> ID of the file or group to look for or create the dimension under
    integer(kind_nf), intent(in) :: file_id
    !> Name of the dimension
    character(len=*), intent(in) :: dim_name
    !> The size of the dimension
    integer(kind_nf), intent(in) :: dim_size
    !> The netCDF ID of the dimension under `file_id`
    integer(kind_nf), intent(out) :: dim_id

    ! Error code of the netCDF calls
    integer(kind_nf) :: status

    status = nf90_inq_dimid(file_id, dim_name, dim_id)
    ! Dimension doesn't exist, so let's create it
    if (status == NF90_EBADDIM) then
      status = nf90_def_dim(file_id, dim_name, dim_size, dim_id)
    end if
    ! Something went wrong with one of the previous two calls
    if (status /= NF90_NOERR) then
      call netcdf_error(status, dim=dim_name, message="(ensure_netcdf_dim_exists)", abort=.true.)
    end if
  end subroutine ensure_netcdf_dim_exists

  !> Wrapper around `nf90_inquire_dimension`. Aborts if any errors are
  !> detected.
  function get_netcdf_dim_length(file_id, dim_id) result(length)
    integer :: length
    !> NetCDF ID of the file
    integer(kind_nf), intent(in) :: file_id
    !> NetCDF ID of the dimension
    integer(kind_nf), intent(in) :: dim_id

    ! Temporary variable for the dimension name
    character(len=NF90_MAX_NAME) :: dim_name
    ! NetCDF error status
    integer :: status

    ! This could be less reliant on module-level variables if we
    ! passed in the dimension name and used nf90_inq_dim to get the ID
    status = nf90_inquire_dimension(file_id, dim_id, dim_name, length)
    if (status /= NF90_NOERR) then
      call netcdf_error(status, &
                        dim=trim(dim_name), &
                        message="(get_netcdf_dim_length)", &
                        abort=.true.)
    end if
  end function get_netcdf_dim_length

  !> Returns the current value of use_compression_default
  logical function get_default_netcdf_compression() result(compression_flag)
    implicit none
    compression_flag = use_compression_default
  end function get_default_netcdf_compression

  !> Set the value of use_compression_default
  subroutine set_default_netcdf_compression(compression_flag)
    implicit none
    logical, intent(in) :: compression_flag
    use_compression_default = compression_flag
  end subroutine set_default_netcdf_compression

  !> Returns the current value of deflate_level_default
  integer(kind_nf) function get_default_netcdf_compression_level() result(compression_level)
    implicit none
    compression_level = deflate_level_default
  end function get_default_netcdf_compression_level

  !> Set the value of deflate_level_default
  subroutine set_default_netcdf_compression_level(compression_level)
    implicit none
    integer(kind_nf), intent(in) :: compression_level
    deflate_level_default = compression_level
  end subroutine set_default_netcdf_compression_level

#endif
end module netcdf_utils