mo_ncwrite.f90

Go to the documentation of this file.

!> \file mo_ncwrite.f90
!> \brief \copybrief mo_ncwrite
!> \details \copydetails mo_ncwrite
 
!> \brief Writing netcdf files
!> \details Subroutines for writing arrays on nc file using the netcdf4 library.
!> \author Stephan Thober, Luis Samaniego, Matthias Cuntz
!> \date Nov 2011
!> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
!! FORCES is released under the LGPLv3+ license \license_note
module mo_ncwrite
 
  use mo_kind, only : i1, i4, sp, dp
  use mo_string_utils, only : nonull
  use mo_utils, only : ne
 
  ! functions and constants of netcdf4 library
  use netcdf, only : nf90_create, nf90_def_dim, nf90_unlimited, nf90_def_var, &
          nf90_char, nf90_put_att, nf90_int, nf90_int, nf90_global, &
          nf90_enddef, nf90_put_var, nf90_float, nf90_double, nf90_byte, &
          nf90_close, nf90_noerr, nf90_strerror, nf90_clobber, &
          nf90_max_name, nf90_write, nf90_inq_varid, nf90_inquire_variable, &
          nf90_inquire_dimension, nf90_open, nf90_64bit_offset, nf90_netcdf4, &
          nf90_inq_varid, nf90_inq_dimid, nf90_inquire, nf90_get_var, nf90_fill_float, &
          nf90_fill_double, nf90_fill_int
 
  ! public routines -------------------------------------------------------------------
  public :: close_netcdf         ! save and close the netcdf file
  public :: create_netcdf        ! create the nc file with variables and their attributes, after they were set
  public :: dump_netcdf          ! simple dump of variable into a netcdf file
  public :: var2nc               ! simple dump of multiple variables with attributes into one netcdf file
  public :: write_dynamic_netcdf ! write dynamically (one record after the other) in the file
  public :: write_static_netcdf  ! write static data in the file
  ! public types -------------------------------------------------------------------
  public :: dims
  public :: variable
  public :: attribute
 
  ! public parameters
  integer(i4), public, parameter :: nmaxdim = 5         !< nr. max dimensions
  integer(i4), public, parameter :: nmaxatt = 20        !< nr. max attributes
  integer(i4), public, parameter :: maxlen = 256       !< nr. string length
  integer(i4), public, parameter :: ngatt = 20        !< nr. global attributes
  integer(i4), public, parameter :: nattdim = 2         !< dim array of attribute values
 
  ! public types -----------------------------------------------------------------
  !> \brief NetCDF dims
  type dims
    character (len = maxLen) :: name                !< dim. name
    integer(i4) :: len                 !< dim. lenght, undefined time => NF90_UNLIMITED
    integer(i4) :: dimid               !< dim. Id
  end type dims
 
  !> \brief NetCDF attribute
  type attribute
    character (len = maxLen) :: name                !< attribute name
    integer(i4) :: xtype               !< attribute of the values
    integer(i4) :: nvalues             !< number of attributes
    character (len = maxLen) :: values              !< numbers or "characters" separed by spaces
  end type attribute
 
  !> \brief NetCDF variable
  type variable
    character (len = maxLen) :: name                !< short name
    integer(i4) :: xtype               !< NF90 var. type
    integer(i4) :: nlvls               !< number of levels
    integer(i4) :: nsubs               !< number of subparts
    logical :: unlimited           !< time limited
    integer(i4) :: varid               !< Id
    integer(i4) :: ndims               !< field dimension
    integer(i4), dimension(nMaxDim) :: dimids              !< passing var. dimensions
    integer(i4), dimension(nMaxDim) :: dimtypes            !< type of dimensions
    integer(i4) :: natt                !< nr. attributes
    type(attribute), dimension(nMaxAtt) :: att                 !< var. attributes
    integer(i4), dimension(nMaxDim) :: start               !< starting indices for netcdf
    integer(i4), dimension(nMaxDim) :: count               !< counter          for netcdf
    logical :: wflag               !< write flag
    integer(i1), pointer :: g0_b              !< array pointing model variables
    integer(i1), dimension(:), pointer :: g1_b              !< array pointing model variables
    integer(i1), dimension(:, :), pointer :: g2_b              !< array pointing model variables
    integer(i1), dimension(:, :, :), pointer :: g3_b              !< array pointing model variables
    integer(i1), dimension(:, :, :, :), pointer :: g4_b              !< array pointing model variables
    integer(i4), pointer :: g0_i              !< array pointing model variables
    integer(i4), dimension(:), pointer :: g1_i              !< array pointing model variables
    integer(i4), dimension(:, :), pointer :: g2_i              !< array pointing model variables
    integer(i4), dimension(:, :, :), pointer :: g3_i              !< array pointing model variables
    integer(i4), dimension(:, :, :, :), pointer :: g4_i              !< array pointing model variables
    real(sp), pointer :: g0_f              !< array pointing model variables
    real(sp), dimension(:), pointer :: g1_f              !< array pointing model variables
    real(sp), dimension(:, :), pointer :: g2_f              !< array pointing model variables
    real(sp), dimension(:, :, :), pointer :: g3_f              !< array pointing model variables
    real(sp), dimension(:, :, :, :), pointer :: g4_f              !< array pointing model variables
    real(dp), pointer :: g0_d              !< array pointing model variables
    real(dp), dimension(:), pointer :: g1_d              !< array pointing model variables
    real(dp), dimension(:, :), pointer :: g2_d              !< array pointing model variables
    real(dp), dimension(:, :, :), pointer :: g3_d              !< array pointing model variables
    real(dp), dimension(:, :, :, :), pointer :: g4_d              !< array pointing model variables
  end type variable
 
  ! public variables -----------------------------------------------------------------
  integer(i4), public :: nvars   !< nr. variables
  integer(i4), public :: ndims   !< nr. dimensions
  type (dims), public, dimension(:), allocatable :: dnc     !< dimensions list
  type(variable), public, dimension(:), allocatable :: v       !< variable list, THIS STRUCTURE WILL BE WRITTEN IN THE FILE
  type(attribute), public, dimension(nGAtt) :: gatt    !< global attributes for netcdf
 
 
  ! ------------------------------------------------------------------
 
  !>    \brief Variable simple write in netcdf.
 
  !>    \details
  !!    Simple write of a variable in a netcdf file.\n
  !!    The variabel can be 1 to 5 dimensional and single or double precision.\n
  !!    1D and 2D are dumped as static variables. From 3 to 5 dimension, the last
  !!    dimension will be defined as time.\n
  !!    The Variable will be called var.
  !!
  !!    \b Example
  !!
  !!         call dump_netcdf('test.nc', myarray)
  !!         call dump_netcdf('test.nc', myarray, netcdf4=.true.)
  !!
  !!    See also example in test!
 
  !>    \param[in]  "character(len=*) :: filename"                  Name of netcdf output file.
  !>    \param[in]  "real(sp/dp)      :: arr(:[,:[,:[,:[,:]]]])"    1D to 5D-array with input numbers.
  !>    \param[in]  "logical, optional :: lfs"                      True: enable netcdf3 large file support, i.e. 64-bit offset.
  !>    \param[in]  "logical, optional :: logical"                  True: use netcdf4 format.
  !>    \param[in]  "integer(i4), optional   :: deflate_level"      Compression level in netcdf4 (default: 1).
 
  !>    \author Matthias Cuntz
  !>    \date Nov 2012
  !>    \author Stephan Thober
  !>    \date Nov 2012
  !!      - added functions for i4 variables
  !>    \author Matthias Cuntz and Juliane Mai
  !>    \date Nov 2012
  !!      - append
  !!      - fake time dimension for 1D and 2D
  !!      - make i4 behave exactly as sp and dp
  !>    \date Mar 2013
  !!      - lfs, netcdf4, deflate_level
  interface dump_netcdf
    module procedure dump_netcdf_1d_sp, dump_netcdf_2d_sp, dump_netcdf_3d_sp, &
            dump_netcdf_4d_sp, dump_netcdf_5d_sp, &
            dump_netcdf_1d_dp, dump_netcdf_2d_dp, dump_netcdf_3d_dp, &
            dump_netcdf_4d_dp, dump_netcdf_5d_dp, &
            dump_netcdf_1d_i4, dump_netcdf_2d_i4, dump_netcdf_3d_i4, &
            dump_netcdf_4d_i4, dump_netcdf_5d_i4
  end interface dump_netcdf
 
 
  ! ------------------------------------------------------------------
 
  !>    \brief Extended dump_netcdf for multiple variables
 
  !>    \details
  !!    Write different variables including attributes to netcdf
  !!    file. The attributes are restricted to long_name, units,
  !!    and missing_value. It is also possible to append variables
  !!    when an unlimited dimension is specified.
  !!
  !!    \b Example
  !!
  !!    Let <field> be some three dimensional array
  !!
  !!    \code{.f90}
  !!    dnames(1) = 'x'
  !!    dnames(2) = 'y'
  !!    dnames(3) = 'time'
  !!    \endcode
  !!
  !!    The simplest call to write <field> to a file is
  !!
  !!    \code{.f90}
  !!    call var2nc('test.nc', field, dnames, 'h')
  !!    \endcode
  !!
  !!    With attributes it looks like
  !!
  !!    \code{.f90}
  !!    call var2nc('test.nc', field, dnames, 'h', &
  !!    long_name = 'height', units = '[m]', missing_value = -9999)
  !!    \endcode
  !!
  !!    or alternatively
  !!
  !!    \code{.f90}
  !!    character(256), dimension(3,2) :: attributes
  !!    attributes(1,1) = 'long_name'
  !!    attributes(1,2) = 'precipitation'
  !!    attributes(2,1) = 'units'
  !!    attributes(2,2) = '[mm/d]'
  !!    attributes(3,1) = 'missing_value'
  !!    attributes(3,2) = '-9999.'
  !!    call var2nc('test.nc', field, dnames, 'h', attributes = attributes, create = .true. )
  !!    \endcode
  !!
  !!    To be able to dynamically write <field>, an unlimited dimension
  !!    needs to be specified (in this example field could also be only two
  !!    dimensional)
  !!
  !!    \code{.f90}
  !!    call var2nc('test.nc', field(:,:,1), dnames, 'h', dim_unlimited=3)
  !!    \endcode
  !!
  !!    Now one can append an arbitrary number of time steps, e.g., the next 9
  !!    and the time has to be added again before
  !!
  !!    \code{.f90}
  !!    call var2nc('test.nc', (/20,...,100/), dnames(3:3), 'time',
  !!                dim_unlimited = 1 )
  !!    call var2nc('test.nc', field(:,:,2:10, dnames, 'h', dim_unlimited=3)
  !!    \endcode
  !!
  !!    You can also write another variable sharing the same dimensions
  !!
  !!    \code{.f90}
  !!    call var2nc('test.nc', field_2, dnames(1:2), 'h_2')
  !!    \endcode
  !!
  !!    The netcdf file can stay open after the first call and subsequent calls can use the file unit
  !!
  !!    \code{.f90}
  !!    ncid = -1_i4
  !!    call var2nc('test.nc', field_1, dnames(1:1), 'h_1', ncid=ncid) ! opens file
  !!    call var2nc('test.nc', field_2, dnames(1:2), 'h_2', ncid=ncid) ! uses ncid from last call
  !!    call close_netcdf(ncid)
  !!    \endcode
  !!
  !!    One can also give the start record number (on the unlimited dimension)
  !!
  !!    \code{.f90}
  !!    ncid = -1_i4
  !!    call var2nc('test.nc', time1, dnames(3:3), 'time', dim_unlimited=1_i4, ncid=ncid, create=.true.)
  !!    do i=1, n
  !!        call var2nc('test.nc', field_2(:,:,i), dnames(1:3), 'h_2', dim_unlimited=3_i4, ncid=ncid, nrec=i)
  !!    end do
  !!    call close_netcdf(ncid)
  !!    \endcode
  !!
  !!    That's it, enjoy!
  !!
  !!    \b Literature
  !!
  !!    The manual of the used netcdf fortran library can be found in
  !!    Robert Pincus & Ross Rew, The netcdf Fortran 90 Interface Guide
 
  !>        \param[in] "character(*) :: f_name"                             filename
  !>        \param[in] "integer(i4)/real(sp,dp) :: arr(:[,:[,:[,:[,:]]]])"  array to write
  !>        \param[in] "character(*) :: dnames(:)"                          dimension names
  !>        \param[in] "character(*) :: v_name"                             variable name
  !>        \param[in] "integer(i4),             optional :: dim_unlimited"     index of unlimited dimension
  !>        \param[in] "character(*),            optional :: long_name"         attribute
  !>        \param[in] "character(*),            optional :: units"             attribute
  !>        \param[in] "integer(i4)/real(sp,dp), optional :: missing_value"     attribute
  !>        \param[in] "character(256), dimension(:,:), optional :: attributes" two dimensional array of attributes
  !!                                                                size of first dimension equals number of attributes
  !!                                                                first entry of second dimension equals attribute name
  !!                                                                (e.g. long_name)\n
  !!                                                                second entry of second dimension equals attribute value
  !!                                                                (e.g. precipitation)\n
  !!                                                                every attribute is written as string with the exception of
  !!                                                                missing_value
  !>        \param[in] "logical,                 optional :: create"            flag - specify whether a
  !>                                                                            output file should be
  !>                                                                            created, default
  !>        \param[inout] "integer(i4)/real(sp,dp), optional :: ncid"   if not given filename will be opened and closed
  !!                                                                    if given and <0 then file will be opened
  !!                                                                    and ncid will return the file unit.
  !!                                                                    if given and >0 then file is assumed open and
  !!                                                                    ncid is used as file unit.
  !>        \param[in] "integer(i4),                optional :: nrec"   if given: start point on unlimited dimension.
  !>        \note It is not allowed to write the following numbers for the indicated type\n
  !!                        number | kind\n
  !!                -2.1474836E+09 | integer(i4)\n
  !!                 9.9692100E+36 | real(sp)\n
  !!        9.9692099683868690E+36 | real(dp)\n
  !!        These numbers are netcdf fortran 90 constants! They are used to determine the
  !!        chunksize of the already written variable. Hence, this routine cannot append
  !!        correctly to variables when these numbers are used. Only five dimensional
  !!        variables can be written, only one unlimited dimension can be defined.
 
  !>    \author Stephan Thober & Matthias Cuntz
  !>    \date May 2014
  !!        - created
  !>    \date Jun 2014
  !!        - added deflate, shuffle, and chunksizes
  !!        - automatically append variable at the end, renamed _FillValue to missing_value
  !>    \date Jul 2014
  !!        - add attributes array, introduced unlimited dimension that is added to the dimensions of the given array
  !>    \date Jan 2015
  !!        - changed chunk_size convention to one chunk per unit in unlimited dimension (typically time)
  !>    \date Feb 2015
  !!        - d_unlimit was not set in 5d cases
  !!        - use ne from mo_utils for fill value comparisons
  !!        - dummy(1) was sp instead of i4 in var2nc_1d_i4
  !>    \date May 2015
  !!        - ncid for opening the file only once
  !!        - nrec for writing a specific record
 
  interface var2nc
    module procedure var2nc_1d_i4, var2nc_1d_sp, var2nc_1d_dp, &
            var2nc_2d_i4, var2nc_2d_sp, var2nc_2d_dp, &
            var2nc_3d_i4, var2nc_3d_sp, var2nc_3d_dp, &
            var2nc_4d_i4, var2nc_4d_sp, var2nc_4d_dp, &
            var2nc_5d_i4, var2nc_5d_sp, var2nc_5d_dp
  end interface var2nc
 
  ! ----------------------------------------------------------------------------
 
  private
 
  ! ----------------------------------------------------------------------------
 
contains
 
  ! ----------------------------------------------------------------------------
 
  !>    \brief  Closes netcdf file stream.
 
  !>    \details
  !!    Closes a stream of an open netcdf file and saves the file.
  !!
  !!    \b Example
  !!
  !!    See test_mo_ncwrite.
  !!
  !!    \b Literature
  !!
  !!    1. http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90.html
 
  !>    \param[in]  "integer(i4) :: ncid"  Stream id of an open netcdf file which shall be closed
 
  !>    \note Closes only an already open stream
 
  !>    \author Luis Samaniego
  !>    \date   Feb 2011
 
  !>    \author Stephan Thober
  !>    \date   Dec 2011
  !!        - added comments and generalized
 
  !>    \author Matthias Cuntz
  !>    \date   Jan 2012
  !!        - Info
  !>    \date   Mar 2013
  !!        - removed Info
 
  subroutine close_netcdf(ncId)
 
    implicit none
 
    integer(i4), intent(in) :: ncid
 
    ! close: save new netcdf dataset
    call check(nf90_close(ncid))
 
  end subroutine close_netcdf
 
  ! ------------------------------------------------------------------------------
 
  !>    \brief Open and write on new netcdf file.
 
  !>    \details
  !!    This subroutine will open a new netcdf file and write the variable
  !!    attributes stored in the structure V in the file. Therefore V
  !!    has to be already set. See the file set_netcdf in test_mo_ncwrite
  !!    for an example.
  !!
  !!    \b Literature
  !!
  !!    1. http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90.html
 
  !>    \param[in]  "character(len=maxLen) :: File"             Filename of file to be written
 
  !>    \param[in]  "logical, optional     :: lfs"              True: enable netcdf3 large file support, i.e. 64-bit offset
  !>    \param[in]  "logical, optional     :: logical"          True: use netcdf4 format
  !>    \param[in]  "integer(i4), optional :: deflate_level"    compression level in netcdf4 (default: 1)
  !>    \param[out] "integer(i4)           :: ncid"             integer value of the stream for the opened file
 
  !>    \note This routine only writes attributes and variables which have been stored in V
  !!     nothing else.
 
  !>    \author Luis Samaniego
  !>    \date  Feb 2011
 
  !>    \author Stephan Thober
  !>    \date   Dec 2011
  !!        - added comments and generalized
  !>    \date   Feb 2013
  !!        - added flag for large file support
  !>    \date   Mar 2013
  !!        - buffersize
 
  !>    \author Matthias Cuntz
  !>    \date   Jan 2012
  !!        - Info
  !>    \date   Mar 2013
  !!        - netcdf4, deflate_level
  !!        - removed Info
 
  subroutine create_netcdf(Filename, ncid, lfs, netcdf4, deflate_level)
 
    implicit none
 
    ! netcdf related variables
    character(len = *), intent(in) :: filename
    integer(i4), intent(out) :: ncid
    logical, intent(in), optional :: lfs           ! netcdf3 Large File Support
    logical, intent(in), optional :: netcdf4       ! netcdf4
    integer(i4), intent(in), optional :: deflate_level ! compression level in netcdf4
 
    integer(i4) :: i, j, k
    integer(i4), dimension(nAttDim) :: att_int
    real(sp), dimension(nAttDim) :: att_float
    real(dp), dimension(nAttDim) :: att_double
    character(len = maxLen), dimension(nAttDim) :: att_char
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
    integer(i4), dimension(:), allocatable :: chunksizes ! Size of chunks in netcdf4 writing
 
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .true.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
    ! 1  Create netcdf dataset: enter define mode     ->  get ncId
    if (inetcdf4) then
      call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
    else
      if (largefile) then
        call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
      else
        ! let the netcdf library chose a buffersize, that results in lesser write system calls
        call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
      end if
    end if
 
    ! 2  Define dimensions                                 -> get dimId
    do i = 1, ndims
      call check(nf90_def_dim(ncid, dnc(i)%name, dnc(i)%len, dnc(i)%dimId))
    end do
 
    ! 3 Define dimids array, which is used to pass the dimids of the dimensions of
    ! the netcdf variables
    do i = 1, nvars
      v(i)%unlimited = .false.
      v(i)%dimids = 0
      v(i)%start = 1
      v(i)%count = 1
      do k = 1, v(i)%nDims
        if (dnc(v(i)%dimTypes(k))%len  == nf90_unlimited) v(i)%unlimited = .true.
        v(i)%dimids(k) = dnc(v(i)%dimTypes(k))%dimId
      end do
      if (v(i)%unlimited) then
        ! set counts for unlimited files (time is always the last dimension)
        if (v(i)%nDims == 1) cycle
        do k = 1, v(i)%nDims - 1
          v(i)%count(k) = dnc(v(i)%dimTypes(k))%len
        end do
      end if
    end do
 
    ! 4 Define the netcdf variables and attributes                            -> get varId
    allocate(chunksizes(maxval(v(1 : nvars)%nDims)))
    do i = 1, nvars
      if (.not. v(i)%wFlag) cycle
      if (inetcdf4) then
        chunksizes(1 : v(i)%nDims) = dnc(v(i)%dimTypes(1 : v(i)%nDims))%len
        chunksizes(v(i)%nDims) = 1
        call check(nf90_def_var(ncid, v(i)%name, v(i)%xtype, v(i)%dimids(1 : v(i)%nDims), v(i)%varId, &
                chunksizes = chunksizes(1 : v(i)%nDims), shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, v(i)%name, v(i)%xtype, v(i)%dimids(1 : v(i)%nDims), v(i)%varId))
      end if
      do k = 1, v(i)%nAtt
        select case (v(i)%att(k)%xType)
        case (nf90_char)
          ! read(V(i)%att(k)%values, *) (att_CHAR(j), j =1, V(i)%att(k)%nValues)
          read(v(i)%att(k)%values, '(a)')  att_char(1)
          call check(nf90_put_att(ncid, v(i)%varId, v(i)%att(k)%name, att_char(1)))
        case (nf90_int)
          read(v(i)%att(k)%values, *) (att_int(j), j = 1, v(i)%att(k)%nValues)
          call check(nf90_put_att(ncid, v(i)%varId, v(i)%att(k)%name, att_int(1 : v(i)%att(k)%nValues)))
        case (nf90_float)
          read(v(i)%att(k)%values, *) (att_float(j), j = 1, v(i)%att(k)%nValues)
          call check(nf90_put_att(ncid, v(i)%varId, v(i)%att(k)%name, att_float(1 : v(i)%att(k)%nValues)))
        case (nf90_double)
          read(v(i)%att(k)%values, *) (att_double(j), j = 1, v(i)%att(k)%nValues)
          call check(nf90_put_att(ncid, v(i)%varId, v(i)%att(k)%name, att_double(1 : v(i)%att(k)%nValues)))
        end select
      end do
    end do
 
    ! 5 Global attributes
    do k = 1, ngatt
      if (nonull(gatt(k)%name)) then
        call check(nf90_put_att(ncid, nf90_global, gatt(k)%name, gatt(k)%values))
      end if
    end do
 
    ! 6 end definitions: leave define mode
    call check(nf90_enddef(ncid))
 
    deallocate(chunksizes)
 
  end subroutine create_netcdf
 
 
  ! ------------------------------------------------------------------
 
  subroutine dump_netcdf_1d_sp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(sp), dimension(:), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 1 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim + 1) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim + 1) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 2) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim + 1) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim + 1) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim + 1) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim + 1)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 2)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 2), ndims = idim, dimids = dimid))
      if (idim /= ndim + 1) stop "dump_netcdf_1d_sp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim + 1
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim + 1) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_1d_sp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_1d_sp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_1d_sp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = dims(ndim + 1) + i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/dims(ndim + 1) + i/), (/dims(ndim + 1) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims(1 : ndim) = shape(arr)
      do i = 1, ndim
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      dims(ndim + 1) = 1
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim + 1)))
 
      ! define dim variables
      do i = 1, ndim
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim) = dims(1 : ndim)
        chunksizes(ndim + 1) = 1
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 2), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 2)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_1d_sp
 
 
  subroutine dump_netcdf_2d_sp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(sp), dimension(:, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 2 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim + 1) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim + 1) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 2) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim + 1) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim + 1) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim + 1) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim + 1)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 2)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 2), ndims = idim, dimids = dimid))
      if (idim /= ndim + 1) stop "dump_netcdf_2d_sp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim + 1
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim + 1) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_2d_sp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_2d_sp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_2d_sp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = dims(ndim + 1) + i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/dims(ndim + 1) + i/), (/dims(ndim + 1) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims(1 : ndim) = shape(arr)
      do i = 1, ndim
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      dims(ndim + 1) = 1
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim + 1)))
 
      ! define dim variables
      do i = 1, ndim
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim) = dims(1 : ndim)
        chunksizes(ndim + 1) = 1
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 2), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 2)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_2d_sp
 
 
  subroutine dump_netcdf_3d_sp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(sp), dimension(:, :, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 3 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 1) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 1)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim /= ndim) stop "dump_netcdf_3d_sp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_3d_sp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_3d_sp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_3d_sp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, size(arr, ndim)
        start(ndim) = dims(ndim) + i
        call check(nf90_put_var(ncid, varid(ndim), (/dims(ndim) + i/), (/dims(ndim) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, i), start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims = shape(arr)
      do i = 1, ndim - 1
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim)))
 
      ! define dim variables
      do i = 1, ndim - 1
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim - 1) = dims(1 : ndim - 1)
        chunksizes(ndim) = 1
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 1), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 1)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim - 1
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, dims(ndim)
        start(ndim) = i
        call check(nf90_put_var(ncid, varid(ndim), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, i), start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_3d_sp
 
 
  subroutine dump_netcdf_4d_sp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(sp), dimension(:, :, :, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 4 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 1) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 1)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim /= ndim) stop "dump_netcdf_4d_sp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_4d_sp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_4d_sp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_4d_sp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, size(arr, ndim)
        start(ndim) = dims(ndim) + i
        call check(nf90_put_var(ncid, varid(ndim), (/dims(ndim) + i/), (/dims(ndim) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, i), start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims = shape(arr)
      do i = 1, ndim - 1
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim)))
 
      ! define dim variables
      do i = 1, ndim - 1
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim - 1) = dims(1 : ndim - 1)
        chunksizes(ndim) = 1
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 1), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 1)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim - 1
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, dims(ndim)
        start(ndim) = i
        call check(nf90_put_var(ncid, varid(ndim), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, i), start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_4d_sp
 
 
  subroutine dump_netcdf_5d_sp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(sp), dimension(:, :, :, :, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 5 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 1) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 1)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim /= ndim) stop "dump_netcdf_5d_sp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_5d_sp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_5d_sp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_5d_sp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, size(arr, ndim)
        start(ndim) = dims(ndim) + i
        call check(nf90_put_var(ncid, varid(ndim), (/dims(ndim) + i/), (/dims(ndim) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, :, i), start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims = shape(arr)
      do i = 1, ndim - 1
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim)))
 
      ! define dim variables
      do i = 1, ndim - 1
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim - 1) = dims(1 : ndim - 1)
        chunksizes(ndim) = 1
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 1), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_float, dimid, varid(ndim + 1)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim - 1
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, dims(ndim)
        start(ndim) = i
        call check(nf90_put_var(ncid, varid(ndim), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, :, i), start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_5d_sp
 
 
  subroutine dump_netcdf_1d_dp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(dp), dimension(:), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 1 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim + 1) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim + 1) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 2) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim + 1) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim + 1) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim + 1) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim + 1)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 2)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 2), ndims = idim, dimids = dimid))
      if (idim /= ndim + 1) stop "dump_netcdf_1d_dp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim + 1
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim + 1) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_1d_dp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_1d_dp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_1d_dp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = dims(ndim + 1) + i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/dims(ndim + 1) + i/), (/dims(ndim + 1) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims(1 : ndim) = shape(arr)
      do i = 1, ndim
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      dims(ndim + 1) = 1
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim + 1)))
 
      ! define dim variables
      do i = 1, ndim
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim) = dims(1 : ndim)
        chunksizes(ndim + 1) = 1
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 2), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 2)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_1d_dp
 
 
  subroutine dump_netcdf_2d_dp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(dp), dimension(:, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 2 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim + 1) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim + 1) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 2) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim + 1) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim + 1) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim + 1) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim + 1)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 2)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 2), ndims = idim, dimids = dimid))
      if (idim /= ndim + 1) stop "dump_netcdf_2d_dp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim + 1
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim + 1) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_2d_dp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_2d_dp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_2d_dp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = dims(ndim + 1) + i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/dims(ndim + 1) + i/), (/dims(ndim + 1) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims(1 : ndim) = shape(arr)
      do i = 1, ndim
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      dims(ndim + 1) = 1
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim + 1)))
 
      ! define dim variables
      do i = 1, ndim
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim) = dims(1 : ndim)
        chunksizes(ndim + 1) = 1
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 2), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 2)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_2d_dp
 
 
  subroutine dump_netcdf_3d_dp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(dp), dimension(:, :, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 3 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 1) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 1)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim /= ndim) stop "dump_netcdf_3d_dp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_3d_dp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_3d_dp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_3d_dp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, size(arr, ndim)
        start(ndim) = dims(ndim) + i
        call check(nf90_put_var(ncid, varid(ndim), (/dims(ndim) + i/), (/dims(ndim) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, i), start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
        ! call check(nf90_set_fill(ncid, NF90_NOFILL, old_fill_mode))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims = shape(arr)
      do i = 1, ndim - 1
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim)))
 
      ! define dim variables
      do i = 1, ndim - 1
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim - 1) = dims(1 : ndim - 1)
        chunksizes(ndim) = 1
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 1), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 1)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim - 1
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, dims(ndim)
        start(ndim) = i
        call check(nf90_put_var(ncid, varid(ndim), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, i), start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_3d_dp
 
 
  subroutine dump_netcdf_4d_dp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(dp), dimension(:, :, :, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 4 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 1) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 1)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim /= ndim) stop "dump_netcdf_4d_dp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_4d_dp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_4d_dp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_4d_dp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, size(arr, ndim)
        start(ndim) = dims(ndim) + i
        call check(nf90_put_var(ncid, varid(ndim), (/dims(ndim) + i/), (/dims(ndim) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, i), start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims = shape(arr)
      do i = 1, ndim - 1
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim)))
 
      ! define dim variables
      do i = 1, ndim - 1
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim - 1) = dims(1 : ndim - 1)
        chunksizes(ndim) = 1
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 1), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 1)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim - 1
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, dims(ndim)
        start(ndim) = i
        call check(nf90_put_var(ncid, varid(ndim), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, i), start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_4d_dp
 
 
  subroutine dump_netcdf_5d_dp(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    real(dp), dimension(:, :, :, :, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 5 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 1) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 1)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim /= ndim) stop "dump_netcdf_5d_dp: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_5d_dp: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_5d_dp: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_5d_dp: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, size(arr, ndim)
        start(ndim) = dims(ndim) + i
        call check(nf90_put_var(ncid, varid(ndim), (/dims(ndim) + i/), (/dims(ndim) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, :, i), start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims = shape(arr)
      do i = 1, ndim - 1
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim)))
 
      ! define dim variables
      do i = 1, ndim - 1
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim - 1) = dims(1 : ndim - 1)
        chunksizes(ndim) = 1
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 1), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_double, dimid, varid(ndim + 1)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim - 1
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, dims(ndim)
        start(ndim) = i
        call check(nf90_put_var(ncid, varid(ndim), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, :, i), start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_5d_dp
 
 
  subroutine dump_netcdf_1d_i4(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    integer(i4), dimension(:), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 1 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim + 1) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim + 1) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 2) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim + 1) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim + 1) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim + 1) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim + 1)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 2)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 2), ndims = idim, dimids = dimid))
      if (idim /= ndim + 1) stop "dump_netcdf_1d_i4: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim + 1
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim + 1) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_1d_i4: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_1d_i4: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_1d_i4: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = dims(ndim + 1) + i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/dims(ndim + 1) + i/), (/dims(ndim + 1) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims(1 : ndim) = shape(arr)
      do i = 1, ndim
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      dims(ndim + 1) = 1
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim + 1)))
 
      ! define dim variables
      do i = 1, ndim
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim) = dims(1 : ndim)
        chunksizes(ndim + 1) = 1
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 2), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 2)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_1d_i4
 
 
  subroutine dump_netcdf_2d_i4(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    integer(i4), dimension(:, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 2 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim + 1) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim + 1) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 2) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim + 1) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim + 1) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim + 1) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim + 1)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 2)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 2), ndims = idim, dimids = dimid))
      if (idim /= ndim + 1) stop "dump_netcdf_2d_i4: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim + 1
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim + 1) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_2d_i4: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_2d_i4: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_2d_i4: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = dims(ndim + 1) + i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/dims(ndim + 1) + i/), (/dims(ndim + 1) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims(1 : ndim) = shape(arr)
      do i = 1, ndim
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      dims(ndim + 1) = 1
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim + 1)))
 
      ! define dim variables
      do i = 1, ndim
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim + 1), varid(ndim + 1)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim) = dims(1 : ndim)
        chunksizes(ndim + 1) = 1
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 2), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 2)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim + 1) = 1
      do i = 1, 1
        start(ndim + 1) = i
        call check(nf90_put_var(ncid, varid(ndim + 1), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 2), arr, start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_2d_i4
 
 
  subroutine dump_netcdf_3d_i4(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    integer(i4), dimension(:, :, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 3 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 1) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 1)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim /= ndim) stop "dump_netcdf_3d_i4: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_3d_i4: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_3d_i4: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_3d_i4: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, size(arr, ndim)
        start(ndim) = dims(ndim) + i
        call check(nf90_put_var(ncid, varid(ndim), (/dims(ndim) + i/), (/dims(ndim) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, i), start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims = shape(arr)
      do i = 1, ndim - 1
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim)))
 
      ! define dim variables
      do i = 1, ndim - 1
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim - 1) = dims(1 : ndim - 1)
        chunksizes(ndim) = 1
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 1), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 1)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim - 1
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, dims(ndim)
        start(ndim) = i
        call check(nf90_put_var(ncid, varid(ndim), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, i), start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_3d_i4
 
 
  subroutine dump_netcdf_4d_i4(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    integer(i4), dimension(:, :, :, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 4 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 1) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 1)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim /= ndim) stop "dump_netcdf_4d_i4: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_4d_i4: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_4d_i4: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_4d_i4: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, size(arr, ndim)
        start(ndim) = dims(ndim) + i
        call check(nf90_put_var(ncid, varid(ndim), (/dims(ndim) + i/), (/dims(ndim) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, i), start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims = shape(arr)
      do i = 1, ndim - 1
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim)))
 
      ! define dim variables
      do i = 1, ndim - 1
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim - 1) = dims(1 : ndim - 1)
        chunksizes(ndim) = 1
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 1), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 1)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim - 1
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, dims(ndim)
        start(ndim) = i
        call check(nf90_put_var(ncid, varid(ndim), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, i), start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_4d_i4
 
 
  subroutine dump_netcdf_5d_i4(filename, arr, append, lfs, netcdf4, deflate_level)
 
    implicit none
 
    character(len = *), intent(in) :: filename ! netcdf file name
    integer(i4), dimension(:, :, :, :, :), intent(in) :: arr      ! input array
    logical, optional, intent(in) :: append   ! append to existing file
    logical, optional, intent(in) :: lfs           ! netcdf3 Large File Support
    logical, optional, intent(in) :: netcdf4       ! netcdf4
    integer(i4), optional, intent(in) :: deflate_level ! compression level in netcdf4
 
    integer(i4), parameter :: ndim = 5 ! Routine for ndim dimensional array
    character(len = 1), dimension(4) :: dnames   ! Common dimension names
    integer(i4), dimension(ndim) :: dims     ! Size of each dimension
    integer(i4), dimension(ndim) :: dimid    ! netcdf IDs of each dimension
    integer(i4), dimension(ndim + 1) :: varid    ! dimension variables and var id
    integer(i4), dimension(ndim) :: start    ! start array for write of each time step
    integer(i4), dimension(ndim) :: counter  ! length array for write of each time step
    integer(i4), dimension(ndim) :: chunksizes ! Size of chunks in netcdf4 writing
    integer(i4) :: ncid                             ! netcdf file id
    integer(i4) :: i, j
    logical :: iappend
    integer(i4) :: idim    ! read dimension on append
    character(NF90_MAX_NAME) :: name    ! name of dimensions from nf90_inquire_dimension
    logical :: largefile
    logical :: inetcdf4
    integer(i4) :: deflate
    integer(i4) :: buffersize
 
    ! append or not
    if (present(append)) then
      if (append) then
        iappend = .true.
      else
        iappend = .false.
      end if
    else
      iappend = .false.
    end if
    largefile = .false.
    if (present(lfs)) largefile = lfs
    inetcdf4 = .false.
    if (present(netcdf4)) inetcdf4 = netcdf4
    deflate = 1
    if (present(deflate_level)) deflate = deflate_level
 
    ! dimension names
    dnames(1 : 4) = (/ 'x', 'y', 'z', 'l' /)
 
    if (iappend) then
      ! open file
      call check(nf90_open(trim(filename), nf90_write, ncid))
 
      ! inquire variables time and var
      call check(nf90_inq_varid(ncid, 'time', varid(ndim)))
      call check(nf90_inq_varid(ncid, 'var', varid(ndim + 1)))
      call check(nf90_inquire_variable(ncid, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim /= ndim) stop "dump_netcdf_5d_i4: number of variable dimensions /= number of file variable dimensions."
 
      ! inquire dimensions
      do i = 1, ndim
        call check(nf90_inquire_dimension(ncid, dimid(i), name, dims(i)))
        if (i < ndim) then
          if (trim(name) /= dnames(i)) stop "dump_netcdf_5d_i4: dimension name problem."
          if (dims(i) /= size(arr, i)) stop "dump_netcdf_5d_i4: variable dimension /= file variable dimension."
        else
          if (trim(name) /= 'time') stop "dump_netcdf_5d_i4: time name problem."
        end if
      enddo
 
      ! append
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, size(arr, ndim)
        start(ndim) = dims(ndim) + i
        call check(nf90_put_var(ncid, varid(ndim), (/dims(ndim) + i/), (/dims(ndim) + i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, :, i), start, counter))
      end do
    else
      ! open file
      if (inetcdf4) then
        call check(nf90_create(trim(filename), nf90_netcdf4, ncid))
      else
        if (largefile) then
          call check(nf90_create(trim(filename), nf90_64bit_offset, ncid, chunksize = buffersize))
        else
          call check(nf90_create(trim(filename), nf90_clobber, ncid, chunksize = buffersize))
        end if
      end if
 
      ! define dims
      dims = shape(arr)
      do i = 1, ndim - 1
        call check(nf90_def_dim(ncid, dnames(i), dims(i), dimid(i)))
      end do
      ! define dim time
      call check(nf90_def_dim(ncid, 'time', nf90_unlimited, dimid(ndim)))
 
      ! define dim variables
      do i = 1, ndim - 1
        if (inetcdf4) then
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        else
          call check(nf90_def_var(ncid, dnames(i), nf90_int, dimid(i), varid(i)))
        end if
      end do
      ! define time variable
      if (inetcdf4) then
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      else
        call check(nf90_def_var(ncid, 'time', nf90_int, dimid(ndim), varid(ndim)))
      end if
 
      ! define variable
      if (inetcdf4) then
        chunksizes(1 : ndim - 1) = dims(1 : ndim - 1)
        chunksizes(ndim) = 1
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 1), &
                chunksizes = chunksizes, shuffle = .true., deflate_level = deflate))
      else
        call check(nf90_def_var(ncid, 'var', nf90_int, dimid, varid(ndim + 1)))
      end if
 
      ! end define mode
      call check(nf90_enddef(ncid))
 
      ! write dimensions
      do i = 1, ndim - 1
        call check(nf90_put_var(ncid, varid(i), (/ (j, j = 1, dims(i)) /)))
      end do
 
      ! write time and variable
      start(:) = 1
      counter(:) = dims
      counter(ndim) = 1
      do i = 1, dims(ndim)
        start(ndim) = i
        call check(nf90_put_var(ncid, varid(ndim), (/i/), (/i/)))
        call check(nf90_put_var(ncid, varid(ndim + 1), arr(:, :, :, :, i), start, counter))
      end do
    end if
 
    ! close netcdf file
    call check(nf90_close(ncid))
 
  end subroutine dump_netcdf_5d_i4
 
 
  ! ------------------------------------------------------------------
 
 
  subroutine var2nc_1d_i4(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 1
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    integer(i4), dimension(:), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! attributes
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    character(256), dimension(:, :), optional, intent(in) :: attributes
    integer(i4), optional, intent(in) :: missing_value
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4) :: f_handle
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    integer(i4), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((size(dnames, 1) .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (size(dnames, 1) .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1) /)
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1
    counter(:) = dims
    dummy = nf90_fill_int
    dummy_count = 1
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim .ne. ndim) stop "var2nc_1d_i4: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_1d_i4: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_1d_i4: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (dummy(1) /= nf90_fill_int) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimension
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_int, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      !
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(I6)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! inquire dimensions
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_1d_i4: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_1d_i4: variable dimension /= file variable dimension."
    enddo
    ! write time and variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_1d_i4
 
  subroutine var2nc_1d_sp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 1
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(sp), dimension(:), intent(in) :: arr
    character(len = *), dimension(:), intent(in) :: dnames
    character(len = *), intent(in) :: v_name
    ! attributes
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    character(256), dimension(:, :), optional, intent(in) :: attributes
    real(sp), optional, intent(in) :: missing_value
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4) :: f_handle
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(sp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((size(dnames, 1) .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (size(dnames, 1) .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1) /)
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1
    counter(:) = dims
    dummy_count = 1
    dummy = nf90_fill_float
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim .ne. ndim) stop "var2nc_1d_sp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_1d_sp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_1d_sp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_float)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimension
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_float, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! inquire dimensions
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_1d_sp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_1d_sp: variable dimension /= file variable dimension."
    enddo
    ! write time and variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_1d_sp
 
  subroutine var2nc_1d_dp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 1
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(dp), dimension(:), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! attributes
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    character(256), dimension(:, :), optional, intent(in) :: attributes
    real(dp), optional, intent(in) :: missing_value
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4) :: f_handle
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(dp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((size(dnames, 1) .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (size(dnames, 1) .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1) /)
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1
    counter(:) = dims
    dummy_count = 1
    dummy = nf90_fill_double
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      if (idim .ne. ndim) stop "var2nc_1d_dp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_1d_dp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_1d_dp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_double)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimension
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_double, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! inquire dimensions
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_1d_dp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_1d_dp: variable dimension /= file variable dimension."
    enddo
    ! write time and variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_1d_dp
 
  subroutine var2nc_2d_i4(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 2
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    integer(i4), dimension(:, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    integer(i4), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    integer(i4), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((ndim .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (ndim .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), size(arr, 2), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1), size(arr, 2) /)
      if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1_i4
    dummy = nf90_fill_int
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_2d_i4: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_2d_i4: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_2d_i4: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (dummy(1) /= nf90_fill_int) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_int, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(I6)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_2d_i4: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_2d_i4: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_2d_i4
 
  subroutine var2nc_2d_sp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 2
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(sp), dimension(:, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    real(sp), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(sp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((ndim .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (ndim .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), size(arr, 2), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1), size(arr, 2) /)
      if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1_i4
    dummy = nf90_fill_float
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_2d_sp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_2d_sp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_2d_sp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_float)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_float, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_2d_sp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_2d_sp: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_2d_sp
 
  subroutine var2nc_2d_dp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 2
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(dp), dimension(:, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    real(dp), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(dp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((ndim .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (ndim .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), size(arr, 2), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1), size(arr, 2) /)
      if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1
    dummy = nf90_fill_double
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_2d_dp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_2d_dp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_2d_dp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_double)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_double, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_2d_dp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_2d_dp: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_2d_dp
 
  subroutine var2nc_3d_i4(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 3
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    integer(i4), dimension(:, :, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    integer(i4), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    integer(i4), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((ndim .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (ndim .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), size(arr, 2), size(arr, 3), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1), size(arr, 2), size(arr, 3) /)
      if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1_i4
    dummy = nf90_fill_int
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_3d_i4: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_3d_i4: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_3d_i4: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (dummy(1) /= nf90_fill_int) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_int, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(I6)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_3d_i4: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_3d_i4: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_3d_i4
 
  subroutine var2nc_3d_sp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 3
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(sp), dimension(:, :, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    real(sp), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(sp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((ndim .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (ndim .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    ! set chunk sizes and dimension names
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), size(arr, 2), size(arr, 3), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1), size(arr, 2), size(arr, 3) /)
      if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1_i4
    dummy = nf90_fill_float
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_3d_sp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_3d_sp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_3d_sp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_float)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_float, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_3d_sp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_3d_sp: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_3d_sp
 
  subroutine var2nc_3d_dp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 3
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(dp), dimension(:, :, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    real(dp), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(dp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((ndim .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (ndim .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), size(arr, 2), size(arr, 3), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1), size(arr, 2), size(arr, 3) /)
      if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1
    dummy = nf90_fill_double
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_3d_dp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_3d_dp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_3d_dp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_double)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_double, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_3d_dp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_3d_dp: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_3d_dp
 
  subroutine var2nc_4d_i4(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 4
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    integer(i4), dimension(:, :, :, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    integer(i4), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    integer(i4), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((ndim .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (ndim .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), size(arr, 2), &
              size(arr, 3), size(arr, 4), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1), size(arr, 2), &
              size(arr, 3), size(arr, 4) /)
      if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1_i4
    dummy = nf90_fill_int
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_4d_i4: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_4d_i4: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_4d_sp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (dummy(1) /= nf90_fill_int) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_int, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(I6)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_4d_i4: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_4d_i4: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_4d_i4
 
  subroutine var2nc_4d_sp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 4
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(sp), dimension(:, :, :, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    real(sp), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(sp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((ndim .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (ndim .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), size(arr, 2), &
              size(arr, 3), size(arr, 4), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1), size(arr, 2), &
              size(arr, 3), size(arr, 4) /)
      if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1_i4
    dummy = nf90_fill_float
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_4d_sp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_4d_sp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_4d_sp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_float)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_float, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_4d_sp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_4d_sp: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_4d_sp
 
  subroutine var2nc_4d_dp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 4
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(dp), dimension(:, :, :, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    real(dp), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(dp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    ! consistency checks
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    if ((ndim .eq. ndim_const + 1) .and. (d_unlimit .ne. ndim_const + 1)) then
      print *, '***ERROR one more dimension name specified than dimension of array, but the last one is not unlimited'
      stop '***ERROR see StdOut'
    end if
    if (ndim .gt. ndim_const + 1) then
      print *, '***ERROR too many dimension name specified, should be atmost ndim_const + 1'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! Initialize
    deflate = 1
    if (ndim .gt. ndim_const) then
      chunksizes = (/ size(arr, 1), size(arr, 2), &
              size(arr, 3), size(arr, 4), 1 /)
      dims(1 : ndim - 1) = shape(arr)
      dims(ndim) = 1
    else
      chunksizes = (/ size(arr, 1), size(arr, 2), &
              size(arr, 3), size(arr, 4) /)
      if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
      dims(1 : ndim_const) = shape(arr)
    end if
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1
    dummy = nf90_fill_double
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_4d_dp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_4d_dp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_4d_dp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_double)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_double, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_4d_dp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_4d_dp: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_4d_dp
 
  subroutine var2nc_5d_i4(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 5
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    integer(i4), dimension(:, :, :, :, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    integer(i4), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    integer(i4), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    ! consistency checks
    if (ndim .gt. ndim_const) then
      print *, '***ERROR more than five dimension names given'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    chunksizes = (/ size(arr, 1), size(arr, 2), size(arr, 3), size(arr, 4), size(arr, 5) /)
    if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
    dims(1 : ndim) = shape(arr)
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1_i4
    dummy = nf90_fill_int
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_5d_i4: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_5d_i4: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_5d_sp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (dummy(1) /= nf90_fill_int) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_int, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(I6)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_5d_i4: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_5d_i4: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_5d_i4
 
  subroutine var2nc_5d_sp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 5
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(sp), dimension(:, :, :, :, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    real(sp), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(sp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    ! consistency checks
    if (ndim .gt. ndim_const) then
      print *, '***ERROR more than five dimension names given'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    chunksizes = (/ size(arr, 1), size(arr, 2), size(arr, 3), size(arr, 4), size(arr, 5) /)
    if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
    dims(1 : ndim) = shape(arr)
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1_i4
    dummy = nf90_fill_float
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_5d_sp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_5d_sp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_5d_sp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_float)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_float, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_5d_sp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_5d_sp: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_5d_sp
 
  subroutine var2nc_5d_dp(f_name, arr, dnames, v_name, dim_unlimited, &
          long_name, units, missing_value, attributes, create, ncid, nrec)
    !
    implicit none
    !
    integer(i4), parameter :: ndim_const = 5
    integer(i4) :: ndim
    ! input variables
    character(len = *), intent(in) :: f_name
    real(dp), dimension(:, :, :, :, :), intent(in) :: arr
    character(len = *), intent(in) :: v_name
    character(len = *), dimension(:), intent(in) :: dnames
    ! optional
    integer(i4), optional, intent(in) :: dim_unlimited
    character(len = *), optional, intent(in) :: long_name
    character(len = *), optional, intent(in) :: units
    real(dp), optional, intent(in) :: missing_value
    character(256), dimension(:, :), optional, intent(in) :: attributes
    logical, optional, intent(in) :: create
    integer(i4), optional, intent(inout) :: ncid
    integer(i4), optional, intent(in) :: nrec
    ! local variables
    logical :: create_loc
    character(256) :: dummy_name
    integer(i4) :: deflate
    integer(i4), dimension(:), allocatable :: chunksizes
    integer(i4), dimension(:), allocatable :: start     ! start array for write
    integer(i4), dimension(:), allocatable :: counter   ! length array for write
    integer(i4) :: idim      ! read dimension on append
    integer(i4) :: f_handle
    integer(i4) :: d_unlimit ! index of unlimited dimension
    integer(i4) :: u_dimid   ! dimid of unlimited dimension
    integer(i4) :: u_len     ! length of unlimited dimension
    integer(i4), dimension(:), allocatable :: dims
    integer(i4), dimension(:), allocatable :: dimid     ! netcdf IDs of each dimension
    integer(i4), dimension(:), allocatable :: varid     ! dimension variables and var id
    integer(i4) :: i         ! loop indices
    integer(i4), dimension(:), allocatable :: dummy_count
    real(dp), dimension(1) :: dummy     ! dummy read
    logical :: openfile  ! tmp logical
    !
    ndim = size(dnames, 1)
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    ! consistency checks
    if (ndim .gt. ndim_const) then
      print *, '***ERROR more than five dimension names given'
      stop '***ERROR see StdOut'
    end if
    if  (((ndim .eq. ndim_const) .and. (d_unlimit .gt. ndim_const)) .or. &
            (d_unlimit .lt. 0_i4)) then
      print*, '***ERROR unlimited dimension out of bounds, must be positive but not greater than number of given dimensions'
      print*, '***Dims: ', ndim, ndim_const, d_unlimit, ndim_const, d_unlimit
      stop '***ERROR see StdOut'
    end if
    !
    allocate(chunksizes(ndim))
    allocate(start(ndim))
    allocate(counter(ndim))
    allocate(dims(ndim))
    allocate(dimid(ndim))
    allocate(varid(1 + ndim))
    allocate(dummy_count(ndim))
    ! initialize
    deflate = 1
    chunksizes = (/ size(arr, 1), size(arr, 2), size(arr, 3), size(arr, 4), size(arr, 5) /)
    if (d_unlimit .gt. 0) chunksizes(d_unlimit) = 1
    dims(1 : ndim) = shape(arr)
    start(:) = 1_i4
    counter(:) = dims
    dummy_count = 1
    dummy = nf90_fill_double
    d_unlimit = 0_i4
    if (present(dim_unlimited)) d_unlimit = dim_unlimited
    ! open the netcdf file
    if (present(ncid)) then
      if (ncid < 0_i4) then
        openfile = .true.
      else
        openfile = .false.
        f_handle = ncid
      end if
    else
      openfile = .true.
    end if
    if (openfile) then
      create_loc = .false.
      if (present(create)) create_loc = create
      f_handle = open_netcdf(f_name, create = create_loc)
    end if
    ! check whether variable exists
    if (nf90_noerr .eq. nf90_inq_varid(f_handle, v_name, varid(ndim + 1))) then
      ! append
      call check(nf90_inquire_variable(f_handle, varid(ndim + 1), ndims = idim, dimids = dimid))
      ! consistency checks
      if (idim .ne. ndim) stop "var2nc_5d_dp: number of variable dimensions /= number of file variable dimensions."
      ! check unlimited dimension
      call check(nf90_inquire(f_handle, unlimiteddimid = u_dimid))
      if (u_dimid .eq. -1) stop 'var2nc_5d_dp: cannot append, no unlimited dimension defined'
      ! check for unlimited dimension
      if (dimid(d_unlimit) .ne. u_dimid) stop 'var2nc_5d_dp: unlimited dimension not specified correctly'
      if (present(nrec)) then
        start(d_unlimit) = nrec
      else
        ! get length of unlimited dimension
        call check(nf90_inquire_dimension(f_handle, u_dimid, len = u_len))
        ! adapt start, that is find last written chunk
        do i = u_len, 1, -1
          if (ne(dummy(1), nf90_fill_double)) exit
          start(d_unlimit) = i
          call check(nf90_get_var(f_handle, varid(ndim + 1), dummy, start, dummy_count))
        end do
        start(d_unlimit) = start(d_unlimit) + 1
      end if
    else
      ! define dimensions
      do i = 1, ndim
        ! check whether dimension exists
        if (nf90_noerr .ne. nf90_inq_dimid(f_handle, dnames(i), dimid(i))) then
          ! create dimension
          if (i .eq. d_unlimit) then
            ! define unlimited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), nf90_unlimited, dimid(i)))
          else
            ! define limited dimension
            call check(nf90_def_dim(f_handle, trim(dnames(i)), dims(i), dimid(i)))
          end if
        end if
      end do
      ! define variable
      call check(nf90_def_var(f_handle, v_name, nf90_double, dimid, varid(ndim + 1), &
              chunksizes = chunksizes(:), shuffle = .true., deflate_level = deflate))
      ! add attributes
      if (present(long_name)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'long_name', long_name))
      if (present(units)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'units', units))
      if (present(missing_value)) call check(nf90_put_att(f_handle, varid(ndim + 1), 'missing_value', missing_value))
      if (present(attributes)) then
        do i = 1, size(attributes, dim = 1)
          if (trim(attributes(i, 1)) .eq. 'missing_value') then
            ! write number
            read(attributes(i, 2), '(F10.2)') dummy(1)
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), dummy(1)))
          else
            ! write string
            call check(nf90_put_att(f_handle, varid(ndim + 1), &
                    trim(attributes(i, 1)), trim(attributes(i, 2))))
          end if
        end do
      end if
      ! end definition
      call check(nf90_enddef(f_handle))
    end if
    ! check dimensions before writing
    do i = 1, ndim_const
      call check(nf90_inquire_dimension(f_handle, dimid(i), dummy_name, dims(i)))
      if (trim(dummy_name) .ne. dnames(i)) &
              stop "var2nc_5d_dp: dimension name problem."
      if ((dims(i) .ne. size(arr, i)) .and. (d_unlimit .ne. i)) &
              stop "var2nc_5d_dp: variable dimension /= file variable dimension."
    end do
    ! write variable
    call check(nf90_put_var(f_handle, varid(ndim + 1), arr, start, counter))
    ! close netcdf_dataset
    if (present(ncid)) then
      if (ncid < 0_i4) ncid = f_handle
    else
      call close_netcdf(f_handle)
    end if
    !
  end subroutine var2nc_5d_dp
 
  ! ----------------------------------------------------------------------------
  !>     \brief dynamic writer
  !>     \details This routine writes data, where one dimension has the unlimited attribute.
  !!     Therefore, the number of the record which should be written has to be
  !!     specified.
  !!     See: http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90.html
 
  !>     \param[in] "integer(i4) :: nc"  stream id of an open netcdf file where data should be written
  !!                                     can be obtained by an create_netcdf call
  !>     \param[in] "integer(i4) :: rec" record id of record which will be written in the file
 
  !>     \note Writes only data, where the data pointers of the structure V are assigned
  !!     and where one dimension has the unlimited attribute. Moreover only one
  !!     record will be written.
  !!     Writes only 1 to 4 dim arrays, integer, single or double precision.
 
  !>     \author Luis Samaniego
  !>     \date Feb 2011
  !>     \author Stephan Thober
  !>     \date Dec 2011
  !!           - added comments and generalized
  !>     \author Matthias Cuntz
  !>     \date Jan 2012
  !!           - Info
  !>     \author Stephan Thober
  !>     \date Jan 2012
  !!           - iRec is not optional
  !>     \author Matthias Cuntz
  !>     \date Mar 2013
  !!           - removed Info
  subroutine write_dynamic_netcdf(ncId, irec)
 
    implicit none
 
    ! netcdf related variables
    integer(i4), intent(in) :: ncid
    integer(i4), intent(in) :: irec
 
    integer(i4) :: i
    ! NOTES: 1) netcdf file must be on *** data mode ***
    !        2) start and end of the data chuck is controled by
    !           V(:)%start and  V(:)%count
 
    ! set values for variables (one scalar or grid at a time)
 
    do i = 1, nvars
      if (.not. v(i)%unlimited) cycle
      if (.not. v(i)%wFlag) cycle
      v(i)%start (v(i)%nDims) = irec
      select case (v(i)%xtype)
      case(nf90_byte)
        select case (v(i)%nDims)
        case (0)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G0_b, v(i)%start))
        case (1)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G1_b, v(i)%start, v(i)%count))
        case (2)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G2_b, v(i)%start, v(i)%count))
        case (3)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G3_b, v(i)%start, v(i)%count))
        case (4)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G4_b, v(i)%start, v(i)%count))
        end select
      case (nf90_int)
        select case (v(i)%nDims - 1)
        case (0)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G0_i, v(i)%start))
        case (1)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G1_i, v(i)%start, v(i)%count))
        case (2)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G2_i, v(i)%start, v(i)%count))
        case (3)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G3_i, v(i)%start, v(i)%count))
        case (4)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G4_i, v(i)%start, v(i)%count))
        end select
      case (nf90_float)
        select case (v(i)%nDims - 1)
        case (0)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G0_f, v(i)%start))
        case (1)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G1_f, v(i)%start, v(i)%count))
        case (2)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G2_f, v(i)%start, v(i)%count))
        case (3)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G3_f, v(i)%start, v(i)%count))
        case (4)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G4_f, v(i)%start, v(i)%count))
        end select
      case (nf90_double)
        select case (v(i)%nDims - 1)
        case (0)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G0_d, v(i)%start))
        case (1)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G1_d, v(i)%start, v(i)%count))
        case (2)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G2_d, v(i)%start, v(i)%count))
        case (3)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G3_d, v(i)%start, v(i)%count))
        case (4)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G4_d, v(i)%start, v(i)%count))
        end select
      end select
    end do
 
  end subroutine write_dynamic_netcdf
 
  ! ----------------------------------------------------------------------------
  !>     \brief static writer
  !>     \details This routines writes static data in the netcdf file that is data
  !!     where no dimension has the unlimited attribute.
  !!     See: http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90.html
 
  !>     \param[in] "integer(i4) :: ncid" stream id of an open netcdf file where data should be written
  !!                                      can be obtained by an create_netcdf call
 
  !>     \note Writes only data, where the data pointers of the structure V are assigned.
  !!     Writes only 1 to 4 dim arrays, integer, single or double precision.
 
  ! HISTORY
  !>     \author Luis Samaniego
  !>     \date Feb 2011
  !>     \author Stephan Thober
  !>     \date Dec 2011
  !!           - added comments and generalized
  !>     \author Matthias Cuntz
  !>     \date Jan 2012
  !!           - Info
  !>     \author Matthias Cuntz
  !>     \date Mar 2013
  !!           - removed Info
 
  subroutine write_static_netcdf(ncId)
 
    implicit none
 
    ! netcdf related variables
    integer(i4), intent(in) :: ncid
 
    integer(i4) :: i
 
    ! write all static variables
    do i = 1, nvars
      if (v(i)%unlimited) cycle
      if (.not. v(i)%wFlag) cycle
      select case (v(i)%xtype)
      case(nf90_byte)
        select case (v(i)%nDims)
        case (0)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G0_b))
        case (1)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G1_b))
        case (2)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G2_b))
        case (3)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G3_b))
        case (4)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G4_b))
        end select
      case (nf90_int)
        select case (v(i)%nDims)
        case (0)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G0_i))
        case (1)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G1_i))
        case (2)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G2_i))
        case (3)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G3_i))
        case (4)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G4_i))
        end select
      case (nf90_float)
        select case (v(i)%nDims)
        case (0)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G0_f))
        case (1)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G1_f))
        case (2)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G2_f))
        case (3)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G3_f))
        case (4)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G4_f))
        end select
      case (nf90_double)
        select case (v(i)%nDims)
        case (0)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G0_d))
        case (1)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G1_d))
        case (2)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G2_d))
        case (3)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G3_d))
        case (4)
          call check(nf90_put_var(ncid, v(i)%varId, v(i)%G4_d))
        end select
      end select
    end do
 
  end subroutine write_static_netcdf
 
  ! -----------------------------------------------------------------------------
  ! PRIVATE PART
  !
 
  ! -----------------------------------------------------------------------------
 
  ! private open netcdf function - returns file handle
  function open_netcdf(f_name, create)
    implicit none
    ! input variables
    character(len = *), intent(in) :: f_name
    logical, intent(in) :: create ! flag indicates that file exists
    ! output variables
    integer(i4) :: open_netcdf
    !
    if (create) then
      ! create file
      call check(nf90_create(trim(f_name), nf90_netcdf4, open_netcdf))
    else
      ! open file
      call check(nf90_open(trim(f_name), nf90_write, open_netcdf))
    end if
  end function open_netcdf
 
 
  ! -----------------------------------------------------------------------------
 
  !  private error checking routine
  subroutine check(status)
 
    implicit none
 
    integer(i4), intent(in) :: status
 
    if (status /= nf90_noerr) then
      write(*, *) 'mo_ncwrite.check error: ', trim(nf90_strerror(status))
      stop
    end if
 
  end subroutine check
 
  ! -----------------------------------------------------------------------------
 
end module mo_ncwrite