netcdf4_support.F90

Go to the documentation of this file.

!> @file
!!  Contain a single module, @ref netcdf4_support, which
!!  provides support for use of netCDF files as input or output.
!!
!! Supports use of netCDF files as input for time-varying,
!! gridded meteorlogic data, or output for any SWB-generated variable.
!!
!! from the C API:
!! The @c nc_get_vars_l type family of functions read a subsampled (strided)
!! array section of values from a netCDF variable of an open netCDF dataset.
!! The subsampled array section is specified by giving a corner,
!! a vector of edge lengths, and a stride vector. The values are read
!! with the last dimension of the netCDF variable varying fastest.
!!          ^^^^                                  ^^^^^^^ ^^^^^^^
!!
!! from the Fortran 90 API:
!! The values to be read are associated with the netCDF variable by
!! assuming that the first dimension of the netCDF variable
!!                   ^^^^^
!! varies fastest in the Fortran 90 interface.
!! ^^^^^^ ^^^^^^^
 
!> Provide support for use of netCDF files as input for time-varying,
!! gridded meteorlogic data, or output for any SWB-generated variable.
module netcdf4_support
 
  use iso_c_binding
  use, intrinsic :: ieee_exceptions
  use constants_and_conversions
  use exceptions
  use logfiles
  use netcdf_c_api_interfaces
  use fstring
  use fstring_list, only           : fstring_list_t
  use proj4_support, only         : create_attributes_from_proj4_string
  use datetime
 
  use version_control, only  : swb_version, git_commit_hash_string,           &
                               git_branch_string, compile_date, compile_time, &
                               compilation_timestamp
 
  use grid
!  use typesizes
  use netcdf_c_api_interfaces
 
  implicit none
 
  private
 
  public :: nc_float, nc_fill_float, nc_fill_int
 
  integer(c_int), public :: nc_readonly          = 0
  integer(c_int), public :: nc_readwrite         = 1
 
  ! from netcdf.h:
 
  ! #define NC_NAT          0
  ! #define NC_BYTE         1
  ! #define NC_CHAR         2
  ! #define NC_SHORT        3
  ! #define NC_INT          4
  ! #define NC_LONG         NC_INT
  ! #define NC_FLOAT        5
  ! #define NC_DOUBLE       6
  ! #define NC_UBYTE        7
  ! #define NC_USHORT       8
  ! #define NC_UINT         9
  ! #define NC_INT64        10
  ! #define NC_UINT64       11
  ! #define NC_STRING       12
 
  integer(c_int), parameter ::  nc_nat    = 0
  integer(c_int), parameter ::  nc_byte   = 1
  integer(c_int), parameter ::  nc_char   = 2
  integer(c_int), parameter ::  nc_short  = 3
  integer(c_int), parameter ::  nc_int    = 4
  integer(c_int), parameter ::  nc_long   = nc_int
  integer(c_int), parameter ::  nc_float  = 5
  integer(c_int), parameter ::  nc_double = 6
  integer(c_int), parameter ::  nc_ubyte  = 7
  integer(c_int), parameter ::  nc_ushort = 8
  integer(c_int), parameter ::  nc_uint   = 9
  integer(c_int), parameter ::  nc_int64  = 10
  integer(c_int), parameter ::  nc_uint64 = 11
  integer(c_int), parameter ::  nc_string = 12
 
  integer(c_int), parameter :: nc_fill_char    = 0
  integer(c_int), parameter :: nc_fill_byte    = -127
  integer(c_int), parameter :: nc_fill_short   = -32767
  integer(c_int), parameter :: nc_fill_int     = -2147483647
  integer(c_int), parameter :: nc_na_int       = -9999
 
  real(c_float),  parameter :: nc_fill_float   = -9999.0_c_float
  real(c_double), parameter :: nc_fill_double  = -9.9e-20_c_double
!  real(c_float),  parameter :: NC_FILL_FLOAT   = -( HUGE( 0_c_float ) - 1.0_c_float )
!  real(c_double), parameter :: NC_FILL_DOUBLE  = -( HUGE( 0_c_double ) - 1.0_c_double )
 
  ! mode flags for opening and creating datasets
  integer(c_int), parameter :: nc_nowrite          = 0
  integer(c_int), parameter :: nc_write            = 1
  integer(c_int), parameter :: nc_clobber          = 0
  integer(c_int), parameter :: nc_noclobber        = 4
  integer(c_int), parameter :: nc_fill             = 0
  integer(c_int), parameter :: nc_nofill           = 256
  integer(c_int), parameter :: nc_lock             = 1024
  integer(c_int), parameter :: nc_share            = 2048
  integer(c_int), parameter :: nc_strict_nc3       = 8
  integer(c_int), parameter :: nc_64bit_offset     = 512
  integer(c_int), parameter :: nc_sizehint_default = 0
  integer(c_int), parameter :: nc_align_chunk      = -1
  integer(c_int), parameter :: nc_format_classic   = 1
  integer(c_int), parameter :: nc_format_64bit     = 2
  integer(c_int), parameter :: nc_format_netcdf4   = 3
  integer(c_int), parameter :: nc_format_netcdf4_classic = 4
 
  ! implementation limits (warning!  should be the same as c interface)
  integer(c_int), parameter :: nc_max_dims     = 1024
  integer(c_int), parameter :: nc_max_attrs    = 8192
  integer(c_int), parameter :: nc_max_vars     = 8192
  integer(c_int), parameter :: nc_max_name     = 256
 
  integer (c_int), parameter :: nc_shuffle_yes = 1
  integer (c_int), parameter :: nc_shuffle_no = 0
  integer (c_int), parameter :: nc_deflate_yes = 1
  integer (c_int), parameter :: nc_deflate_no = 0
 
  integer(c_int), parameter :: nc_netcdf4        = 4096
  integer(c_int), parameter :: nc_classic_model  = 256
 
 
  integer(c_int),  parameter :: nc_unlimited = 0
  integer(c_int),  parameter :: nc_global    = -1
 
  !> @TODO: implement a more flexible way of tracking
  !!        variable IDs; presently the code can break if
  !!        lat and lon are omitted, but time_bnds is included;
  !!        this is because if the lat and lon variables are
  !!        not defined, NC_TIME_BNDS may be 4 or 5 or some other value.
  integer (c_int), public, parameter :: nc_time      = 0
  integer (c_int), public, parameter :: nc_y         = 1
  integer (c_int), public, parameter :: nc_x         = 2
  integer (c_int), public, parameter :: nc_z         = 3
  integer (c_int), public, parameter :: nc_aux       = 3
  integer (c_int), public, parameter :: nc_crs       = 4
  integer (c_int), public, parameter :: nc_lat       = 5
  integer (c_int), public, parameter :: nc_lon       = 6
  integer (c_int), public            :: nc_time_bnds = 7
 
  integer (c_int), parameter :: nc_first = 0
  integer (c_int), parameter :: nc_last  = 1
  integer (c_int), parameter :: nc_by    = 2
 
  integer (c_int), public, parameter :: nc_left  = 0
  integer (c_int), public, parameter :: nc_right = 1
  integer (c_int), public, parameter :: nc_top    = 0
  integer (c_int), public, parameter :: nc_bottom = 1
 
  integer (c_int), parameter :: column = 1
  integer (c_int), parameter :: row = 2
 
  integer (c_int), parameter :: leap_year = 0
  integer (c_int), parameter :: noleap_year = 1
  integer (c_int), parameter :: year_is_360_days = 2 
 
  character (len=25), dimension(4), parameter :: netcdf_format_string = &
    ["NC_FORMAT_CLASSIC        ", &
     "NC_FORMAT_64BIT          ", &
     "NC_FORMAT_NETCDF4        ", &
     "NC_FORMAT_NETCDF4_CLASSIC" ]
 
  character (len=6), dimension(0:6), parameter :: netcdf_data_type = &
    ["nat   ", &
     "byte  ", &
     "char  ", &
     "short ", &
     "int   ", &
     "float ", &
     "double" ]
 
  type t_netcdf_dimension
    character (len=64)      :: sdimensionname
    integer (c_int)    :: inc_dimid = nc_na_int
    integer (c_size_t) :: inc_dimsize
    logical (c_bool)   :: lunlimited = false
  end type t_netcdf_dimension
 
  type t_netcdf_attribute
    character (len=64) :: sattributename
    character (len=256), dimension(:), allocatable    :: sattvalue
    integer (c_short), dimension(:), allocatable :: i2attvalue
    integer (c_int), dimension(:), allocatable   :: iattvalue
    real (c_float), dimension(:), allocatable    :: rattvalue
    real (c_double), dimension(:), allocatable   :: dpattvalue
    integer (c_int) :: inc_atttype
    integer (c_size_t) :: inc_attsize
  end type t_netcdf_attribute
 
  type t_netcdf_variable
    character (len=64) :: svariablename
    integer (c_int) :: inc_varid = nc_na_int
    integer (c_int) :: inc_vartype
    integer (c_int) :: inumberofdimensions
    integer (c_int), dimension(0:3) :: inc_dimid = nc_na_int
    integer (c_int) :: inumberofattributes
    type (t_netcdf_attribute), dimension(:), pointer :: pnc_att => null()
  end type t_netcdf_variable
 
  type t_netcdf4_file
    integer (c_int) :: incid
    character (len=256)  :: sfilename = ''
    integer (c_int) :: ifileformat = 0
    integer (c_int) :: inumberofdimensions = 0
    integer (c_int) :: inumberofvariables = 0
    integer (c_int) :: inumberofattributes = 0
    integer (c_int) :: inc3_unlimiteddimensionnumber = 0
    integer (c_int) :: ioriginjd   = nc_na_int
    integer (c_int) :: ifirstdayjd = nc_na_int
    integer (c_int) :: ilastdayjd  = nc_na_int
    integer (c_int) :: ioriginmonth = 0
    integer (c_int) :: ioriginday = 0
    integer (c_int) :: ioriginyear = 0
    integer (c_int) :: ioriginhh = 0
    integer (c_int) :: ioriginmm = 0
    integer (c_int) :: ioriginss = 0
    integer (c_int) :: lleapyeartreatment = leap_year
    integer (c_size_t), dimension(0:3) :: istart = 0
    integer (c_size_t), dimension(0:3) :: icount = 0
    integer (c_size_t), dimension(0:3) :: istride = 1
    integer (c_size_t), dimension(0:1) :: icolbounds = 0
    integer (c_size_t), dimension(0:1) :: irowbounds = 0
    integer (c_int) :: inx = 0
    integer (c_int) :: iny = 0
    character (len=3) :: svariableorder = "tyx"
    real (c_double), dimension(0:1) :: rx
    real (c_double), dimension(0:1) :: ry
    real (c_double)   :: rcoordinatetolerance = 0.0_c_double  ! set this to be > 0.0 to allow some 'slop' when comparing coordinates
    logical (c_bool)  :: lx_increaseswithindex = true
    logical (c_bool)  :: ly_increaseswithindex = false
    logical (c_bool)  :: lallowautomaticdataflipping = true
 
    real (c_double), dimension(0:1) :: dpfirstandlasttimevalues = 0.
    character (len=64), dimension(0:3) :: svarname = ["time","y   ","x   ","z   "]
    integer (c_int), dimension(0:3) :: ivarid = nc_na_int
    integer (c_int), dimension(0:3) :: ivarindex = nc_na_int
    integer (c_int), dimension(0:3) :: ivartype = nc_na_int
    character (len=64), dimension(0:3) :: svarunits = "NA"
    integer (c_int), dimension(0:3, 0:3) :: ivar_dimid = nc_na_int
    real (c_double), dimension(0:3) :: rscalefactor = 1.0_c_double
    real (c_double), dimension(0:3) :: raddoffset = 0.0_c_double
    integer (c_int), dimension(0:2) :: irowiter = 0
    integer (c_int), dimension(0:2) :: icoliter = 0
    logical (c_bool) :: lfliphorizontal = false
    logical (c_bool) :: lflipvertical = false
 
    real (c_double), allocatable, dimension(:) :: rx_coords
    real (c_double), allocatable, dimension(:) :: ry_coords
    real (c_double)                            :: rx_coord_addoffset = 0.0_c_double
    real (c_double)                            :: ry_coord_addoffset = 0.0_c_double
    real (c_double), allocatable, dimension(:) :: rdatetimevalues
    real (c_double) :: rgridcellsizex
    real (c_double) :: rgridcellsizey
 
    type (t_netcdf_dimension), dimension(:), pointer :: pnc_dim => null()
    type (t_netcdf_variable), dimension(:), pointer  :: pnc_var => null()
    type (t_netcdf_attribute), dimension(:), pointer :: pnc_att => null()
  end type t_netcdf4_file
 
!
! For a netCDF file to be read successfully, the reader must have a way to
! tie coordinate values to index values. By default the code assumes that coordinates
! decrease while the indices increase (i.e. read the file in as rows from top to bottom).
!
! If this is not the case, the data/coordinates must be 'flipped' relative to one another.
!
! DEFAULT ASSUMPTION REGARDING VERTICAL COORDINATES
!
!                  /\
!  coordinates     |   |
!  increase        |   |
!  in upward       |   |
!  direction       |   |
!                  |   |
!                  |   | column index increases in downward direction
!                  |   |
!                  |   |
!                      \/
!
 
  public :: t_netcdf_dimension, t_netcdf_variable, t_netcdf_attribute
  public :: t_netcdf4_file
 
  public :: netcdf_get_attribute_list_for_variable
  public :: netcdf_open_and_prepare_for_merging
  public :: netcdf_open_and_prepare_as_input
  public :: netcdf_open_and_prepare_as_output_archive
  public :: netcdf_open_and_prepare_as_output
  public :: netcdf_date_within_range
  public :: netcdf_deallocate_data_struct
  public :: netcdf_nullify_data_struct
  public :: netcdf_dump_cdl
  public :: netcdf_open_file
  public :: netcdf_close_file
  public :: netcdf_get_variable_list
  public :: netcdf_get_variable_slice
  public :: netcdf_update_time_starting_index
  public :: netcdf_put_variable_array
  public :: netcdf_put_packed_variable_array
  public :: netcdf_put_variable_vector
  public :: netcdf_coord_to_col_row
  public :: netcdf_set_coordinate_tolerance
  public :: netcdf_get_variable_id_for_variable
  public :: netcdf_rewrite_attribute
 
contains
 
!----------------------------------------------------------------------
 
function netcdf_date_within_range( NCFILE, iJulianDay)  result( lWithinRange )
 
  type (t_netcdf4_file ) :: ncfile
  integer (c_int) :: ijulianday
  logical (c_bool) :: lwithinrange
 
  if ( ijulianday >= ncfile%iFirstDayJD &
      .and. ijulianday <= ncfile%iLastDayJD ) then
 
    lwithinrange = true
 
  else
 
    lwithinrange = false
 
  endif
 
end function netcdf_date_within_range
 
!----------------------------------------------------------------------
 
!> We need two functions to convert from index to timeval, and timeval to JD;
!> note that timeval refers to the number of days from the origin
!> of the netCDF file
 
!> return the day value (number of days since origin
function nf_julian_day_to_index(NCFILE, rJulianDay)  result (iIndex)
 
  type (t_netcdf4_file) :: ncfile
  real (c_double) :: rjulianday
  integer (c_int) :: iindex
 
  iindex = int(aint(rjulianday), c_int) - ncfile%iFirstDayJD
 
end function nf_julian_day_to_index
 
!----------------------------------------------------------------------
 
function nf_index_to_dayvalue(NCFILE, iIndex)   result(rDayValue)
 
  type (t_netcdf4_file) :: ncfile
  integer (c_int) :: iindex
  real (c_double) :: rdayvalue
 
  call assert(iindex >= lbound(ncfile%rDateTimeValues, 1) &
    .and. iindex <= ubound(ncfile%rDateTimeValues, 1), &
    "Dimension out of bounds", __file__, __line__)
  rdayvalue = ncfile%rDateTimeValues(iindex)
 
end function nf_index_to_dayvalue
 
!----------------------------------------------------------------------
 
function nf_dayvalue_to_julian_day(NCFILE, rDayValue)   result(rJulianDay)
 
  type (t_netcdf4_file) :: ncfile
  real (c_double) :: rdayvalue
  real (c_double) :: rjulianday
 
  rjulianday = real(ncfile%iOriginJD, c_double) &
    + real(ncfile%iOriginHH, c_double) / 24.0_c_double &
    + real(ncfile%iOriginMM, c_double) / 1440.0_c_double &
    + real(ncfile%iOriginSS, c_double) / 86400.0_c_double &
    + rdayvalue
 
end function nf_dayvalue_to_julian_day
 
!----------------------------------------------------------------------
 
function nf_julian_day_to_index_adj( NCFILE, rJulianDay )  result(iStart)
 
  type (t_netcdf4_file ) :: ncfile
  real (c_double) :: rjulianday
  integer (c_size_t) :: istart
 
  ! [ LOCALS ]
  integer (c_int) :: imindiff, idiff
  integer (c_int) :: icandidateindex
  integer (c_int) :: iinitialcandidateindex
  integer (c_int) :: itestindex
  real (c_double) :: rtestjd
  integer (c_int) :: iindexlower, iindexupper, iindex
  logical (c_bool) :: lchanged
 
  istart = -9999
  imindiff = ibigval
  !> First guess at what the appropriate index value should be.
  !> Current JD minus the Origin JD is a good guess.
  icandidateindex = nf_julian_day_to_index(ncfile, rjulianday)
 
  call assert(icandidateindex >=0, "Problem finding the index number of the time " &
    //"variable in netCDF file "//dquote(ncfile%sFilename), __file__, __line__)
 
  iinitialcandidateindex = icandidateindex
 
  do
 
    !> calculate the range of *INDEX* values to search over
    iindexlower = max( lbound(ncfile%rDateTimeValues, 1), icandidateindex - 1)
    iindexupper = min( ubound(ncfile%rDateTimeValues, 1), icandidateindex + 1)
 
    lchanged = false
 
    do iindex=iindexlower,iindexupper
 
      rtestjd = nf_dayvalue_to_julian_day(ncfile=ncfile, &
          rdayvalue=ncfile%rDateTimeValues(iindex))
 
      itestindex = int(aint(rtestjd), c_int) - ncfile%iFirstDayJD
      idiff = abs(itestindex - iinitialcandidateindex)
 
      if (idiff < imindiff ) then
 
        imindiff = idiff
        icandidateindex = iindex
        lchanged = true
 
      endif
 
    enddo
 
    if (.not. lchanged ) exit
 
  enddo
 
  if (imindiff == 0) istart = int(icandidateindex, c_size_t)
 
end function nf_julian_day_to_index_adj
 
!----------------------------------------------------------------------
 
function nf_return_varid( NCFILE, iVarIndex)   result(iVarID)
 
  type (t_netcdf4_file ) :: ncfile
   integer (c_int) :: ivarindex
   integer (c_int) :: ivarid
 
   type (t_netcdf_variable), pointer :: pnc_var
 
   pnc_var => ncfile%pNC_VAR(ivarindex)
 
   ivarid = pnc_var%iNC_VarID
 
end function nf_return_varid
 
!----------------------------------------------------------------------
 
function nf_return_dimid( NCFILE, iDimIndex)   result(iDimID)
 
  type (t_netcdf4_file ) :: ncfile
   integer (c_int) :: idimindex
   integer (c_int) :: idimid
 
   type (t_netcdf_dimension), pointer :: pnc_dim
 
   pnc_dim => ncfile%pNC_DIM(idimindex)
 
   idimid = pnc_dim%iNC_DimID
 
end function nf_return_dimid
 
!----------------------------------------------------------------------
 
function nf_return_varindex( NCFILE, iVarID)   result(iVarIndex)
 
  type (t_netcdf4_file ) :: ncfile
   integer (c_int) :: ivarid
   integer (c_int) :: ivarindex
 
   type (t_netcdf_variable), pointer :: pnc_var
   integer (c_int) :: iindex
   logical (c_bool) :: lfound
 
  lfound = false
 
  do iindex=0, ncfile%iNumberOfVariables - 1
 
    pnc_var => ncfile%pNC_VAR(iindex)
 
    if (pnc_var%iNC_VarID == ivarid) then
      lfound = true
      exit
    endif
 
  enddo
 
  call assert(lfound, "INTERNAL PROGRAMMING ERROR - No matching variable " &
    //"ID found: was looking for Variable ID: "//trim(ascharacter(ivarid)), &
    __file__, __line__)
 
  ivarindex = iindex
 
end function nf_return_varindex
 
!----------------------------------------------------------------------
 
function nf_return_attvalue( NCFILE, iVarIndex, sAttName)   result(sAttValue)
 
  type (t_netcdf4_file ) :: ncfile
   integer (c_int) :: ivarindex
   character (len=*) :: sattname
   character (len=256) :: sattvalue
 
   type (t_netcdf_attribute), dimension(:), pointer :: pnc_att
   integer (c_int) :: iindex, iindex2
   logical (c_bool) :: lfound
 
  if (ivarindex < 0) then
 
    pnc_att => ncfile%pNC_ATT
 
  else
 
    call assert(ivarindex >= lbound(ncfile%pNC_VAR,1) &
      .and. ivarindex <= ubound(ncfile%pNC_VAR,1), &
      "Index out of bounds referencing NCFILE%pNC_VAR" &
      //"~Offending index value: "//trim(ascharacter(ivarindex)), &
      __file__, __line__)
 
    pnc_att => ncfile%pNC_VAR(ivarindex)%pNC_ATT
 
  endif
 
  lfound = false
 
  do iindex=lbound(pnc_att,1), ubound(pnc_att,1)
 
    if ( sattname .strequal. pnc_att(iindex)%sAttributeName ) then
      lfound = true
      exit
    endif
 
  enddo
 
  call assert(lfound, "INTERNAL PROGRAMMING ERROR - No matching attribute " &
    //"name found: was looking for attribute with name: "//dquote(sattname), &
    __file__, __line__)
 
  sattvalue = ""
  do iindex2=0, ubound(pnc_att(iindex)%sAttValue,1)
    sattvalue = sattvalue//" "//trim(pnc_att(iindex)%sAttValue(iindex))
  enddo
 
 
end function nf_return_attvalue
 
!----------------------------------------------------------------------
 
function nf_return_dimindex( NCFILE, iDimID)   result(iDimIndex)
 
  type (t_netcdf4_file ) :: ncfile
   integer (c_int) :: idimid
   integer (c_int) :: idimindex
 
   type (t_netcdf_dimension), pointer :: pnc_dim
   integer (c_int) :: iindex
   logical (c_bool) :: lfound
 
  lfound = false
 
  do iindex=0, ncfile%iNumberOfDimensions - 1
 
    pnc_dim => ncfile%pNC_DIM(iindex)
 
    if (pnc_dim%iNC_DimID == idimid) then
      lfound = true
      exit
    endif
 
  enddo
 
  call assert(lfound, "INTERNAL PROGRAMMING ERROR - No matching dimension " &
    //"ID found: was looking for Dimension ID: "//trim(ascharacter(idimid)), &
    __file__, __line__)
 
  idimindex = iindex
 
end function nf_return_dimindex
 
!----------------------------------------------------------------------
 
function nf_return_dimsize( NCFILE, iDimID)   result(iDimSize)
 
  type (t_netcdf4_file ) :: ncfile
   integer (c_int) :: idimid
   integer (c_size_t) :: idimsize
 
   type (t_netcdf_dimension), pointer :: pnc_dim
   integer (c_int) :: iindex
   logical (c_bool) :: lfound
 
  pnc_dim => null()
  lfound = false
 
  do iindex=0, ncfile%iNumberOfDimensions - 1
 
    pnc_dim => ncfile%pNC_DIM(iindex)
 
    if (pnc_dim%iNC_DimID == idimid) then
      lfound = true
      exit
    endif
 
  enddo
 
  call assert(lfound, "INTERNAL PROGRAMMING ERROR - No matching dimension " &
    //"ID found: was looking for Dimension ID: "//trim(ascharacter(idimid)), &
    __file__, __line__)
 
  idimsize = pnc_dim%iNC_DimSize
 
end function nf_return_dimsize
 
!--------------------------------------------------------------------------------------------------
 
subroutine nf_guess_z_variable_name(NCFILE)
 
  type (T_NETCDF4_FILE ) :: NCFILE
 
  ! [ LOCALS ]
  integer (c_int) :: iIndex
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
 
  do iindex=lbound(ncfile%pNC_VAR,1),ubound(ncfile%pNC_VAR,1)
 
    pnc_var => ncfile%pNC_VAR(iindex)
 
    if (       ( pnc_var%sVariableName .strapprox. 'time')             &
         .or.  ( pnc_var%sVariableName .strapprox. 'x')                &
         .or.  ( pnc_var%sVariableName .strapprox. 'y')                &
         .or.  ( pnc_var%sVariableName .strapprox. 'time_bnds')        &
         .or.  ( pnc_var%sVariableName .strapprox. 'lat')              &
         .or.  ( pnc_var%sVariableName .strapprox. 'lon') )  cycle
 
    ncfile%sVarName(nc_z) = pnc_var%sVariableName
    exit
 
  enddo
 
end subroutine nf_guess_z_variable_name
 
!--------------------------------------------------------------------------------------------------
 
subroutine netcdf_open_and_prepare_for_merging( NCFILE, sFilename, guess_z_var_name )
 
  type (t_netcdf4_file ), intent(inout)       :: ncfile
  character (len=*), intent(in)               :: sfilename
  logical (c_bool), intent(in), optional :: guess_z_var_name
 
  ! [ LOCALS ]
  logical (c_bool) :: guess_z_var_name_l
 
  if (present( guess_z_var_name) ) then
    guess_z_var_name_l = guess_z_var_name
  else
    guess_z_var_name_l = false
  endif
 
  call nf_open_file(ncfile=ncfile, sfilename=sfilename)
 
  call nf_populate_dimension_struct( ncfile )
 
  call nf_populate_variable_struct( ncfile )
 
  if (guess_z_var_name_l)  call nf_guess_z_variable_name(ncfile)
 
  call nf_get_variable_id_and_type( ncfile, strict_asserts=true )
 
  ! OK. We only want to attempt to call functions that
  ! process the time variable if a time variable actually exists!!
  if ( ncfile%iVarID(nc_time) >= 0 ) then
 
    ncfile%dpFirstAndLastTimeValues = nf_get_first_and_last(ncfile=ncfile, &
        ivarindex=ncfile%iVarIndex(nc_time) )
 
    !> look for and process the "days since MM-DD-YYYY" attribute
    call nf_get_time_units(ncfile=ncfile)
 
    call nf_calculate_time_range(ncfile)
 
    !> retrieve the time values as included in the netCDF file
    call nf_get_time_vals(ncfile)
 
  endif
 
  !> retrieve the X and Y coordinates from the netCDF file...
  call nf_get_x_and_y(ncfile)
 
  !> define the entire grid area as the AOI
  ncfile%iColBounds(nc_left) = lbound(ncfile%rX_Coords,1)
  ncfile%iColBounds(nc_right) = ubound(ncfile%rX_Coords,1)
 
  ncfile%iRowBounds(nc_top) = lbound(ncfile%rY_Coords,1)
  ncfile%iRowBounds(nc_bottom) = ubound(ncfile%rY_Coords,1)
 
  !> based on the subset of the netCDF file as determined above, set the
  !> start, count, and stride parameters for use in all further data
  !> retrievals
  call nf_set_start_count_stride(ncfile)
 
  !> establish the bounds to iterate over; this can enable horiz or vert flipping
  call nf_set_iteration_bounds(ncfile)
 
  !> now that we have (possibly) created a subset, need to get the
  !> **NATIVE** coordinate bounds so that the intermediate grid file
  !> can be created
  call nf_return_native_coord_bounds(ncfile)
 
end subroutine netcdf_open_and_prepare_for_merging
 
!----------------------------------------------------------------------
 
subroutine netcdf_open_and_prepare_as_input(NCFILE, sFilename, &
    lFlipHorizontal, lFlipVertical,                            &
    lAllowAutomaticDataFlipping,                               &
    rX_Coord_AddOffset, rY_Coord_AddOffset,                    &
    sVariableOrder, sVarName_x,                                &
    sVarName_y, sVarName_z, sVarName_time,                     &
    rCoordinateTolerance,                                      &
    tGridBounds, iLU)
 
  type (t_netcdf4_file ) :: ncfile
  character (len=*) :: sfilename
  logical (c_bool), optional :: lfliphorizontal
  logical (c_bool), optional :: lflipvertical
  logical (c_bool), optional :: lallowautomaticdataflipping
  character (len=*), optional :: svariableorder
  real (c_double), optional    :: rx_coord_addoffset
  real (c_double), optional    :: ry_coord_addoffset
  character (len=*), optional :: svarname_x
  character (len=*), optional :: svarname_y
  character (len=*), optional :: svarname_z
  character (len=*), optional :: svarname_time
  real (c_double), optional   :: rcoordinatetolerance
  type (grid_bounds_t), optional :: tgridbounds
  integer (c_int), optional :: ilu
 
  ! [ LOCALS ]
  logical :: lfileopen
  integer (c_size_t), dimension(2) :: icolrow_ll, icolrow_ur, icolrow_lr, icolrow_ul
 
  call nf_open_file(ncfile=ncfile, sfilename=sfilename)
 
  call nf_populate_dimension_struct( ncfile )
  call nf_populate_variable_struct( ncfile )
 
  if (present(lfliphorizontal) ) ncfile%lFlipHorizontal = lfliphorizontal
  if (present(lflipvertical) ) ncfile%lFlipVertical = lflipvertical
  if (present( lallowautomaticdataflipping) )                                 &
      ncfile%lAllowAutomaticDataFlipping = lallowautomaticdataflipping
  if (present(rx_coord_addoffset))  ncfile%rX_Coord_AddOffset = rx_coord_addoffset
  if (present(ry_coord_addoffset))  ncfile%rY_Coord_AddOffset = ry_coord_addoffset
 
  if (present(rcoordinatetolerance))  ncfile%rCoordinateTolerance = rcoordinatetolerance
 
  if (present(svariableorder) )  ncfile%sVariableOrder = svariableorder
 
  if( present(ilu) ) then
    inquire (unit=ilu, opened=lfileopen)
    if ( lfileopen )  call netcdf_dump_cdl( ncfile, ilu)
  endif
 
  if (present(svarname_x) ) then
    ncfile%sVarName(nc_x) = svarname_x
  else
    ncfile%sVarName(nc_x) = "x"
  endif
 
  if (present(svarname_y) ) then
    ncfile%sVarName(nc_y) = svarname_y
  else
    ncfile%sVarName(nc_y) = "y"
  endif
 
  if (present(svarname_z) ) then
    ncfile%sVarName(nc_z) = svarname_z
  else
    ncfile%sVarName(nc_z) = "prcp"
  endif
 
  if (present(svarname_time) ) then
    ncfile%sVarName(nc_time) = svarname_time
  else
    ncfile%sVarName(nc_time) = "time"
  endif
 
  call nf_get_variable_id_and_type( ncfile )
 
  ! OK. We only want to attempt to call functions that
  ! process the time variable if a time variable actually exists!!
  if ( ncfile%iVarID(nc_time) >= 0 ) then
 
    ncfile%dpFirstAndLastTimeValues = nf_get_first_and_last(ncfile=ncfile, &
        ivarindex=ncfile%iVarIndex(nc_time) )
   
    !> look for and process the "days since MM-D-YYYY" attribute
    call nf_get_time_units(ncfile=ncfile)
 
    call nf_calculate_time_range(ncfile)
 
    !> retrieve the time values as included in the netCDF file
    call nf_get_time_vals(ncfile)
 
  endif
 
  call nf_get_xyz_units(ncfile=ncfile)
 
  !> establish scale_factor and add_offset values, if present
  call nf_get_scale_and_offset(ncfile=ncfile)
 
  !> retrieve the X and Y coordinates from the netCDF file...
  call nf_get_x_and_y(ncfile)
 
  if (present(tgridbounds) ) then
 
    !> define a subset of the grid as the AOI
    !> need all four corner points since it is likely that
    !> the AOI rectangle is rotated relative to the base
    !> projection
 
    icolrow_ll = netcdf_coord_to_col_row(ncfile=ncfile, &
                                     rx=tgridbounds%rXll, &
                                     ry=tgridbounds%rYll)
  
    icolrow_lr = netcdf_coord_to_col_row(ncfile=ncfile, &
                                     rx=tgridbounds%rXlr, &
                                     ry=tgridbounds%rYlr)
 
    icolrow_ul = netcdf_coord_to_col_row(ncfile=ncfile, &
                                     rx=tgridbounds%rXul, &
                                     ry=tgridbounds%rYul)
 
    icolrow_ur = netcdf_coord_to_col_row(ncfile=ncfile, &
                                     rx=tgridbounds%rXur, &
                                     ry=tgridbounds%rYur)
 
#ifdef DEBUG_PRINT
    write(*, fmt="(a)") "subroutine 'netcdf_open_and_prepare_as_input'"
    write(*, fmt="(a,a,i6)") "Find correspondence between project bounds (in native projection) and row, col of dataset | ", &
      trim(__file__), __line__
    write(*, fmt="(a)") "      column     row              X              Y"
    write(*, fmt="(a,i6,i6,a,f14.3,f14.3)") "LL: ", icolrow_ll(column), icolrow_ll(row), " <==> ", tgridbounds%rXll, tgridbounds%rYll
    write(*, fmt="(a,i6,i6,a,f14.3,f14.3)") "LR: ", icolrow_lr(column), icolrow_lr(row), " <==> ", tgridbounds%rXlr, tgridbounds%rYlr
    write(*, fmt="(a,i6,i6,a,f14.3,f14.3)") "UL: ", icolrow_ul(column), icolrow_ul(row), " <==> ", tgridbounds%rXul, tgridbounds%rYul
    write(*, fmt="(a,i6,i6,a,f14.3,f14.3)") "UR: ", icolrow_ur(column), icolrow_ur(row), " <==> ", tgridbounds%rXur, tgridbounds%rYur
#endif
 
    ncfile%iColBounds(nc_left) = &
      max( min( icolrow_ul(column), icolrow_ur(column), icolrow_ll(column), icolrow_lr(column) ) - 4_c_size_t, &
                int(lbound(ncfile%rX_Coords,1), c_size_t) )
 
    ncfile%iColBounds(nc_right) = &
      min( max( icolrow_ul(column), icolrow_ur(column), icolrow_ll(column), icolrow_lr(column) ) + 4_c_size_t, &
                int(ubound(ncfile%rX_Coords,1), c_size_t) )
 
 
      ncfile%iRowBounds(nc_top) = &
        max( min( icolrow_ul(row), icolrow_ur(row), icolrow_ll(row), icolrow_lr(row) ) - 4_c_size_t, &
                  int(lbound(ncfile%rY_Coords,1), c_size_t) )
 
      ncfile%iRowBounds(nc_bottom) = &
        min( max( icolrow_ul(row), icolrow_ur(row), icolrow_ll(row), icolrow_lr(row) ) + 4_c_size_t, &
                  int(ubound(ncfile%rY_Coords,1), c_size_t) )
 
  else
 
    !> define the entire grid area as the AOI
    ncfile%iColBounds(nc_left) = lbound(ncfile%rX_Coords,1)
    ncfile%iColBounds(nc_right) = ubound(ncfile%rX_Coords,1)
 
    ncfile%iRowBounds(nc_top) = lbound(ncfile%rY_Coords,1)
    ncfile%iRowBounds(nc_bottom) = ubound(ncfile%rY_Coords,1)
 
  endif
 
  !> based on the subset of the netCDF file as determined above, set the
  !> start, count, and stride parameters for use in all further data
  !> retrievals
  call nf_set_start_count_stride(ncfile)
 
  !> establish the bounds to iterate over; this can enable horiz or vert flipping
  call nf_set_iteration_bounds(ncfile)
 
  !> now that we have (possibly) created a subset, need to get the
  !> **NATIVE** coordinate bounds so that the intermediate grid file
  !> can be created
  call nf_return_native_coord_bounds(ncfile)
 
end subroutine netcdf_open_and_prepare_as_input
 
!----------------------------------------------------------------------
 
subroutine netcdf_open_and_prepare_as_output_archive(NCFILE, NCFILE_ARCHIVE, &
   iOriginMonth, iOriginDay, iOriginYear, iStartYear, iEndYear)
 
  type (t_netcdf4_file ) :: ncfile
  type (t_netcdf4_file ) :: ncfile_archive
  integer (c_int) :: ioriginmonth
  integer (c_int) :: ioriginday
  integer (c_int) :: ioriginyear
  integer (c_int) :: istartyear
  integer (c_int) :: iendyear
 
  ! [ LOCALS ]
  integer (c_int) :: inumcols, inumrows
  integer (c_int) :: imincol, imaxcol
  integer (c_int) :: iminrow, imaxrow
  real (c_double), dimension(:), allocatable :: rx, ry
  character (len=10) :: sorigintext
  character (len=256) :: sfilename
 
  write(sorigintext, fmt="(i4.4,'-',i2.2,'-',i2.2)") ioriginyear, &
    ioriginmonth, ioriginday
 
  imaxrow = int(maxval(ncfile%iRowBounds), c_int)
  iminrow = int(minval(ncfile%iRowBounds), c_int)
  imaxcol = int(maxval(ncfile%iColBounds), c_int)
  imincol = int(minval(ncfile%iColBounds), c_int)
 
  inumrows = imaxrow - iminrow + 1
  inumcols = imaxcol - imincol + 1
 
  allocate(rx(inumcols))
  allocate(ry(inumrows))
  rx = ncfile%rX_Coords(imincol:imaxcol)
  ry = ncfile%rY_Coords(iminrow:imaxrow)
 
  sfilename = trim(ncfile%sVarName(nc_z))//"_"//trim(ascharacter(istartyear)) &
    //"_"//trim(ascharacter(iendyear))//"__" &
    //trim(ascharacter(inumrows)) &
     //"_by_"//trim(ascharacter(inumcols))//".nc"
 
  call nf_create(ncfile=ncfile_archive, sfilename=trim(sfilename) )
 
  !> set dimension values in the NCFILE struct
  call nf_set_standard_dimensions(ncfile=ncfile_archive, &
                       inx=inumcols, &
                       iny=inumrows)
 
  ncfile_archive%sVarUnits(nc_x) =   ncfile%sVarUnits(nc_x)
  ncfile_archive%sVarUnits(nc_y) =   ncfile%sVarUnits(nc_y)
  ncfile_archive%sVarUnits(nc_z) =   ncfile%sVarUnits(nc_z)
 
  !> transfer dimension values to netCDF file
  call nf_define_dimensions( ncfile=ncfile_archive )
 
  !> set variable values in the NCFILE struct
  call nf_set_standard_variables(ncfile=ncfile_archive, &
       svarname_z = trim(ncfile%sVarName(nc_z)) )
 
  !> transfer variable values to netCDF file
  call nf_define_variables(ncfile=ncfile_archive)
 
  call nf_get_variable_id_and_type( ncfile=ncfile_archive )
 
  call nf_set_standard_attributes(ncfile=ncfile_archive, &
    sorigintext=sorigintext)
 
  call nf_set_global_attributes(ncfile=ncfile_archive, &
     sdatatype=trim(ncfile%sVarName(nc_z)), &
     ssourcefile=trim(ncfile%sFilename))
 
  call nf_put_attributes(ncfile=ncfile_archive)
 
  !> enable a low level of data compression for the
  !> variable of interest
  call nf_define_deflate(ncfile=ncfile_archive, &
     ivarid=ncfile_archive%iVarID(nc_z), &
     ishuffle=nc_shuffle_yes, &
     ideflate=nc_deflate_yes, &
     ideflate_level=2 )
 
  call nf_enddef(ncfile=ncfile_archive)
 
  call nf_put_x_and_y(ncfile=ncfile_archive, &
       dpx=ncfile%rX_Coords(imincol:imaxcol), &
       dpy=ncfile%rY_Coords(iminrow:imaxrow) )
!       dpX=rX, &
!       dpY=rY )
 
!  call netcdf_close_file(NCFILE_ARCHIVE)
 
end subroutine netcdf_open_and_prepare_as_output_archive
 
 
 
subroutine netcdf_open_and_prepare_as_output( NCFILE, sVariableName, sVariableUnits,    &
   iNX, iNY, fX, fY, StartDate, EndDate, PROJ4_string, history_list, executable_name,   &
   dpLat, dpLon, fValidMin, fValidMax, write_time_bounds, filename_prefix,              &
   filename_modifier)
 
  type (t_netcdf4_file ), pointer, intent(inout)       :: ncfile
  character (len=*), intent(in)                        :: svariablename
  character (len=*), intent(in)                        :: svariableunits
  integer (c_int), intent(in)                          :: iny
  real (c_double), intent(in)                          :: fx(:)
  integer (c_int), intent(in)                          :: inx
  real (c_double), intent(in)                          :: fy(:)
  type (datetime_t), intent(in)                        :: startdate
  type (datetime_t), intent(in)                        :: enddate
  character (len=*), intent(in)                        :: proj4_string
  type (fstring_list_t), intent(in), pointer, optional :: history_list
  character (len=*), intent(in), optional              :: executable_name
  real (c_double), intent(in), optional                :: dplat(:,:)
  real (c_double), intent(in), optional                :: dplon(:,:)
  real (c_float), intent(in), optional                 :: fvalidmin
  real (c_float), intent(in), optional                 :: fvalidmax
  logical (c_bool), intent(in), optional               :: write_time_bounds
  character (len=*), intent(in), optional              :: filename_prefix
  character (len=*), intent(in), optional              :: filename_modifier
 
!   ! [ LOCALS ]
  character (len=10)                                :: sorigintext
  character (len=:), allocatable                    :: sfilename
  type (fstring_list_t), pointer                    :: history_list_l
  logical (c_bool)                                  :: write_time_bounds_l
  character (len=:), allocatable                    :: executable_name_l
  real (c_float)                                    :: valid_minimum
  real (c_float)                                    :: valid_maximum
  logical (c_bool)                                  :: include_latlon
  type (datetime_t)                                 :: dt
  character (len=:), allocatable                    :: date_time_text
  character (len=:), allocatable                    :: filename_prefix_l
  character (len=:), allocatable                    :: filename_modifier_l
  call dt%systime()
  date_time_text = dt%prettydatetime()
 
 
  if (present( fvalidmin ) ) then
    valid_minimum = fvalidmin
  else
    valid_minimum = -1.0e+3
  endif
 
  if (present( fvalidmax ) ) then
    valid_maximum = fvalidmax
  else
    valid_maximum = 1.0e+10
  endif
 
  if ( present( executable_name ) ) then
    executable_name_l = trim( executable_name )
  else
    executable_name_l = "SWB"
  endif
 
  if ( present( filename_modifier ) ) then
    filename_modifier_l = "_"//trim( filename_modifier )//"_"
  else
    filename_modifier_l = "_"
  endif
 
  if ( present( filename_prefix ) ) then
    if (len_trim(filename_prefix) > 0) then
      filename_prefix_l = trim( filename_prefix )//"_"
    else 
      filename_prefix_l = ""
    endif
  else
    filename_prefix_l = ""
  endif
 
  if ( present( history_list) ) then
    history_list_l => history_list
  else
    allocate( history_list_l )
    call history_list_l%append(date_time_text//": Soil-Water-Balance run started.")
  endif
 
  if ( present( write_time_bounds ) ) then
    write_time_bounds_l = write_time_bounds
  else
    write_time_bounds_l = false
  endif
 
  include_latlon = logical( present( dpLat ) .and. present( dpLon ), c_bool )
 
  write(sorigintext, fmt="(i4.4,'-',i2.2,'-',i2.2)") startdate%iYear, startdate%iMonth, startdate%iDay
 
  ! if a filename_prefix argument has been supplied, override the variable specified in
  ! the 'output' module
  if (len_trim(filename_prefix_l) > 0 ) then
 
  sfilename = trim(output_directory_name)//trim(filename_prefix_l)              &
    //trim(svariablename)//"_"//filename_modifier_l                             &
    //startdate%prettydate()//"_to_"//enddate%prettydate()//"__"                &
    //trim(ascharacter(iny))//"_by_"//trim(ascharacter(inx))//".nc"
 
  else
 
    sfilename = trim(output_directory_name)//trim(output_prefix_name)           &
    //trim(svariablename)//"_"//filename_modifier_l                             &
    //startdate%prettydate()//"_to_"//enddate%prettydate()//"__"                &
    //trim(ascharacter(iny))//"_by_"//trim(ascharacter(inx))//".nc"
 
  endif
 
  call logs%write("Attempting to open netCDF file for writing with filename "//dquote(sfilename))
 
  call nf_create(ncfile=ncfile, sfilename=trim(sfilename) )
 
  !> set dimension values in the NCFILE struct
  call nf_set_standard_dimensions(ncfile=ncfile, inx=inx, iny=iny,            &
                                  write_time_bounds=write_time_bounds_l )
 
  !> @todo implement more flexible method of assigning units
!   NCFILE%sVarUnits(NC_X)    = sXY_units
!   NCFILE%sVarUnits(NC_Y)    = sXY_units
  ncfile%sVarUnits(nc_z)    = svariableunits
 
  !> transfer dimension values to netCDF file
  call nf_define_dimensions( ncfile=ncfile )
 
  !> set variable values in the NCFILE struct
  call nf_set_standard_variables(ncfile=ncfile, svarname_z = svariablename,   &
    llatlon=include_latlon, write_time_bounds=write_time_bounds_l )
 
  !> transfer variable values to netCDF file
  call nf_define_variables(ncfile=ncfile)
 
  call nf_get_variable_id_and_type( ncfile=ncfile )
 
  call nf_set_standard_attributes(ncfile=ncfile, sorigintext=sorigintext,     &
      proj4_string=proj4_string, llatlon=include_latlon, fvalidmin=valid_minimum, &
      fvalidmax=valid_maximum, write_time_bounds=write_time_bounds_l )
 
  call nf_set_global_attributes(ncfile=ncfile, &
     sdatatype=trim(ncfile%sVarName(nc_z)), history_list=history_list_l, &
     executable_name=executable_name_l )
 
  call nf_put_attributes(ncfile=ncfile)
 
  !> enable a low level of data compression for the variable of interest
   call nf_define_deflate(ncfile=ncfile, &
      ivarid=ncfile%iVarID(nc_z), &
      ishuffle=nc_shuffle_yes, &
      ideflate=nc_deflate_yes, &
      ideflate_level=2 )
 
  call nf_enddef(ncfile=ncfile)
 
  ! we are only supplying a vector of X and Y on the assumption that the base projection
  ! results in a uniform grid (in other words, we have the same X value for all coluns of a given row)
  call nf_put_x_and_y(ncfile=ncfile, &
       dpx=fx, &
       dpy=fy )
 
!  allocate( NCFILE%rX_Coords( size( fX) ), stat=iStat )
!  allocate( NCFILE%rY_Coords( size( fY) ), stat=iStat )
  ncfile%rX_Coords = fx
  ncfile%rY_Coords = fy
 
  if (present( dplat ) .and. present( dplon ) ) then
 
    call nf_put_lat_and_lon(ncfile=ncfile, &
       dplat=dplat, &
       dplon=dplon )
 
  endif
 
end subroutine netcdf_open_and_prepare_as_output
 
 
!----------------------------------------------------------------------
 
subroutine nf_set_z_variable_name(NCFILE, sVarName_z)
 
  type (T_NETCDF4_FILE ) :: NCFILE
  character (len=*) :: sVarName_z
 
  ncfile%sVarName(nc_z) = svarname_z
 
end subroutine nf_set_z_variable_name
 
!----------------------------------------------------------------------
 
subroutine nf_set_iteration_bounds(NCFILE)
 
  type (T_NETCDF4_FILE ) :: NCFILE
 
!  if (NCFILE%lFlipVertical) then
!    NCFILE%iRowIter(NC_FIRST) = NCFILE%iNY
!    NCFILE%iRowIter(NC_LAST) = 1
!    NCFILE%iRowIter(NC_BY) = -1
! else
    ncfile%iRowIter(nc_first) = 1
    ncfile%iRowIter(nc_last) = ncfile%iNY
    ncfile%iRowIter(nc_by) = 1
!  endif
 
  ! if (NCFILE%lFlipHorizontal) then
  !   NCFILE%iColIter(NC_FIRST) = NCFILE%iNX
  !   NCFILE%iColIter(NC_LAST) = 1
  !   NCFILE%iColIter(NC_BY) = -1
  ! else
    ncfile%iColIter(nc_first) = 1
    ncfile%iColIter(nc_last) = ncfile%iNX
    ncfile%iColIter(nc_by) = 1
  ! endif
 
end subroutine nf_set_iteration_bounds
 
!----------------------------------------------------------------------
 
subroutine nf_set_start_count_stride(NCFILE)
 
  type (T_NETCDF4_FILE ) :: NCFILE
 
  ! [ LOCALS ]
  integer (c_int) :: iIndex
 
  ! loop over the three (assumed) dimensions of the "Z" variable;
  ! assign appropriate bounds to each
  do iindex = 0,3
 
    select case (iindex)
 
      case (nc_x)
 
        !> need to subtract 1 from the start index: we are using the
        !> netCDF C API, in which index values are relative to zero
        ncfile%iStart(iindex) = minval(ncfile%iColBounds) - 1_c_size_t
        ncfile%iNX = int(maxval(ncfile%iColBounds) - minval(ncfile%iColBounds) + 1_c_size_t, c_int)
        ncfile%iCount(iindex) = int(ncfile%iNX, c_size_t)
!        NCFILE%iCount(iIndex) = maxval(NCFILE%iColBounds) - minval(NCFILE%iColBounds)
        ncfile%iStride(iindex) = 1_c_size_t
 
      case (nc_y)
 
        !> note: this assumes that the row numbers increase from top to bottom,
        !>       while the Y coordinates decrease top to bottom
 
        ncfile%iStart(iindex) = minval(ncfile%iRowBounds) - 1_c_size_t
        ncfile%iNY = int(maxval(ncfile%iRowBounds) - minval(ncfile%iRowBounds) + 1_c_size_t, c_int)
        ncfile%iCount(iindex) = int(ncfile%iNY, c_size_t)
        !>
        !> count must be set to the number of values! maxval minus minval results
        !> in a diagonal pattern in the input as we read in the incorrect number
        !> of results
!        NCFILE%iCount(iIndex) = maxval(NCFILE%iRowBounds) - minval(NCFILE%iRowBounds)
        ncfile%iStride(iindex) = 1_c_size_t
 
      case (nc_time)
 
        ncfile%iStart(iindex) = 0_c_size_t
        ncfile%iCount(iindex) = 1_c_size_t
        ncfile%iStride(iindex) = 1_c_size_t
 
      case default
 
    end select
 
  enddo
 
end subroutine nf_set_start_count_stride
 
!----------------------------------------------------------------------
 
subroutine nf_return_native_coord_bounds(NCFILE)
 
  type (T_NETCDF4_FILE ) :: NCFILE
 
  ! [ LOCALS ]
  real (c_double) :: rXmin, rXmax
  real (c_double) :: rYmin, rYmax
 
  !> find the (x,y) associated with the column and row number bounds
  rxmin = minval(ncfile%rX_Coords(ncfile%iColBounds(nc_left):ncfile%iColBounds(nc_right)) )
  rxmax = maxval(ncfile%rX_Coords(ncfile%iColBounds(nc_left):ncfile%iColBounds(nc_right)) )
  rymin = minval(ncfile%rY_Coords(ncfile%iRowBounds(nc_top):ncfile%iRowBounds(nc_bottom)) )
  rymax = maxval(ncfile%rY_Coords(ncfile%iRowBounds(nc_top):ncfile%iRowBounds(nc_bottom)) )
 
  ncfile%rX(nc_left) = rxmin - ncfile%rGridCellSizeX * 0.5_c_double
  ncfile%rX(nc_right) = rxmax + ncfile%rGridCellSizeX * 0.5_c_double
  ncfile%rY(nc_top) = rymax + ncfile%rGridCellSizeY * 0.5_c_double
  ncfile%rY(nc_bottom) = rymin - ncfile%rGridCellSizeY * 0.5_c_double
 
#ifdef DEBUG_PRINT
  print *, ""
  print *, repeat("-", 80)
  print *, "Filename: ", ncfile%sFilename
  print *, "Grid cell size (X): ", ncfile%rGridCellSizeX
  print *, "Grid cell size (Y): ", ncfile%rGridCellSizeY
 
  print *, "Bounds of data subset area, in native coordinates"
  print *, "X (left): ", ncfile%rX(nc_left)
  print *, "X (right): ", ncfile%rX(nc_right)
  print *, "Y (top): ", ncfile%rY(nc_top)
  print *, "Y (bottom): ", ncfile%rY(nc_bottom)
  print *, ""
#endif
 
end subroutine nf_return_native_coord_bounds
 
!----------------------------------------------------------------------
 
subroutine nf_get_time_vals(NCFILE)
 
  type (T_NETCDF4_FILE ) :: NCFILE
 
  integer (c_int) :: iVarIndex_time
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR_time
  type (T_NETCDF_DIMENSION), pointer :: pNC_DIM_time
  integer (c_int) :: iStat
 
  istat = 0
 
  ivarindex_time = ncfile%iVarIndex(nc_time)
 
  call assert(ivarindex_time >= lbound(ncfile%pNC_VAR,1) &
    .and. ivarindex_time <= ubound(ncfile%pNC_VAR,1), &
    "INTERNAL PROGRAMMING ERROR - Index out of bounds", __file__, __line__)
 
  pnc_var_time => ncfile%pNC_VAR(ivarindex_time)
  pnc_dim_time => ncfile%pNC_DIM( pnc_var_time%iNC_DimID(0) )
 
  if (allocated(ncfile%rDateTimeValues) ) deallocate(ncfile%rDateTimeValues, stat=istat)
  call assert(istat==0, "Failed to deallocate memory for time values", &
    __file__, __line__)
 
  allocate( ncfile%rDateTimeValues(0 : pnc_dim_time%iNC_DimSize-1_c_size_t ), stat=istat )
  call assert(istat==0, "Failed to allocate memory for time values", &
    __file__, __line__)
 
  !> @todo allow time to be read in as float, short, or int as well
 
  call nf_get_variable_vector_double(ncfile=ncfile, &
       inc_varid=pnc_var_time%iNC_VarID, &
       inc_start=0_c_size_t, &
       inc_count=pnc_dim_time%iNC_DimSize, &
       inc_stride=1_c_size_t, &
       dpnc_vars=ncfile%rDateTimeValues)
 
end subroutine nf_get_time_vals
 
!----------------------------------------------------------------------
 
subroutine nf_get_x_and_y(NCFILE)
 
  type (T_NETCDF4_FILE ) :: NCFILE
 
  integer (c_int) :: iVarIndex_x, iVarIndex_y
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR_x, pNC_VAR_y
  type (T_NETCDF_DIMENSION), pointer :: pNC_DIM_x, pNC_DIM_y
  integer (c_int) :: iLowerBound, iUpperBound
  integer (c_int) :: iStat
 
  ivarindex_x = ncfile%iVarIndex(nc_x)
  ivarindex_y = ncfile%iVarIndex(nc_y)
 
  call assert(ivarindex_x >= lbound(ncfile%pNC_VAR,1) &
    .and. ivarindex_x <= ubound(ncfile%pNC_VAR,1), &
    "INTERNAL PROGRAMMING ERROR - Index out of bounds", __file__, __line__)
 
  call assert(ivarindex_y >= lbound(ncfile%pNC_VAR,1) &
    .and. ivarindex_y <= ubound(ncfile%pNC_VAR,1), &
    "INTERNAL PROGRAMMING ERROR - Index out of bounds", __file__, __line__)
 
  pnc_var_x => ncfile%pNC_VAR(ivarindex_x)
  pnc_var_y => ncfile%pNC_VAR(ivarindex_y)
 
  call assert( pnc_var_x%iNumberOfDimensions == 1, &
    "Dimensions other than one for the x-coordinate variable are currently unsupported.", &
    __file__, __line__)
 
  call assert( pnc_var_y%iNumberOfDimensions == 1, &
    "Dimensions other than one for the y-coordinate variable are currently unsupported.", &
    __file__, __line__)
 
  pnc_dim_x => ncfile%pNC_DIM( pnc_var_x%iNC_DimID(0) )
  pnc_dim_y => ncfile%pNC_DIM( pnc_var_y%iNC_DimID(0) )
 
  allocate( ncfile%rX_Coords( pnc_dim_x%iNC_DimSize ), stat=istat )
  call assert(istat==0, "Failed to allocate memory for X-coordinate values", &
    __file__, __line__)
 
  allocate (ncfile%rY_Coords( pnc_dim_y%iNC_DimSize  ), stat=istat )
  call assert(istat==0, "Failed to allocate memory for Y-coordinate values", &
    __file__, __line__)
 
  call nf_get_variable_vector_double(ncfile=ncfile, &
       inc_varid=pnc_var_x%iNC_VarID, &
       inc_start=0_c_size_t, &
       inc_count=pnc_dim_x%iNC_DimSize, &
       inc_stride=1_c_size_t, &
       dpnc_vars=ncfile%rX_Coords)
 
  call nf_get_variable_vector_double(ncfile=ncfile, &
       inc_varid=pnc_var_y%iNC_VarID, &
       inc_start=0_c_size_t, &
       inc_count=pnc_dim_y%iNC_DimSize, &
       inc_stride=1_c_size_t, &
       dpnc_vars=ncfile%rY_Coords)
 
  ncfile%rX_Coords = ncfile%rX_Coords + ncfile%rX_Coord_AddOffset
  ncfile%rY_Coords = ncfile%rY_Coords + ncfile%rY_Coord_AddOffset
 
  ilowerbound = lbound(ncfile%rX_Coords, 1)
  iupperbound = ubound(ncfile%rX_Coords, 1)
 
  if (ncfile%rX_Coords(iupperbound) > ncfile%rX_Coords(ilowerbound) ) then
    ncfile%lX_IncreasesWithIndex = true
    if ( ncfile%lAllowAutomaticDataFlipping ) ncfile%lFlipHorizontal = false
  else
    ncfile%lX_IncreasesWithIndex = false
    if ( ncfile%lAllowAutomaticDataFlipping ) then
      ncfile%lFlipHorizontal       = true
      call logs%write( "** Horizontal coordinates decrease with index values **", &
                       ilinesbefore=1,                                            &
                       iloglevel=log_all )
      call logs%write( "  ==> flipping grid horizontally",                        &
                       ilinesafter=1,                                             &
                       iloglevel=log_all )
    endif
  endif
 
  ilowerbound = lbound(ncfile%rY_Coords, 1)
  iupperbound = ubound(ncfile%rY_Coords, 1)
 
  if (ncfile%rY_Coords(iupperbound) > ncfile%rY_Coords(ilowerbound) ) then
    ncfile%lY_IncreasesWithIndex = true
    if ( ncfile%lAllowAutomaticDataFlipping ) then
      ncfile%lFlipVertical         = true
      call logs%write( "** Vertical coordinates increase with index values **", &
                       ilinesbefore=1,                                          &
                       iloglevel=log_all )
      call logs%write( "  ==> flipping grid vertically",                        &
                       ilinesafter=1,                                           &
                       iloglevel=log_all )
    endif
  else
    ncfile%lY_IncreasesWithIndex = false
    if ( ncfile%lAllowAutomaticDataFlipping ) ncfile%lFlipVertical = false
  endif
 
  call assert(pnc_dim_x%iNC_DimSize > 2, "INTERNAL PROGRAMMING ERROR - " &
    //"netCDF X dimension size must be greater than 2.", __file__, __line__)
 
  call assert(pnc_dim_y%iNC_DimSize > 2, "INTERNAL PROGRAMMING ERROR - " &
    //"netCDF Y dimension size must be greater than 2.", __file__, __line__)
 
  ncfile%rGridCellSizeX = ( maxval(ncfile%rX_Coords) &
                                - minval(ncfile%rX_Coords) ) &
                                / real(pnc_dim_x%iNC_DimSize - 1, c_double)
 
  ncfile%rGridCellSizeY = ( maxval(ncfile%rY_Coords) &
                                - minval(ncfile%rY_Coords) ) &
                                / real(pnc_dim_y%iNC_DimSize - 1, c_double)
 
  ! print *, '*** netCDF x- and y- coord read ***'
  ! print *, trim(__FILE__), ': ', __LINE__
  ! print *, 'file:', trim(NCFILE%sFilename)
  ! print *, 'x-coords: ', NCFILE%rX_Coords
  ! print *, 'y-coords: ', NCFILE%rY_Coords
  ! print *, '***********************************'
 
end subroutine nf_get_x_and_y
 
!----------------------------------------------------------------------
 
subroutine nf_open_file(NCFILE, sFilename, iLU)
 
  type (T_NETCDF4_FILE ) :: NCFILE
  character (len=*) :: sFilename
  integer (c_int), optional :: iLU
 
  ! [ LOCALS ]
  logical :: lFileOpen
 
  call logs%write("Attempting to open READONLY netCDF file: " &
    //dquote(sfilename))
 
  call nf_trap( nc_open(trim(sfilename)//c_null_char, &
                nc_readonly, ncfile%iNCID), __file__, __line__ )
 
  call nf_trap( nc_inq_format(ncid=ncfile%iNCID, formatp=ncfile%iFileFormat), &
               __file__, __line__)
 
  call logs%write("   Succeeded.  ncid: "//trim(ascharacter(ncfile%iNCID)) &
         //"  format: "//trim(netcdf_format_string(ncfile%iFileFormat) ) )
 
  ncfile%sFilename = sfilename
 
  if( present(ilu) ) then
    inquire (unit=ilu, opened=lfileopen)
    if ( lfileopen )  call netcdf_dump_cdl( ncfile, ilu)
  endif
 
end subroutine nf_open_file
 
!----------------------------------------------------------------------
 
subroutine netcdf_open_file(NCFILE, sFilename, iLU)
 
  type (t_netcdf4_file ) :: ncfile
  character (len=*) :: sfilename
  integer (c_int), optional :: ilu
 
  if (present(ilu) ) then
 
    call nf_open_file(ncfile=ncfile, &
                    sfilename=sfilename, &
                    ilu=ilu)
 
  else
 
    call nf_open_file(ncfile=ncfile, &
                    sfilename=sfilename)
 
  endif
 
  !> Similarly, the structure of the file may be slightly different from the
  !> previous file
  call nf_populate_dimension_struct( ncfile )
  call nf_populate_variable_struct( ncfile )
 
  !> CANNOT ASSUME THAT THIS WILL REMAIN CONSTANT ACROSS FILES FROM THE
  !> SAME PROVIDER!! MUST UPDATE TO ENSURE THAT THE INDICES ARE STILL RELEVANT
  call nf_get_variable_id_and_type( ncfile )
 
  ! OK. We only want to attempt to call functions that
  ! process the time variable if a time variable actually exists!!
  !
  ! NOTE: this was previously coded as 'if ( NCFILE%iVarID(NC_TIME) > 0 ) then', which meant that the code
  !       failed to properly initialize the time units in the event that the TIME variable happened to be
  !       variable number 0
  !
  if ( ncfile%iVarID(nc_time) >= 0 ) then
 
    ncfile%dpFirstAndLastTimeValues = nf_get_first_and_last(ncfile=ncfile, &
        ivarindex=ncfile%iVarIndex(nc_time) )
 
    !> retrieve the origin for the time units associated with this file
    call nf_get_time_units(ncfile=ncfile)
 
    !> retrieve the time value specific to this file
    call nf_get_time_vals(ncfile)
 
    call nf_calculate_time_range(ncfile)
 
  endif
 
  !> establish scale_factor and add_offset values, if present
  call nf_get_scale_and_offset(ncfile=ncfile)
 
end subroutine netcdf_open_file
 
!----------------------------------------------------------------------
 
subroutine nf_trap( iResultCode, sFilename, iLineNumber, netcdf_filename )
 
  integer (c_int) :: iResultCode
  character (len=*), optional :: sFilename
  integer (c_int), optional :: iLineNumber
  character (len=*), optional :: netcdf_filename
 
  ! [ LOCALS ]
  type(c_ptr) :: cpResult
  character (len=256) :: sTextString
  character (len=256) :: sFile
  integer (c_int) :: iLine
 
  if (iresultcode /= 0) then
 
    if (present(sfilename)) then
      sfile = trim(sfilename)
    else
      sfile = __file__
    endif
 
    if (present(ilinenumber)) then
      iline = ilinenumber
    else
      iline = __line__
    endif
 
    cpresult = nc_strerror(iresultcode)
    stextstring = char_ptr_to_fortran_string( cpresult )
 
    if (present( netcdf_filename ) )                                           &
      call logs%write("netCDF filename: "//dquote( netcdf_filename  ) )
 
    call logs%write("netCDF ERROR: "//dquote( stextstring  )//" | error code was: " &
      //trim(ascharacter(iresultcode)) )
 
    call assert(false, "SWB is stopping due to a problem reading or writing" &
      //" a netCDF file", trim(sfile), iline)
 
  endif
 
end subroutine nf_trap
 
!----------------------------------------------------------------------
 
subroutine netcdf_close_file( NCFILE)
 
  type (t_netcdf4_file ) :: ncfile
 
  call logs%write("Closing netCDF file with name: "//dquote(ncfile%sFilename))
  call nf_trap( nc_close(ncfile%iNCID), __file__, __line__ )
 
!  call nf_deallocate_data_struct( NCFILE=NCFILE )
 
end subroutine netcdf_close_file
 
!----------------------------------------------------------------------
 
subroutine netcdf_deallocate_data_struct( NCFILE )
 
  type (t_netcdf4_file ) :: ncfile
 
  ! [ LOCALS ]
  type (t_netcdf_variable), pointer :: pnc_var
  integer (c_int) :: iindex
 
  do iindex=0, ncfile%iNumberOfVariables - 1
 
    pnc_var => ncfile%pNC_VAR(iindex)
 
    if (pnc_var%iNumberOfAttributes == 0 ) cycle
 
    if (associated( pnc_var%pNC_ATT ))  deallocate( pnc_var%pNC_ATT )
    pnc_var%pNC_ATT => null()
 
  enddo
 
  if (associated( ncfile%pNC_VAR ))  deallocate( ncfile%pNC_VAR )
  if (associated( ncfile%pNC_ATT ))  deallocate( ncfile%pNC_ATT )
  if (associated( ncfile%pNC_DIM ))  deallocate( ncfile%pNC_DIM )
 
  ncfile%pNC_VAR => null()
  ncfile%pNC_ATT => null()
  ncfile%pNC_DIM => null()
 
end subroutine netcdf_deallocate_data_struct
 
!----------------------------------------------------------------------
 
subroutine netcdf_nullify_data_struct( NCFILE )
 
  type (t_netcdf4_file ) :: ncfile
 
  ! [ LOCALS ]
 
  ncfile%pNC_VAR => null()
  ncfile%pNC_ATT => null()
  ncfile%pNC_DIM => null()
 
end subroutine netcdf_nullify_data_struct
 
!----------------------------------------------------------------------
 
subroutine nf_populate_dimension_struct( NCFILE )
 
  type (T_NETCDF4_FILE) :: NCFILE
  integer (c_int) :: iStat
  integer (c_int) :: iIndex
  character (len=256) :: sDimName
 
  call nf_trap( nc_inq_ndims(ncid=ncfile%iNCID, ndimsp=ncfile%iNumberOfDimensions), &
                __file__, __line__ )
 
  istat = 0
  if (associated(ncfile%pNC_DIM) ) deallocate(ncfile%pNC_DIM, stat=istat)
  call assert(istat == 0, "Could not deallocate memory for NC_DIM member in NC_FILE defined type", &
    __file__, __line__)
 
  allocate(ncfile%pNC_DIM( 0 : ncfile%iNumberOfDimensions-1), stat=istat )
  call assert(istat == 0, "Could not allocate memory for NC_DIM member in NC_FILE defined type", &
    __file__, __line__)
 
  ! netCDF 3 function
  call nf_trap( nc_inq_unlimdim(ncid=ncfile%iNCID, unlimdimidp=ncfile%iNC3_UnlimitedDimensionNumber), &
               __file__, __line__ )
 
  do iindex = 0, ncfile%iNumberOfDimensions-1
 
    call nf_trap(nc_inq_dim(ncid=ncfile%iNCID, dimid=iindex, &
      name=sdimname, &
      lenp=ncfile%pNC_DIM(iindex)%iNC_DimSize), __file__, __line__ )
 
    ncfile%pNC_DIM(iindex)%iNC_DimID = iindex
    ncfile%pNC_DIM(iindex)%sDimensionName = c_to_fortran_string(sdimname)
 
  enddo
 
end subroutine nf_populate_dimension_struct
 
!----------------------------------------------------------------------
 
subroutine nf_populate_variable_struct( NCFILE )
 
  type (T_NETCDF4_FILE) :: NCFILE
 
  type (T_NETCDF_ATTRIBUTE), pointer :: pNC_ATT
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
  integer (c_int) :: iStat
  integer (c_int) :: iIndex, iIndex2
  character (len=256) :: sVarName
 
  call nf_trap( nc_inq_nvars(ncid=ncfile%iNCID, nvarsp=ncfile%iNumberOfVariables), &
       __file__, __line__ )
 
  istat = 0
  if (associated(ncfile%pNC_VAR) ) deallocate(ncfile%pNC_VAR, stat=istat)
  call assert(istat == 0, "Could not deallocate memory for NC_VAR member in NC_FILE defined type", &
    __file__, __line__)
 
  allocate(ncfile%pNC_VAR( 0 : ncfile%iNumberOfVariables-1), stat=istat )
  call assert(istat == 0, "Could not allocate memory for NC_VAR member in NC_FILE defined type", &
    __file__, __line__)
 
  do iindex = 0, ncfile%iNumberOfVariables-1
 
    pnc_var => ncfile%pNC_VAR(iindex)
 
    call nf_trap(nc_inq_var(ncid=ncfile%iNCID, &
        varid=iindex, &
        name=svarname, &
        xtypep=pnc_var%iNC_VarType, &
        ndimsp=pnc_var%iNumberOfDimensions, &
        dimidsp=pnc_var%iNC_DimID, &
        nattsp=pnc_var%iNumberOfAttributes ), __file__, __line__ )
 
    pnc_var%iNC_VarID = iindex
    pnc_var%sVariableName = c_to_fortran_string(svarname)
 
    if( pnc_var%iNumberOfAttributes > 0 ) then
 
      if (associated(pnc_var%pNC_ATT) ) deallocate(pnc_var%pNC_ATT, stat=istat)
      call assert(istat == 0, "Could not deallocate memory for NC_ATT member within NC_VAR in NC_FILE defined type", &
        __file__, __line__)
 
      allocate( pnc_var%pNC_ATT( 0:pnc_var%iNumberOfAttributes - 1 ), stat = istat)
      call assert(istat == 0, "Could not allocate memory for NC_ATT member within NC_VAR in NC_FILE defined type", &
        __file__, __line__)
 
      do iindex2=0, pnc_var%iNumberOfAttributes - 1
 
        pnc_att => pnc_var%pNC_ATT(iindex2)
 
        call nf_populate_attribute_struct( ncfile=ncfile, pnc_att=pnc_att, &
          inc_varid=iindex, iattnum=iindex2 )
 
      enddo
 
    endif
 
  enddo
 
  call nf_trap( nc_inq_natts(ncid=ncfile%iNCID, ngattsp=ncfile%iNumberOfAttributes), &
       __file__, __line__ )
 
 
  if (associated(ncfile%pNC_ATT) )  deallocate(ncfile%pNC_ATT, stat=istat)
  call assert(istat == 0, "Could not deallocate memory for NC_ATT member within NC_FILE defined type", &
    __file__, __line__)
 
  allocate(ncfile%pNC_ATT(0:ncfile%iNumberOfAttributes - 1), stat=istat )
  call assert(istat == 0, "Could not allocate memory for NC_ATT member within NC_FILE defined type", &
    __file__, __line__)
 
  do iindex=0, ncfile%iNumberOfAttributes - 1
    pnc_att => ncfile%pNC_ATT(iindex)
 
    call nf_populate_attribute_struct( ncfile=ncfile, pnc_att=pnc_att, &
      inc_varid=nc_global, iattnum=iindex )
 
  enddo
 
end subroutine nf_populate_variable_struct
 
!----------------------------------------------------------------------
 
subroutine nf_populate_attribute_struct( NCFILE, pNC_ATT, iNC_VarID, iAttNum )
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  type (T_NETCDF_ATTRIBUTE), pointer :: pNC_ATT
  integer (c_int) :: iNC_VarID
  integer (c_int) :: iAttNum
 
  ![ LOCALS ]
  integer (c_int) :: iStat
  character (len=256) :: sAttName
  integer (c_int) :: iIndex
  integer (c_int) :: iLength
 
  call nf_trap( nc_inq_attname(ncid=ncfile%iNCID, &
    varid=inc_varid, &
    attnum=iattnum, &
    name=sattname), __file__, __line__ )
 
  pnc_att%sAttributeName = c_to_fortran_string(sattname)
 
  call nf_trap( nc_inq_att(ncid=ncfile%iNCID, &
    varid=inc_varid, &
    name=sattname, &
    xtypep=pnc_att%iNC_AttType, &
    lenp=pnc_att%iNC_AttSize), __file__, __line__ )
 
  ilength = int(pnc_att%iNC_AttSize, c_int)
 
  istat = 0
  allocate(pnc_att%sAttValue(0:ilength-1), stat=istat )
  call assert(istat==0, "INTERNAL PROGRAMMING ERROR - problem allocating memory", &
    __file__, __line__)
  pnc_att%sAttValue = ""
 
  select case(pnc_att%iNC_AttType)
 
    case (nc_char)
 
      call nf_trap( nc_get_att_text(ncid=ncfile%iNCID, &
        varid=inc_varid, &
        name=sattname, &
        ip=pnc_att%sAttValue), __file__, __line__ )
 
        pnc_att%sAttValue = c_to_fortran_string(pnc_att%sAttValue)
 
    case (nc_short)
 
      allocate(pnc_att%i2AttValue(0:ilength-1), stat=istat )
      call assert(istat==0, "INTERNAL PROGRAMMING ERROR - problem allocating memory", &
        __file__, __line__)
 
 
      call nf_trap( nc_get_att_short(ncid=ncfile%iNCID, &
        varid=inc_varid, &
        name=sattname, &
        ip=pnc_att%i2AttValue), __file__, __line__ )
 
      do iindex=1,ubound(pnc_att%i2AttValue,1)
        pnc_att%sAttValue(iindex) = ascharacter(pnc_att%i2AttValue(iindex))
      enddo
 
    case (nc_int)
 
      allocate(pnc_att%iAttValue(0:ilength-1), stat=istat )
      call assert(istat==0, "INTERNAL PROGRAMMING ERROR - problem allocating memory", &
        __file__, __line__)
 
      call nf_trap( nc_get_att_int(ncid=ncfile%iNCID, &
        varid=inc_varid, &
        name=sattname, &
        ip=pnc_att%iAttValue), __file__, __line__ )
 
      do iindex=1,ubound(pnc_att%iAttValue,1)
        pnc_att%sAttValue(iindex) = ascharacter(pnc_att%iAttValue(iindex))
      enddo
 
    case (nc_float)
 
      allocate(pnc_att%rAttValue(0:ilength-1), stat=istat )
      call assert(istat==0, "INTERNAL PROGRAMMING ERROR - problem allocating memory", &
        __file__, __line__)
 
 
      call nf_trap( nc_get_att_float(ncid=ncfile%iNCID, &
        varid=inc_varid, &
        name=sattname, &
        ip=pnc_att%rAttValue), __file__, __line__ )
 
      do iindex=1,ubound(pnc_att%rAttValue,1)
        pnc_att%sAttValue(iindex) = ascharacter(pnc_att%rAttValue(iindex))
      enddo
 
    case (nc_double)
 
      allocate(pnc_att%dpAttValue(0:ilength-1), stat=istat )
      call assert(istat==0, "INTERNAL PROGRAMMING ERROR - problem allocating memory", &
        __file__, __line__)
 
      call nf_trap( nc_get_att_double(ncid=ncfile%iNCID, &
        varid=inc_varid, &
        name=sattname, &
        ip=pnc_att%dpAttValue), __file__, __line__ )
 
      do iindex=1,ubound(pnc_att%dpAttValue,1)
        pnc_att%sAttValue(iindex) = ascharacter(pnc_att%dpAttValue(iindex))
      enddo
 
    case default
 
  end select
 
end subroutine nf_populate_attribute_struct
 
!----------------------------------------------------------------------
 
 
subroutine netcdf_get_variable_id_for_variable( NCFILE, variable_name, &
                                                variable_id      )
 
  type (t_netcdf4_file), intent(inout) :: ncfile
  character (len=*), intent(in)        :: variable_name
  integer (c_int), intent(out)    :: variable_id
 
  ! [ LOCALS ]
  integer (c_int) :: indx
  type (t_netcdf_variable), pointer  :: pnc_var
  logical (c_bool)              :: variable_was_found
 
  variable_was_found = false
 
  do indx=0, ncfile%iNumberOfVariables-1
    pnc_var => ncfile%pNC_VAR( indx )
    if ( associated( pnc_var ) ) then
      if ( pnc_var%sVariableName .strequal. variable_name ) then
        variable_was_found = true
        exit
      endif
    endif
  enddo
 
  if ( variable_was_found ) then
 
    variable_id = pnc_var%iNC_VarID
 
  else
 
    variable_id = -9999
 
  endif
 
end subroutine netcdf_get_variable_id_for_variable
 
!--------------------------------------------------------------------------------------------------
 
subroutine netcdf_get_attribute_list_for_variable( NCFILE, variable_name, &
                                                   attribute_name_list,   &
                                                   attribute_value_list )
 
  type (t_netcdf4_file), intent(inout) :: ncfile
  character (len=*), intent(in)        :: variable_name
  type (fstring_list_t), intent(out)    :: attribute_name_list
  type (fstring_list_t), intent(out)    :: attribute_value_list
 
  ! [ LOCALS ]
  integer (c_int) :: indx
  type (t_netcdf_attribute), pointer :: pnc_att
  type (t_netcdf_variable), pointer  :: pnc_var
  logical (c_bool)              :: variable_was_found
  character (len=256)                :: tempstring
 
  variable_was_found = false
 
  call attribute_name_list%clear
  call attribute_value_list%clear
 
  do indx=0, ncfile%iNumberOfVariables-1
    pnc_var => ncfile%pNC_VAR( indx )
    if ( associated( pnc_var ) ) then
      if ( pnc_var%sVariableName .strequal. variable_name ) then
        variable_was_found = true
        exit
      endif
    endif
  enddo
 
  if ( variable_was_found ) then
 
    do indx=0, pnc_var%iNumberOfAttributes-1
      pnc_att => pnc_var%pNC_ATT( indx )
      if (associated( pnc_att) ) then
 
        tempstring = pnc_att%sAttributeName
        call attribute_name_list%append( trim( tempstring ) )
 
        select case ( pnc_att%iNC_AttType )
 
          case ( nc_char )
 
            tempstring = pnc_att%sAttValue(0)
            call attribute_value_list%append( trim( tempstring ) )
 
          case ( nc_short )
 
            tempstring = ascharacter( pnc_att%i2AttValue(0) )
            call attribute_value_list%append( trim( tempstring ) )
 
          case ( nc_int )
 
            tempstring = ascharacter( pnc_att%iAttValue(0) )
            call attribute_value_list%append( trim( tempstring ) )
 
          case ( nc_float )
 
            tempstring = ascharacter( pnc_att%rAttValue(0) )
            call attribute_value_list%append( trim( tempstring ) )
 
          case ( nc_double )
 
            tempstring = ascharacter( pnc_att%dpAttValue(0) )
            call attribute_value_list%append( trim( tempstring ) )
 
        end select
 
      endif
 
    enddo
 
  else
 
    call attribute_name_list%append("<NA>")
    call attribute_value_list%append("<NA>")
 
  endif
 
end subroutine netcdf_get_attribute_list_for_variable
 
!----------------------------------------------------------------------
 
subroutine netcdf_set_coordinate_tolerance(NCFILE, tolerance)
 
  type (t_netcdf4_file)        :: ncfile
  real (c_double), intent(in)  :: tolerance
 
  ncfile%rCoordinateTolerance = tolerance
 
end subroutine netcdf_set_coordinate_tolerance
 
!----------------------------------------------------------------------
 
subroutine netcdf_get_variable_list( NCFILE, variable_list )
 
  type (t_netcdf4_file) :: ncfile
  type (fstring_list_t)  :: variable_list
 
  ! [ LOCALS ]
  integer (c_int) :: indx
  type (t_netcdf_variable), pointer :: pnc_var
 
  call variable_list%clear
 
  do indx=0, ncfile%iNumberOfVariables-1
    pnc_var => ncfile%pNC_VAR( indx )
    if ( associated( pnc_var ) )  call variable_list%append( pnc_var%sVariableName )
  enddo
 
end subroutine netcdf_get_variable_list
 
!----------------------------------------------------------------------
 
function netcdf_update_time_starting_index(NCFILE, iJulianDay)  result(lDateTimeFound)
 
  type (t_netcdf4_file) :: ncfile
  integer (c_int) :: ijulianday
  logical (c_bool) :: ldatetimefound
 
  ! [ LOCALS ]
 
  ncfile%iStart(nc_time) = nf_julian_day_to_index_adj( ncfile=ncfile, &
                                     rjulianday=real(ijulianday, c_double ) )
 
  if (ncfile%iStart(nc_time) < 0) then
    ncfile%iStart(nc_time) = 0
    ldatetimefound = false
  else
    ldatetimefound = true
  endif
 
end function netcdf_update_time_starting_index
 
!----------------------------------------------------------------------
 
subroutine netcdf_get_variable_slice(NCFILE, rValues, dpValues, iValues)
 
  type (t_netcdf4_file), intent(inout)       :: ncfile
  real (c_float), dimension(:,:), optional   :: rvalues
  real (c_double), dimension(:,:), optional  :: dpvalues
  integer (c_int), dimension(:,:), optional  :: ivalues
 
  real (c_float), allocatable     :: rtempvals(:,:)
  integer (c_short), allocatable  :: i2tempvals(:,:)
  integer (c_int), allocatable    :: itempvals(:,:)
  real (c_double), allocatable    :: dptempvals(:,:)
  integer (c_int)                 :: nrow, ncol
 
  nrow = 0
  ncol = 0
 
  if (present( rvalues) ) then
    nrow = size(rvalues,2)
    ncol = size(rvalues,1)
  elseif (present( dpvalues) ) then
    nrow = size(dpvalues,2)
    ncol = size(dpvalues,1)
  elseif (present( ivalues) ) then
    nrow = size(ivalues,2)
    ncol = size(ivalues,1)
  else
    call die("Internal programming error: unhandled data type",                &
      __file__, __line__)
  endif
 
  if (   ncfile%iVarType(nc_z) == nc_short                                     &
    .or. ncfile%iVarType(nc_z) == nc_ushort) then
 
    allocate( i2tempvals(ncol, nrow))
 
    call nf_get_variable_slice_short(ncfile, i2tempvals)
 
    if (present(rvalues))  rvalues = asfloat(i2tempvals)
    if (present(ivalues))  ivalues = i2tempvals
    if (present(dpvalues))  dpvalues = asdouble(i2tempvals)
 
  elseif (   ncfile%iVarType(nc_z) == nc_int                                   &
    .or. ncfile%iVarType(nc_z) == nc_uint) then
 
    allocate( itempvals(ncol, nrow))
    call nf_get_variable_slice_int(ncfile, itempvals)
 
    if (present(rvalues))  rvalues = asfloat(itempvals)
    if (present(ivalues))  ivalues = itempvals
    if (present(dpvalues))  dpvalues = asdouble(itempvals)
 
  elseif (ncfile%iVarType(nc_z) == nc_float) then
 
    allocate( rtempvals(ncol, nrow))
    call nf_get_variable_slice_float(ncfile, rtempvals)
 
    if (present(rvalues))  rvalues = rtempvals
    if (present(ivalues))  ivalues = asint(rtempvals)
    if (present(dpvalues))  dpvalues = asdouble(rtempvals)
 
  elseif (ncfile%iVarType(nc_z) == nc_double) then
 
    allocate( dptempvals(ncol, nrow))
    call nf_get_variable_slice_double(ncfile, dptempvals)
 
    if (present(rvalues))  rvalues = asfloat(dptempvals)
    if (present(ivalues))  ivalues = asint(dptempvals)
    if (present(dpvalues))  dpvalues = dptempvals
 
  else
 
    call warn("Failed to find a method to retrieve data of the given type.", __file__, __line__)
 
    print *, 'var type: ',ncfile%iVarType(nc_z)
 
  endif
 
end subroutine netcdf_get_variable_slice
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_slice_short(NCFILE, i2Values)
 
  type (T_NETCDF4_FILE) :: NCFILE
  integer (c_short), intent(inout) :: i2Values(:,:)
 
  ! [ LOCALS ]
  !! dimension #1 = column (iNX)
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
 
  integer (c_short), dimension(size(i2Values,2) * size(i2Values,1)) :: iTemp
 
  integer (c_int) :: iRow, iCol, iIndex
  integer (c_int) :: iFromRow, iToRow, iByRow
  integer (c_int) :: iFromCol, iToCol, iByCol
 
  ifromrow = ncfile%iRowIter(nc_first)
  itorow = ncfile%iRowIter(nc_last)
  ibyrow = ncfile%iRowIter(nc_by)
 
  ifromcol = ncfile%iColIter(nc_first)
  itocol = ncfile%iColIter(nc_last)
  ibycol = ncfile%iColIter(nc_by)
 
  pnc_var => ncfile%pNC_VAR(nf_return_varindex( ncfile, ncfile%iVarID(nc_z)) )
 
  select case (ncfile%sVariableOrder)
 
    case ("txy")    ! time, col, row
 
      call nf_get_variable_array_as_vector_short(ncfile=ncfile, &
        inc_varid=ncfile%iVarID(nc_z), &
        inc_start=[ncfile%iStart(nc_time), ncfile%iStart(nc_x), ncfile%iStart(nc_y)], &
        inc_count=[ncfile%iCount(nc_time), ncfile%iCount(nc_x), ncfile%iCount(nc_y)], &
        inc_stride=[ncfile%iStride(nc_time), ncfile%iStride(nc_x), ncfile%iStride(nc_y)], &
        inc_vars=itemp)
 
        iindex = 0
        do icol=ifromcol, itocol, ibycol
          do irow=ifromrow, itorow, ibyrow
            iindex = iindex + 1
            i2values(icol,irow) = itemp(iindex)
          enddo
        enddo
 
    case ("tyx")    ! time, row, col
 
      call nf_get_variable_array_as_vector_short(ncfile=ncfile, &
        inc_varid=ncfile%iVarID(nc_z), &
        inc_start=[ncfile%iStart(nc_time), ncfile%iStart(nc_y), ncfile%iStart(nc_x)], &
        inc_count=[ncfile%iCount(nc_time), ncfile%iCount(nc_y), ncfile%iCount(nc_x)], &
        inc_stride=[ncfile%iStride(nc_time), ncfile%iStride(nc_y), ncfile%iStride(nc_x)], &
        inc_vars=itemp)
 
        iindex = 0
        do irow=ifromrow, itorow, ibyrow
          do icol=ifromcol, itocol, ibycol
            iindex = iindex + 1
            i2values(icol,irow) = itemp(iindex)
          enddo
        enddo
 
    case default
 
      call warn("INTERNAL PROGRAMMING ERROR: Unhandled select case. Program will probably fail.", __file__, __line__)
 
  end select
 
end subroutine nf_get_variable_slice_short
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_slice_int(NCFILE, iValues)
 
  type (T_NETCDF4_FILE) :: NCFILE
  integer (c_int), dimension(:,:) :: iValues
 
  ! [ LOCALS ]
  !! dimension #1 = column (iNX)
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
 
  integer (c_int), dimension(size(iValues,2) * size(iValues,1)) :: iTemp
 
  integer (c_int) :: iRow, iCol, iIndex
  integer (c_int) :: iFromRow, iToRow, iByRow
  integer (c_int) :: iFromCol, iToCol, iByCol
 
  ifromrow = ncfile%iRowIter(nc_first)
  itorow = ncfile%iRowIter(nc_last)
  ibyrow = ncfile%iRowIter(nc_by)
 
  ifromcol = ncfile%iColIter(nc_first)
  itocol = ncfile%iColIter(nc_last)
  ibycol = ncfile%iColIter(nc_by)
 
  pnc_var => ncfile%pNC_VAR(nf_return_varindex( ncfile, ncfile%iVarID(nc_z)) )
 
  select case (ncfile%sVariableOrder)
 
    case ("txy")    ! time, col, row
 
      call nf_get_variable_array_as_vector_int(ncfile=ncfile, &
        inc_varid=ncfile%iVarID(nc_z), &
        inc_start=[ncfile%iStart(nc_time), ncfile%iStart(nc_x), ncfile%iStart(nc_y)], &
        inc_count=[ncfile%iCount(nc_time), ncfile%iCount(nc_x), ncfile%iCount(nc_y)], &
        inc_stride=[ncfile%iStride(nc_time), ncfile%iStride(nc_x), ncfile%iStride(nc_y)], &
        inc_vars=itemp)
 
        iindex = 0
        do icol=ifromcol, itocol, ibycol
          do irow=ifromrow, itorow, ibyrow
            iindex = iindex + 1
            ivalues(icol,irow) = itemp(iindex)
          enddo
        enddo
 
    case ("tyx")    ! time, row, col
 
      call nf_get_variable_array_as_vector_int(ncfile=ncfile, &
        inc_varid=ncfile%iVarID(nc_z), &
        inc_start=[ncfile%iStart(nc_time), ncfile%iStart(nc_y), ncfile%iStart(nc_x)], &
        inc_count=[ncfile%iCount(nc_time), ncfile%iCount(nc_y), ncfile%iCount(nc_x)], &
        inc_stride=[ncfile%iStride(nc_time), ncfile%iStride(nc_y), ncfile%iStride(nc_x)], &
        inc_vars=itemp)
 
        iindex = 0
        do irow=ifromrow, itorow, ibyrow
          do icol=ifromcol, itocol, ibycol
            iindex = iindex + 1
            ivalues(icol,irow) = itemp(iindex)
          enddo
        enddo
 
    case default
 
      call warn("INTERNAL PROGRAMMING ERROR: Unhandled select case. Program will probably fail.", __file__, __line__)
 
 
  end select
 
end subroutine nf_get_variable_slice_int
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_slice_float(NCFILE, rValues)
 
  type (T_NETCDF4_FILE) :: NCFILE
  real (c_float), dimension(:,:) :: rValues
 
  ! [ LOCALS ]
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
  real (c_float), dimension(size(rValues,2) * size(rValues,1)) :: rTemp
  integer (c_int) :: iRow, iCol, iIndex
  integer (c_int) :: iFromRow, iToRow, iByRow
  integer (c_int) :: iFromCol, iToCol, iByCol
 
  ifromrow = ncfile%iRowIter(nc_first)
  itorow = ncfile%iRowIter(nc_last)
  ibyrow = ncfile%iRowIter(nc_by)
 
  ifromcol = ncfile%iColIter(nc_first)
  itocol = ncfile%iColIter(nc_last)
  ibycol = ncfile%iColIter(nc_by)
 
  pnc_var => ncfile%pNC_VAR(nf_return_varindex( ncfile, ncfile%iVarID(nc_z)) )
 
  select case (ncfile%sVariableOrder)
 
    case ("txy")    ! time, col, row
 
      call nf_get_variable_array_as_vector_float(ncfile=ncfile, &
        inc_varid=ncfile%iVarID(nc_z), &
        inc_start=[ncfile%iStart(nc_time), ncfile%iStart(nc_x), ncfile%iStart(nc_y)], &
        inc_count=[ncfile%iCount(nc_time), ncfile%iCount(nc_x), ncfile%iCount(nc_y)], &
        inc_stride=[ncfile%iStride(nc_time), ncfile%iStride(nc_x), ncfile%iStride(nc_y)], &
        rnc_vars=rtemp)
 
      iindex = 0
      do icol=ifromcol, itocol, ibycol
        do irow=ifromrow, itorow, ibyrow
          iindex = iindex + 1
          rvalues(icol,irow) = rtemp(iindex)
        enddo
      enddo
 
    case ("tyx")    ! time, row, col
 
      call nf_get_variable_array_as_vector_float(ncfile=ncfile, &
        inc_varid=ncfile%iVarID(nc_z), &
        inc_start=[ncfile%iStart(nc_time), ncfile%iStart(nc_y), ncfile%iStart(nc_x)], &
        inc_count=[ncfile%iCount(nc_time), ncfile%iCount(nc_y), ncfile%iCount(nc_x)], &
        inc_stride=[ncfile%iStride(nc_time), ncfile%iStride(nc_y), ncfile%iStride(nc_x)], &
        rnc_vars=rtemp)
 
      iindex = 0
      do irow=ifromrow, itorow, ibyrow
        do icol=ifromcol, itocol, ibycol
          iindex = iindex + 1
          rvalues(icol,irow) = rtemp(iindex)
        enddo
      enddo
 
    case default
 
      call warn("INTERNAL PROGRAMMING ERROR: Unhandled select case. Program will probably fail.", __file__, __line__)
 
  end select
 
end subroutine nf_get_variable_slice_float
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_slice_double(NCFILE, dpValues)
 
  type (T_NETCDF4_FILE) :: NCFILE
  real (c_double), dimension(:,:) :: dpValues
 
  ! [ LOCALS ]
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
  real (c_double), dimension(size(dpValues,2) * size(dpValues,1)) :: dpTemp
  integer (c_int) :: iRow, iCol, iIndex
  integer (c_int) :: iFromRow, iToRow, iByRow
  integer (c_int) :: iFromCol, iToCol, iByCol
 
  ifromrow = ncfile%iRowIter(nc_first)
  itorow = ncfile%iRowIter(nc_last)
  ibyrow = ncfile%iRowIter(nc_by)
 
  ifromcol = ncfile%iColIter(nc_first)
  itocol = ncfile%iColIter(nc_last)
  ibycol = ncfile%iColIter(nc_by)
 
  pnc_var => ncfile%pNC_VAR(nf_return_varindex( ncfile, ncfile%iVarID(nc_z)) )
 
  select case (ncfile%sVariableOrder)
 
    case ("txy")    ! time, col, row
 
      call nf_get_variable_array_as_vector_double(ncfile=ncfile, &
        inc_varid=ncfile%iVarID(nc_z), &
        inc_start=[ncfile%iStart(nc_time), ncfile%iStart(nc_x), ncfile%iStart(nc_y)], &
        inc_count=[ncfile%iCount(nc_time), ncfile%iCount(nc_x), ncfile%iCount(nc_y)], &
        inc_stride=[ncfile%iStride(nc_time), ncfile%iStride(nc_x), ncfile%iStride(nc_y)], &
        dpnc_vars=dptemp)
 
      iindex = 0
      do icol=ifromcol, itocol, ibycol
        do irow=ifromrow, itorow, ibyrow
          iindex = iindex + 1
          dpvalues(icol,irow) = dptemp(iindex)
        enddo
      enddo
 
    case ("tyx")    ! time, row, col
 
      call nf_get_variable_array_as_vector_double(ncfile=ncfile, &
        inc_varid=ncfile%iVarID(nc_z), &
        inc_start=[ncfile%iStart(nc_time), ncfile%iStart(nc_y), ncfile%iStart(nc_x)], &
        inc_count=[ncfile%iCount(nc_time), ncfile%iCount(nc_y), ncfile%iCount(nc_x)], &
        inc_stride=[ncfile%iStride(nc_time), ncfile%iStride(nc_y), ncfile%iStride(nc_x)], &
        dpnc_vars=dptemp)
 
      iindex = 0
      do irow=ifromrow, itorow, ibyrow
        do icol=ifromcol, itocol, ibycol
          iindex = iindex + 1
          dpvalues(icol,irow) = dptemp(iindex)
        enddo
      enddo
 
    case default
 
      call warn("INTERNAL PROGRAMMING ERROR: Unhandled select case. Program will probably fail.", __file__, __line__)
 
  end select
 
end subroutine nf_get_variable_slice_double
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_vector_short(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, iNC_Vars)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t) :: iNC_Start
  integer (c_size_t) :: iNC_Count
  integer (c_size_t) :: iNC_Stride
  integer (c_short), dimension(:) :: iNC_Vars
 
  call nf_trap(nc_get_vars_short(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=inc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_vector_short
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_array_short(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, iNC_Vars)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t), dimension(:) :: iNC_Start
  integer (c_size_t), dimension(:) :: iNC_Count
  integer (c_size_t), dimension(:) :: iNC_Stride
  integer (c_short), dimension(:,:) :: iNC_Vars
 
  call nf_trap(nc_get_vars_short(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=inc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_array_short
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_array_as_vector_short(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, iNC_Vars)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t), dimension(:) :: iNC_Start
  integer (c_size_t), dimension(:) :: iNC_Count
  integer (c_size_t), dimension(:) :: iNC_Stride
  integer (c_short), dimension(:) :: iNC_Vars
 
  call nf_trap(nc_get_vars_short(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=inc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_array_as_vector_short
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_array_as_vector_int(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, iNC_Vars)
 
  type (T_NETCDF4_FILE), intent(inout)     :: NCFILE
  integer (c_int)                     :: iNC_VarID
  integer (c_size_t), dimension(:)    :: iNC_Start
  integer (c_size_t), dimension(:)    :: iNC_Count
  integer (c_size_t), dimension(:) :: iNC_Stride
  integer (c_int), dimension(:)       :: iNC_Vars
 
  call nf_trap(nc_get_vars_int(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=inc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_array_as_vector_int
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_vector_int(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, iNC_Vars)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t) :: iNC_Start
  integer (c_size_t) :: iNC_Count
  integer (c_size_t) :: iNC_Stride
  integer (c_int), dimension(:) :: iNC_Vars
 
  call nf_trap(nc_get_vars_int(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=inc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_vector_int
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_vector_double(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, dpNC_Vars)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t) :: iNC_Start
  integer (c_size_t) :: iNC_Count
  integer (c_size_t) :: iNC_Stride
  real (c_double), dimension(:) :: dpNC_Vars
 
  call nf_trap(nc_get_vars_double(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=dpnc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_vector_double
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_array_double(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, dpNC_Vars)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t), dimension(:) :: iNC_Start
  integer (c_size_t), dimension(:) :: iNC_Count
  integer (c_size_t), dimension(:) :: iNC_Stride
  real (c_double), dimension(:,:) :: dpNC_Vars
 
  call nf_trap(nc_get_vars_double(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=dpnc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_array_double
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_array_as_vector_double(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, dpNC_Vars)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t), dimension(:) :: iNC_Start
  integer (c_size_t), dimension(:) :: iNC_Count
  integer (c_size_t), dimension(:) :: iNC_Stride
  real (c_double), dimension(:) :: dpNC_Vars
 
  call nf_trap(nc_get_vars_double(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=dpnc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_array_as_vector_double
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_vector_float(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, rNC_Vars )
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t) :: iNC_Start
  integer (c_size_t) :: iNC_Count
  integer (c_size_t) :: iNC_Stride
  real (c_float), dimension(:) :: rNC_Vars
 
  call nf_trap(nc_get_vars_float(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=rnc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_vector_float
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_array_float(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, rNC_Vars )
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t), dimension(:) :: iNC_Start
  integer (c_size_t), dimension(:) :: iNC_Count
  integer (c_size_t), dimension(:) :: iNC_Stride
  real (c_float), dimension(:,:) :: rNC_Vars
 
  call nf_trap(nc_get_vars_float(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=rnc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_array_float
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_array_as_vector_float(NCFILE, iNC_VarID, iNC_Start, iNC_Count, &
   iNC_Stride, rNC_Vars )
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
  integer (c_int) :: iNC_VarID
  integer (c_size_t), dimension(:) :: iNC_Start
  integer (c_size_t), dimension(:) :: iNC_Count
  integer (c_size_t), dimension(:) :: iNC_Stride
  real (c_float), dimension(:) :: rNC_Vars
 
  call nf_trap(nc_get_vars_float(ncid=ncfile%iNCID, &
       varid=inc_varid, &
       startp=[inc_start], &
       countp=[inc_count], &
       stridep=[inc_stride], &
       vars=rnc_vars), __file__, __line__ )
 
end subroutine nf_get_variable_array_as_vector_float
 
!----------------------------------------------------------------------
 
subroutine netcdf_dump_cdl(NCFILE, iLU)
 
  type (t_netcdf4_file ) :: ncfile
  type (t_netcdf_attribute), pointer :: pnc_att
  type (t_netcdf_variable), pointer :: pnc_var
  type (t_netcdf_dimension), pointer :: pnc_dim
  integer :: ilu
  character (len=256) :: sbuf, sbuf2
  character (len=256) :: sdimname
  integer (c_int) :: idimid
  integer (c_int) :: iubound
 
  integer :: iindex, iindex3, iindex4
 
  sbuf=""; sbuf2=""
 
  write(unit=ilu, fmt="(a)") "netcdf "//trim(ncfile%sFilename)//" {"
  write(unit=ilu, fmt="(a)") "  dimensions:"
 
  do iindex = 0, ncfile%iNumberOfDimensions - 1
    write(unit=ilu, fmt="(4x,a, ' = ', i0, ';')") trim(ncfile%pNC_DIM(iindex)%sDimensionName), &
      ncfile%pNC_DIM(iindex)%iNC_DimSize
  enddo
 
  do iindex = 0, ncfile%iNumberOfVariables - 1
 
    pnc_var => ncfile%pNC_VAR(iindex)
 
    if(pnc_var%iNumberOfDimensions > 0) then
 
      sbuf = ' ('
 
      iubound = pnc_var%iNumberOfDimensions - 1
      do iindex3 = 0, iubound
 
        idimid = pnc_var%iNC_DimID(iindex3)
 
        call assert(idimid >=0 .and. &
          idimid <= ubound( ncfile%pNC_DIM, 1 ), &
          "INTERNAL PROGRAMMING ERROR -- iDimID out of bounds", &
          __file__, __line__)
 
        pnc_dim => ncfile%pNC_DIM(idimid)
        sdimname = pnc_dim%sDimensionName
 
        write(sbuf2, fmt="(i12)") pnc_dim%iNC_DimSize
        sbuf = trim(sbuf)//trim(pnc_dim%sDimensionName)//"=" &
           //trim(adjustl(sbuf2))
 
        if (iindex3 /= iubound) sbuf = trim(sbuf)//", "
 
      enddo
 
      sbuf = trim(sbuf)//')'
 
    else
 
      sbuf = ""
 
    endif
 
    sbuf = trim(netcdf_data_type(pnc_var%iNC_VarType)) &
       //" "//trim(pnc_var%sVariableName)//sbuf//";"
 
    write(unit=ilu, fmt="(2x,a)") trim(sbuf)
 
    iubound = pnc_var%iNumberOfAttributes - 1
    do iindex3 = 0, iubound
 
      pnc_att => ncfile%pNC_VAR(iindex)%pNC_ATT(iindex3)
 
      sbuf = trim(pnc_var%sVariableName)//":"//trim(pnc_att%sAttributeName )//" ="
 
      do iindex4=0, ubound(pnc_att%sAttValue, 1)
 
          sbuf = trim(sbuf)//" "//trim(pnc_att%sAttValue(iindex4))
 
      enddo
 
      sbuf=trim(sbuf)//"; // "//trim(netcdf_data_type(pnc_att%iNC_AttType) )
 
      write(unit=ilu, fmt="(4x,a)") trim(sbuf)
 
    enddo
 
  enddo
 
  do iindex = 0, ncfile%iNumberOfAttributes - 1
 
    pnc_att => ncfile%pNC_ATT(iindex)
 
    sbuf = ":"//trim(pnc_att%sAttributeName )//" ="
 
    do iindex4=0, ubound(pnc_att%sAttValue, 1)
 
        sbuf = trim(sbuf)//" "//trim(pnc_att%sAttValue(iindex4))
 
    enddo
 
    sbuf=trim(sbuf)//"; // "//trim(netcdf_data_type(pnc_att%iNC_AttType) )
 
    write(unit=ilu, fmt="(a)") trim(sbuf)
 
  enddo
 
  write(unit=ilu, fmt="(a,/,/)") "}"
 
 
end subroutine netcdf_dump_cdl
 
!----------------------------------------------------------------------
 
function nf_get_first_and_last(NCFILE, iVarIndex)  result(dpValues)
 
  type (t_netcdf4_file ) :: ncfile
  integer (c_int) :: ivarindex
  real (c_double), dimension(0:1) :: dpvalues
 
  ! [ LOCALS ]
  type (t_netcdf_variable), pointer :: pnc_var
  integer (c_int) :: idimsize
  integer (c_size_t) :: istride
  integer (c_size_t) :: icount
  integer (c_short), dimension(0:1) :: spvalues
  integer (c_int), dimension(0:1) :: ipvalues
  real (c_float), dimension(0:1) :: rpvalues
 
  call assert (ivarindex >= lbound(ncfile%pNC_VAR,1) &
    .and. ivarindex <= ubound(ncfile%pNC_VAR,1), &
      "INTERNAL PROGRAMMING ERROR - index out of bounds NC_FILE%pNC_VAR" &
      //"~Offending index value: "//trim(ascharacter(ivarindex)), &
      __file__, __line__)
 
  pnc_var => ncfile%pNC_VAR(ivarindex)
  idimsize = int(nf_return_dimsize(ncfile, pnc_var%iNC_DimID(0) ), c_int)
 
  if (idimsize > 1) then
    icount = 2_c_size_t
    istride = int(idimsize, c_size_t) - 1_c_size_t
  else
    icount = 1_c_size_t
    istride = 1_c_size_t
  endif
 
  select case (pnc_var%iNC_VarType )
 
    case (nc_short)
 
      call nf_get_variable_vector_short(ncfile=ncfile, &
        inc_varid=pnc_var%iNC_VarID, &
        inc_start=0_c_size_t, &
        inc_count=icount, &
        inc_stride=istride, &
        inc_vars=spvalues)
 
      dpvalues = real(spvalues, c_double)
 
    case (nc_int)
 
      call nf_get_variable_vector_int(ncfile=ncfile, &
        inc_varid=pnc_var%iNC_VarID, &
        inc_start=0_c_size_t, &
        inc_count=icount, &
        inc_stride=istride, &
        inc_vars=ipvalues)
 
      dpvalues = real(ipvalues, c_double)
 
    case (nc_float)
 
      call nf_get_variable_vector_float(ncfile=ncfile, &
        inc_varid=pnc_var%iNC_VarID, &
        inc_start=0_c_size_t, &
        inc_count=icount, &
        inc_stride=istride, &
        rnc_vars=rpvalues)
 
      dpvalues = real(rpvalues, c_double)
 
    case (nc_double)
 
      call nf_get_variable_vector_double(ncfile=ncfile, &
        inc_varid=pnc_var%iNC_VarID, &
        inc_start=0_c_size_t, &
        inc_count=icount, &
        inc_stride=istride, &
        dpnc_vars=dpvalues)
 
    case default
 
      call warn("INTERNAL PROGRAMMING ERROR: Unhandled select case. Program will probably fail.", __file__, __line__)
 
  end select
 
  !> if there is only one day of data in this netCDF file, the
  !> first day equals the last day
  if (icount == 1) dpvalues(nc_last) = dpvalues(nc_first)
 
end function nf_get_first_and_last
 
!----------------------------------------------------------------------
 
subroutine nf_calculate_time_range(NCFILE)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
 
  ncfile%iOriginJD = julian_day(ncfile%iOriginYear, &
    ncfile%iOriginMonth, ncfile%iOriginDay)
 
  ncfile%iFirstDayJD = ncfile%iOriginJD + int(ncfile%dpFirstAndLastTimeValues(nc_first), c_int)
  ncfile%iLastDayJD = ncfile%iOriginJD + int(ncfile%dpFirstAndLastTimeValues(nc_last), c_int)
 
end subroutine nf_calculate_time_range
 
!----------------------------------------------------------------------
 
subroutine nf_get_time_units(NCFILE)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
 
  ! [ LOCALS ]
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
  character (len=256)   :: sDateTime
  character (len=256)   :: sItem
  integer (c_int)  :: iIndex
  logical (c_bool) :: lFound
  integer (c_int)  :: iStat
  real (c_float)   :: fTempVal
 
  call assert(ncfile%iVarID(nc_time) >= 0, "INTERNAL PROGRAMMING ERROR -- " &
    //"nf_get_time_units must be called only after a call is made to ~" &
    //"netcdf_get_variable_ids", __file__, __line__)
 
  pnc_var => ncfile%pNC_VAR(ncfile%iVarID(nc_time) )
 
  lfound = false
 
  do iindex=0, pnc_var%iNumberOfAttributes - 1
 
    if ( pnc_var%pNC_ATT(iindex)%sAttributeName .strequal. "units"   ) then
      lfound = true
      exit
    endif
 
  enddo
 
  call assert (lfound, "Failed to find the 'units' attribute associated " &
    //"with time variable "//dquote(pnc_var%sVariableName), &
    __file__, __line__)
 
  sdatetime = pnc_var%pNC_ATT(iindex)%sAttValue(0)
 
  call chomp(sdatetime, sitem)    ! should be "days"
  call chomp(sdatetime, sitem)    ! should be "since"
 
  call chomp(sdatetime, sitem, "/-")
  read(sitem, *) ncfile%iOriginYear
 
  call chomp(sdatetime, sitem, "/-")
  read(sitem, *) ncfile%iOriginMonth
 
  !> @todo this does not appear to have the fix that was applied to the master swb branch to
  !!       deal with cases where no time values are given at all
 
  read(sdatetime, *) ncfile%iOriginDay
 
  call chomp(sdatetime, sitem, ":")
  read(sitem, *, iostat=istat) ncfile%iOriginHH
  if (istat /=0) ncfile%iOriginHH = 0
 
  call chomp(sdatetime, sitem, ":")
  read(sitem, *, iostat=istat) ncfile%iOriginMM
  if (istat /=0) ncfile%iOriginMM = 0
 
  ! changed this to a real value, since some data providers encode the date and
  ! time as YYYY-MM-DD HH-MM-SS.S
  read(sdatetime, *, iostat=istat) ftempval
  if (istat ==0) then
    ncfile%iOriginSS = int(ftempval, c_int)
  else
    ncfile%iOriginSS = 0
  endif
 
end subroutine nf_get_time_units
 
!----------------------------------------------------------------------
 
subroutine nf_get_xyz_units(NCFILE)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
 
  ! [ LOCALS ]
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
  integer (c_int) :: iIndex, iIndex2
  logical (c_bool) :: lFound
 
  do iindex = nc_y, nc_z
 
    call assert(ncfile%iVarID(iindex) >= 0, "INTERNAL PROGRAMMING ERROR -- " &
    //"nc_get_XYZ_units must be called only after a call is made to ~" &
    //"netcdf_get_variable_ids", __file__, __line__)
 
    pnc_var => ncfile%pNC_VAR(ncfile%iVarID(iindex) )
 
    lfound = false
 
    do iindex2=0, pnc_var%iNumberOfAttributes - 1
 
      if ( pnc_var%pNC_ATT(iindex2)%sAttributeName .strequal. "units" ) then
        lfound = true
        exit
      endif
 
    enddo
 
    if (lfound) then
      ncfile%sVarUnits(iindex) = trim(pnc_var%pNC_ATT(iindex2)%sAttValue(0))
    endif
 
  enddo
 
end subroutine nf_get_xyz_units
 
!----------------------------------------------------------------------
 
subroutine nf_get_scale_and_offset(NCFILE)
 
  type (T_NETCDF4_FILE), intent(inout) :: NCFILE
 
  ! [ LOCALS ]
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
  integer (c_int) :: iIndex
  logical (c_bool) :: lFound
  character (len=32) :: sBuf
 
  pnc_var => ncfile%pNC_VAR(ncfile%iVarID(nc_z) )
 
  lfound = false
 
  do iindex=0, pnc_var%iNumberOfAttributes - 1
 
    if ( pnc_var%pNC_ATT(iindex)%sAttributeName .strequal. "scale_factor" ) then
      lfound = true
      exit
    endif
 
  enddo
 
  if (lfound) then
    if (allocated( pnc_var%pNC_ATT(iindex)%i2AttValue )) then
      ncfile%rScaleFactor(nc_z) = asfloat( pnc_var%pNC_ATT(iindex)%i2AttValue(0) )
    elseif (allocated( pnc_var%pNC_ATT(iindex)%iAttValue )) then
      ncfile%rScaleFactor(nc_z) = asfloat( pnc_var%pNC_ATT(iindex)%iAttValue(0) )
    elseif (allocated( pnc_var%pNC_ATT(iindex)%rAttValue )) then
      ncfile%rScaleFactor(nc_z) = pnc_var%pNC_ATT(iindex)%rAttValue(0)
    elseif (allocated( pnc_var%pNC_ATT(iindex)%dpAttValue )) then
      ncfile%rScaleFactor(nc_z) = asfloat( pnc_var%pNC_ATT(iindex)%dpAttValue(0) )
    elseif (len_trim(pnc_var%pNC_ATT(iindex)%sAttValue(0)) > 0) then
      sbuf = trim(pnc_var%pNC_ATT(iindex)%sAttValue(0) )
      read(sbuf,*) ncfile%rScaleFactor(nc_z)
    else
      call die("Error reading the 'scale_factor' attribute from the netCDF file" &
        //dquote(ncfile%sFilename))
    endif
  endif
 
  !> Now repeat the process for "add_offset" attribute
  lfound = false
 
  do iindex=0, pnc_var%iNumberOfAttributes - 1
 
    if ( pnc_var%pNC_ATT(iindex)%sAttributeName .strequal. "add_offset" ) then
      lfound = true
      exit
    endif
 
  enddo
 
  if (lfound) then
    if (allocated( pnc_var%pNC_ATT(iindex)%i2AttValue )) then
      ncfile%rAddOffset(nc_z) = asfloat( pnc_var%pNC_ATT(iindex)%i2AttValue(0) )
    elseif (allocated( pnc_var%pNC_ATT(iindex)%iAttValue )) then
      ncfile%rAddOffset(nc_z) = asfloat( pnc_var%pNC_ATT(iindex)%iAttValue(0) )
    elseif (allocated( pnc_var%pNC_ATT(iindex)%rAttValue )) then
      ncfile%rAddOffset(nc_z) = pnc_var%pNC_ATT(iindex)%rAttValue(0)
    elseif (allocated( pnc_var%pNC_ATT(iindex)%dpAttValue )) then
      ncfile%rAddOffset(nc_z) = asfloat( pnc_var%pNC_ATT(iindex)%dpAttValue(0) )
    elseif (len_trim(pnc_var%pNC_ATT(iindex)%sAttValue(0)) > 0) then
      sbuf = trim(pnc_var%pNC_ATT(iindex)%sAttValue(0) )
      read(sbuf,*) ncfile%rAddOffset(nc_z)
    else
      call die("Error reading the 'add_offset' attribute from the netCDF file" &
        //dquote(ncfile%sFilename))
    endif
  endif
 
end subroutine nf_get_scale_and_offset
 
!----------------------------------------------------------------------
 
subroutine nf_get_variable_id_and_type( NCFILE, strict_asserts )
  type (T_NETCDF4_FILE), intent(inout)          :: NCFILE
  logical (c_bool), intent(in), optional   :: strict_asserts
 
  ! [ LOCALS ]
  integer (c_int)              :: iIndex
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
  logical (c_bool)             :: strict_asserts_l
 
  if ( present( strict_asserts ) ) then
    strict_asserts_l = strict_asserts
  else
    strict_asserts_l = true
  endif
 
  ncfile%iVarID = -9999
 
  do iindex=0, ncfile%iNumberOfVariables - 1
 
    pnc_var => ncfile%pNC_VAR(iindex)
 
    if ( pnc_var%sVariableName .strequal. ncfile%sVarName(nc_x)  ) then
     ncfile%iVarIndex(nc_x) = iindex
     ncfile%iVarID(nc_x) = pnc_var%iNC_VarID
     ncfile%iVarType(nc_x) = pnc_var%iNC_VarType
     ncfile%iVar_DimID(nc_x,:) = pnc_var%iNC_DimID
 
    elseif ( pnc_var%sVariableName .strequal. ncfile%sVarName(nc_y) ) then
     ncfile%iVarIndex(nc_y) = iindex
     ncfile%iVarID(nc_y) = pnc_var%iNC_VarID
     ncfile%iVarType(nc_y) = pnc_var%iNC_VarType
     ncfile%iVar_DimID(nc_y,:) = pnc_var%iNC_DimID
 
   elseif ( pnc_var%sVariableName .strequal. ncfile%sVarName(nc_z) ) then
     ncfile%iVarIndex(nc_z) = iindex
     ncfile%iVarID(nc_z) = pnc_var%iNC_VarID
     ncfile%iVarType(nc_z) = pnc_var%iNC_VarType
     ncfile%iVar_DimID(nc_z,:) = pnc_var%iNC_DimID
 
    elseif ( pnc_var%sVariableName .strequal. ncfile%sVarName(nc_time) ) then
     ncfile%iVarIndex(nc_time) = iindex
     ncfile%iVarID(nc_time) = pnc_var%iNC_VarID
     ncfile%iVarType(nc_time) = pnc_var%iNC_VarType
     ncfile%iVar_DimID(nc_time,:) = pnc_var%iNC_DimID
    endif
 
  enddo
 
  if ( strict_asserts_l ) then
 
    call assert(ncfile%iVarID(nc_x) >= 0, &
     "Unable to find the variable named "//dquote(ncfile%sVarName(nc_x) )//" in " &
     //"file "//dquote(ncfile%sFilename), __file__, __line__)
 
    call assert(ncfile%iVarID(nc_y) >= 0, &
     "Unable to find the variable named "//dquote(ncfile%sVarName(nc_y))//" in " &
     //"file "//dquote(ncfile%sFilename), __file__, __line__)
 
    call assert(ncfile%iVarID(nc_z) >= 0, &
     "Unable to find the variable named "//dquote(ncfile%sVarName(nc_z))//" in " &
     //"file "//dquote(ncfile%sFilename), __file__, __line__)
 
    if ( ncfile%iVarID(nc_time) < 0 )  &
     call warn("Unable to find the variable named "//dquote(ncfile%sVarName(nc_time))//" in " &
       //"file "//dquote(ncfile%sFilename) )
 
  endif
 
end subroutine nf_get_variable_id_and_type
 
!----------------------------------------------------------------------
 
function nf_return_index_double(rValues, rTargetValue, rOffsetValue)  result(iIndex)
 
  real (c_double), dimension(:) :: rvalues
  real (c_double) :: rtargetvalue
  real (c_double) :: roffsetvalue
  integer (c_int) :: iindex
 
  ! [ LOCALS ]
  integer (c_int) :: icount
  real (c_double) :: rdiff, rdiffmin
 
  iindex = -1
 
  ! attempting to account for the fact that coordinates are specified
  ! initially relative toi cell edges, but are stored internally within netCDF
  ! relative to cell centers
  if ( .not. (rtargetvalue >= (minval(rvalues) - roffsetvalue) )             &
         .and. (rtargetvalue <= (maxval(rvalues) + roffsetvalue) ) ) then
    call logs%write("rTargetValue (" &
    //trim(ascharacter(rtargetvalue))//") is not within the range " &
    //trim(ascharacter(minval(rvalues)))//" to "//trim(ascharacter(maxval(rvalues))), lecho=true )
 
    call assert(false, "INTERNAL PROGRAMMING ERROR", __file__, __line__)
  endif
 
  rdiffmin = 1.e+20
 
  do icount=lbound(rvalues,1), ubound(rvalues,1)
 
    rdiff = abs(rvalues(icount) - rtargetvalue)
 
    if ( rdiff < rdiffmin ) then
      iindex = icount
      rdiffmin =rdiff
    endif
 
  enddo
 
!  print *, trim(__FILE__), ": ", __LINE__
!  print *, "index: ", iIndex, "  value: ", rValues(iIndex),           &
!           "  target value: ", rTargetValue
 
end function nf_return_index_double
 
!----------------------------------------------------------------------
 
function netcdf_coord_to_col_row(NCFILE, rX, rY)  result(iColRow)
 
  type (t_netcdf4_file ) :: ncfile
  real (c_double) :: rx
  real (c_double) :: ry
  integer (c_size_t), dimension(2) :: icolrow
 
  ! [ LOCALS ]
  integer (c_int) :: icolnum, irownum
  real (c_double) :: x_offset
  real (c_double) :: y_offset
 
  x_offset = ncfile%rGridCellSizeX / 2.0_c_double + ncfile%rCoordinateTolerance
  y_offset = ncfile%rGridCellSizeY / 2.0_c_double + ncfile%rCoordinateTolerance
 
  call assert( allocated( ncfile%rX_Coords ), "Internal programming error -- attempt " &
  //"to access unallocated array rX_Coords.", __file__, __line__ )
 
  if (rx < (minval(ncfile%rX_Coords) - x_offset) ) &
    call die( "X coordinate value "//ascharacter(rx)//" is less than the minimum X coordinate " &
      //"value ("//ascharacter(minval(ncfile%rX_Coords)-x_offset)//") contained in the netCDF file " &
      //dquote(ncfile%sFilename), trim(__file__), __line__ )
 
  if (rx > (maxval(ncfile%rX_Coords) + x_offset) ) &
    call die( "X coordinate value "//ascharacter(rx)//" is greater than the maximum X coordinate " &
      //"value ("//ascharacter(maxval(ncfile%rX_Coords)+x_offset)//") contained in the netCDF file " &
      //dquote(ncfile%sFilename), trim(__file__), __line__  )
 
  if (ry < (minval(ncfile%rY_Coords) - y_offset) ) &
    call die( "Y coordinate value "//ascharacter(ry)//" is less than the minimum Y coordinate " &
      //"value ("//ascharacter(minval(ncfile%rY_Coords)-y_offset)//") contained in the netCDF file " &
      //dquote(ncfile%sFilename), trim(__file__), __line__  )
 
  if (ry > (maxval(ncfile%rY_Coords) + y_offset) ) &
    call die( "Y coordinate value "//ascharacter(ry)//" is greater than the maximum Y coordinate " &
      //"value ("//ascharacter(maxval(ncfile%rY_Coords)+y_offset)//") contained in the netCDF file " &
      //dquote(ncfile%sFilename), trim(__file__), __line__  )
 
  icolnum = nf_return_index_double(ncfile%rX_Coords, rx, x_offset)
  irownum = nf_return_index_double(ncfile%rY_Coords, ry, y_offset)
 
  icolrow(column) = int(icolnum, c_size_t)
  icolrow(row) = int(irownum, c_size_t)
 
end function netcdf_coord_to_col_row
 
!----------------------------------------------------------------------
 
function nf_get_varid(NCFILE, sVariableName)  result(iNC_VarID)
 
  type (t_netcdf4_file ) :: ncfile
  character (len=*) :: svariablename
  integer (c_int) :: inc_varid
 
  integer (c_int) :: iindex
  type (t_netcdf_variable), pointer :: pnc_var
 
  inc_varid = -9999
 
  do iindex=0, ncfile%iNumberOfVariables - 1
 
    pnc_var => ncfile%pNC_VAR(iindex)
 
    if(trim(svariablename) .eq. trim(pnc_var%sVariableName) ) then
 
      inc_varid = iindex
      exit
 
    endif
 
  enddo
 
end function nf_get_varid
 
!----------------------------------------------------------------------
 
subroutine nf_create(NCFILE, sFilename, iLU)
 
  type (T_NETCDF4_FILE ) :: NCFILE
  character (len=*) :: sFilename
  integer (c_int), optional :: iLU
 
  logical :: old_invalid, old_div0, old_overflow
 
  ncfile%sFilename = trim(sfilename)
  ncfile%iFileFormat = nc_format_netcdf4
 
  ! attempting to keep code from choking on NaN comparisons during 
  ! file initialization
  call ieee_get_halting_mode(ieee_invalid,        old_invalid)
  call ieee_get_halting_mode(ieee_divide_by_zero, old_div0)
  call ieee_get_halting_mode(ieee_overflow,       old_overflow)
  
  ! Mask traps before calling into NetCDF/HDF5:
  call ieee_set_halting_mode(ieee_invalid,        .false.)
  call ieee_set_halting_mode(ieee_divide_by_zero, .false.)
  call ieee_set_halting_mode(ieee_overflow,       .false.)
 
  call nf_trap(nc_create(path=trim(sfilename)//c_null_char, &
                         cmode=nc_netcdf4, &
                         ncidp=ncfile%iNCID), &
                         __file__, __line__)
 
  ! Restore previous behavior:
  call ieee_set_halting_mode(ieee_invalid,        old_invalid)
  call ieee_set_halting_mode(ieee_divide_by_zero, old_div0)
  call ieee_set_halting_mode(ieee_overflow,       old_overflow)
 
!
! had read somewhere that the interface:
! character (c_char)   :: varname(*)
!
! would properly pass fortran string to c; however,
! this does not appear ot be the case.
!
 
!   call nf_trap(nc_create(path=trim(sFilename), &
!                  cmode=NC_NETCDF4, &
!                  ncidp=NCFILE%iNCID), &
!                  __FILE__, __LINE__)
 
 
  ! NCFILE%sFilename = trim(sFilename)
  ! NCFILE%iFileFormat = NC_FORMAT_NETCDF4
 
  if (present(ilu) ) then
    call logs%write("Created netCDF file for output. Filename: " &
      //dquote(ncfile%sFilename)//"; NCID="//trim(ascharacter(ncfile%iNCID) ) )
  endif
 
end subroutine nf_create
 
!----------------------------------------------------------------------
 
subroutine nf_define_deflate(NCFILE, iVarID, iShuffle, iDeflate, iDeflate_level)
 
  type (T_NETCDF4_FILE ) :: NCFILE
  integer (c_int) :: iVarID
  integer (c_int) :: iShuffle
  integer (c_int) :: iDeflate
  integer (c_int) :: iDeflate_level
 
  call nf_trap(nc_def_var_deflate(ncid=ncfile%iNCID, &
          varid=ivarid, &
          shuffle=ishuffle, &
          deflate=ideflate, &
          deflate_level=ideflate_level), &
          __file__, __line__)
 
end subroutine nf_define_deflate
 
!----------------------------------------------------------------------
 
subroutine nf_enddef(NCFILE)
 
  type (T_NETCDF4_FILE ) :: NCFILE
 
  call nf_trap(nc_enddef(ncid=ncfile%iNCID), &
       __file__, __line__)
 
end subroutine nf_enddef
 
!----------------------------------------------------------------------
 
subroutine nf_redef(NCFILE)
 
  type (T_NETCDF4_FILE ) :: NCFILE
 
  call nf_trap(nc_redef(ncid=ncfile%iNCID), &
       __file__, __line__)
 
end subroutine nf_redef
 
!----------------------------------------------------------------------
 
subroutine nf_define_dimension(NCFILE, sDimensionName, iDimensionSize)
 
  type (T_NETCDF4_FILE ) :: NCFILE
  character (len=*) :: sDimensionName
  integer (c_int) :: iDimensionSize
  integer (c_int) :: iDimID
 
  integer (c_size_t) :: iDimSize
 
  idimsize = int(idimensionsize, c_size_t)
 
  call nf_trap(nc_def_dim(ncid=ncfile%iNCID, &
                          name=trim(sdimensionname)//c_null_char, &
                          lenv=idimsize, &
                          dimidp=idimid), &
                          __file__, __line__)
 
end subroutine nf_define_dimension
 
!----------------------------------------------------------------------
 
subroutine nf_delete_attribute(NCFILE, sVariableName, sAttributeName)
 
  type (T_NETCDF4_FILE ) :: NCFILE
  character (len=*)      :: sVariableName
  character (len=*)      :: sAttributeName
 
  integer (c_int) :: iVarID
 
  ivarid = nf_get_varid(ncfile, svariablename//c_null_char)
 
  call nf_trap(nc_del_att(ncid=ncfile%iNCID,                       &
                          varid=ivarid,                            &
                          name=trim(sattributename)//c_null_char), &
                          __file__, __line__)
 
end subroutine nf_delete_attribute
 
!----------------------------------------------------------------------
 
subroutine nf_define_dimensions( NCFILE )
 
  type (T_NETCDF4_FILE) :: NCFILE
 
  ! [ LOCALS ]
  integer (c_int) :: iIndex
  type (T_NETCDF_DIMENSION), pointer :: pNC_DIM
 
  do iindex = 0, ncfile%iNumberOfDimensions-1
 
    pnc_dim => ncfile%pNC_DIM(iindex)
 
    call nf_trap(nc_def_dim(ncid=ncfile%iNCID, &
      name=trim(pnc_dim%sDimensionName)//c_null_char, &
      lenv=pnc_dim%iNC_DimSize, &
      dimidp=pnc_dim%iNC_DimID), &
      __file__, __line__ )
 
  enddo
 
end subroutine nf_define_dimensions
 
!----------------------------------------------------------------------
 
subroutine nf_set_standard_dimensions(NCFILE, iNX, iNY, write_time_bounds )
 
  type (T_NETCDF4_FILE )           :: NCFILE
  integer (c_int)             :: iNX
  integer (c_int)             :: iNY
  logical (c_bool), optional  :: write_time_bounds
 
  ! [ LOCALS ]
  integer (c_int)   :: iStat
  logical (c_bool)  :: write_time_bounds_l
 
  istat = 0
 
  ncfile%iNumberOfDimensions = 3
 
  if ( present( write_time_bounds ) ) then
 
    write_time_bounds_l = write_time_bounds
 
  else
 
    write_time_bounds_l = false
 
  endif
 
  if ( write_time_bounds_l ) ncfile%iNumberOfDimensions = 4
 
  if (associated(ncfile%pNC_DIM) ) deallocate(ncfile%pNC_DIM, stat=istat)
  call assert(istat == 0, "Could not deallocate memory for NC_DIM member in NC_FILE defined type", &
    __file__, __line__)
 
  allocate(ncfile%pNC_DIM( 0 : ncfile%iNumberOfDimensions-1), stat=istat )
  call assert(istat == 0, "Could not allocate memory for NC_DIM member in NC_FILE defined type", &
    __file__, __line__)
 
  !> define the time dimension;
  ncfile%pNC_DIM(nc_time)%sDimensionName = "time"
  ncfile%pNC_DIM(nc_time)%iNC_DimSize = nc_unlimited
 
  !> define the y dimension;
  ncfile%pNC_DIM(nc_y)%sDimensionName = "y"
  ncfile%pNC_DIM(nc_y)%iNC_DimSize = int(iny, c_size_t)
 
  !> define the x dimension;
  ncfile%pNC_DIM(nc_x)%sDimensionName = "x"
  ncfile%pNC_DIM(nc_x)%iNC_DimSize = int(inx, c_size_t)
 
  if ( write_time_bounds_l ) then
    !> define the auxiliary dimension;
    ncfile%pNC_DIM(nc_aux)%sDimensionName = "nv"
    ncfile%pNC_DIM(nc_aux)%iNC_DimSize = 2
  endif
 
end subroutine nf_set_standard_dimensions
 
!----------------------------------------------------------------------
 
subroutine nf_set_standard_variables(NCFILE, sVarName_z, lLatLon, write_time_bounds )
 
  type (T_NETCDF4_FILE )            :: NCFILE
  character (len=*)                 :: sVarName_z
  logical (c_bool), optional   :: lLatLon
  logical (c_bool), optional   :: write_time_bounds
 
  ! [ LOCALS ]
  integer (c_int) :: iStat
  logical (c_bool) :: lLatLon_l
  logical (c_bool) :: write_time_bounds_l
 
  if (present( llatlon) ) then
    llatlon_l = llatlon
  else
    llatlon_l = false
  endif
 
  if (present(write_time_bounds) ) then
    write_time_bounds_l = write_time_bounds
  else
    write_time_bounds_l = false
  endif
 
  istat = 0
 
  ncfile%iNumberOfVariables = 5
 
  if ( llatlon_l ) ncfile%iNumberOfVariables = 7
  if ( write_time_bounds_l )  ncfile%iNumberOfVariables =                       &
                               ncfile%iNumberOfVariables + 1
 
  ! reset the ID for TIME BNDS to the last varid
  nc_time_bnds = ncfile%iNumberOfVariables - 1
 
  if (associated(ncfile%pNC_VAR) ) deallocate(ncfile%pNC_VAR, stat=istat)
  call assert(istat == 0, "Could not deallocate memory for NC_VAR member in NC_FILE defined type", &
    __file__, __line__)
 
  allocate(ncfile%pNC_VAR( 0 : ncfile%iNumberOfVariables-1), stat=istat )
  call assert(istat == 0, "Could not allocate memory for NC_VAR member in NC_FILE defined type", &
    __file__, __line__)
 
  ncfile%pNC_VAR(nc_time)%sVariableName = "time"
  ncfile%pNC_VAR(nc_time)%iNC_VarType = nc_float
  ncfile%pNC_VAR(nc_time)%iNumberOfDimensions = 1
  ncfile%pNC_VAR(nc_time)%iNC_DimID(0) = ncfile%pNC_DIM(nc_time)%iNC_DimID
 
  ncfile%pNC_VAR(nc_y)%sVariableName = "y"
  ncfile%pNC_VAR(nc_y)%iNC_VarType = nc_double
  ncfile%pNC_VAR(nc_y)%iNumberOfDimensions = 1
  ncfile%pNC_VAR(nc_y)%iNC_DimID = ncfile%pNC_DIM(nc_y)%iNC_DimID
 
  ncfile%pNC_VAR(nc_x)%sVariableName = "x"
  ncfile%pNC_VAR(nc_x)%iNC_VarType = nc_double
  ncfile%pNC_VAR(nc_x)%iNumberOfDimensions = 1
  ncfile%pNC_VAR(nc_x)%iNC_DimID = ncfile%pNC_DIM(nc_x)%iNC_DimID
 
  ncfile%pNC_VAR(nc_crs)%sVariableName = "crs"
  ncfile%pNC_VAR(nc_crs)%iNC_VarType = nc_int
  ncfile%pNC_VAR(nc_crs)%iNumberOfDimensions = 0
 
  ncfile%pNC_VAR(nc_z)%sVariableName = trim(svarname_z)
  ncfile%pNC_VAR(nc_z)%iNC_VarType = nc_float
  ncfile%pNC_VAR(nc_z)%iNumberOfDimensions = 3
  ncfile%pNC_VAR(nc_z)%iNC_DimID = [ncfile%pNC_DIM(nc_time)%iNC_DimID, &
                                 ncfile%pNC_DIM(nc_y)%iNC_DimID, &
                                 ncfile%pNC_DIM(nc_x)%iNC_DimID,0]
 
  ncfile%sVarName(nc_z) = trim(svarname_z)
 
  if ( llatlon_l ) then
 
    ncfile%pNC_VAR(nc_lat)%sVariableName = "lat"
    ncfile%pNC_VAR(nc_lat)%iNC_VarType = nc_double
    ncfile%pNC_VAR(nc_lat)%iNumberOfDimensions = 2
    ncfile%pNC_VAR(nc_lat)%iNC_DimID = [ncfile%pNC_DIM(nc_y)%iNC_DimID, &
                                        ncfile%pNC_DIM(nc_x)%iNC_DimID,0,0]
 
    ncfile%pNC_VAR(nc_lon)%sVariableName = "lon"
    ncfile%pNC_VAR(nc_lon)%iNC_VarType = nc_double
    ncfile%pNC_VAR(nc_lon)%iNumberOfDimensions = 2
    ncfile%pNC_VAR(nc_lon)%iNC_DimID = [ncfile%pNC_DIM(nc_y)%iNC_DimID, &
                                        ncfile%pNC_DIM(nc_x)%iNC_DimID,0,0]
  endif
 
  if ( write_time_bounds_l ) then
 
    ncfile%pNC_VAR(nc_time_bnds)%sVariableName = "time_bnds"
    ncfile%pNC_VAR(nc_time_bnds)%iNC_VarType = nc_double
    ncfile%pNC_VAR(nc_time_bnds)%iNumberOfDimensions = 2
    ncfile%pNC_VAR(nc_time_bnds)%iNC_DimID = [ncfile%pNC_DIM(nc_time)%iNC_DimID,       &
                                              ncfile%pNC_DIM(nc_aux)%iNC_DimID,0,0]
 
  endif
 
end subroutine nf_set_standard_variables
 
!----------------------------------------------------------------------
 
subroutine nf_set_global_attributes(NCFILE, sDataType, executable_name, &
                                    history_list, sSourceFile )
 
  type (T_NETCDF4_FILE ) :: NCFILE
  character (len=*), intent(in) :: sDataType
  character (len=*), intent(in), optional    :: executable_name
  type (FSTRING_LIST_T), intent(in), pointer, optional :: history_list
  character (len=*), intent(in), optional    :: sSourceFile
 
  ! [ LOCALS ]
  integer (c_int)           :: iStat
  type (DATETIME_T)              :: DT
  character (len=20)             :: sDateTime
  character (len=:), allocatable :: executable_name_l
  type (FSTRING_LIST_T), pointer           :: history_list_l
  integer (c_int)           :: indx, jndx
  integer (c_int)  :: records
 
  if (present( executable_name) ) then
    executable_name_l = trim( executable_name )
  else
    executable_name_l = "SWB"
  endif
 
  if ( present( history_list ) ) then
    history_list_l => history_list
    ncfile%iNumberOfAttributes = 3 + history_list_l%count
  else
    allocate( history_list_l )
    call history_list_l%append(trim(sdatetime)//": Soil-Water-Balance model run started.")
    ncfile%iNumberOfAttributes = 4
  endif
 
  call dt%systime()
  sdatetime = dt%prettydatetime()
 
  allocate( ncfile%pNC_ATT(0:ncfile%iNumberOfAttributes-1), stat=istat)
  call assert(istat == 0, "Could not allocate memory for NC_ATT member of NC_FILE", &
    __file__, __line__)
 
  block
 
    if (present(ssourcefile) ) then
 
      ncfile%pNC_ATT(0)%sAttributeName = "source"
      allocate(ncfile%pNC_ATT(0)%sAttValue(0:0))
      ncfile%pNC_ATT(0)%sAttValue(0) = trim(sdatatype)//" data from file "//dquote(ssourcefile)
      ncfile%pNC_ATT(0)%iNC_AttType = nc_char
      ncfile%pNC_ATT(0)%iNC_AttSize = 1_c_size_t
 
    else
 
      ncfile%pNC_ATT(0)%sAttributeName = "source"
      allocate(ncfile%pNC_ATT(0)%sAttValue(0:0))
      ncfile%pNC_ATT(0)%sAttValue(0) = trim(sdatatype)//" output from " &
        //executable_name_l//" run "   &
        //"started on "//trim(sdatetime)//"."
      ncfile%pNC_ATT(0)%iNC_AttType = nc_char
      ncfile%pNC_ATT(0)%iNC_AttSize = 1_c_size_t
 
    endif
 
    ncfile%pNC_ATT(1)%sAttributeName = "executable_version"
    allocate(ncfile%pNC_ATT(1)%sAttValue(0:0))
    ncfile%pNC_ATT(1)%sAttValue(0) = "version "//trim(swb_version)          &
      //", Git branch: "//trim(git_branch_string)//", Git commit hash string: "   &
      //trim(git_commit_hash_string)//", compiled on: "//trim(compile_date)       &
      //" "//trim(compile_time)//"."
    ncfile%pNC_ATT(1)%iNC_AttType = nc_char
    ncfile%pNC_ATT(1)%iNC_AttSize = 1_c_size_t
 
    ncfile%pNC_ATT(2)%sAttributeName = "Conventions"
    allocate(ncfile%pNC_ATT(2)%sAttValue(0:0))
    ncfile%pNC_ATT(2)%sAttValue(0) = "CF-1.6"
    ncfile%pNC_ATT(2)%iNC_AttType = nc_char
    ncfile%pNC_ATT(2)%iNC_AttSize = 1_c_size_t
 
    ! special case: history may have meny records
    records = history_list_l%count
 
    do indx=1, records
 
      jndx = indx + 2
      ncfile%pNC_ATT( jndx )%sAttributeName = "history"
      allocate(ncfile%pNC_ATT(jndx)%sAttValue( 0:0 ) )
      ncfile%pNC_ATT(jndx)%iNC_AttType = nc_char
      ncfile%pNC_ATT(jndx)%iNC_AttSize = int( records, c_size_t )
      ncfile%pNC_ATT(jndx)%sAttValue( 0 ) = trim(history_list_l%get( indx ))//c_null_char
 
    enddo
 
  end block
 
end subroutine nf_set_global_attributes
 
!----------------------------------------------------------------------
 
subroutine nf_set_standard_attributes(NCFILE, sOriginText, PROJ4_string,    &
                                      lLatLon, fValidMin, fValidMax,        &
                                      write_time_bounds )
 
  type (T_NETCDF4_FILE )             :: NCFILE
  character (len=*)                  :: sOriginText
  character (len=*), optional        :: PROJ4_string
  logical (c_bool), optional    :: lLatLon
  real (c_float), optional      :: fValidMin
  real (c_float), optional      :: fValidMax
  logical (c_bool), optional    :: write_time_bounds
 
  ! [ LOCALS ]
  integer (c_int)                             :: iStat
  integer (c_int)                             :: iNumAttributes
  type (T_NETCDF_ATTRIBUTE), dimension(:), pointer :: pNC_ATT
  logical (c_bool)                            :: lLatLon_l
  logical (c_bool)                            :: write_time_bounds_l
  type (FSTRING_LIST_T)                       :: attribute_name_list
  type (FSTRING_LIST_T)                       :: attribute_value_list
  character (len=:), allocatable              :: tempstring
  character (len=:), allocatable              :: value_string
  character (len=:), allocatable              :: value_string1
  character (len=:), allocatable              :: value_string2
  integer (c_int)                             :: indx
 
  if (present( llatlon ) ) then
    llatlon_l = llatlon
  else
    llatlon_l = false
  endif
 
  if ( present( write_time_bounds ) ) then
    write_time_bounds_l = write_time_bounds
  else
    write_time_bounds_l = false
  endif
 
  if (present( proj4_string ) ) then
    call create_attributes_from_proj4_string( proj4_string, attribute_name_list,   &
                                              attribute_value_list )
 
    ! Define attributes for the coordinate reference system (CRS)
    inumattributes = attribute_name_list%count + 1
 
    allocate( ncfile%pNC_VAR(nc_crs)%pNC_ATT(0:inumattributes-1), stat=istat)
    call assert(istat == 0, "Could not allocate memory for NC_ATT member in NC_VAR struct of NC_FILE", &
      __file__, __line__)
    ncfile%pNC_VAR(nc_crs)%iNumberOfAttributes = inumattributes
 
    pnc_att => ncfile%pNC_VAR(nc_crs)%pNC_ATT
 
    ! the following block of code parses the PROJ4 string and rips standard
    ! CF attributes from the string
    do indx=0, inumattributes-2
 
      tempstring = attribute_name_list%get( indx + 1 )
      pnc_att(indx)%sAttributeName = tempstring
 
      select case ( tempstring )
 
        case ( "units" )
 
          value_string = attribute_value_list%get( indx + 1 )
          ncfile%sVarUnits(nc_x) = value_string
          ncfile%sVarUnits(nc_y) = value_string
 
          allocate(pnc_att(indx)%sAttValue(0:0))
          pnc_att(indx)%sAttValue(0) = value_string
          pnc_att(indx)%iNC_AttType = nc_char
          pnc_att(indx)%iNC_AttSize = 1_c_size_t
 
        case ( "datum", "spheroid", "grid_mapping_name" )
 
          allocate(pnc_att(indx)%sAttValue(0:0))
          pnc_att(indx)%sAttValue(0) = attribute_value_list%get( indx + 1 )
          pnc_att(indx)%iNC_AttType = nc_char
          pnc_att(indx)%iNC_AttSize = 1_c_size_t
 
        case ( "standard_parallel" )
 
          allocate(pnc_att(indx)%dpAttValue(0:1))
 
          value_string = attribute_value_list%get( indx + 1 )
          value_string1 = left(value_string, substring=",")
          value_string2 = right(value_string, substring=",")
          
          call assert(len_trim(value_string1) > 0, "standard_parallel requires valid values for '+lat_1' and '+lat_2'.", &
            shints="Are '+lat_1' or '+lat_2' missing or out of order in the control file PROJ string?")
 
          pnc_att(indx)%dpAttValue(0) = asdouble( value_string1 )
          pnc_att(indx)%dpAttValue(1) = asdouble( value_string2 )
          pnc_att(indx)%iNC_AttType = nc_double
          pnc_att(indx)%iNC_AttSize = 2_c_size_t
 
        case ( "UTM_zone" )
 
          allocate(pnc_att(indx)%iAttValue(0:0))
          pnc_att(indx)%iAttValue(0) = asint( attribute_value_list%get( indx + 1 ) )
          pnc_att(indx)%iNC_AttType = nc_int
          pnc_att(indx)%iNC_AttSize = 1_c_size_t
 
        case default
 
          allocate(pnc_att(indx)%dpAttValue(0:0))
          pnc_att(indx)%dpAttValue(0) = asdouble( attribute_value_list%get( indx + 1 ) )
          pnc_att(indx)%iNC_AttType = nc_double
          pnc_att(indx)%iNC_AttSize = 1_c_size_t
 
      end select
 
    enddo
 
    ! last, store the actual PROJ4 string
    pnc_att(indx)%sAttributeName = "proj4_string"
    allocate(pnc_att(indx)%sAttValue(0:0))
    pnc_att(indx)%sAttValue(0) = proj4_string
    pnc_att(indx)%iNC_AttType = nc_char
    pnc_att(indx)%iNC_AttSize = 1_c_size_t
 
    call attribute_name_list%clear()
    call attribute_value_list%clear()
 
  endif
 
  !! define attributes associated with TIME variable
  inumattributes = 3
  if ( write_time_bounds_l ) inumattributes=4
  allocate( ncfile%pNC_VAR(nc_time)%pNC_ATT(0:inumattributes-1), stat=istat)
  call assert(istat == 0, "Could not allocate memory for NC_ATT member in NC_VAR struct of NC_FILE", &
    __file__, __line__)
  ncfile%pNC_VAR(nc_time)%iNumberOfAttributes = inumattributes
 
  pnc_att => ncfile%pNC_VAR(nc_time)%pNC_ATT
 
  pnc_att(0)%sAttributeName = "units"
  allocate(pnc_att(0)%sAttValue(0:0))
  pnc_att(0)%sAttValue(0) = "days since "//trim(sorigintext)//" 00:00:00"
  pnc_att(0)%iNC_AttType = nc_char
  pnc_att(0)%iNC_AttSize = 1_c_size_t
 
  pnc_att(1)%sAttributeName = "calendar"
  allocate(pnc_att(1)%sAttValue(0:0))
  pnc_att(1)%sAttValue(0) = "standard"
  pnc_att(1)%iNC_AttType = nc_char
  pnc_att(1)%iNC_AttSize = 1_c_size_t
 
  pnc_att(2)%sAttributeName = "long_name"
  allocate(pnc_att(2)%sAttValue(0:0))
  pnc_att(2)%sAttValue(0) = "time"
  pnc_att(2)%iNC_AttType = nc_char
  pnc_att(2)%iNC_AttSize = 1_c_size_t
 
  if ( write_time_bounds_l ) then
    pnc_att(3)%sAttributeName = "bounds"
    allocate(pnc_att(3)%sAttValue(0:0))
    pnc_att(3)%sAttValue(0) = "time_bounds"
    pnc_att(3)%iNC_AttType = nc_char
    pnc_att(3)%iNC_AttSize = 1_c_size_t
  endif
 
  if (present( fvalidmin ) .and. present( fvalidmax) ) then
 
    inumattributes = 7
    allocate( ncfile%pNC_VAR(nc_z)%pNC_ATT(0:inumattributes-1), stat=istat)
    call assert(istat == 0, "Could not allocate memory for NC_ATT member in NC_VAR struct of NC_FILE", &
      __file__, __line__)
    ncfile%pNC_VAR(nc_z)%iNumberOfAttributes = inumattributes
 
    pnc_att => ncfile%pNC_VAR(nc_z)%pNC_ATT
 
    pnc_att(0)%sAttributeName = "units"
    allocate(pnc_att(0)%sAttValue(0:0))
    pnc_att(0)%sAttValue(0) = ncfile%sVarUnits(nc_z)
    pnc_att(0)%iNC_AttType = nc_char
    pnc_att(0)%iNC_AttSize = 1_c_size_t
 
    pnc_att(1)%sAttributeName = "valid_min"
    allocate(pnc_att(1)%rAttValue(0:0))
    pnc_att(1)%rAttValue(0) = fvalidmin
    pnc_att(1)%iNC_AttType = nc_float
    pnc_att(1)%iNC_AttSize = 1_c_size_t
 
    pnc_att(2)%sAttributeName = "valid_max"
    allocate(pnc_att(2)%rAttValue(0:0))
    pnc_att(2)%rAttValue(0) = fvalidmax
    pnc_att(2)%iNC_AttType = nc_float
    pnc_att(2)%iNC_AttSize = 1_c_size_t
 
    pnc_att(3)%sAttributeName = "valid_range"
    allocate(pnc_att(3)%rAttValue(0:1))
    pnc_att(3)%rAttValue(0) = fvalidmin
    pnc_att(3)%rAttValue(1) = fvalidmax
    pnc_att(3)%iNC_AttType = nc_float
    pnc_att(3)%iNC_AttSize = 2_c_size_t
 
    pnc_att(4)%sAttributeName = "_FillValue"
    allocate(pnc_att(4)%rAttValue(0:0))
    pnc_att(4)%rAttValue(0) = nc_fill_float
    pnc_att(4)%iNC_AttType = nc_float
    pnc_att(4)%iNC_AttSize = 1_c_size_t
 
    pnc_att(5)%sAttributeName = "coordinates"
    allocate(pnc_att(5)%sAttValue(0:0))
!    pNC_ATT(5)%sAttValue(0) = "lat lon"
    pnc_att(5)%sAttValue(0) = "crs"
    pnc_att(5)%iNC_AttType = nc_char
    pnc_att(5)%iNC_AttSize = 1_c_size_t
 
    pnc_att(6)%sAttributeName = "grid_mapping"
    allocate(pnc_att(6)%sAttValue(0:0))
    pnc_att(6)%sAttValue(0) = "crs"
    pnc_att(6)%iNC_AttType = nc_char
    pnc_att(6)%iNC_AttSize = 1_c_size_t
 
  else
 
    inumattributes = 3
    allocate( ncfile%pNC_VAR(nc_z)%pNC_ATT(0:inumattributes-1), stat=istat)
    call assert(istat == 0, "Could not allocate memory for NC_ATT member in NC_VAR struct of NC_FILE", &
      __file__, __line__)
    ncfile%pNC_VAR(nc_z)%iNumberOfAttributes = inumattributes
 
    pnc_att => ncfile%pNC_VAR(nc_z)%pNC_ATT
 
    pnc_att(0)%sAttributeName = "units"
    allocate(pnc_att(0)%sAttValue(0:0))
    pnc_att(0)%sAttValue(0) = ncfile%sVarUnits(nc_z)
    pnc_att(0)%iNC_AttType = nc_char
    pnc_att(0)%iNC_AttSize = 1_c_size_t
 
    pnc_att(1)%sAttributeName = "coordinates"
    allocate(pnc_att(1)%sAttValue(0:0))
    pnc_att(1)%sAttValue(0) = "lat lon"
    pnc_att(1)%iNC_AttType = nc_char
    pnc_att(1)%iNC_AttSize = 1_c_size_t
 
    pnc_att(2)%sAttributeName = "grid_mapping"
    allocate(pnc_att(2)%sAttValue(0:0))
    pnc_att(2)%sAttValue(0) = "crs"
    pnc_att(2)%iNC_AttType = nc_char
    pnc_att(2)%iNC_AttSize = 1_c_size_t
 
  endif
 
  inumattributes = 3
  allocate( ncfile%pNC_VAR(nc_y)%pNC_ATT(0:inumattributes-1), stat=istat)
  call assert(istat == 0, "Could not allocate memory for NC_ATT member in NC_VAR struct of NC_FILE", &
    __file__, __line__)
  ncfile%pNC_VAR(nc_y)%iNumberOfAttributes = inumattributes
 
  pnc_att => ncfile%pNC_VAR(nc_y)%pNC_ATT
 
  block
 
    pnc_att(0)%sAttributeName = "units"
    allocate(pnc_att(0)%sAttValue(0:0))
    pnc_att(0)%sAttValue(0) = ncfile%sVarUnits(nc_y)
    pnc_att(0)%iNC_AttType = nc_char
    pnc_att(0)%iNC_AttSize = 1_c_size_t
 
    pnc_att(1)%sAttributeName = "long_name"
    allocate(pnc_att(1)%sAttValue(0:0))
    pnc_att(1)%sAttValue(0) = "y coordinate of projection"
    pnc_att(1)%iNC_AttType = nc_char
    pnc_att(1)%iNC_AttSize = 1_c_size_t
 
    pnc_att(2)%sAttributeName = "standard_name"
    allocate(pnc_att(2)%sAttValue(0:0))
    pnc_att(2)%sAttValue(0) = "projection_y_coordinate"
    pnc_att(2)%iNC_AttType = nc_char
    pnc_att(2)%iNC_AttSize = 1_c_size_t
 
 
  end block
 
  inumattributes = 3
  allocate( ncfile%pNC_VAR(nc_x)%pNC_ATT(0:inumattributes-1), stat=istat)
  call assert(istat == 0, "Could not allocate memory for NC_ATT member in NC_VAR struct of NC_FILE", &
    __file__, __line__)
  ncfile%pNC_VAR(nc_x)%iNumberOfAttributes = inumattributes
 
  pnc_att => ncfile%pNC_VAR(nc_x)%pNC_ATT
 
  block
 
    pnc_att(0)%sAttributeName = "units"
    allocate(pnc_att(0)%sAttValue(0:0))
    pnc_att(0)%sAttValue(0) = ncfile%sVarUnits(nc_x)
    pnc_att(0)%iNC_AttType = nc_char
    pnc_att(0)%iNC_AttSize = 1_c_size_t
 
    pnc_att(1)%sAttributeName = "long_name"
    allocate(pnc_att(1)%sAttValue(0:0))
    pnc_att(1)%sAttValue(0) = "x coordinate of projection"
    pnc_att(1)%iNC_AttType = nc_char
    pnc_att(1)%iNC_AttSize = 1_c_size_t
 
    pnc_att(2)%sAttributeName = "standard_name"
    allocate(pnc_att(2)%sAttValue(0:0))
    pnc_att(2)%sAttValue(0) = "projection_x_coordinate"
    pnc_att(2)%iNC_AttType = nc_char
    pnc_att(2)%iNC_AttSize = 1_c_size_t
 
 
  end block
 
  if ( llatlon_l ) then
 
    inumattributes = 3
    allocate( ncfile%pNC_VAR(nc_lat)%pNC_ATT(0:inumattributes-1), stat=istat)
    call assert(istat == 0, "Could not allocate memory for NC_ATT member in NC_VAR struct of NC_FILE", &
      __file__, __line__)
    ncfile%pNC_VAR(nc_lat)%iNumberOfAttributes = inumattributes
 
    pnc_att => ncfile%pNC_VAR(nc_lat)%pNC_ATT
 
    block
 
      pnc_att(0)%sAttributeName = "units"
      allocate(pnc_att(0)%sAttValue(0:0))
      pnc_att(0)%sAttValue(0) = "degrees_north"
      pnc_att(0)%iNC_AttType = nc_char
      pnc_att(0)%iNC_AttSize = 1_c_size_t
 
      pnc_att(1)%sAttributeName = "long_name"
      allocate(pnc_att(1)%sAttValue(0:0))
      pnc_att(1)%sAttValue(0) = "latitude"
      pnc_att(1)%iNC_AttType = nc_char
      pnc_att(1)%iNC_AttSize = 1_c_size_t
 
      pnc_att(2)%sAttributeName = "standard_name"
      allocate(pnc_att(2)%sAttValue(0:0))
      pnc_att(2)%sAttValue(0) = "latitude"
      pnc_att(2)%iNC_AttType = nc_char
      pnc_att(2)%iNC_AttSize = 1_c_size_t
 
 
    end block
 
 
    inumattributes = 3
    allocate( ncfile%pNC_VAR(nc_lon)%pNC_ATT(0:inumattributes-1), stat=istat)
    call assert(istat == 0, "Could not allocate memory for NC_ATT member in NC_VAR struct of NC_FILE", &
      __file__, __line__)
    ncfile%pNC_VAR(nc_lon)%iNumberOfAttributes = inumattributes
 
    pnc_att => ncfile%pNC_VAR(nc_lon)%pNC_ATT
 
    block
 
      pnc_att(0)%sAttributeName = "units"
      allocate(pnc_att(0)%sAttValue(0:0))
      pnc_att(0)%sAttValue(0) = "degrees_east"
      pnc_att(0)%iNC_AttType = nc_char
      pnc_att(0)%iNC_AttSize = 1_c_size_t
 
      pnc_att(1)%sAttributeName = "long_name"
      allocate(pnc_att(1)%sAttValue(0:0))
      pnc_att(1)%sAttValue(0) = "longitude"
      pnc_att(1)%iNC_AttType = nc_char
      pnc_att(1)%iNC_AttSize = 1_c_size_t
 
      pnc_att(2)%sAttributeName = "standard_name"
      allocate(pnc_att(2)%sAttValue(0:0))
      pnc_att(2)%sAttValue(0) = "longitude"
      pnc_att(2)%iNC_AttType = nc_char
      pnc_att(2)%iNC_AttSize = 1_c_size_t
 
 
    end block
 
  endif
 
end subroutine nf_set_standard_attributes
 
!----------------------------------------------------------------------
 
subroutine nf_put_x_and_y(NCFILE, dpX, dpY)
 
  type (T_NETCDF4_FILE) :: NCFILE
  real (c_double), dimension(:) :: dpX
  real (c_double), dimension(:) :: dpY
 
  ! [ LOCALS ]
  integer (c_size_t) :: iLength
 
  ilength = int(size(dpx, 1), c_size_t)
 
  call netcdf_put_variable_vector(ncfile=ncfile, &
                   ivarid=ncfile%pNC_VAR(nc_x)%iNC_VarID, &
                   istart=[0_c_size_t], &
                   icount=[ilength], &
                   istride=[1_c_size_t], &
                   dpvalues=dpx)
 
  ilength = int(size(dpy, 1), c_size_t)
 
  call netcdf_put_variable_vector(ncfile=ncfile, &
                   ivarid=ncfile%pNC_VAR(nc_y)%iNC_VarID, &
                   istart=[0_c_size_t], &
                   icount=[ilength], &
                   istride=[1_c_size_t], &
                   dpvalues=dpy)
 
end subroutine nf_put_x_and_y
 
!----------------------------------------------------------------------
 
subroutine nf_put_lat_and_lon(NCFILE, dpLat, dpLon)
 
  type (T_NETCDF4_FILE) :: NCFILE
  real (c_double), dimension(:,:) :: dpLat
  real (c_double), dimension(:,:) :: dpLon
 
  ! [ LOCALS ]
  integer (c_size_t) :: iNX, iNY
 
  inx = int( size(dplat, 1), c_size_t)
  iny = int( size(dplat, 2), c_size_t)
 
  call netcdf_put_variable_array(ncfile=ncfile, &
                   ivarid=ncfile%pNC_VAR(nc_lat)%iNC_VarID, &
                   istart=[0_c_size_t, 0_c_size_t], &
                   icount=[ iny, inx ], &
                   istride=[1_c_size_t,1_c_size_t], &
                   dpvalues=dplat)
 
 
  call netcdf_put_variable_array(ncfile=ncfile, &
                   ivarid=ncfile%pNC_VAR(nc_lon)%iNC_VarID, &
                   istart=[0_c_size_t, 0_c_size_t], &
                   icount=[ iny, inx ],&
                   istride=[1_c_size_t,1_c_size_t], &
                   dpvalues=dplon)
 
end subroutine nf_put_lat_and_lon
 
!----------------------------------------------------------------------
 
function nf_define_variable(NCFILE, sVariableName, iVariableType, &
   iNumberOfDimensions, iDimIDs)    result(iVarID)
 
  type (t_netcdf4_file ) :: ncfile
  character (len=*) :: svariablename
  integer (c_int) :: ivariabletype
  integer (c_int) :: inumberofdimensions
  integer (c_int), dimension(:) :: idimids
  integer (c_int) :: ivarid
 
  call nf_trap( nc_def_var(ncid=ncfile%iNCID,&
                           name=trim(fortran_to_c_string(svariablename)), &
                           xtype=ivariabletype, &
                           ndims=inumberofdimensions, &
                           dimidsp=idimids, &
                           varidp=ivarid), &
                           __file__, __line__)
 
end function nf_define_variable
 
!----------------------------------------------------------------------
 
!! before this function is called, the values associated with NCFILE must be defined.
 
subroutine nf_define_variables( NCFILE )
 
  type (T_NETCDF4_FILE) :: NCFILE
 
  ! [ LOCALS ]
  integer (c_int) :: iIndex
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
 
  !! note: the default number of variables for a simple archive file is 4:
  !! 0) time, 1) Y, 2) X, 3) variable of interest
 
  do iindex = 0, ncfile%iNumberOfVariables-1
 
    pnc_var => ncfile%pNC_VAR(iindex)
 
    call nf_trap( nc_def_var(ncid=ncfile%iNCID,&
                             name=trim(fortran_to_c_string(pnc_var%sVariableName)), &
                             xtype=pnc_var%iNC_VarType, &
                             ndims=pnc_var%iNumberOfDimensions, &
                             dimidsp=pnc_var%iNC_DimID, &
                             varidp=pnc_var%iNC_VarID), &
                             __file__, __line__)
 
  enddo
 
end subroutine nf_define_variables
 
!----------------------------------------------------------------------
 
subroutine netcdf_rewrite_attribute(NCFILE, sVariableName, sAttributeName,       &
  sAttributeValue, iAttributeValue, rAttributeValue, dpAttributeValue)
 
  type (t_netcdf4_file ) :: ncfile
  character (len=*) :: svariablename
  character (len=*) :: sattributename
  character (len=*), dimension(:), optional :: sattributevalue
  integer (c_int), dimension(:), optional :: iattributevalue
  real (c_float), dimension(:), optional :: rattributevalue
  real (c_double), dimension(:), optional :: dpattributevalue
 
  integer (c_int) :: ivarid
 
  ! put netCDF file into define mode again before attempting to redefine the attribute
  call nf_redef(ncfile)
 
  ivarid = nf_get_varid(ncfile, trim(svariablename))
  
  if (present(sattributevalue))     &
    call nf_put_attribute( ncfile, ivarid, trim(sattributename)//c_null_char, sattributevalue )
  
  if (present(iattributevalue))     &
    call nf_put_attribute( ncfile, ivarid, trim(sattributename)//c_null_char, iattributevalue=iattributevalue )
  
  if (present(rattributevalue))     &
    call nf_put_attribute( ncfile, ivarid, trim(sattributename)//c_null_char, rattributevalue=rattributevalue )
  
  if (present(dpattributevalue))    &
    call nf_put_attribute( ncfile, ivarid, trim(sattributename)//c_null_char, dpattributevalue=dpattributevalue )
 
end subroutine netcdf_rewrite_attribute
 
!----------------------------------------------------------------------
 
subroutine nf_put_attribute(NCFILE, iVarID, sAttributeName, &
  sAttributeValue, iAttributeValue, rAttributeValue, dpAttributeValue)
 
  type (T_NETCDF4_FILE ) :: NCFILE
  integer (c_int) :: iVarID
  character (len=*) :: sAttributeName
  character (len=*), dimension(:), optional :: sAttributeValue
  integer (c_int), dimension(:), optional :: iAttributeValue
  real (c_float), dimension(:), optional :: rAttributeValue
  real (c_double), dimension(:), optional :: dpAttributeValue
 
  ! [ LOCALS ]
  integer (c_size_t) :: iNumberOfAttributes
 
  if (present(sattributevalue) ) then
 
    inumberofattributes = int(size( sattributevalue, 1), c_size_t)
    inumberofattributes = int(len_trim(sattributevalue(1)), c_size_t)
 
    call nf_trap( nc_put_att_text(ncid=ncfile%iNCID, &
                    varid=ivarid, &
                    name=trim(sattributename), &
                    nlen=inumberofattributes, &
                    tp=trim(sattributevalue(1))), &
                    __file__, __line__)
 
  elseif (present(iattributevalue) ) then
 
    inumberofattributes = int(size( iattributevalue, 1), c_size_t)
 
    call nf_trap( nc_put_att_int(ncid=ncfile%iNCID, &
                    varid=ivarid, &
                    name=trim(sattributename), &
                    xtype=nc_int, &
                    nlen=inumberofattributes, &
                     ip=iattributevalue), &
                     __file__, __line__)
 
  elseif (present(rattributevalue) ) then
 
    inumberofattributes = int(size( rattributevalue, 1), c_size_t)
 
    call nf_trap( nc_put_att_float(ncid=ncfile%iNCID, &
                    varid=ivarid, &
                    name=trim(sattributename), &
                    xtype=nc_float, &
                    nlen=inumberofattributes, &
                    fp=rattributevalue), &
                    __file__, __line__)
 
  elseif (present(dpattributevalue) ) then
 
    inumberofattributes = int(size( dpattributevalue, 1), c_size_t)
 
    call nf_trap( nc_put_att_double(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       name=trim(sattributename), &
                       xtype=nc_double, &
                       nlen=inumberofattributes, &
                       dp=dpattributevalue), &
                       __file__, __line__)
 
  endif
 
 
end subroutine nf_put_attribute
 
!----------------------------------------------------------------------
 
subroutine nf_put_attributes(NCFILE)
 
  type (T_NETCDF4_FILE ) :: NCFILE
 
  ! [ LOCALS ]
  type (T_NETCDF_VARIABLE), pointer :: pNC_VAR
  type (T_NETCDF_ATTRIBUTE), pointer :: pNC_ATT
  integer (c_int) :: iIndex
  integer (c_int) :: iIndex2
  integer (c_int) :: indx
 
  ! loop over variables
  do iindex = 0, ncfile%iNumberOfVariables-1
 
    pnc_var => ncfile%pNC_VAR(iindex)
 
    ! for each variable, loop over the associated attributes
    do iindex2 = 0, pnc_var%iNumberOfAttributes-1
 
      pnc_att => pnc_var%pNC_ATT(iindex2)
 
        select case (pnc_att%iNC_AttType)
 
          case (nc_double)
 
            if (.not. allocated(pnc_att%dpAttValue) ) &
              call die("INTERNAL PROGRAMMING ERROR - attempt to use unallocated variable; " &
              //"attribute name: "//dquote(pnc_att%sAttributeName), &
              __file__, __line__)
 
            call nf_put_attribute(ncfile=ncfile, &
                ivarid=pnc_var%iNC_VarID, &
                sattributename=trim(pnc_att%sAttributeName)//c_null_char, &
                dpattributevalue=pnc_att%dpAttValue)
 
          case (nc_int)
 
            if (.not. allocated(pnc_att%iAttValue) ) &
              call die("INTERNAL PROGRAMMING ERROR - attempt to use unallocated variable; " &
              //"attribute name: "//dquote(pnc_att%sAttributeName), &
              __file__, __line__)
 
            call nf_put_attribute(ncfile=ncfile, &
                ivarid=pnc_var%iNC_VarID, &
                sattributename=trim(pnc_att%sAttributeName)//c_null_char, &
                iattributevalue=pnc_att%iAttValue)
 
          case (nc_float)
 
            if (.not. allocated(pnc_att%rAttValue) ) &
              call die("INTERNAL PROGRAMMING ERROR - attempt to use unallocated variable; " &
              //"attribute name: "//dquote(pnc_att%sAttributeName), &
              __file__, __line__)
 
            call nf_put_attribute(ncfile=ncfile, &
                ivarid=pnc_var%iNC_VarID, &
                sattributename=trim(pnc_att%sAttributeName)//c_null_char, &
                rattributevalue=pnc_att%rAttValue)
 
          case (nc_char)
 
            if (.not. allocated(pnc_att%sAttValue) ) &
              call die("INTERNAL PROGRAMMING ERROR - attempt to use unallocated variable; " &
              //"attribute name: "//dquote(pnc_att%sAttributeName), &
              __file__, __line__)
 
            do indx=0,ubound(pnc_att%sAttValue, 1)
 
              call nf_put_attribute(ncfile=ncfile, &
                  ivarid=pnc_var%iNC_VarID, &
                  sattributename=trim(pnc_att%sAttributeName)//c_null_char, &
                  sattributevalue=[trim(pnc_att%sAttValue(indx))//c_null_char])
 
            enddo
 
        end select
 
    enddo
 
  enddo
 
  ! now loop over global attributes
  do iindex2 = 0, ncfile%iNumberOfAttributes-1
 
    pnc_att => ncfile%pNC_ATT(iindex2)
 
    select case (pnc_att%iNC_AttType)
 
      case (nc_double)
 
        if (.not. allocated(pnc_att%sAttValue) ) &
          call die("INTERNAL PROGRAMMING ERROR - attempt to use unallocated variable; " &
          //"attribute name: "//dquote(pnc_att%sAttributeName), &
          __file__, __line__)
 
        call nf_put_attribute(ncfile=ncfile, &
            ivarid=nc_global, &
            sattributename=trim(pnc_att%sAttributeName)//c_null_char, &
            dpattributevalue=pnc_att%dpAttValue)
 
      case (nc_int)
 
        if (.not. allocated(pnc_att%sAttValue) ) &
          call die("INTERNAL PROGRAMMING ERROR - attempt to use unallocated variable; " &
          //"attribute name: "//dquote(pnc_att%sAttributeName), &
          __file__, __line__)
 
        call nf_put_attribute(ncfile=ncfile, &
            ivarid=nc_global, &
            sattributename=trim(pnc_att%sAttributeName)//c_null_char, &
            iattributevalue=pnc_att%iAttValue)
 
      case (nc_float)
 
        if (.not. allocated(pnc_att%sAttValue) ) &
          call die("INTERNAL PROGRAMMING ERROR - attempt to use unallocated variable; " &
          //"attribute name: "//dquote(pnc_att%sAttributeName), &
          __file__, __line__)
 
        call nf_put_attribute(ncfile=ncfile, &
            ivarid=nc_global, &
            sattributename=trim(pnc_att%sAttributeName)//c_null_char, &
            rattributevalue=pnc_att%rAttValue)
 
      case (nc_char)
 
        if (.not. allocated(pnc_att%sAttValue) ) &
          call die("INTERNAL PROGRAMMING ERROR - attempt to use unallocated variable; " &
          //"attribute name: "//dquote(pnc_att%sAttributeName), &
          __file__, __line__)
 
        do indx=0, ubound(pnc_att%sAttValue, 1)
 
          call nf_put_attribute(ncfile=ncfile,                               &
              ivarid=nc_global,                                              &
              sattributename=trim(pnc_att%sAttributeName)//c_null_char,      &
              sattributevalue=[trim(pnc_att%sAttValue(indx))//c_null_char])
 
        enddo
 
    end select
 
  enddo
 
end subroutine nf_put_attributes
 
!----------------------------------------------------------------------
 
subroutine netcdf_put_variable_array(NCFILE, iVarID, iStart, iCount, iStride, &
   iValues, i2Values, rValues, dpValues)
 
  type (t_netcdf4_file ) :: ncfile
  integer (c_int) :: ivarid
  integer (c_size_t), dimension(:) :: istart
  integer (c_size_t), dimension(:) :: icount
  integer (c_size_t), dimension(:) :: istride
  integer (c_int), dimension(:,:), optional :: ivalues
  integer (c_short), dimension(:,:), optional :: i2values
  real (c_float), dimension(:,:), optional :: rvalues
  real (c_double), dimension(:,:), optional :: dpvalues
 
  if (present(ivalues) ) then
 
    call nf_trap(nc_put_vars_int(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=ivalues), &
                       __file__, __line__)
 
  elseif (present(i2values) ) then
 
    call nf_trap(nc_put_vars_short(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=i2values), &
                       __file__, __line__)
 
  elseif (present(rvalues) ) then
 
   call nf_trap(nc_put_vars_float(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=rvalues), &
                       __file__, __line__)
 
  elseif (present(dpvalues) ) then
 
 
 
    call nf_trap(nc_put_vars_double(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=dpvalues), &
                       __file__, __line__)
 
  endif
 
end subroutine netcdf_put_variable_array
 
subroutine netcdf_put_packed_variable_array(NCFILE, iVarID, iStart, iCount, iStride, &
   lMask, iValues, iField, i2Values, i2Field, rValues, rField, dpValues, dpField)
 
  type (t_netcdf4_file ) :: ncfile
  integer (c_int) :: ivarid
  integer (c_size_t), dimension(:)            :: istart
  integer (c_size_t), dimension(:)            :: icount
  integer (c_size_t), dimension(:)         :: istride
  logical (c_bool), dimension(:,:)            :: lmask
  integer (c_int), dimension(:), optional     :: ivalues
  integer (c_int), dimension(:,:), optional   :: ifield
  integer (c_short), dimension(:), optional   :: i2values
  integer (c_short), dimension(:,:), optional :: i2field
  real (c_float), dimension(:), optional      :: rvalues
  real (c_float), dimension(:,:), optional    :: rfield
  real (c_double), dimension(:), optional     :: dpvalues
  real (c_double), dimension(:,:), optional   :: dpfield
 
  if (present(ivalues) ) then
 
    call nf_trap(nc_put_vars_int(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=unpack(ivalues, lmask, ifield)), &
                       __file__, __line__)
 
  elseif (present(i2values) ) then
 
    call nf_trap(nc_put_vars_short(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=unpack(i2values, lmask, i2field)), &
                       __file__, __line__)
 
  elseif (present(rvalues) ) then
 
   call nf_trap(nc_put_vars_float(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=unpack(rvalues, lmask, rfield)), &
                       __file__, __line__)
 
  elseif (present(dpvalues) ) then
 
    call nf_trap(nc_put_vars_double(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=unpack(dpvalues, lmask, dpfield)), &
                       __file__, __line__)
 
  endif
 
end subroutine netcdf_put_packed_variable_array
 
 
subroutine netcdf_put_variable_vector(NCFILE, iVarID, iStart, iCount, iStride, &
   iValues, i2Values, rValues, dpValues)
 
  type (t_netcdf4_file ) :: ncfile
  integer (c_int) :: ivarid
  integer (c_size_t), dimension(:) :: istart
  integer (c_size_t), dimension(:) :: icount
  integer (c_size_t), dimension(:) :: istride
  integer (c_int), dimension(:), optional :: ivalues
  integer (c_short), dimension(:), optional :: i2values
  real (c_float), dimension(:), optional :: rvalues
  real (c_double), dimension(:), optional :: dpvalues
 
  if (present(ivalues) ) then
 
    call nf_trap(nc_put_vars_int(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=ivalues), &
                       __file__, __line__)
 
  elseif (present(i2values) ) then
 
    call nf_trap(nc_put_vars_short(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=i2values), &
                       __file__, __line__)
 
  elseif (present(rvalues) ) then
 
    call nf_trap(nc_put_vars_float(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=rvalues), &
                       __file__, __line__)
 
  elseif (present(dpvalues) ) then
 
    call nf_trap(nc_put_vars_double(ncid=ncfile%iNCID, &
                       varid=ivarid, &
                       startp=istart, &
                       countp=icount, &
                       stridep=istride, &
                       vars=dpvalues), &
                       __file__, __line__, ncfile%sFilename )
 
  endif
 
end subroutine netcdf_put_variable_vector
 
 
end module netcdf4_support