DisNCStructured.f90

Go to the documentation of this file.

!> @brief This module contains the DisNCStructuredModule
!!
!! This module defines a STRUCTURED (non-ugrid) netcdf
!! export type for DIS models. It is dependent on netcdf
!! libraries.
!!
!<
module disncstructuredmodule
 
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: linelength, lenbigline, lencomponentname, &
                             lenmempath, dnodata, dzero, dhalf
  use simvariablesmodule, only: errmsg, warnmsg
  use simmodule, only: store_error, store_warning, store_error_filename
  use memorymanagermodule, only: mem_setptr
  use inputdefinitionmodule, only: inputparamdefinitiontype
  use characterstringmodule, only: characterstringtype
  use ncmodelexportmodule, only: ncbasemodelexporttype, export_varname, &
                                 export_longname, wkt_to_cf_gridmapping, &
                                 wrap_rotated_crs, &
                                 wkt_transverse_mercator_params, &
                                 wkt_lambert_conformal_conic_params, &
                                 wkt_albers_conical_equal_area_params
  use dismodule, only: distype
  use netcdfcommonmodule, only: nf_verify
  use netcdf
 
  implicit none
  private
  public :: disncstructuredtype
 
  type :: structuredncdimidtype
    integer(I4B) :: x !< number of columns
    integer(I4B) :: y !< number of rows
    integer(I4B) :: layer !< number of layers
    integer(I4B) :: time !< number of steps
    integer(I4B) :: bnd !< number in boundary
  contains
  end type structuredncdimidtype
 
  type :: structuredncvaridtype
    integer(I4B) :: x !< x coordinate variable
    integer(I4B) :: y !< y coordinate variable
    integer(I4B) :: layer !< layer coordinate variable (discrete layer index, not a position)
    integer(I4B) :: elevation !< z auxiliary coordinate variable (cell center elevation, a real vertical position)
    integer(I4B) :: time !< time coordinate variable
    integer(I4B) :: dependent !< dependent variable
    integer(I4B) :: x_bnds !< x boundaries 2D array
    integer(I4B) :: y_bnds !< y boundaries 2D array
    integer(I4B) :: latitude !< latitude 2D array
    integer(I4B) :: longitude !< longitude 2D array
    integer(I4B) :: export !< in scope export
  contains
  end type structuredncvaridtype
 
  type, extends(ncbasemodelexporttype) :: disncstructuredtype
    type(structuredncdimidtype) :: dim_ids !< structured dimension ids type
    type(structuredncvaridtype) :: var_ids !< structured variable ids type
    type(distype), pointer :: dis => null() !< pointer to model dis package
    integer(I4B) :: nlay !< number of layers
    real(dp), dimension(:), pointer, contiguous :: latitude => null() !< lat input array pointer
    real(dp), dimension(:), pointer, contiguous :: longitude => null() !< lon input array pointer
    integer(I4B), pointer :: chunk_z !< chunking parameter for z dimension
    integer(I4B), pointer :: chunk_y !< chunking parameter for y dimension
    integer(I4B), pointer :: chunk_x !< chunking parameter for x dimension
    integer(I4B), dimension(:), allocatable :: layers !< layers array
    logical(LGP) :: latlon !< are lat and lon arrays to be written to netcdf file
  contains
    procedure :: init => dis_export_init
    procedure :: destroy => dis_export_destroy
    procedure :: df
    procedure :: df_export
    procedure :: step
    procedure :: export_input_array
    procedure :: export_df
    procedure :: create_timeseries
    procedure :: export_input_arrays
    procedure :: package_step
    procedure :: add_pkg_data
    procedure :: add_global_att
    procedure :: define_dim
    procedure :: define_dependent
    procedure :: define_gridmap
    procedure :: define_geocoords
    procedure :: add_proj_data
    procedure :: add_grid_data
  end type disncstructuredtype
 
  interface nc_export_array
    module procedure nc_export_int1d, nc_export_int2d, &
      nc_export_int3d, nc_export_dbl1d, &
      nc_export_dbl2d, nc_export_dbl3d
  end interface nc_export_array
 
contains
 
  !> @brief netcdf export dis init
  !<
  subroutine dis_export_init(this, modelname, modeltype, modelfname, nc_fname, &
                             disenum, nctype, iout)
    use memorymanagermodule, only: get_isize
    use memorymanagerextmodule, only: mem_set_value
    class(disncstructuredtype), intent(inout) :: this
    character(len=*), intent(in) :: modelname
    character(len=*), intent(in) :: modeltype
    character(len=*), intent(in) :: modelfname
    character(len=*), intent(in) :: nc_fname
    integer(I4B), intent(in) :: disenum
    integer(I4B), intent(in) :: nctype
    integer(I4B), intent(in) :: iout
    integer(I4B) :: k, latsz, lonsz
    logical(LGP) :: found
 
    ! set nlay
    this%nlay = this%dis%nlay
 
    ! allocate
    allocate (this%chunk_z)
    allocate (this%chunk_y)
    allocate (this%chunk_x)
    allocate (this%layers(this%nlay))
 
    ! initialize
    this%chunk_z = -1
    this%chunk_y = -1
    this%chunk_x = -1
    do k = 1, this%nlay
      this%layers(k) = k
    end do
 
    this%latlon = .false.
 
    ! initialize base class
    call this%NCModelExportType%init(modelname, modeltype, modelfname, nc_fname, &
                                     disenum, nctype, iout)
 
    ! update values from input context
    if (this%ncf_mempath /= '') then
      call mem_set_value(this%chunk_z, 'CHUNK_Z', this%ncf_mempath, found)
      call mem_set_value(this%chunk_y, 'CHUNK_Y', this%ncf_mempath, found)
      call mem_set_value(this%chunk_x, 'CHUNK_X', this%ncf_mempath, found)
 
      if (this%chunk_time > 0 .and. this%chunk_z > 0 .and. &
          this%chunk_y > 0 .and. this%chunk_x > 0) then
        this%chunking_active = .true.
      else if (this%chunk_time > 0 .or. this%chunk_z > 0 .or. &
               this%chunk_y > 0 .or. this%chunk_x > 0) then
        this%chunk_time = -1
        this%chunk_z = -1
        this%chunk_y = -1
        this%chunk_x = -1
        write (warnmsg, '(a)') 'Ignoring user provided NetCDF chunking &
          &parameters. Define chunk_time, chunk_x, chunk_y and chunk_z input &
          &parameters to see an effect in file "'//trim(nc_fname)//'".'
        call store_warning(warnmsg)
      end if
 
      call get_isize('LATITUDE', this%ncf_mempath, latsz)
      call get_isize('LONGITUDE', this%ncf_mempath, lonsz)
 
      if (latsz > 0 .and. lonsz > 0) then
        this%latlon = .true.
        if (this%wkt /= '' .or. this%crs_wkt /= '') then
          write (warnmsg, '(a)') 'Ignoring user provided NetCDF wkt/crs_wkt &
         &parameter(s) as longitude and latitude arrays have been provided. &
            &Applies to file "'//trim(nc_fname)//'".'
          call store_warning(warnmsg)
          this%wkt = ''
          this%crs_wkt = ''
          this%gridmap_name = ''
        end if
        call mem_setptr(this%latitude, 'LATITUDE', this%ncf_mempath)
        call mem_setptr(this%longitude, 'LONGITUDE', this%ncf_mempath)
        if (this%dis%angrot /= dzero) then
          write (warnmsg, '(a)') 'Structured rotated grid with &
            &user-provided latitude/longitude: MF6 cannot verify these &
            &values correctly account for grid rotation (angrot). They &
            &must be computed from each cell''s true rotated position, &
            &not from grid-local x/y. Applies to file "'// &
            trim(nc_fname)//'".'
          call store_warning(warnmsg)
        end if
      end if
 
      if (this%gridmap_name /= '') then
        if (this%dis%angrot /= dzero) then
          write (warnmsg, '(a)') 'CRS parameter set with structured rotated &
            &grid. The x/y coordinate variables have grid-local, not &
            &real-world, values. A CRS-aware consumer must apply the full &
            &crs_wkt (which encodes rotation when a valid WKT2 PROJCRS was &
            &supplied) via proper reprojection to resolve the true &
            &position; deriving real longitude/latitude directly from x/y &
            &and the flat/numeric CF grid_mapping parameters alone will be &
            &incorrect, since those parameters never encode grid rotation. &
            &Applies to file "'//trim(nc_fname)//'".'
          call store_warning(warnmsg)
        end if
      end if
    end if
 
    if (this%dis%lenuni == 1) then
      this%lenunits = 'ft'
    else if (this%dis%lenuni == 3) then
      this%lenunits = 'cm'
    else
      this%lenunits = 'm'
    end if
 
    ! create the netcdf file
    call nf_verify(nf90_create(this%nc_fname, &
                               ior(nf90_clobber, nf90_netcdf4), this%ncid), &
                   this%nc_fname)
  end subroutine dis_export_init
 
  !> @brief netcdf export dis destroy
  !<
  subroutine dis_export_destroy(this)
    class(disncstructuredtype), intent(inout) :: this
    call nf_verify(nf90_close(this%ncid), this%nc_fname)
    deallocate (this%chunk_z)
    deallocate (this%chunk_y)
    deallocate (this%chunk_x)
    deallocate (this%layers)
    nullify (this%chunk_z)
    nullify (this%chunk_y)
    nullify (this%chunk_x)
    ! destroy base class
    call this%NCModelExportType%destroy()
  end subroutine dis_export_destroy
 
  !> @brief netcdf export define
  !<
  subroutine df(this)
    use constantsmodule, only: mvalidate
    use simvariablesmodule, only: isim_mode
    class(disncstructuredtype), intent(inout) :: this
    ! put root group file scope attributes
    call this%add_global_att()
    ! define root group dimensions and coordinate variables
    call this%define_dim()
    ! define grid projection variables
    call this%define_geocoords()
    if (isim_mode == mvalidate) then
      ! define period input arrays
      call this%df_export()
    else
      ! define the dependent variable
      call this%define_dependent()
    end if
    ! exit define mode
    call nf_verify(nf90_enddef(this%ncid), this%nc_fname)
    ! add data locations
    call this%add_grid_data()
    ! add projection data
    call this%add_proj_data()
    if (isim_mode == mvalidate) then
      ! define and set package input griddata
      call this%add_pkg_data()
    end if
    ! define and set gridmap variable
    call this%define_gridmap()
    ! synchronize file
    call nf_verify(nf90_sync(this%ncid), this%nc_fname)
  end subroutine df
 
  !> @brief define timeseries input variables
  !<
  subroutine df_export(this)
    use ncmodelexportmodule, only: exportpackagetype
    class(disncstructuredtype), intent(inout) :: this
    class(exportpackagetype), pointer :: export_pkg
    integer(I4B) :: idx
    do idx = 1, this%pkglist%Count()
      export_pkg => this%get(idx)
      call this%export_df(export_pkg)
    end do
  end subroutine df_export
 
  !> @brief netcdf export step
  !<
  subroutine step(this)
    use constantsmodule, only: dhnoflo
    use tdismodule, only: totim
    class(disncstructuredtype), intent(inout) :: this
    real(DP), dimension(:), pointer, contiguous :: dbl1d
    integer(I4B) :: n, istp
 
    ! set global step index
    istp = this%istp()
 
    if (size(this%dis%nodeuser) < &
        size(this%dis%nodereduced)) then
      allocate (dbl1d(size(this%dis%nodereduced)))
      dbl1d = dhnoflo
      do n = 1, size(this%dis%nodereduced)
        if (this%dis%nodereduced(n) > 0) then
          dbl1d(n) = this%x(this%dis%nodereduced(n))
        end if
      end do
      ! write step data to dependent variable
      call nf_verify(nf90_put_var(this%ncid, &
                                  this%var_ids%dependent, dbl1d, &
                                  start=(/1, 1, 1, istp/), &
                                  count=(/this%dis%ncol, &
                                          this%dis%nrow, &
                                          this%dis%nlay, 1/)), &
                     this%nc_fname)
      deallocate (dbl1d)
    else
      ! write step data to dependent variable
      call nf_verify(nf90_put_var(this%ncid, &
                                  this%var_ids%dependent, this%x, &
                                  start=(/1, 1, 1, istp/), &
                                  count=(/this%dis%ncol, &
                                          this%dis%nrow, &
                                          this%dis%nlay, 1/)), &
                     this%nc_fname)
    end if
 
    ! write to time coordinate variable
    call nf_verify(nf90_put_var(this%ncid, this%var_ids%time, &
                                totim, start=(/istp/)), &
                   this%nc_fname)
 
    ! synchronize file
    call nf_verify(nf90_sync(this%ncid), this%nc_fname)
  end subroutine step
 
  !> @brief netcdf export an input array
  !<
  subroutine export_input_array(this, pkgtype, pkgname, mempath, idt)
    class(disncstructuredtype), intent(inout) :: this
    character(len=*), intent(in) :: pkgtype
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: mempath
    type(inputparamdefinitiontype), pointer, intent(in) :: idt
    integer(I4B), dimension(:), pointer, contiguous :: int1d
    integer(I4B), dimension(:, :), pointer, contiguous :: int2d
    integer(I4B), dimension(:, :, :), pointer, contiguous :: int3d
    real(DP), dimension(:), pointer, contiguous :: dbl1d
    real(DP), dimension(:, :), pointer, contiguous :: dbl2d
    real(DP), dimension(:, :, :), pointer, contiguous :: dbl3d
    character(len=LINELENGTH) :: nc_tag
    integer(I4B) :: iper, iaux
 
    ! initialize
    iper = 0
    iaux = 0
 
    ! set variable name and input attribute string
    nc_tag = this%input_attribute(pkgname, idt)
 
    select case (idt%datatype)
    case ('INTEGER1D')
      call mem_setptr(int1d, idt%mf6varname, mempath)
      call nc_export_array(int1d, this%ncid, this%dim_ids, this%var_ids, &
                           this%dis, idt, mempath, nc_tag, pkgname, &
                           this%gridmap_name, this%latlon, this%deflate, &
                           this%shuffle, this%chunk_z, this%chunk_y, &
                           this%chunk_x, iper, this%nc_fname)
    case ('INTEGER2D')
      call mem_setptr(int2d, idt%mf6varname, mempath)
      call nc_export_array(int2d, this%ncid, this%dim_ids, this%var_ids, &
                           this%dis, idt, mempath, nc_tag, pkgname, &
                           this%gridmap_name, this%latlon, this%deflate, &
                           this%shuffle, this%chunk_z, this%chunk_y, &
                           this%chunk_x, this%nc_fname)
    case ('INTEGER3D')
      call mem_setptr(int3d, idt%mf6varname, mempath)
      call nc_export_array(int3d, this%ncid, this%dim_ids, this%var_ids, &
                           this%dis, idt, mempath, nc_tag, pkgname, &
                           this%gridmap_name, this%latlon, this%deflate, &
                           this%shuffle, this%chunk_z, this%chunk_y, &
                           this%chunk_x, this%nc_fname)
    case ('DOUBLE1D')
      call mem_setptr(dbl1d, idt%mf6varname, mempath)
      call nc_export_array(dbl1d, this%ncid, this%dim_ids, this%var_ids, &
                           this%dis, idt, mempath, nc_tag, pkgname, &
                           this%gridmap_name, this%latlon, this%deflate, &
                           this%shuffle, this%chunk_z, this%chunk_y, &
                           this%chunk_x, iper, iaux, this%nc_fname)
    case ('DOUBLE2D')
      call mem_setptr(dbl2d, idt%mf6varname, mempath)
      call nc_export_array(dbl2d, this%ncid, this%dim_ids, this%var_ids, &
                           this%dis, idt, mempath, nc_tag, pkgname, &
                           this%gridmap_name, this%latlon, this%deflate, &
                           this%shuffle, this%chunk_z, this%chunk_y, &
                           this%chunk_x, this%nc_fname)
    case ('DOUBLE3D')
      call mem_setptr(dbl3d, idt%mf6varname, mempath)
      call nc_export_array(dbl3d, this%ncid, this%dim_ids, this%var_ids, &
                           this%dis, idt, mempath, nc_tag, pkgname, &
                           this%gridmap_name, this%latlon, this%deflate, &
                           this%shuffle, this%chunk_z, this%chunk_y, &
                           this%chunk_x, this%nc_fname)
    case default
      ! no-op, no other datatypes exported
    end select
  end subroutine export_input_array
 
  !> @brief define export package
  !<
  subroutine export_df(this, export_pkg)
    use ncmodelexportmodule, only: exportpackagetype
    use definitionselectmodule, only: get_param_definition_type
    class(disncstructuredtype), intent(inout) :: this
    class(exportpackagetype), pointer, intent(in) :: export_pkg
    type(inputparamdefinitiontype), pointer :: idt
    integer(I4B) :: iparam, iaux
 
    ! export defined period input
    do iparam = 1, export_pkg%nparam
      ! initialize
      iaux = 0
      ! set input definition
      idt => &
        get_param_definition_type(export_pkg%mf6_input%param_dfns, &
                                  export_pkg%mf6_input%component_type, &
                                  export_pkg%mf6_input%subcomponent_type, &
                                  'PERIOD', export_pkg%param_names(iparam), '')
      select case (idt%shape)
      case ('NCPL', 'NODES')
        call this%create_timeseries(idt, iparam, iaux, export_pkg)
      case ('NAUX NCPL', 'NAUX NODES')
        do iaux = 1, export_pkg%naux
          call this%create_timeseries(idt, iparam, iaux, export_pkg)
        end do
      case default
      end select
    end do
  end subroutine export_df
 
  !> @brief create timeseries export variable
  !<
  subroutine create_timeseries(this, idt, iparam, iaux, export_pkg)
    use ncmodelexportmodule, only: exportpackagetype
    class(disncstructuredtype), intent(inout) :: this
    type(inputparamdefinitiontype), pointer, intent(in) :: idt
    integer(I4B), intent(in) :: iparam
    integer(I4B), intent(in) :: iaux
    class(exportpackagetype), pointer, intent(in) :: export_pkg
    character(len=LINELENGTH) :: varname, longname, nc_tag
    integer(I4B) :: varid
    logical(LGP) :: has_layer
 
    has_layer = .not. (idt%shape == 'NCPL' .or. idt%shape == 'NAUX NCPL')
 
    ! set variable input tag
    nc_tag = this%input_attribute(export_pkg%mf6_input%subcomponent_name, &
                                  idt)
 
    ! set names
    varname = export_varname(export_pkg%mf6_input%subcomponent_name, &
                             idt%tagname, export_pkg%mf6_input%mempath, &
                             iaux=iaux)
    longname = export_longname(idt%longname, &
                               export_pkg%mf6_input%subcomponent_name, &
                               idt%tagname, export_pkg%mf6_input%mempath, &
                               iaux=iaux, &
                               component_type=idt%component_type, &
                               subcomponent_type=idt%subcomponent_type)
 
    ! create the netcdf timeseries variable
    select case (idt%datatype)
    case ('DOUBLE1D', 'DOUBLE2D')
      if (idt%shape == 'NCPL' .or. &
          idt%shape == 'NAUX NCPL') then
        call nf_verify(nf90_def_var(this%ncid, varname, nf90_double, &
                                    (/this%dim_ids%x, &
                                      this%dim_ids%y, &
                                      this%dim_ids%time/), varid), &
                       this%nc_fname)
      else
        call nf_verify(nf90_def_var(this%ncid, varname, nf90_double, &
                                    (/this%dim_ids%x, &
                                      this%dim_ids%y, &
                                      this%dim_ids%layer, &
                                      this%dim_ids%time/), varid), &
                       this%nc_fname)
      end if
      call nf_verify(nf90_put_att(this%ncid, varid, &
                                  '_FillValue', (/dnodata/)), &
                     this%nc_fname)
    case ('INTEGER1D')
      if (idt%shape == 'NCPL' .or. &
          idt%shape == 'NAUX NCPL') then
        call nf_verify(nf90_def_var(this%ncid, varname, nf90_int, &
                                    (/this%dim_ids%x, &
                                      this%dim_ids%y, &
                                      this%dim_ids%time/), varid), &
                       this%nc_fname)
      else
        call nf_verify(nf90_def_var(this%ncid, varname, nf90_int, &
                                    (/this%dim_ids%x, &
                                      this%dim_ids%y, &
                                      this%dim_ids%layer, &
                                      this%dim_ids%time/), varid), &
                       this%nc_fname)
      end if
      call nf_verify(nf90_put_att(this%ncid, varid, &
                                  '_FillValue', (/nf90_fill_int/)), &
                     this%nc_fname)
    end select
 
    ! apply chunking parameters
    if (this%chunking_active) then
      call nf_verify(nf90_def_var_chunking(this%ncid, &
                                           varid, &
                                           nf90_chunked, &
                                           (/this%chunk_x, this%chunk_y, &
                                             this%chunk_z, this%chunk_time/)), &
                     this%nc_fname)
    end if
 
    ! deflate and shuffle
    call ncvar_deflate(this%ncid, varid, this%deflate, &
                       this%shuffle, this%nc_fname)
 
    ! variable attributes
    call nf_verify(nf90_put_att(this%ncid, varid, &
                                'units', this%lenunits), this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, varid, &
                                'long_name', longname), this%nc_fname)
 
    ! add grid mapping and mf6 attr
    call ncvar_gridmap(this%ncid, varid, this%gridmap_name, this%latlon, &
                       has_layer, this%nc_fname)
    call ncvar_mf6attr(this%ncid, varid, iaux, nc_tag, this%nc_fname)
 
    ! store variable id
    if (idt%tagname == 'AUX') then
      export_pkg%varids_aux(iaux, 1) = varid
    else
      export_pkg%varids_param(iparam, 1) = varid
    end if
  end subroutine create_timeseries
 
  !> @brief write package gridded input data
  !<
  subroutine export_input_arrays(this, pkgtype, pkgname, mempath, param_dfns)
    use memorymanagermodule, only: get_isize
    class(disncstructuredtype), intent(inout) :: this
    character(len=*), intent(in) :: pkgtype
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: mempath
    type(inputparamdefinitiontype), dimension(:), pointer, &
      intent(in) :: param_dfns
    type(inputparamdefinitiontype), pointer :: idt
    integer(I4B) :: iparam, isize
    do iparam = 1, size(param_dfns)
      ! assign param definition pointer
      idt => param_dfns(iparam)
      ! for now only griddata is exported
      if (idt%blockname == 'GRIDDATA') then
        ! check if variable is already allocated
        call get_isize(idt%mf6varname, mempath, isize)
        if (isize > 0) then
          call this%export_input_array(pkgtype, pkgname, mempath, idt)
        end if
      end if
    end do
  end subroutine export_input_arrays
 
  !> @brief netcdf export package dynamic input
  !<
  subroutine package_step(this, export_pkg)
    use tdismodule, only: totim, kper
    use definitionselectmodule, only: get_param_definition_type
    use ncmodelexportmodule, only: exportpackagetype
    class(disncstructuredtype), intent(inout) :: this
    class(exportpackagetype), pointer, intent(in) :: export_pkg
    type(inputparamdefinitiontype), pointer :: idt
    integer(I4B), dimension(:), pointer, contiguous :: int1d
    real(DP), dimension(:), pointer, contiguous :: dbl1d, nodes
    real(DP), dimension(:, :), pointer, contiguous :: dbl2d
    character(len=LINELENGTH) :: nc_tag
    integer(I4B) :: iaux, iparam, nvals, n
    integer(I4B), pointer :: nbound
 
    ! initialize
    iaux = 0
 
    ! export defined period input
    do iparam = 1, export_pkg%nparam
      if (export_pkg%iper /= kper) cycle
      ! check if variable was read this period
      if (export_pkg%param_reads(iparam)%invar < 1) cycle
 
      ! set input definition
      idt => &
        get_param_definition_type(export_pkg%mf6_input%param_dfns, &
                                  export_pkg%mf6_input%component_type, &
                                  export_pkg%mf6_input%subcomponent_type, &
                                  'PERIOD', export_pkg%param_names(iparam), '')
 
      ! set variable input tag
      nc_tag = this%input_attribute(export_pkg%mf6_input%subcomponent_name, &
                                    idt)
 
      ! export arrays
      select case (idt%datatype)
      case ('INTEGER1D')
        call mem_setptr(int1d, idt%mf6varname, export_pkg%mf6_input%mempath)
        this%var_ids%export = export_pkg%varids_param(iparam, 1)
        call nc_export_array(int1d, this%ncid, this%dim_ids, this%var_ids, &
                             this%dis, idt, export_pkg%mf6_input%mempath, &
                             nc_tag, export_pkg%mf6_input%subcomponent_name, &
                             this%gridmap_name, this%latlon, this%deflate, &
                             this%shuffle, this%chunk_z, this%chunk_y, &
                             this%chunk_x, kper, this%nc_fname)
      case ('DOUBLE1D')
        call mem_setptr(dbl1d, idt%mf6varname, export_pkg%mf6_input%mempath)
        this%var_ids%export = export_pkg%varids_param(iparam, 1)
        select case (idt%shape)
        case ('NCPL')
          call nc_export_array(dbl1d, this%ncid, this%dim_ids, this%var_ids, &
                               this%dis, idt, export_pkg%mf6_input%mempath, &
                               nc_tag, export_pkg%mf6_input%subcomponent_name, &
                               this%gridmap_name, this%latlon, this%deflate, &
                               this%shuffle, this%chunk_z, this%chunk_y, &
                               this%chunk_x, kper, iaux, this%nc_fname)
        case ('NODES')
          nvals = this%dis%nodesuser
          allocate (nodes(nvals))
          nodes = dnodata
          call mem_setptr(dbl1d, idt%mf6varname, export_pkg%mf6_input%mempath)
          call mem_setptr(int1d, 'NODEULIST', export_pkg%mf6_input%mempath)
          call mem_setptr(nbound, 'NBOUND', export_pkg%mf6_input%mempath)
          do n = 1, nbound
            nodes(int1d(n)) = dbl1d(n)
          end do
          call nc_export_array(nodes, this%ncid, this%dim_ids, this%var_ids, &
                               this%dis, idt, export_pkg%mf6_input%mempath, &
                               nc_tag, export_pkg%mf6_input%subcomponent_name, &
                               this%gridmap_name, this%latlon, this%deflate, &
                               this%shuffle, this%chunk_z, this%chunk_y, &
                               this%chunk_x, kper, iaux, this%nc_fname)
          deallocate (nodes)
        case default
        end select
      case ('DOUBLE2D')
        call mem_setptr(dbl2d, idt%mf6varname, export_pkg%mf6_input%mempath)
        select case (idt%shape)
        case ('NAUX NCPL')
          nvals = this%dis%nrow * this%dis%ncol
          allocate (nodes(nvals))
          do iaux = 1, size(dbl2d, dim=1) !naux
            this%var_ids%export = export_pkg%varids_aux(iaux, 1)
            do n = 1, nvals
              nodes(n) = dbl2d(iaux, n)
            end do
            call nc_export_array(nodes, this%ncid, this%dim_ids, this%var_ids, &
                                 this%dis, idt, export_pkg%mf6_input%mempath, &
                                 nc_tag, export_pkg%mf6_input%subcomponent_name, &
                                 this%gridmap_name, this%latlon, this%deflate, &
                                 this%shuffle, this%chunk_z, this%chunk_y, &
                                 this%chunk_x, kper, iaux, this%nc_fname)
          end do
          deallocate (nodes)
        case ('NAUX NODES')
          nvals = this%dis%nodesuser
          allocate (nodes(nvals))
          call mem_setptr(int1d, 'NODEULIST', export_pkg%mf6_input%mempath)
          call mem_setptr(nbound, 'NBOUND', export_pkg%mf6_input%mempath)
          do iaux = 1, size(dbl2d, dim=1) ! naux
            nodes = dnodata
            this%var_ids%export = export_pkg%varids_aux(iaux, 1)
            do n = 1, nbound
              nodes(int1d(n)) = dbl2d(iaux, n)
            end do
            call nc_export_array(nodes, this%ncid, this%dim_ids, this%var_ids, &
                                 this%dis, idt, export_pkg%mf6_input%mempath, &
                                 nc_tag, export_pkg%mf6_input%subcomponent_name, &
                                 this%gridmap_name, this%latlon, this%deflate, &
                                 this%shuffle, this%chunk_z, this%chunk_y, &
                                 this%chunk_x, kper, iaux, this%nc_fname)
 
          end do
          deallocate (nodes)
        case default
        end select
      case default
        ! no-op, no other datatypes exported
      end select
    end do
 
    ! write to time coordinate variable
    call nf_verify(nf90_put_var(this%ncid, this%var_ids%time, &
                                totim, start=(/kper/)), &
                   this%nc_fname)
 
    ! synchronize file
    call nf_verify(nf90_sync(this%ncid), this%nc_fname)
  end subroutine package_step
 
  !> @brief determine packages to write gridded input
  !<
  subroutine add_pkg_data(this)
    use simvariablesmodule, only: idm_context
    use memorymanagerextmodule, only: mem_set_value
    use memoryhelpermodule, only: create_mem_path
    use idmdfnselectormodule, only: param_definitions
    use sourcecommonmodule, only: idm_subcomponent_type
    class(disncstructuredtype), intent(inout) :: this
    character(LENCOMPONENTNAME) :: ptype, pname, pkgtype
    type(characterstringtype), dimension(:), contiguous, &
      pointer :: pkgtypes => null()
    type(characterstringtype), dimension(:), contiguous, &
      pointer :: pkgnames => null()
    type(characterstringtype), dimension(:), contiguous, &
      pointer :: mempaths => null()
    type(inputparamdefinitiontype), dimension(:), pointer :: param_dfns
    character(len=LENMEMPATH) :: input_mempath, mempath
    integer(I4B) :: n
    integer(I4B), pointer :: export_arrays
    logical(LGP) :: found
 
    input_mempath = create_mem_path(component=this%modelname, context=idm_context)
 
    ! set pointers to model path package info
    call mem_setptr(pkgtypes, 'PKGTYPES', input_mempath)
    call mem_setptr(pkgnames, 'PKGNAMES', input_mempath)
    call mem_setptr(mempaths, 'MEMPATHS', input_mempath)
 
    do n = 1, size(mempaths)
      ! allocate export_arrays
      allocate (export_arrays)
      export_arrays = 0
 
      ! set package attributes
      mempath = mempaths(n)
      pname = pkgnames(n)
      ptype = pkgtypes(n)
 
      ! export input arrays
      if (mempath /= '') then
        ! update export
        call mem_set_value(export_arrays, 'EXPORT_NC', mempath, found, &
                           release=.false.)
 
        if (export_arrays > 0) then
          pkgtype = idm_subcomponent_type(this%modeltype, ptype)
          param_dfns => param_definitions(this%modeltype, pkgtype)
          call this%export_input_arrays(ptype, pname, mempath, param_dfns)
        end if
      end if
 
      ! cleanup
      deallocate (export_arrays)
    end do
  end subroutine add_pkg_data
 
  !> @brief create file (group) attributes
  !<
  subroutine add_global_att(this)
    class(disncstructuredtype), intent(inout) :: this
    ! file scoped title
    call nf_verify(nf90_put_att(this%ncid, nf90_global, 'title', &
                                this%annotation%title), this%nc_fname)
    ! source (MODFLOW 6)
    call nf_verify(nf90_put_att(this%ncid, nf90_global, 'source', &
                                this%annotation%source), this%nc_fname)
    ! grid type (MODFLOW 6)
    call nf_verify(nf90_put_att(this%ncid, nf90_global, 'modflow_grid', &
                                this%annotation%grid), this%nc_fname)
    ! MODFLOW 6 model type
    call nf_verify(nf90_put_att(this%ncid, nf90_global, 'modflow_model', &
                                this%annotation%model), this%nc_fname)
    ! generation datetime
    call nf_verify(nf90_put_att(this%ncid, nf90_global, 'history', &
                                this%annotation%history), this%nc_fname)
    ! supported conventions
    call nf_verify(nf90_put_att(this%ncid, nf90_global, 'Conventions', &
                                this%annotation%conventions), &
                   this%nc_fname)
  end subroutine add_global_att
 
  !> @brief netcdf export define dimensions
  !<
  subroutine define_dim(this)
    use constantsmodule, only: mvalidate
    use simvariablesmodule, only: isim_mode
    class(disncstructuredtype), intent(inout) :: this
 
    ! bound dim
    call nf_verify(nf90_def_dim(this%ncid, 'bnd', 2, this%dim_ids%bnd), &
                   this%nc_fname)
 
    if (isim_mode /= mvalidate .or. this%pkglist%Count() > 0) then
      ! Time
      call nf_verify(nf90_def_dim(this%ncid, 'time', this%totnstp, &
                                  this%dim_ids%time), this%nc_fname)
      call nf_verify(nf90_def_var(this%ncid, 'time', nf90_double, &
                                  this%dim_ids%time, this%var_ids%time), &
                     this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%time, 'calendar', &
                                  'standard'), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%time, 'units', &
                                  this%datetime), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%time, 'axis', 'T'), &
                     this%nc_fname)
      !call nf_verify(nf90_put_att(ncid, var_ids%time, 'bounds', 'time_bnds'), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%time, 'standard_name', &
                                  'time'), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%time, 'long_name', &
                                  'time'), this%nc_fname)
    end if
 
    ! Z dimension
    call nf_verify(nf90_def_dim(this%ncid, 'layer', this%dis%nlay, &
                                this%dim_ids%layer), this%nc_fname)
    call nf_verify(nf90_def_var(this%ncid, 'layer', nf90_int, &
                                this%dim_ids%layer, this%var_ids%layer), &
                   this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%layer, 'units', '1'), &
                   this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%layer, 'axis', 'Z'), &
                   this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%layer, 'positive', &
                                'down'), this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%layer, 'long_name', &
                                'model layer'), this%nc_fname)
 
    ! Y dimension
    call nf_verify(nf90_def_dim(this%ncid, 'y', this%dis%nrow, this%dim_ids%y), &
                   this%nc_fname)
    call nf_verify(nf90_def_var(this%ncid, 'y', nf90_double, this%dim_ids%y, &
                                this%var_ids%y), this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'units', &
                                this%lenunits), this%nc_fname)
    if (this%dis%angrot == dzero .or. &
        wrap_rotated_crs(this%crs_wkt, this%dis%xorigin, this%dis%yorigin, &
                         this%dis%angrot) /= '') then
      ! axis/standard_name assert that y holds the true projected
      ! position. For a rotated grid this only holds once crs_wkt
      ! encodes the rotation (see wrap_rotated_crs in define_gridmap);
      ! otherwise y is grid-local (see add_grid_data) and asserting it
      ! would cause GDAL's netCDF driver to derive an incorrect,
      ! unrotated position directly from y's grid-local values.
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'axis', 'Y'), &
                     this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'standard_name', &
                                  'projection_y_coordinate'), this%nc_fname)
    end if
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'long_name', &
                                'Northing'), this%nc_fname)
    if (this%gridmap_name /= '') then
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'grid_mapping', &
                                  this%gridmap_name), this%nc_fname)
    end if
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%y, 'bounds', 'y_bnds'), &
                   this%nc_fname)
    call nf_verify(nf90_def_var(this%ncid, 'y_bnds', nf90_double, &
                                (/this%dim_ids%bnd, this%dim_ids%y/), &
                                this%var_ids%y_bnds), this%nc_fname)
 
    ! X dimension
    call nf_verify(nf90_def_dim(this%ncid, 'x', this%dis%ncol, this%dim_ids%x), &
                   this%nc_fname)
    call nf_verify(nf90_def_var(this%ncid, 'x', nf90_double, this%dim_ids%x, &
                                this%var_ids%x), this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'units', &
                                this%lenunits), this%nc_fname)
    if (this%dis%angrot == dzero .or. &
        wrap_rotated_crs(this%crs_wkt, this%dis%xorigin, this%dis%yorigin, &
                         this%dis%angrot) /= '') then
      ! see matching comment on the Y dimension above
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'axis', 'X'), &
                     this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'standard_name', &
                                  'projection_x_coordinate'), this%nc_fname)
    end if
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'long_name', &
                                'Easting'), this%nc_fname)
    if (this%gridmap_name /= '') then
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'grid_mapping', &
                                  this%gridmap_name), this%nc_fname)
    end if
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%x, 'bounds', 'x_bnds'), &
                   this%nc_fname)
    call nf_verify(nf90_def_var(this%ncid, 'x_bnds', nf90_double, &
                                (/this%dim_ids%bnd, this%dim_ids%x/), &
                                this%var_ids%x_bnds), this%nc_fname)
 
    ! z: cell center elevation auxiliary coordinate -- a real vertical
    ! position, unlike the discrete layer index.  Always written,
    ! independent of CRS/NCF configuration.
    call nf_verify(nf90_def_var(this%ncid, 'z', nf90_double, &
                                (/this%dim_ids%x, this%dim_ids%y, &
                                  this%dim_ids%layer/), &
                                this%var_ids%elevation), this%nc_fname)
    call ncvar_chunk3d(this%ncid, this%var_ids%elevation, this%chunk_x, &
                       this%chunk_y, this%chunk_z, this%nc_fname)
    call ncvar_deflate(this%ncid, this%var_ids%elevation, this%deflate, &
                       this%shuffle, this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%elevation, 'units', &
                                this%lenunits), this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%elevation, &
                                'standard_name', 'altitude'), this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%elevation, &
                                'positive', 'up'), this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%elevation, &
                                'long_name', 'cell center elevation'), &
                   this%nc_fname)
  end subroutine define_dim
 
  !> @brief create the model layer dependent variables
  !<
  subroutine define_dependent(this)
    use constantsmodule, only: dhnoflo
    class(disncstructuredtype), intent(inout) :: this
 
    call nf_verify(nf90_def_var(this%ncid, this%xname, nf90_double, &
                                (/this%dim_ids%x, this%dim_ids%y, &
                                  this%dim_ids%layer, this%dim_ids%time/), &
                                this%var_ids%dependent), &
                   this%nc_fname)
 
    ! apply chunking parameters
    if (this%chunking_active) then
      call nf_verify(nf90_def_var_chunking(this%ncid, &
                                           this%var_ids%dependent, &
                                           nf90_chunked, &
                                           (/this%chunk_x, this%chunk_y, &
                                             this%chunk_z, this%chunk_time/)), &
                     this%nc_fname)
    end if
 
    ! deflate and shuffle
    call ncvar_deflate(this%ncid, this%var_ids%dependent, this%deflate, &
                       this%shuffle, this%nc_fname)
 
    ! put attr
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, &
                                'units', this%lenunits), this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, 'long_name', &
                                this%annotation%longname), this%nc_fname)
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, '_FillValue', &
                                (/dhnoflo/)), this%nc_fname)
    ! cell_methods (CF-1.13 7.3): dependent variable values are the
    ! instantaneous simulated state at each output time, not a time
    ! mean/accumulation over the interval -- an intensive quantity whose
    ! default method is already "point" per CF's own Appendix E, stated
    ! explicitly here per CF's recommendation to do so for every
    ! spatio-temporal dimension.
    call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, &
                                'cell_methods', 'time: point'), this%nc_fname)
 
    ! add grid mapping
    call ncvar_gridmap(this%ncid, this%var_ids%dependent, this%gridmap_name, &
                       this%latlon, .true., this%nc_fname)
  end subroutine define_dependent
 
  !> @brief create the file grid mapping container variable
  !<
  subroutine define_gridmap(this)
    class(disncstructuredtype), intent(inout) :: this
    integer(I4B) :: var_id
    character(len=LINELENGTH) :: gmname
    character(len=LENBIGLINE) :: effective_crs_wkt, derived_wkt
    character(len=100) :: rotinfo
    real(DP) :: longitude_of_central_meridian, latitude_of_projection_origin
    real(DP) :: scale_factor_at_central_meridian, false_easting, false_northing
    real(DP) :: semi_major_axis, inverse_flattening
    real(DP) :: standard_parallel_1, standard_parallel_2
    logical(LGP) :: params_found, is_2sp
 
    if (this%wkt /= '' .or. this%crs_wkt /= '') then
      call nf_verify(nf90_redef(this%ncid), this%nc_fname)
      call nf_verify(nf90_def_var(this%ncid, this%gridmap_name, nf90_int, &
                                  var_id), this%nc_fname)
      ! wkt (WKT1, OGC 01-009) -- for legacy/QGIS compatibility.  WKT1
      ! has no derived-CRS syntax, so it cannot encode grid rotation and
      ! is always written unmodified.
      if (this%wkt /= '') then
        call nf_verify(nf90_put_att(this%ncid, var_id, 'wkt', this%wkt), &
                       this%nc_fname)
      end if
      ! crs_wkt (WKT2, ISO 19162:2019) -- required by CF-1.13
      if (this%crs_wkt /= '') then
        effective_crs_wkt = this%crs_wkt
      else
        effective_crs_wkt = this%wkt
      end if
      ! For a rotated structured grid, wrap crs_wkt in a derived CRS
      ! encoding the rotation (see wrap_rotated_crs) so a CRS-aware
      ! consumer resolves the correct real-world position and shape
      ! via ordinary reprojection.  Falls back to an unmodified
      ! crs_wkt, with a warning, if CRS_WKT was not provided or is not
      ! a recognized WKT2 PROJCRS.
      if (this%dis%angrot /= dzero) then
        if (this%crs_wkt == '') then
          write (warnmsg, '(a)') 'Structured rotated grid with no CRS_WKT &
            &input: no rotation-positioning information will be written &
            &to the output CRS; x/y coordinate values will not reflect &
            &the true rotated position. Applies to file "'// &
            trim(this%nc_fname)//'".'
          call store_warning(warnmsg)
        else
          derived_wkt = wrap_rotated_crs(this%crs_wkt, this%dis%xorigin, &
                                         this%dis%yorigin, this%dis%angrot)
          if (derived_wkt /= '') then
            write (rotinfo, '(a,g0,a,g0,a,g0,a)') 'xorigin=', &
              this%dis%xorigin, ', yorigin=', this%dis%yorigin, &
              ', angrot=', this%dis%angrot, ' deg'
            write (warnmsg, '(a)') 'Structured rotated grid: wrapped &
              &CRS_WKT in a derived CRS to encode grid rotation ('// &
              trim(rotinfo)//'). Applies to file "'// &
              trim(this%nc_fname)//'".'
            call store_warning(warnmsg)
            effective_crs_wkt = derived_wkt
          else
            write (warnmsg, '(a)') 'Structured rotated grid: CRS_WKT is &
              &not a recognized WKT2 PROJCRS; no rotation-positioning &
              &information will be written to the output CRS; x/y &
              &coordinate values will not reflect the true rotated &
              &position. Applies to file "'//trim(this%nc_fname)//'".'
            call store_warning(warnmsg)
          end if
        end if
      end if
      call nf_verify(nf90_put_att(this%ncid, var_id, 'crs_wkt', &
                                  trim(effective_crs_wkt)), this%nc_fname)
      ! grid_mapping_name derived from WKT projection keyword
      gmname = wkt_to_cf_gridmapping(trim(effective_crs_wkt))
      if (gmname /= '') then
        call nf_verify(nf90_put_att(this%ncid, var_id, 'grid_mapping_name', &
                                    trim(gmname)), this%nc_fname)
      end if
      ! CF-standard numeric grid_mapping parameters (in addition to
      ! wkt/crs_wkt): some netCDF readers -- notably ArcGIS's classic
      ! netCDF connector -- position a projected grid using these
      ! individual attributes rather than parsing wkt/crs_wkt, and
      ! default every parameter to 0 if they are absent. Covers
      ! transverse_mercator, lambert_conformal_conic (2SP only -- CF has
      ! no scale_factor attribute for this method, so a 1SP/scale-factor
      ! WKT cannot be exactly represented), and albers_conical_equal_area;
      ! other grid_mapping_name values are unaffected. Best-effort
      ! extraction from effective_crs_wkt -- MF6 has no CRS parsing
      ! library, so this is skipped if any parameter cannot be located.
      if (trim(gmname) == 'transverse_mercator') then
        call wkt_transverse_mercator_params(trim(effective_crs_wkt), &
                                            longitude_of_central_meridian, &
                                            latitude_of_projection_origin, &
                                            scale_factor_at_central_meridian, &
                                            false_easting, false_northing, &
                                            semi_major_axis, inverse_flattening, &
                                            params_found)
        if (params_found) then
          call nf_verify(nf90_put_att(this%ncid, var_id, &
                                      'longitude_of_central_meridian', &
                                      longitude_of_central_meridian), &
                         this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, &
                                      'latitude_of_projection_origin', &
                                      latitude_of_projection_origin), &
                         this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, &
                                      'scale_factor_at_central_meridian', &
                                      scale_factor_at_central_meridian), &
                         this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'false_easting', &
                                      false_easting), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'false_northing', &
                                      false_northing), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'semi_major_axis', &
                                      semi_major_axis), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'inverse_flattening', &
                                      inverse_flattening), this%nc_fname)
        else
          write (warnmsg, '(a)') 'Structured grid: unable to extract CF &
            &transverse_mercator grid_mapping parameters from CRS_WKT; &
            &some netCDF readers may not correctly position this grid. &
            &Applies to file "'//trim(this%nc_fname)//'".'
          call store_warning(warnmsg)
        end if
      else if (trim(gmname) == 'lambert_conformal_conic') then
        call wkt_lambert_conformal_conic_params(trim(effective_crs_wkt), &
                                                standard_parallel_1, &
                                                standard_parallel_2, &
                                                longitude_of_central_meridian, &
                                                latitude_of_projection_origin, &
                                                false_easting, false_northing, &
                                                semi_major_axis, &
                                                inverse_flattening, is_2sp, &
                                                params_found)
        if (params_found) then
          call nf_verify(nf90_put_att(this%ncid, var_id, 'standard_parallel', &
                                      (/standard_parallel_1, &
                                        standard_parallel_2/)), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, &
                                      'longitude_of_central_meridian', &
                                      longitude_of_central_meridian), &
                         this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, &
                                      'latitude_of_projection_origin', &
                                      latitude_of_projection_origin), &
                         this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'false_easting', &
                                      false_easting), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'false_northing', &
                                      false_northing), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'semi_major_axis', &
                                      semi_major_axis), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'inverse_flattening', &
                                      inverse_flattening), this%nc_fname)
        else if (.not. is_2sp) then
          write (warnmsg, '(a)') 'Structured grid: 1SP lambert_conformal_conic &
            &variants are not supported for CF grid_mapping parameter &
            &extraction (CF has no scale_factor attribute for this method); &
            &some netCDF readers may not correctly position this grid. &
            &Applies to file "'//trim(this%nc_fname)//'".'
          call store_warning(warnmsg)
        else
          write (warnmsg, '(a)') 'Structured grid: unable to extract CF &
            &lambert_conformal_conic grid_mapping parameters from CRS_WKT; &
            &some netCDF readers may not correctly position this grid. &
            &Applies to file "'//trim(this%nc_fname)//'".'
          call store_warning(warnmsg)
        end if
      else if (trim(gmname) == 'albers_conical_equal_area') then
        call wkt_albers_conical_equal_area_params(trim(effective_crs_wkt), &
                                                  standard_parallel_1, &
                                                  standard_parallel_2, &
                                                  longitude_of_central_meridian, &
                                                  latitude_of_projection_origin, &
                                                  false_easting, false_northing, &
                                                  semi_major_axis, &
                                                  inverse_flattening, &
                                                  params_found)
        if (params_found) then
          call nf_verify(nf90_put_att(this%ncid, var_id, 'standard_parallel', &
                                      (/standard_parallel_1, &
                                        standard_parallel_2/)), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, &
                                      'longitude_of_central_meridian', &
                                      longitude_of_central_meridian), &
                         this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, &
                                      'latitude_of_projection_origin', &
                                      latitude_of_projection_origin), &
                         this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'false_easting', &
                                      false_easting), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'false_northing', &
                                      false_northing), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'semi_major_axis', &
                                      semi_major_axis), this%nc_fname)
          call nf_verify(nf90_put_att(this%ncid, var_id, 'inverse_flattening', &
                                      inverse_flattening), this%nc_fname)
        else
          write (warnmsg, '(a)') 'Structured grid: unable to extract CF &
            &albers_conical_equal_area grid_mapping parameters from &
            &CRS_WKT; some netCDF readers may not correctly position this &
            &grid. Applies to file "'//trim(this%nc_fname)//'".'
          call store_warning(warnmsg)
        end if
      end if
      ! GeoTransform/spatial_ref intentionally not written: a GDAL-only,
      ! non-CF extension ArcGIS does not read, superseded by crs_wkt's
      ! derived-CRS encoding above, which any CRS-aware consumer resolves
      ! via ordinary reprojection.
      call nf_verify(nf90_enddef(this%ncid), this%nc_fname)
      call nf_verify(nf90_put_var(this%ncid, var_id, 1), &
                     this%nc_fname)
    end if
  end subroutine define_gridmap
 
  !> @brief define grid projection variables
  !<
  subroutine define_geocoords(this)
    class(disncstructuredtype), intent(inout) :: this
    if (this%latlon) then
      ! lat
      call nf_verify(nf90_def_var(this%ncid, 'lat', nf90_double, &
                                  (/this%dim_ids%x, this%dim_ids%y/), &
                                  this%var_ids%latitude), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%latitude, &
                                  'units', 'degrees_north'), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%latitude, &
                                  'standard_name', 'latitude'), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%latitude, &
                                  'long_name', 'latitude'), this%nc_fname)
 
      ! lon
      call nf_verify(nf90_def_var(this%ncid, 'lon', nf90_double, &
                                  (/this%dim_ids%x, this%dim_ids%y/), &
                                  this%var_ids%longitude), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%longitude, &
                                  'units', 'degrees_east'), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%longitude, &
                                  'standard_name', 'longitude'), this%nc_fname)
      call nf_verify(nf90_put_att(this%ncid, this%var_ids%longitude, &
                                  'long_name', 'longitude'), this%nc_fname)
    end if
  end subroutine define_geocoords
 
  !> @brief add grid projection data
  !<
  subroutine add_proj_data(this)
    class(disncstructuredtype), intent(inout) :: this
    if (this%latlon) then
      ! lat
      call nf_verify(nf90_put_var(this%ncid, this%var_ids%latitude, &
                                  this%latitude, start=(/1, 1/), &
                                  count=(/this%dis%ncol, this%dis%nrow/)), &
                     this%nc_fname)
 
      ! lon
      call nf_verify(nf90_put_var(this%ncid, this%var_ids%longitude, &
                                  this%longitude, start=(/1, 1/), &
                                  count=(/this%dis%ncol, this%dis%nrow/)), &
                     this%nc_fname)
    end if
  end subroutine add_proj_data
 
  !> @brief add grid coordinates
  !<
  subroutine add_grid_data(this)
    class(disncstructuredtype), intent(inout) :: this
    integer(I4B) :: ibnd, n, k
    real(DP), dimension(:, :), pointer, contiguous :: dbl2d
    real(DP), dimension(:), allocatable :: x, y
    real(DP), dimension(:, :, :), allocatable :: elev3d
    real(DP) :: xoff, yoff
 
    if (this%dis%angrot /= dzero) then
      xoff = dzero
      yoff = dzero
    else
      xoff = this%dis%xorigin
      yoff = this%dis%yorigin
    end if
 
    allocate (x(size(this%dis%cellx)))
    allocate (y(size(this%dis%celly)))
 
    do n = 1, size(this%dis%cellx)
      x(n) = this%dis%cellx(n) + xoff
    end do
 
    do n = 1, size(this%dis%celly)
      y(n) = this%dis%celly(n) + yoff
    end do
 
    call nf_verify(nf90_put_var(this%ncid, this%var_ids%x, x), &
                   this%nc_fname)
    call nf_verify(nf90_put_var(this%ncid, this%var_ids%y, y), &
                   this%nc_fname)
    ! TODO see cf-conventions 4.3.3. Parametric Vertical Coordinate
    call nf_verify(nf90_put_var(this%ncid, this%var_ids%layer, this%layers), &
                   this%nc_fname)
 
    deallocate (x)
    deallocate (y)
 
    ! z: cell center elevation, layer 1 = midpoint of top/bot3d(1),
    ! layer k>1 = midpoint of bot3d(k-1)/bot3d(k)
    allocate (elev3d(this%dis%ncol, this%dis%nrow, this%nlay))
    do k = 1, this%nlay
      if (k == 1) then
        elev3d(:, :, k) = (this%dis%top2d(:, :) + this%dis%bot3d(:, :, k)) * &
                          dhalf
      else
        elev3d(:, :, k) = (this%dis%bot3d(:, :, k - 1) + &
                           this%dis%bot3d(:, :, k)) * dhalf
      end if
    end do
    call nf_verify(nf90_put_var(this%ncid, this%var_ids%elevation, elev3d), &
                   this%nc_fname)
    deallocate (elev3d)
 
    ! bounds x
    allocate (dbl2d(2, size(this%dis%cellx)))
    ibnd = 1
    do n = 1, size(this%dis%cellx)
      if (ibnd == 1) then
        dbl2d(1, ibnd) = xoff
        dbl2d(2, ibnd) = xoff + this%dis%delr(ibnd)
      else
        dbl2d(1, ibnd) = dbl2d(1, ibnd - 1) + this%dis%delr(ibnd)
        dbl2d(2, ibnd) = dbl2d(2, ibnd - 1) + this%dis%delr(ibnd)
      end if
      ibnd = ibnd + 1
    end do
    call nf_verify(nf90_put_var(this%ncid, this%var_ids%x_bnds, dbl2d), &
                   this%nc_fname)
    deallocate (dbl2d)
 
    ! bounds y
    allocate (dbl2d(2, size(this%dis%celly)))
    ibnd = 1
    do n = size(this%dis%celly), 1, -1
      if (ibnd == 1) then
        dbl2d(1, ibnd) = yoff + sum(this%dis%delc) - this%dis%delc(n)
        dbl2d(2, ibnd) = yoff + sum(this%dis%delc)
      else
        dbl2d(1, ibnd) = dbl2d(1, ibnd - 1) - this%dis%delc(n)
        dbl2d(2, ibnd) = dbl2d(2, ibnd - 1) - this%dis%delc(n)
      end if
      ibnd = ibnd + 1
    end do
    call nf_verify(nf90_put_var(this%ncid, this%var_ids%y_bnds, dbl2d), &
                   this%nc_fname)
    deallocate (dbl2d)
  end subroutine add_grid_data
 
  !> @brief define 2d variable chunking
  !<
  subroutine ncvar_chunk2d(ncid, varid, chunk_x, chunk_y, nc_fname)
    integer(I4B), intent(in) :: ncid
    integer(I4B), intent(in) :: varid
    integer(I4B), intent(in) :: chunk_x
    integer(I4B), intent(in) :: chunk_y
    character(len=*), intent(in) :: nc_fname
    if (chunk_y > 0 .and. chunk_x > 0) then
      call nf_verify(nf90_def_var_chunking(ncid, varid, nf90_chunked, &
                                           (/chunk_x, chunk_y/)), nc_fname)
    end if
  end subroutine ncvar_chunk2d
 
  !> @brief define 3d variable chunking
  !<
  subroutine ncvar_chunk3d(ncid, varid, chunk_x, chunk_y, chunk_z, nc_fname)
    integer(I4B), intent(in) :: ncid
    integer(I4B), intent(in) :: varid
    integer(I4B), intent(in) :: chunk_x
    integer(I4B), intent(in) :: chunk_y
    integer(I4B), intent(in) :: chunk_z
    character(len=*), intent(in) :: nc_fname
    if (chunk_z > 0 .and. chunk_y > 0 .and. chunk_x > 0) then
      call nf_verify(nf90_def_var_chunking(ncid, varid, nf90_chunked, &
                                           (/chunk_x, chunk_y, chunk_z/)), &
                     nc_fname)
    end if
  end subroutine ncvar_chunk3d
 
  !> @brief define variable compression
  !<
  subroutine ncvar_deflate(ncid, varid, deflate, shuffle, nc_fname)
    integer(I4B), intent(in) :: ncid
    integer(I4B), intent(in) :: varid
    integer(I4B), intent(in) :: deflate
    integer(I4B), intent(in) :: shuffle
    character(len=*), intent(in) :: nc_fname
    ! deflate and shuffle
    if (deflate >= 0) then
      call nf_verify(nf90_def_var_deflate(ncid, varid, shuffle=shuffle, &
                                          deflate=1, deflate_level=deflate), &
                     nc_fname)
    end if
  end subroutine ncvar_deflate
 
  !> @brief put variable gridmap attributes
  !!
  !! has_layer must be true only for variables carrying the layer, y, and x
  !! dimensions -- z(layer, y, x) can only be legally referenced via
  !! `coordinates` from a variable whose dimensions are a superset of its
  !! own (CF-1.13 5.2). 2D (y, x)-only variables such as dis_top must pass
  !! has_layer=.false. z is written unconditionally, independent of
  !! gridmap_name/latlon.
  !<
  subroutine ncvar_gridmap(ncid, varid, gridmap_name, latlon, has_layer, &
                           nc_fname)
    integer(I4B), intent(in) :: ncid
    integer(I4B), intent(in) :: varid
    character(len=*), intent(in) :: gridmap_name
    logical(LGP), intent(in) :: latlon
    logical(LGP), intent(in) :: has_layer
    character(len=*), intent(in) :: nc_fname
    character(len=LINELENGTH) :: coords
    coords = ''
    ! NOTE: x/y are not listed here even when gridmap_name is set -- unlike
    ! lon/lat (real 2D auxiliary coordinates), x/y are true CF dimension
    ! coordinates already discoverable by dimension-name matching alone, so
    ! listing them in 'coordinates' is redundant (CF-1.13 Ch.5 preamble).
    if (latlon) then
      coords = 'lon lat'
    end if
    if (has_layer) then
      if (len_trim(coords) > 0) then
        coords = trim(coords)//' z'
      else
        coords = 'z'
      end if
    end if
    if (len_trim(coords) > 0) then
      call nf_verify(nf90_put_att(ncid, varid, 'coordinates', trim(coords)), &
                     nc_fname)
    end if
    if (gridmap_name /= '') then
      call nf_verify(nf90_put_att(ncid, varid, 'grid_mapping', gridmap_name), &
                     nc_fname)
    end if
  end subroutine ncvar_gridmap
 
  !> @brief put variable internal modflow6 attributes
  !<
  subroutine ncvar_mf6attr(ncid, varid, iaux, nc_tag, nc_fname)
    integer(I4B), intent(in) :: ncid
    integer(I4B), intent(in) :: varid
    integer(I4B), intent(in) :: iaux
    character(len=*), intent(in) :: nc_tag
    character(len=*), intent(in) :: nc_fname
    if (nc_tag /= '') then
      call nf_verify(nf90_put_att(ncid, varid, 'modflow_input', &
                                  nc_tag), nc_fname)
      if (iaux > 0) then
        call nf_verify(nf90_put_att(ncid, varid, 'modflow_iaux', &
                                    iaux), nc_fname)
      end if
    end if
  end subroutine ncvar_mf6attr
 
  !> @brief netcdf export 1D integer
  !<
  subroutine nc_export_int1d(p_mem, ncid, dim_ids, var_ids, dis, idt, mempath, &
                             nc_tag, pkgname, gridmap_name, latlon, deflate, &
                             shuffle, chunk_z, chunk_y, chunk_x, iper, nc_fname)
    use tdismodule, only: kper
    integer(I4B), dimension(:), pointer, contiguous, intent(in) :: p_mem
    integer(I4B), intent(in) :: ncid
    type(structuredncdimidtype), intent(inout) :: dim_ids
    type(structuredncvaridtype), intent(inout) :: var_ids
    type(distype), pointer, intent(in) :: dis
    type(inputparamdefinitiontype), pointer, intent(in) :: idt
    character(len=*), intent(in) :: mempath
    character(len=*), intent(in) :: nc_tag
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: gridmap_name
    logical(LGP), intent(in) :: latlon
    integer(I4B), intent(in) :: deflate
    integer(I4B), intent(in) :: shuffle
    integer(I4B), intent(in) :: chunk_z
    integer(I4B), intent(in) :: chunk_y
    integer(I4B), intent(in) :: chunk_x
    integer(I4B), intent(in) :: iper
    character(len=*), intent(in) :: nc_fname
    integer(I4B) :: var_id, axis_sz
    character(len=LINELENGTH) :: varname, longname
 
    varname = export_varname(pkgname, idt%tagname, mempath)
 
    if (idt%shape == 'NROW' .or. &
        idt%shape == 'NCOL' .or. &
        idt%shape == 'NCPL' .or. &
        idt%shape == 'NAUX NCPL') then
 
      if (iper == 0) then
        select case (idt%shape)
        case ('NROW')
          axis_sz = dim_ids%y
        case ('NCOL')
          axis_sz = dim_ids%x
        end select
 
        longname = export_longname(idt%longname, pkgname, idt%tagname, mempath, &
                                   component_type=idt%component_type, &
                                   subcomponent_type=idt%subcomponent_type)
 
        ! reenter define mode and create variable
        call nf_verify(nf90_redef(ncid), nc_fname)
        call nf_verify(nf90_def_var(ncid, varname, nf90_int, &
                                    (/axis_sz/), var_id), &
                       nc_fname)
 
        ! NROW/NCOL shapes use default chunking
        call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname)
 
        ! put attr
        call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', &
                                    (/nf90_fill_int/)), nc_fname)
        call nf_verify(nf90_put_att(ncid, var_id, 'long_name', &
                                    longname), nc_fname)
 
        ! add mf6 attr
        call ncvar_mf6attr(ncid, var_id, 0, nc_tag, nc_fname)
 
        ! exit define mode and write data
        call nf_verify(nf90_enddef(ncid), nc_fname)
        call nf_verify(nf90_put_var(ncid, var_id, p_mem), &
                       nc_fname)
      else
        ! timeseries
        call nf_verify(nf90_put_var(ncid, &
                                    var_ids%export, p_mem, &
                                    start=(/1, 1, kper/), &
                                    count=(/dis%ncol, dis%nrow, 1/)), nc_fname)
      end if
 
    else
 
      if (iper == 0) then
        ! reenter define mode and create variable
        call nf_verify(nf90_redef(ncid), nc_fname)
        call nf_verify(nf90_def_var(ncid, varname, nf90_int, &
                                    (/dim_ids%x, dim_ids%y, dim_ids%layer/), &
                                    var_id), nc_fname)
 
        ! apply chunking parameters
        call ncvar_chunk3d(ncid, var_id, chunk_x, chunk_y, chunk_z, nc_fname)
        ! deflate and shuffle
        call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname)
 
        ! put attr
        call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', &
                                    (/nf90_fill_int/)), nc_fname)
        call nf_verify(nf90_put_att(ncid, var_id, 'long_name', &
                                    idt%longname), nc_fname)
 
        ! add grid mapping and mf6 attr
        call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, .true., &
                           nc_fname)
        call ncvar_mf6attr(ncid, var_id, 0, nc_tag, nc_fname)
 
        ! exit define mode and write data
        call nf_verify(nf90_enddef(ncid), nc_fname)
        call nf_verify(nf90_put_var(ncid, var_id, p_mem, start=(/1, 1, 1/), &
                                    count=(/dis%ncol, dis%nrow, dis%nlay/)), &
                       nc_fname)
      else
        ! timeseries
        call nf_verify(nf90_put_var(ncid, &
                                    var_ids%export, p_mem, &
                                    start=(/1, 1, 1, kper/), &
                                    count=(/dis%ncol, dis%nrow, dis%nlay, 1/)), &
                       nc_fname)
      end if
    end if
  end subroutine nc_export_int1d
 
  !> @brief netcdf export 2D integer
  !<
  subroutine nc_export_int2d(p_mem, ncid, dim_ids, var_ids, dis, idt, mempath, &
                             nc_tag, pkgname, gridmap_name, latlon, deflate, &
                             shuffle, chunk_z, chunk_y, chunk_x, nc_fname)
    integer(I4B), dimension(:, :), pointer, contiguous, intent(in) :: p_mem
    integer(I4B), intent(in) :: ncid
    type(structuredncdimidtype), intent(inout) :: dim_ids
    type(structuredncvaridtype), intent(inout) :: var_ids
    type(distype), pointer, intent(in) :: dis
    type(inputparamdefinitiontype), pointer, intent(in) :: idt
    character(len=*), intent(in) :: mempath
    character(len=*), intent(in) :: nc_tag
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: gridmap_name
    logical(LGP), intent(in) :: latlon
    integer(I4B), intent(in) :: deflate
    integer(I4B), intent(in) :: shuffle
    integer(I4B), intent(in) :: chunk_z
    integer(I4B), intent(in) :: chunk_y
    integer(I4B), intent(in) :: chunk_x
    character(len=*), intent(in) :: nc_fname
    character(len=LINELENGTH) :: varname
    integer(I4B) :: var_id
 
    varname = export_varname(pkgname, idt%tagname, mempath)
 
    ! reenter define mode and create variable
    call nf_verify(nf90_redef(ncid), nc_fname)
    call nf_verify(nf90_def_var(ncid, varname, nf90_int, &
                                (/dim_ids%x, dim_ids%y/), var_id), &
                   nc_fname)
 
    ! apply chunking parameters
    call ncvar_chunk2d(ncid, var_id, chunk_x, chunk_y, nc_fname)
    ! deflate and shuffle
    call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname)
 
    ! put attr
    call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', &
                                (/nf90_fill_int/)), nc_fname)
    call nf_verify(nf90_put_att(ncid, var_id, 'long_name', &
                                idt%longname), nc_fname)
 
    ! add grid mapping and mf6 attr
    call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, .false., &
                       nc_fname)
    call ncvar_mf6attr(ncid, var_id, 0, nc_tag, nc_fname)
 
    ! exit define mode and write data
    call nf_verify(nf90_enddef(ncid), nc_fname)
    call nf_verify(nf90_put_var(ncid, var_id, p_mem, start=(/1, 1/), &
                                count=(/dis%ncol, dis%nrow/)), &
                   nc_fname)
  end subroutine nc_export_int2d
 
  !> @brief netcdf export 3D integer
  !<
  subroutine nc_export_int3d(p_mem, ncid, dim_ids, var_ids, dis, idt, mempath, &
                             nc_tag, pkgname, gridmap_name, latlon, deflate, &
                             shuffle, chunk_z, chunk_y, chunk_x, nc_fname)
    integer(I4B), dimension(:, :, :), pointer, contiguous, intent(in) :: p_mem
    integer(I4B), intent(in) :: ncid
    type(structuredncdimidtype), intent(inout) :: dim_ids
    type(structuredncvaridtype), intent(inout) :: var_ids
    type(distype), pointer, intent(in) :: dis
    type(inputparamdefinitiontype), pointer, intent(in) :: idt
    character(len=*), intent(in) :: mempath
    character(len=*), intent(in) :: nc_tag
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: gridmap_name
    logical(LGP), intent(in) :: latlon
    integer(I4B), intent(in) :: deflate
    integer(I4B), intent(in) :: shuffle
    integer(I4B), intent(in) :: chunk_z
    integer(I4B), intent(in) :: chunk_y
    integer(I4B), intent(in) :: chunk_x
    character(len=*), intent(in) :: nc_fname
    character(len=LINELENGTH) :: varname
    integer(I4B) :: var_id
 
    varname = export_varname(pkgname, idt%tagname, mempath)
 
    ! reenter define mode and create variable
    call nf_verify(nf90_redef(ncid), nc_fname)
    call nf_verify(nf90_def_var(ncid, varname, nf90_int, &
                                (/dim_ids%x, dim_ids%y, dim_ids%layer/), &
                                var_id), nc_fname)
 
    ! apply chunking parameters
    call ncvar_chunk3d(ncid, var_id, chunk_x, chunk_y, chunk_z, nc_fname)
    ! deflate and shuffle
    call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname)
 
    ! put attr
    call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', &
                                (/nf90_fill_int/)), nc_fname)
    call nf_verify(nf90_put_att(ncid, var_id, 'long_name', &
                                idt%longname), nc_fname)
 
    ! add grid mapping and mf6 attr
    call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, .true., &
                       nc_fname)
    call ncvar_mf6attr(ncid, var_id, 0, nc_tag, nc_fname)
 
    ! exit define mode and write data
    call nf_verify(nf90_enddef(ncid), nc_fname)
    call nf_verify(nf90_put_var(ncid, var_id, p_mem, start=(/1, 1, 1/), &
                                count=(/dis%ncol, dis%nrow, dis%nlay/)), &
                   nc_fname)
  end subroutine nc_export_int3d
 
  !> @brief netcdf export 1D double
  !<
  subroutine nc_export_dbl1d(p_mem, ncid, dim_ids, var_ids, dis, idt, mempath, &
                             nc_tag, pkgname, gridmap_name, latlon, deflate, &
                             shuffle, chunk_z, chunk_y, chunk_x, iper, iaux, &
                             nc_fname)
    use tdismodule, only: kper
    real(DP), dimension(:), pointer, contiguous, intent(in) :: p_mem
    integer(I4B), intent(in) :: ncid
    type(structuredncdimidtype), intent(inout) :: dim_ids
    type(structuredncvaridtype), intent(inout) :: var_ids
    type(distype), pointer, intent(in) :: dis
    type(inputparamdefinitiontype), pointer, intent(in) :: idt
    character(len=*), intent(in) :: mempath
    character(len=*), intent(in) :: nc_tag
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: gridmap_name
    logical(LGP), intent(in) :: latlon
    integer(I4B), intent(in) :: deflate
    integer(I4B), intent(in) :: shuffle
    integer(I4B), intent(in) :: chunk_z
    integer(I4B), intent(in) :: chunk_y
    integer(I4B), intent(in) :: chunk_x
    integer(I4B), intent(in) :: iper
    integer(I4B), intent(in) :: iaux
    character(len=*), intent(in) :: nc_fname
    integer(I4B) :: var_id, axis_sz
    character(len=LINELENGTH) :: varname, longname
 
    if (idt%shape == 'NROW' .or. &
        idt%shape == 'NCOL' .or. &
        idt%shape == 'NCPL' .or. &
        idt%shape == 'NAUX NCPL') then
 
      if (iper == 0) then
        select case (idt%shape)
        case ('NROW')
          axis_sz = dim_ids%y
        case ('NCOL')
          axis_sz = dim_ids%x
        end select
 
        varname = export_varname(pkgname, idt%tagname, mempath)
        longname = export_longname(idt%longname, pkgname, idt%tagname, mempath, &
                                   iaux=iaux, &
                                   component_type=idt%component_type, &
                                   subcomponent_type=idt%subcomponent_type)
 
        ! reenter define mode and create variable
        call nf_verify(nf90_redef(ncid), nc_fname)
        call nf_verify(nf90_def_var(ncid, varname, nf90_double, &
                                    (/axis_sz/), var_id), &
                       nc_fname)
 
        ! NROW/NCOL shapes use default chunking
        call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname)
 
        ! put attr
        call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', &
                                    (/nf90_fill_double/)), nc_fname)
        call nf_verify(nf90_put_att(ncid, var_id, 'long_name', &
                                    longname), nc_fname)
 
        ! add mf6 attr
        call ncvar_mf6attr(ncid, var_id, iaux, nc_tag, nc_fname)
 
        ! exit define mode and write data
        call nf_verify(nf90_enddef(ncid), nc_fname)
        call nf_verify(nf90_put_var(ncid, var_id, p_mem), &
                       nc_fname)
      else
        ! timeseries
        call nf_verify(nf90_put_var(ncid, &
                                    var_ids%export, p_mem, &
                                    start=(/1, 1, kper/), &
                                    count=(/dis%ncol, dis%nrow, 1/)), nc_fname)
      end if
 
    else
 
      if (iper == 0) then
        varname = export_varname(pkgname, idt%tagname, mempath, iaux=iaux)
        longname = export_longname(idt%longname, pkgname, idt%tagname, mempath, &
                                   iaux=iaux, &
                                   component_type=idt%component_type, &
                                   subcomponent_type=idt%subcomponent_type)
 
        ! reenter define mode and create variable
        call nf_verify(nf90_redef(ncid), nc_fname)
        call nf_verify(nf90_def_var(ncid, varname, nf90_double, &
                                    (/dim_ids%x, dim_ids%y, dim_ids%layer/), &
                                    var_id), nc_fname)
 
        ! apply chunking parameters
        call ncvar_chunk3d(ncid, var_id, chunk_x, chunk_y, chunk_z, nc_fname)
        ! deflate and shuffle
        call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname)
 
        ! put attr
        call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', &
                                    (/nf90_fill_double/)), nc_fname)
        call nf_verify(nf90_put_att(ncid, var_id, 'long_name', &
                                    longname), nc_fname)
 
        ! add grid mapping and mf6 attr
        call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, .true., &
                           nc_fname)
        call ncvar_mf6attr(ncid, var_id, iaux, nc_tag, nc_fname)
 
        ! exit define mode and write data
        call nf_verify(nf90_enddef(ncid), nc_fname)
        call nf_verify(nf90_put_var(ncid, var_id, p_mem, start=(/1, 1, 1/), &
                                    count=(/dis%ncol, dis%nrow, dis%nlay/)), &
                       nc_fname)
      else
        ! timeseries
        call nf_verify(nf90_put_var(ncid, &
                                    var_ids%export, p_mem, &
                                    start=(/1, 1, 1, kper/), &
                                    count=(/dis%ncol, dis%nrow, dis%nlay, 1/)), &
                       nc_fname)
      end if
    end if
  end subroutine nc_export_dbl1d
 
  !> @brief netcdf export 2D double
  !<
  subroutine nc_export_dbl2d(p_mem, ncid, dim_ids, var_ids, dis, idt, mempath, &
                             nc_tag, pkgname, gridmap_name, latlon, deflate, &
                             shuffle, chunk_z, chunk_y, chunk_x, nc_fname)
    real(DP), dimension(:, :), pointer, contiguous, intent(in) :: p_mem
    integer(I4B), intent(in) :: ncid
    type(structuredncdimidtype), intent(inout) :: dim_ids
    type(structuredncvaridtype), intent(inout) :: var_ids
    type(distype), pointer, intent(in) :: dis
    type(inputparamdefinitiontype), pointer, intent(in) :: idt
    character(len=*), intent(in) :: mempath
    character(len=*), intent(in) :: nc_tag
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: gridmap_name
    logical(LGP), intent(in) :: latlon
    integer(I4B), intent(in) :: deflate
    integer(I4B), intent(in) :: shuffle
    integer(I4B), intent(in) :: chunk_z
    integer(I4B), intent(in) :: chunk_y
    integer(I4B), intent(in) :: chunk_x
    character(len=*), intent(in) :: nc_fname
    character(len=LINELENGTH) :: varname
    integer(I4B) :: var_id
 
    varname = export_varname(pkgname, idt%tagname, mempath)
 
    ! reenter define mode and create variable
    call nf_verify(nf90_redef(ncid), nc_fname)
    call nf_verify(nf90_def_var(ncid, varname, nf90_double, &
                                (/dim_ids%x, dim_ids%y/), var_id), &
                   nc_fname)
 
    ! apply chunking parameters
    call ncvar_chunk2d(ncid, var_id, chunk_x, chunk_y, nc_fname)
    ! deflate and shuffle
    call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname)
 
    ! put attr
    call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', &
                                (/nf90_fill_double/)), nc_fname)
    call nf_verify(nf90_put_att(ncid, var_id, 'long_name', &
                                idt%longname), nc_fname)
 
    ! add grid mapping and mf6 attr
    call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, .false., &
                       nc_fname)
    call ncvar_mf6attr(ncid, var_id, 0, nc_tag, nc_fname)
 
    ! exit define mode and write data
    call nf_verify(nf90_enddef(ncid), nc_fname)
    call nf_verify(nf90_put_var(ncid, var_id, p_mem, start=(/1, 1/), &
                                count=(/dis%ncol, dis%nrow/)), &
                   nc_fname)
  end subroutine nc_export_dbl2d
 
  !> @brief netcdf export 3D double
  !<
  subroutine nc_export_dbl3d(p_mem, ncid, dim_ids, var_ids, dis, idt, mempath, &
                             nc_tag, pkgname, gridmap_name, latlon, deflate, &
                             shuffle, chunk_z, chunk_y, chunk_x, nc_fname)
    real(DP), dimension(:, :, :), pointer, contiguous, intent(in) :: p_mem
    integer(I4B), intent(in) :: ncid
    type(structuredncdimidtype), intent(inout) :: dim_ids
    type(structuredncvaridtype), intent(inout) :: var_ids
    type(distype), pointer, intent(in) :: dis
    type(inputparamdefinitiontype), pointer, intent(in) :: idt
    character(len=*), intent(in) :: mempath
    character(len=*), intent(in) :: nc_tag
    character(len=*), intent(in) :: pkgname
    character(len=*), intent(in) :: gridmap_name
    logical(LGP), intent(in) :: latlon
    integer(I4B), intent(in) :: deflate
    integer(I4B), intent(in) :: shuffle
    integer(I4B), intent(in) :: chunk_z
    integer(I4B), intent(in) :: chunk_y
    integer(I4B), intent(in) :: chunk_x
    character(len=*), intent(in) :: nc_fname
    integer(I4B) :: var_id
    character(len=LINELENGTH) :: varname, longname
 
    varname = export_varname(pkgname, idt%tagname, mempath)
    longname = export_longname(idt%longname, pkgname, idt%tagname, mempath, &
                               component_type=idt%component_type, &
                               subcomponent_type=idt%subcomponent_type)
 
    ! reenter define mode and create variable
    call nf_verify(nf90_redef(ncid), nc_fname)
    call nf_verify(nf90_def_var(ncid, varname, nf90_double, &
                                (/dim_ids%x, dim_ids%y, dim_ids%layer/), &
                                var_id), nc_fname)
 
    ! apply chunking parameters
    call ncvar_chunk3d(ncid, var_id, chunk_x, chunk_y, chunk_z, nc_fname)
    ! deflate and shuffle
    call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname)
 
    ! put attr
    call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', &
                                (/nf90_fill_double/)), nc_fname)
    call nf_verify(nf90_put_att(ncid, var_id, 'long_name', &
                                longname), nc_fname)
 
    ! add grid mapping and mf6 attr
    call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, .true., &
                       nc_fname)
    call ncvar_mf6attr(ncid, var_id, 0, nc_tag, nc_fname)
 
    ! exit define mode and write data
    call nf_verify(nf90_enddef(ncid), nc_fname)
    call nf_verify(nf90_put_var(ncid, var_id, p_mem, start=(/1, 1, 1/), &
                                count=(/dis%ncol, dis%nrow, dis%nlay/)), &
                   nc_fname)
  end subroutine nc_export_dbl3d
 
end module disncstructuredmodule