LoadMf6File.f90

Go to the documentation of this file.

!> @brief This module contains the LoadMf6FileModule
!!
!! This module contains the input data model routines for
!! loading static data from a MODFLOW 6 input file using the
!! block parser.
!!
!<
module loadmf6filemodule
 
  use kindmodule, only: dp, i4b, lgp
  use simvariablesmodule, only: errmsg
  use simmodule, only: store_error
  use constantsmodule, only: linelength, lenvarname
  use blockparsermodule, only: blockparsertype
  use layeredarrayreadermodule, only: read_dbl1d_layered, &
                                      read_dbl2d_layered, &
                                      read_dbl3d_layered, &
                                      read_int1d_layered, &
                                      read_int2d_layered, &
                                      read_int3d_layered
  use double1dreadermodule, only: read_dbl1d
  use double2dreadermodule, only: read_dbl2d
  use integer1dreadermodule, only: read_int1d
  use integer2dreadermodule, only: read_int2d
  use inputoutputmodule, only: parseline
  use inputdefinitionmodule, only: inputparamdefinitiontype
  use definitionselectmodule, only: get_param_definition_type, &
                                    get_aggregate_definition_type
  use modflowinputmodule, only: modflowinputtype
  use memorymanagermodule, only: mem_allocate, mem_setptr, get_isize, get_mem_rank
  use structarraymodule, only: structarraytype
  use structvectormodule, only: structvectortype
  use ncfilevarsmodule, only: ncpackagevarstype
  use idmloggermodule, only: idm_log_var, idm_log_header, idm_log_close, &
                             idm_export
 
  implicit none
  private
  public :: loadmf6filetype
  public :: read_control_record
 
  !> @brief Fortran workaround for allocatable arrays of pointers; wraps a StructArray pointer for deferred TS linking.
  !<
  type :: staticsatype
    type(structarraytype), pointer :: sa => null()
  end type staticsatype
 
  !> @brief Static parser based input loader
  !!
  !! This type defines a static input context loader
  !! for traditional mf6 ascii input files.
  !!
  !<
  type :: loadmf6filetype
    type(blockparsertype), pointer :: parser !< ascii block parser
    integer(I4B), dimension(:), pointer, contiguous :: mshape => null() !< model shape
    type(structarraytype), pointer :: structarray => null() !< structarray for loading list input
    type(staticsatype), allocatable :: ts_sas(:) !< saved structarrays for deferred TS linking
    type(modflowinputtype) :: mf6_input !< description of input
    type(ncpackagevarstype), pointer :: nc_vars => null()
    character(len=LINELENGTH) :: filename !< name of ascii input file
    character(len=LINELENGTH), dimension(:), allocatable :: block_tags !< read block tags
    logical(LGP) :: ts_active !< is timeseries active
    logical(LGP) :: export !< is array export active
    logical(LGP) :: readasarrays
    logical(LGP) :: readarraygrid
    integer(I4B) :: inamedbound
    integer(I4B) :: iauxiliary
    integer(I4B) :: iout !< inunit for list log
  contains
    procedure :: load
    procedure :: init
    procedure :: load_block
    procedure :: finalize
    procedure :: parse_block
    procedure :: block_post_process
    procedure :: parse_io_tag
    procedure :: parse_record_tag
    procedure :: load_tag
    procedure :: block_index_dfn
    procedure :: parse_structarray_block
    procedure :: save_ts_sa
    procedure :: ts_sa_count
    procedure :: get_ts_sa
    procedure :: cleanup
  end type loadmf6filetype
 
contains
 
  !> @brief load all static input blocks
  !!
  !! Invoke this routine to load all static input blocks
  !! in single call.
  !!
  !<
  subroutine load(this, parser, mf6_input, nc_vars, filename, iout)
    class(loadmf6filetype) :: this
    type(blockparsertype), target, intent(inout) :: parser
    type(modflowinputtype), intent(in) :: mf6_input
    type(ncpackagevarstype), pointer, intent(in) :: nc_vars
    character(len=*), intent(in) :: filename
    integer(I4B), intent(in) :: iout
    integer(I4B) :: iblk
 
    ! initialize static load
    call this%init(parser, mf6_input, filename, iout)
 
    ! set netcdf vars
    this%nc_vars => nc_vars
 
    ! process blocks
    do iblk = 1, size(this%mf6_input%block_dfns)
      ! don't load dynamic input data
      if (this%mf6_input%block_dfns(iblk)%blockname == 'PERIOD') exit
      ! load the block
      call this%load_block(iblk)
    end do
 
    ! finalize static load
    call this%finalize()
  end subroutine load
 
  !> @brief init
  !!
  !! init / finalize are only used when load_block() will be called
  !!
  !<
  subroutine init(this, parser, mf6_input, filename, iout)
    class(loadmf6filetype) :: this
    type(blockparsertype), target, intent(inout) :: parser
    type(modflowinputtype), intent(in) :: mf6_input
    character(len=*), intent(in) :: filename
    integer(I4B), intent(in) :: iout
    integer(I4B) :: isize
 
    this%parser => parser
    this%mf6_input = mf6_input
    this%filename = filename
    this%ts_active = .false.
    this%export = .false.
    this%readasarrays = .false.
    this%readarraygrid = .false.
    this%inamedbound = 0
    this%iauxiliary = 0
    this%iout = iout
 
    call get_isize('MODEL_SHAPE', mf6_input%component_mempath, isize)
    if (isize > 0) then
      call mem_setptr(this%mshape, 'MODEL_SHAPE', mf6_input%component_mempath)
    end if
 
    ! init ts stuctarray list
    allocate (this%ts_sas(0))
 
    ! log lst file header
    call idm_log_header(this%mf6_input%component_name, &
                        this%mf6_input%subcomponent_name, this%iout)
  end subroutine init
 
  !> @brief load a single block
  !!
  !! Assumed in order load of single (next) block. If a
  !! StructArray object is allocated to load this block
  !! it persists until this routine (or finalize) is
  !! called again.
  !!
  !<
  subroutine load_block(this, iblk)
    use structarraymodule, only: destructstructarray
    class(loadmf6filetype) :: this
    integer(I4B), intent(in) :: iblk
 
    ! reset structarray if it was created for previous block
    if (associated(this%structarray)) then
      ! destroy the structured array reader
      call destructstructarray(this%structarray)
    end if
 
    allocate (this%block_tags(0))
    ! load the block
    call this%parse_block(iblk, .false.)
    ! post process block
    call this%block_post_process(iblk)
    ! cleanup
    deallocate (this%block_tags)
  end subroutine load_block
 
  !> @brief finalize
  !!
  !! init / finalize are only used when load_block() will be called
  !!
  !<
  subroutine finalize(this)
    use structarraymodule, only: destructstructarray
    class(loadmf6filetype) :: this
    ! cleanup
    if (associated(this%structarray)) then
      ! destroy the structured array reader
      call destructstructarray(this%structarray)
    end if
    ! close logging block
    call idm_log_close(this%mf6_input%component_name, &
                       this%mf6_input%subcomponent_name, this%iout)
  end subroutine finalize
 
  !> @brief Post parse block handling
  !!
  !<
  subroutine block_post_process(this, iblk)
    use sourcecommonmodule, only: set_model_shape
    class(loadmf6filetype) :: this
    integer(I4B), intent(in) :: iblk
    type(inputparamdefinitiontype), pointer :: idt
    integer(I4B) :: iparam
    integer(I4B), pointer :: intptr
 
    ! update state based on read tags
    do iparam = 1, size(this%block_tags)
      select case (this%mf6_input%block_dfns(iblk)%blockname)
      case ('OPTIONS')
        if (this%block_tags(iparam) == 'AUXILIARY') then
          this%iauxiliary = 1
        else if (this%block_tags(iparam) == 'BOUNDNAMES') then
          this%inamedbound = 1
        else if (this%block_tags(iparam) == 'READASARRAYS') then
          this%readasarrays = .true.
        else if (this%block_tags(iparam) == 'READARRAYGRID') then
          this%readarraygrid = .true.
        else if (this%block_tags(iparam) == 'TS6_FILENAME') then
          this%ts_active = .true.
        else if (this%block_tags(iparam) == 'EXPORT_ARRAY_ASCII') then
          this%export = .true.
        end if
      case default
      end select
    end do
 
    ! update input context allocations based on dfn set and input
    select case (this%mf6_input%block_dfns(iblk)%blockname)
    case ('OPTIONS')
      ! allocate naux and set to 0 if not allocated
      do iparam = 1, size(this%mf6_input%param_dfns)
        idt => this%mf6_input%param_dfns(iparam)
        if (idt%blockname == 'OPTIONS' .and. &
            idt%tagname == 'AUXILIARY') then
          if (this%iauxiliary == 0) then
            call mem_allocate(intptr, 'NAUX', this%mf6_input%mempath)
            intptr = 0
          end if
          exit
        end if
      end do
    case ('DIMENSIONS')
      ! set model shape if discretization dimensions have been read
      if (this%mf6_input%pkgtype(1:3) == 'DIS') then
        call set_model_shape(this%mf6_input%pkgtype, this%filename, &
                             this%mf6_input%component_mempath, &
                             this%mf6_input%mempath, this%mshape)
      end if
    case default
    end select
  end subroutine block_post_process
 
  !> @brief parse block
  !!
  !<
  recursive subroutine parse_block(this, iblk, recursive_call)
    use memorytypemodule, only: memorytype
    use memorymanagermodule, only: get_from_memorystore
    class(loadmf6filetype) :: this
    integer(I4B), intent(in) :: iblk
    logical(LGP), intent(in) :: recursive_call !< true if recursive call
    logical(LGP) :: isblockfound
    logical(LGP) :: endofblock
    logical(LGP) :: supportopenclose
    integer(I4B) :: ierr
    logical(LGP) :: found, required
    type(memorytype), pointer :: mt
    character(len=LINELENGTH) :: tag
    type(inputparamdefinitiontype), pointer :: idt
 
    ! disu vertices/cell2d blocks are contingent on NVERT dimension
    if (this%mf6_input%pkgtype == 'DISU6' .or. &
        this%mf6_input%pkgtype == 'DISV1D6' .or. &
        this%mf6_input%pkgtype == 'DISV2D6') then
      if (this%mf6_input%block_dfns(iblk)%blockname == 'VERTICES' .or. &
          this%mf6_input%block_dfns(iblk)%blockname == 'CELL2D') then
        call get_from_memorystore('NVERT', this%mf6_input%mempath, mt, found, &
                                  .false.)
        if (.not. found) return
        if (mt%intsclr == 0) return
      end if
    end if
 
    ! block open/close support
    supportopenclose = (this%mf6_input%block_dfns(iblk)%blockname /= 'GRIDDATA')
 
    ! parser search for block
    required = this%mf6_input%block_dfns(iblk)%required .and. .not. recursive_call
    call this%parser%GetBlock(this%mf6_input%block_dfns(iblk)%blockname, &
                              isblockfound, ierr, &
                              supportopenclose=supportopenclose, &
                              blockrequired=required)
    ! process block
    if (isblockfound) then
      if (this%mf6_input%block_dfns(iblk)%aggregate) then
        ! process block recarray type, set of variable 1d/2d types
        call this%parse_structarray_block(iblk)
      else
        do
          ! process each line in block
          call this%parser%GetNextLine(endofblock)
          if (endofblock) exit
          ! process line as tag(s)
          call this%parser%GetStringCaps(tag)
          idt => get_param_definition_type( &
                 this%mf6_input%param_dfns, &
                 this%mf6_input%component_type, &
                 this%mf6_input%subcomponent_type, &
                 this%mf6_input%block_dfns(iblk)%blockname, &
                 tag, this%filename)
          if (idt%in_record) then
            call this%parse_record_tag(iblk, idt, .false.)
          else
            call this%load_tag(iblk, idt)
          end if
        end do
      end if
    end if
 
    ! recurse if block is reloadable and was just read
    if (this%mf6_input%block_dfns(iblk)%block_variable) then
      if (isblockfound) then
        call this%parse_block(iblk, .true.)
      end if
    end if
  end subroutine parse_block
 
  subroutine parse_io_tag(this, iblk, pkgtype, which, tag)
    use arrayhandlersmodule, only: expandarray
    class(loadmf6filetype) :: this
    integer(I4B), intent(in) :: iblk
    character(len=*), intent(in) :: pkgtype
    character(len=*), intent(in) :: which
    character(len=*), intent(in) :: tag
    type(inputparamdefinitiontype), pointer :: idt !< input data type object describing this record
    ! matches, read and load file name
    idt => &
      get_param_definition_type(this%mf6_input%param_dfns, &
                                this%mf6_input%component_type, &
                                this%mf6_input%subcomponent_type, &
                                this%mf6_input%block_dfns(iblk)%blockname, &
                                tag, this%filename)
    ! load io tag
    call load_io_tag(this%parser, idt, this%mf6_input%mempath, which, this%iout)
    call expandarray(this%block_tags)
    this%block_tags(size(this%block_tags)) = trim(idt%tagname)
  end subroutine parse_io_tag
 
  recursive subroutine parse_record_tag(this, iblk, inidt, recursive_call)
    use definitionselectmodule, only: split_record_dfn_tag1, &
                                      split_record_dfn_tag2, &
                                      idt_datatype
    class(loadmf6filetype) :: this
    integer(I4B), intent(in) :: iblk
    type(inputparamdefinitiontype), pointer, intent(in) :: inidt
    logical(LGP), intent(in) :: recursive_call !< true if recursive call
    type(inputparamdefinitiontype), pointer :: idt
    character(len=40), dimension(:), allocatable :: words
    integer(I4B) :: n, istart, nwords
    character(len=LINELENGTH) :: tag
 
    nullify (idt)
    istart = 1
 
    if (recursive_call) then
      call split_record_dfn_tag1(this%mf6_input%param_dfns, &
                                 this%mf6_input%component_type, &
                                 this%mf6_input%subcomponent_type, &
                                 inidt%tagname, nwords, words)
      call this%load_tag(iblk, inidt)
      istart = 3
    else
      call this%parser%GetStringCaps(tag)
      if (tag /= '') then
        call split_record_dfn_tag2(this%mf6_input%param_dfns, &
                                   this%mf6_input%component_type, &
                                   this%mf6_input%subcomponent_type, &
                                   inidt%tagname, tag, nwords, words)
        if (nwords == 4 .and. &
            (tag == 'FILEIN' .or. &
             tag == 'FILEOUT')) then
          call this%parse_io_tag(iblk, words(2), words(3), words(4))
          nwords = 0
        else
          idt => get_param_definition_type( &
                 this%mf6_input%param_dfns, &
                 this%mf6_input%component_type, &
                 this%mf6_input%subcomponent_type, &
                 this%mf6_input%block_dfns(iblk)%blockname, &
                 tag, this%filename)
          ! avoid namespace collisions (CIM)
          if (tag /= 'PRINT_FORMAT') call this%load_tag(iblk, inidt)
          call this%load_tag(iblk, idt)
          istart = 4
        end if
      else
        call this%load_tag(iblk, inidt)
        nwords = 0
      end if
    end if
 
    if (istart > 1 .and. nwords == 0) then
      write (errmsg, '(5a)') &
        '"', trim(this%mf6_input%block_dfns(iblk)%blockname), &
        '" block input record that includes keyword "', trim(inidt%tagname), &
        '" is not properly formed.'
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
 
    do n = istart, nwords
      idt => get_param_definition_type( &
             this%mf6_input%param_dfns, &
             this%mf6_input%component_type, &
             this%mf6_input%subcomponent_type, &
             this%mf6_input%block_dfns(iblk)%blockname, &
             words(n), this%filename)
      if (idt_datatype(idt) == 'RECORD') then
        call this%parser%GetStringCaps(tag)
        idt => get_param_definition_type( &
               this%mf6_input%param_dfns, &
               this%mf6_input%component_type, &
               this%mf6_input%subcomponent_type, &
               this%mf6_input%block_dfns(iblk)%blockname, &
               tag, this%filename)
        call this%parse_record_tag(iblk, idt, .true.)
        exit
      else
        if (idt%tagname /= 'FORMAT') then
          call this%parser%GetStringCaps(tag)
          if (tag == '') then
            exit
          else if (idt%tagname /= tag) then
            write (errmsg, '(5a)') 'Expecting record input tag "', &
              trim(idt%tagname), '" but instead found "', trim(tag), '".'
            call store_error(errmsg)
            call this%parser%StoreErrorUnit()
          end if
        end if
        call this%load_tag(iblk, idt)
      end if
    end do
 
    if (allocated(words)) deallocate (words)
  end subroutine parse_record_tag
 
  !> @brief load input keyword
  !! Load input associated with tag key into the memory manager.
  !<
  subroutine load_tag(this, iblk, idt)
    use featureflagsmodule, only: developmode
    use arrayhandlersmodule, only: expandarray
    class(loadmf6filetype) :: this
    integer(I4B), intent(in) :: iblk
    type(inputparamdefinitiontype), pointer, intent(in) :: idt !< input data type object describing this record
    character(len=LINELENGTH) :: dev_msg
 
    ! check if input param is developmode
    if (idt%developmode) then
      dev_msg = 'Input tag "'//trim(idt%tagname)// &
        &'" read from file "'//trim(this%filename)// &
        &'" is still under development. Install the &
        &nightly build or compile from source with IDEVELOPMODE = 1.'
      call developmode(dev_msg, this%iout)
    end if
 
    ! allocate and load data type
    select case (idt%datatype)
    case ('KEYWORD')
      call load_keyword_type(this%parser, idt, this%mf6_input%mempath, this%iout)
      ! check/set as dev option
      if (idt%tagname(1:4) == 'DEV_' .and. &
          this%mf6_input%block_dfns(iblk)%blockname == 'OPTIONS') then
        call this%parser%DevOpt()
      end if
    case ('STRING')
      if (idt%shape == 'NAUX') then
        call load_auxvar_names(this%parser, idt, this%mf6_input%mempath, &
                               this%iout)
      else
        call load_string_type(this%parser, idt, this%mf6_input%mempath, this%iout)
      end if
    case ('INTEGER')
      call load_integer_type(this%parser, idt, this%mf6_input%mempath, this%iout)
    case ('INTEGER1D')
      call load_integer1d_type(this%parser, idt, this%mf6_input, this%mshape, &
                               this%export, this%nc_vars, this%filename, &
                               this%iout)
    case ('INTEGER2D')
      call load_integer2d_type(this%parser, idt, this%mf6_input, this%mshape, &
                               this%export, this%nc_vars, this%filename, &
                               this%iout)
    case ('INTEGER3D')
      call load_integer3d_type(this%parser, idt, this%mf6_input, this%mshape, &
                               this%export, this%nc_vars, this%filename, &
                               this%iout)
    case ('DOUBLE')
      call load_double_type(this%parser, idt, this%mf6_input%mempath, this%iout)
    case ('DOUBLE1D')
      call load_double1d_type(this%parser, idt, this%mf6_input, this%mshape, &
                              this%export, this%nc_vars, this%filename, this%iout)
    case ('DOUBLE2D')
      call load_double2d_type(this%parser, idt, this%mf6_input, this%mshape, &
                              this%export, this%nc_vars, this%filename, this%iout)
    case ('DOUBLE3D')
      call load_double3d_type(this%parser, idt, this%mf6_input, this%mshape, &
                              this%export, this%nc_vars, this%filename, this%iout)
    case default
      write (errmsg, '(a,a)') 'Failure reading data for tag: ', trim(idt%tagname)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end select
 
    call expandarray(this%block_tags)
    this%block_tags(size(this%block_tags)) = trim(idt%tagname)
  end subroutine load_tag
 
  function block_index_dfn(this, iblk) result(idt)
    class(loadmf6filetype) :: this
    integer(I4B), intent(in) :: iblk
    type(inputparamdefinitiontype) :: idt !< input data type object describing this record
    character(len=LENVARNAME) :: varname
    integer(I4B) :: ilen
    character(len=3) :: block_suffix = 'NUM'
 
    ! assign first column as the block number
    ilen = len_trim(this%mf6_input%block_dfns(iblk)%blockname)
 
    if (ilen > (lenvarname - len(block_suffix))) then
      varname = &
        this%mf6_input%block_dfns(iblk)% &
        blockname(1:(lenvarname - len(block_suffix)))//block_suffix
    else
      varname = trim(this%mf6_input%block_dfns(iblk)%blockname)//block_suffix
    end if
 
    idt%component_type = trim(this%mf6_input%component_type)
    idt%subcomponent_type = trim(this%mf6_input%subcomponent_type)
    idt%blockname = trim(this%mf6_input%block_dfns(iblk)%blockname)
    idt%tagname = varname
    idt%mf6varname = varname
    idt%datatype = 'INTEGER'
  end function block_index_dfn
 
  !> @brief parse a structured array record into memory manager
  !!
  !! A structarray is similar to a numpy recarray.  It it used to
  !! load a list of data in which each column in the list may be a
  !! different type.  Each column in the list is stored as a 1d
  !! vector.
  !!
  !<
  subroutine parse_structarray_block(this, iblk)
    use structarraymodule, only: structarraytype, constructstructarray
    use loadcontextmodule, only: loadcontexttype, is_feature_keystring
    class(loadmf6filetype) :: this
    integer(I4B), intent(in) :: iblk
    type(loadcontexttype) :: ctx
    character(len=LINELENGTH), dimension(:), allocatable :: param_names
    type(inputparamdefinitiontype), pointer :: idt !< input data type object describing this record
    type(inputparamdefinitiontype), pointer :: shape_idt
    type(inputparamdefinitiontype), target :: blockvar_idt
    integer(I4B) :: blocknum
    integer(I4B), pointer :: nrow
    integer(I4B), dimension(:), pointer, contiguous :: int1d
    integer(I4B) :: nrows, nrowsread
    integer(I4B) :: mem_rank, isize
    integer(I4B) :: ibinary, oc_inunit
    integer(I4B) :: icol, iparam
    integer(I4B) :: ncol, nparam
    logical(LGP) :: shape_found
    logical(LGP) :: is_sum_shape
    character(len=LINELENGTH) :: shape_name
    integer(I4B), pointer :: pkgdata_maxbound
 
    nrows = -1
 
    ! initialize load context
    call ctx%init(this%mf6_input, blockname= &
                  this%mf6_input%block_dfns(iblk)%blockname)
    ! set in-scope params directly from context
    param_names = ctx%params
    nparam = size(ctx%params)
    call ctx%check_developmode(this%filename)
    ! set input definition for this block
    idt => &
      get_aggregate_definition_type(this%mf6_input%aggregate_dfns, &
                                    this%mf6_input%component_type, &
                                    this%mf6_input%subcomponent_type, &
                                    this%mf6_input%block_dfns(iblk)%blockname)
    ! if block is reloadable read the block number
    if (this%mf6_input%block_dfns(iblk)%block_variable) then
      blocknum = this%parser%GetInteger()
    else
      blocknum = 0
    end if
 
    ! set ncol
    ncol = nparam
    ! add col if block is reloadable
    if (blocknum > 0) ncol = ncol + 1
    ! use shape to set the max num of rows
    if (idt%shape /= '') then
      ! "SUM(name)" means total rows = sum of a per-feature array (e.g.
      ! LAK's CONNECTIONDATA: NLAKECONN summed across all lakes); this
      ! notation is declared explicitly in the dfn, not inferred
      is_sum_shape = .false.
      shape_name = idt%shape
      if (idt%shape(1:4) == 'SUM(') then
        is_sum_shape = .true.
        shape_name = idt%shape(5:len_trim(idt%shape) - 1)
      end if
      shape_idt => &
        get_param_definition_type(this%mf6_input%param_dfns, &
                                  this%mf6_input%component_type, &
                                  this%mf6_input%subcomponent_type, &
                                  'DIMENSIONS', shape_name, this%filename, &
                                  found=shape_found)
      if (.not. shape_found) then
        ! also allow a per-row PACKAGEDATA field (e.g. LAK's NLAKECONN)
        shape_idt => &
          get_param_definition_type(this%mf6_input%param_dfns, &
                                    this%mf6_input%component_type, &
                                    this%mf6_input%subcomponent_type, &
                                    'PACKAGEDATA', shape_name, this%filename, &
                                    found=shape_found)
      end if
      if (shape_found) then
        call get_isize(shape_idt%mf6varname, this%mf6_input%mempath, isize)
        if (isize < 0) then
          if (shape_idt%required) then
            write (errmsg, '(3a)') 'Required dimension "', &
              trim(shape_name), '" not found.'
            call store_error(errmsg)
            call this%parser%StoreErrorUnit()
          end if
        else
          ! cross-check the dfn's declared SUM(...)/plain shape against the
          ! variable's actual rank, so a mismatch errors instead of
          ! silently mis-using a scalar/array pointer
          call get_mem_rank(shape_idt%mf6varname, this%mf6_input%mempath, &
                            mem_rank)
          if (is_sum_shape .and. mem_rank == 1) then
            call mem_setptr(int1d, shape_idt%mf6varname, this%mf6_input%mempath)
            nrows = sum(int1d)
            nullify (int1d)
          else if (.not. is_sum_shape .and. mem_rank == 0) then
            ! scalar shape variable (e.g. NLAKES, MAXBOUND) — use value directly
            call mem_setptr(nrow, shape_idt%mf6varname, this%mf6_input%mempath)
            nrows = nrow
          else
            write (errmsg, '(5a)') 'Shape variable "', trim(shape_name), &
              '" rank does not match its declared shape "', &
              trim(idt%shape), '".'
            call store_error(errmsg)
            call this%parser%StoreErrorUnit()
          end if
        end if
      end if
      ! -- else: shape variable not found; nrows stays at its initial -1, so
      !    the block falls back to deferred sizing
    end if
 
    ! create a structured array; use a larger deferred init for blocks with no
    ! explicit shape, which include APT based advanced packages.
    if (nrows < 0) then
      this%structarray => constructstructarray(this%mf6_input, ncol, nrows, &
                                               blocknum, this%mf6_input%mempath, &
                                               this%mf6_input%component_mempath, &
                                               size_init=64)
    else
      this%structarray => constructstructarray(this%mf6_input, ncol, nrows, &
                                               blocknum, this%mf6_input%mempath, &
                                               this%mf6_input%component_mempath)
    end if
    ! create structarray vectors for each column
    do icol = 1, ncol
      ! if block is reloadable, block number is first column
      if (blocknum > 0) then
        if (icol == 1) then
          blockvar_idt = this%block_index_dfn(iblk)
          idt => blockvar_idt
          call this%structarray%mem_create_vector(icol, idt)
          ! continue as this column managed by internally SA object
          cycle
        end if
        ! set indexes (where first column is blocknum)
        iparam = icol - 1
      else
        ! set indexes (no blocknum column)
        iparam = icol
      end if
      ! set pointer to input definition for this 1d vector
      idt => &
        get_param_definition_type(this%mf6_input%param_dfns, &
                                  this%mf6_input%component_type, &
                                  this%mf6_input%subcomponent_type, &
                                  this%mf6_input%block_dfns(iblk)%blockname, &
                                  param_names(iparam), this%filename)
      ! allocate variable in memory manager
      call this%structarray%mem_create_vector(icol, idt)
    end do
 
    ! finish context setup after allocating vectors
    call ctx%allocate_arrays()
 
    ! read the block control record
    ibinary = read_control_record(this%parser, oc_inunit, this%iout)
 
    if (ibinary == 1) then
      ! read from binary
      nrowsread = this%structarray%read_from_binary(oc_inunit, this%iout)
      call this%parser%terminateblock()
      close (oc_inunit)
    else
      ! read from ascii
      nrowsread = this%structarray%read_from_parser(this%parser, this%ts_active, &
                                                    this%iout, this%filename)
      ! save structarray for deferred TS linking in df() if any strlocs were stored
      if (this%ts_active) call this%save_ts_sa()
    end if
 
    ! an advanced package's PACKAGEDATA row count is published as
    ! MAXBOUND, the feature count the PERIOD block's load context needs
    if (this%mf6_input%block_dfns(iblk)%blockname == 'PACKAGEDATA' .and. &
        is_feature_keystring(this%mf6_input) .and. ctx%is_advanced) then
      call get_isize('MAXBOUND', this%mf6_input%mempath, isize)
      if (isize < 0) then
        call mem_allocate(pkgdata_maxbound, 'MAXBOUND', this%mf6_input%mempath)
        pkgdata_maxbound = nrowsread
      end if
    end if
 
    ! clean up
    call ctx%destroy()
  end subroutine parse_structarray_block
 
  !> @brief Return number of saved static StructArrays with deferred TS strlocs
  !<
  function ts_sa_count(this) result(n)
    class(loadmf6filetype), intent(in) :: this
    integer(I4B) :: n
    n = size(this%ts_sas)
  end function ts_sa_count
 
  !> @brief Return the n-th saved static StructArray pointer
  !<
  function get_ts_sa(this, n) result(sa)
    class(loadmf6filetype), intent(in) :: this
    integer(I4B), intent(in) :: n
    type(structarraytype), pointer :: sa
    sa => this%ts_sas(n)%sa
  end function get_ts_sa
 
  !> @brief load type keyword
  !<
  subroutine load_keyword_type(parser, idt, memoryPath, iout)
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    character(len=*), intent(in) :: memoryPath !< memorypath to put loaded information
    integer(I4B), intent(in) :: iout !< unit number for output
    integer(I4B), pointer :: intvar
    call mem_allocate(intvar, idt%mf6varname, memorypath)
    intvar = 1
    call idm_log_var(intvar, idt%tagname, memorypath, idt%datatype, iout)
  end subroutine load_keyword_type
 
  !> @brief load type string
  !<
  subroutine load_string_type(parser, idt, memoryPath, iout)
    use constantsmodule, only: lenbigline
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    character(len=*), intent(in) :: memoryPath !< memorypath to put loaded information
    integer(I4B), intent(in) :: iout !< unit number for output
    character(len=LINELENGTH), pointer :: cstr
    character(len=LENBIGLINE), pointer :: bigcstr
    integer(I4B) :: ilen
    select case (idt%shape)
    case ('LENBIGLINE')
      ilen = lenbigline
      call mem_allocate(bigcstr, ilen, idt%mf6varname, memorypath)
      call parser%GetString(bigcstr, (.not. idt%preserve_case))
      call idm_log_var(bigcstr, idt%tagname, memorypath, iout)
    case default
      ilen = linelength
      call mem_allocate(cstr, ilen, idt%mf6varname, memorypath)
      call parser%GetString(cstr, (.not. idt%preserve_case))
      call idm_log_var(cstr, idt%tagname, memorypath, iout)
    end select
  end subroutine load_string_type
 
  !> @brief load io tag
  !<
  subroutine load_io_tag(parser, idt, memoryPath, which, iout)
    use memorymanagermodule, only: mem_allocate, mem_reallocate
    use characterstringmodule, only: characterstringtype
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    character(len=*), intent(in) :: memoryPath !< memorypath to put loaded information
    character(len=*), intent(in) :: which
    integer(I4B), intent(in) :: iout !< unit number for output
    character(len=LINELENGTH) :: cstr
    type(characterstringtype), dimension(:), pointer, contiguous :: charstr1d
    integer(I4B) :: ilen, isize, idx
    ilen = linelength
    if (which == 'FILEIN') then
      call get_isize(idt%mf6varname, memorypath, isize)
      if (isize < 0) then
        call mem_allocate(charstr1d, ilen, 1, idt%mf6varname, memorypath)
        idx = 1
      else
        call mem_setptr(charstr1d, idt%mf6varname, memorypath)
        call mem_reallocate(charstr1d, ilen, isize + 1, idt%mf6varname, &
                            memorypath)
        idx = isize + 1
      end if
      call parser%GetString(cstr, (.not. idt%preserve_case))
      charstr1d(idx) = cstr
    else if (which == 'FILEOUT') then
      call load_string_type(parser, idt, memorypath, iout)
    end if
  end subroutine load_io_tag
 
  !> @brief load aux variable names
  !!
  !<
  subroutine load_auxvar_names(parser, idt, memoryPath, iout)
    use constantsmodule, only: lenauxname, linelength, lenpackagename
    use inputoutputmodule, only: urdaux
    use characterstringmodule, only: characterstringtype
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    character(len=*), intent(in) :: memoryPath !< memorypath to put loaded information
    integer(I4B), intent(in) :: iout !< unit number for output
    character(len=:), allocatable :: line
    character(len=LENAUXNAME), dimension(:), allocatable :: caux
    integer(I4B) :: lloc
    integer(I4B) :: istart
    integer(I4B) :: istop
    integer(I4B) :: i
    character(len=LENPACKAGENAME) :: text = ''
    integer(I4B), pointer :: intvar
    type(characterstringtype), dimension(:), &
      pointer, contiguous :: acharstr1d !< variable for allocation
    call mem_allocate(intvar, idt%shape, memorypath)
    intvar = 0
    call parser%GetRemainingLine(line)
    lloc = 1
    call urdaux(intvar, parser%iuactive, iout, lloc, &
                istart, istop, caux, line, text)
    call mem_allocate(acharstr1d, lenauxname, intvar, idt%mf6varname, memorypath)
    do i = 1, intvar
      acharstr1d(i) = caux(i)
    end do
    deallocate (line)
    deallocate (caux)
  end subroutine load_auxvar_names
 
  !> @brief load type integer
  !<
  subroutine load_integer_type(parser, idt, memoryPath, iout)
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    character(len=*), intent(in) :: memoryPath !< memorypath to put loaded information
    integer(I4B), intent(in) :: iout !< unit number for output
    integer(I4B), pointer :: intvar
    call mem_allocate(intvar, idt%mf6varname, memorypath)
    intvar = parser%GetInteger()
    call idm_log_var(intvar, idt%tagname, memorypath, idt%datatype, iout)
  end subroutine load_integer_type
 
  !> @brief load type 1d integer
  !<
  subroutine load_integer1d_type(parser, idt, mf6_input, mshape, export, &
                                 nc_vars, input_fname, iout)
    use sourcecommonmodule, only: get_shape_from_string, get_layered_shape
    use loadncinputmodule, only: netcdf_read_array
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    type(modflowinputtype), intent(in) :: mf6_input !< description of input
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape
    logical(LGP), intent(in) :: export !< export to ascii layer files
    type(ncpackagevarstype), pointer, intent(in) :: nc_vars
    character(len=*), intent(in) :: input_fname !< ascii input file name
    integer(I4B), intent(in) :: iout !< unit number for output
    integer(I4B), dimension(:), pointer, contiguous :: int1d
    integer(I4B) :: nlay
    integer(I4B) :: nvals
    integer(I4B), dimension(:), allocatable :: array_shape
    integer(I4B), dimension(:), allocatable :: layer_shape
    character(len=LINELENGTH) :: keyword
 
    ! Check if it is a full grid sized array (NODES), otherwise use
    ! idt%shape to construct shape from variables in memoryPath
    if (idt%shape == 'NODES') then
      nvals = product(mshape)
    else
      call get_shape_from_string(idt%shape, array_shape, mf6_input%mempath)
      nvals = array_shape(1)
    end if
 
    ! allocate memory for the array
    call mem_allocate(int1d, nvals, idt%mf6varname, mf6_input%mempath)
 
    ! read keyword
    keyword = ''
    call parser%GetStringCaps(keyword)
 
    ! check for "NETCDF" and "LAYERED"
    if (keyword == 'NETCDF') then
      call netcdf_read_array(int1d, mshape, idt, mf6_input, nc_vars, &
                             input_fname, iout)
    else if (keyword == 'LAYERED' .and. idt%layered) then
      call get_layered_shape(mshape, nlay, layer_shape)
      call read_int1d_layered(parser, int1d, idt%mf6varname, nlay, layer_shape)
    else
      call read_int1d(parser, int1d, idt%mf6varname)
    end if
 
    ! log information on the loaded array to the list file
    call idm_log_var(int1d, idt%tagname, mf6_input%mempath, iout)
 
    ! create export file for griddata parameters if optioned
    if (export) then
      if (idt%blockname == 'GRIDDATA') then
        call idm_export(int1d, idt%tagname, mf6_input%mempath, idt%shape, iout)
      end if
    end if
  end subroutine load_integer1d_type
 
  !> @brief load type 2d integer
  !<
  subroutine load_integer2d_type(parser, idt, mf6_input, mshape, export, &
                                 nc_vars, input_fname, iout)
    use sourcecommonmodule, only: get_shape_from_string, get_layered_shape
    use loadncinputmodule, only: netcdf_read_array
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    type(modflowinputtype), intent(in) :: mf6_input !< description of input
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape
    logical(LGP), intent(in) :: export !< export to ascii layer files
    type(ncpackagevarstype), pointer, intent(in) :: nc_vars
    character(len=*), intent(in) :: input_fname !< ascii input file name
    integer(I4B), intent(in) :: iout !< unit number for output
    integer(I4B), dimension(:, :), pointer, contiguous :: int2d
    integer(I4B) :: nlay
    integer(I4B) :: nsize1, nsize2
    integer(I4B), dimension(:), allocatable :: array_shape
    integer(I4B), dimension(:), allocatable :: layer_shape
    character(len=LINELENGTH) :: keyword
 
    ! determine the array shape from the input data definition (idt%shape),
    ! which looks like "NCOL, NROW, NLAY"
    call get_shape_from_string(idt%shape, array_shape, mf6_input%mempath)
    nsize1 = array_shape(1)
    nsize2 = array_shape(2)
 
    ! create a new 3d memory managed variable
    call mem_allocate(int2d, nsize1, nsize2, idt%mf6varname, mf6_input%mempath)
 
    ! read keyword
    keyword = ''
    call parser%GetStringCaps(keyword)
 
    ! check for "NETCDF" and "LAYERED"
    if (keyword == 'NETCDF') then
      call netcdf_read_array(int2d, mshape, idt, mf6_input, nc_vars, &
                             input_fname, iout)
    else if (keyword == 'LAYERED' .and. idt%layered) then
      call get_layered_shape(mshape, nlay, layer_shape)
      call read_int2d_layered(parser, int2d, idt%mf6varname, nlay, layer_shape)
    else
      call read_int2d(parser, int2d, idt%mf6varname)
    end if
 
    ! log information on the loaded array to the list file
    call idm_log_var(int2d, idt%tagname, mf6_input%mempath, iout)
 
    ! create export file for griddata parameters if optioned
    if (export) then
      if (idt%blockname == 'GRIDDATA') then
        call idm_export(int2d, idt%tagname, mf6_input%mempath, idt%shape, iout)
      end if
    end if
  end subroutine load_integer2d_type
 
  !> @brief load type 3d integer
  !<
  subroutine load_integer3d_type(parser, idt, mf6_input, mshape, export, &
                                 nc_vars, input_fname, iout)
    use sourcecommonmodule, only: get_shape_from_string, get_layered_shape
    use loadncinputmodule, only: netcdf_read_array
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    type(modflowinputtype), intent(in) :: mf6_input !< description of input
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape
    logical(LGP), intent(in) :: export !< export to ascii layer files
    type(ncpackagevarstype), pointer, intent(in) :: nc_vars
    character(len=*), intent(in) :: input_fname !< ascii input file name
    integer(I4B), intent(in) :: iout !< unit number for output
    integer(I4B), dimension(:, :, :), pointer, contiguous :: int3d
    integer(I4B) :: nlay
    integer(I4B) :: nsize1, nsize2, nsize3
    integer(I4B), dimension(:), allocatable :: array_shape
    integer(I4B), dimension(:), allocatable :: layer_shape
    integer(I4B), dimension(:), pointer, contiguous :: int1d_ptr
    character(len=LINELENGTH) :: keyword
 
    ! determine the array shape from the input data definition (idt%shape),
    ! which looks like "NCOL, NROW, NLAY"
    call get_shape_from_string(idt%shape, array_shape, mf6_input%mempath)
    nsize1 = array_shape(1)
    nsize2 = array_shape(2)
    nsize3 = array_shape(3)
 
    ! create a new 3d memory managed variable
    call mem_allocate(int3d, nsize1, nsize2, nsize3, idt%mf6varname, &
                      mf6_input%mempath)
 
    ! read keyword
    keyword = ''
    call parser%GetStringCaps(keyword)
 
    ! check for "NETCDF" and "LAYERED"
    if (keyword == 'NETCDF') then
      call netcdf_read_array(int3d, mshape, idt, mf6_input, nc_vars, &
                             input_fname, iout)
    else if (keyword == 'LAYERED' .and. idt%layered) then
      call get_layered_shape(mshape, nlay, layer_shape)
      call read_int3d_layered(parser, int3d, idt%mf6varname, nlay, &
                              layer_shape)
    else
      int1d_ptr(1:nsize1 * nsize2 * nsize3) => int3d(:, :, :)
      call read_int1d(parser, int1d_ptr, idt%mf6varname)
    end if
 
    ! log information on the loaded array to the list file
    call idm_log_var(int3d, idt%tagname, mf6_input%mempath, iout)
 
    ! create export file for griddata parameters if optioned
    if (export) then
      if (idt%blockname == 'GRIDDATA') then
        call idm_export(int3d, idt%tagname, mf6_input%mempath, idt%shape, iout)
      end if
    end if
  end subroutine load_integer3d_type
 
  !> @brief load type double
  !<
  subroutine load_double_type(parser, idt, memoryPath, iout)
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    character(len=*), intent(in) :: memoryPath !< memorypath to put loaded information
    integer(I4B), intent(in) :: iout !< unit number for output
    real(DP), pointer :: dblvar
    call mem_allocate(dblvar, idt%mf6varname, memorypath)
    dblvar = parser%GetDouble()
    call idm_log_var(dblvar, idt%tagname, memorypath, iout)
  end subroutine load_double_type
 
  !> @brief load type 1d double
  !<
  subroutine load_double1d_type(parser, idt, mf6_input, mshape, export, &
                                nc_vars, input_fname, iout)
    use sourcecommonmodule, only: get_shape_from_string, get_layered_shape
    use loadncinputmodule, only: netcdf_read_array
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    type(modflowinputtype), intent(in) :: mf6_input !< description of input
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape
    logical(LGP), intent(in) :: export !< export to ascii layer files
    type(ncpackagevarstype), pointer, intent(in) :: nc_vars
    character(len=*), intent(in) :: input_fname !< ascii input file name
    integer(I4B), intent(in) :: iout !< unit number for output
    real(DP), dimension(:), pointer, contiguous :: dbl1d
    integer(I4B) :: nlay
    integer(I4B) :: nvals
    integer(I4B), dimension(:), allocatable :: array_shape
    integer(I4B), dimension(:), allocatable :: layer_shape
    character(len=LINELENGTH) :: keyword
 
    ! Check if it is a full grid sized array (NODES)
    if (idt%shape == 'NODES') then
      nvals = product(mshape)
    else
      call get_shape_from_string(idt%shape, array_shape, mf6_input%mempath)
      nvals = array_shape(1)
    end if
 
    ! allocate memory for the array
    call mem_allocate(dbl1d, nvals, idt%mf6varname, mf6_input%mempath)
 
    ! read keyword
    keyword = ''
    call parser%GetStringCaps(keyword)
 
    ! check for "NETCDF" and "LAYERED"
    if (keyword == 'NETCDF') then
      call netcdf_read_array(dbl1d, mshape, idt, mf6_input, nc_vars, &
                             input_fname, iout)
    else if (keyword == 'LAYERED' .and. idt%layered) then
      call get_layered_shape(mshape, nlay, layer_shape)
      call read_dbl1d_layered(parser, dbl1d, idt%mf6varname, nlay, layer_shape)
    else
      call read_dbl1d(parser, dbl1d, idt%mf6varname)
    end if
 
    ! log information on the loaded array to the list file
    call idm_log_var(dbl1d, idt%tagname, mf6_input%mempath, iout)
 
    ! create export file for griddata parameters if optioned
    if (export) then
      if (idt%blockname == 'GRIDDATA') then
        call idm_export(dbl1d, idt%tagname, mf6_input%mempath, idt%shape, iout)
      end if
    end if
  end subroutine load_double1d_type
 
  !> @brief load type 2d double
  !<
  subroutine load_double2d_type(parser, idt, mf6_input, mshape, export, &
                                nc_vars, input_fname, iout)
    use sourcecommonmodule, only: get_shape_from_string, get_layered_shape
    use loadncinputmodule, only: netcdf_read_array
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    type(modflowinputtype), intent(in) :: mf6_input !< description of input
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape
    logical(LGP), intent(in) :: export !< export to ascii layer files
    type(ncpackagevarstype), pointer, intent(in) :: nc_vars
    character(len=*), intent(in) :: input_fname !< ascii input file name
    integer(I4B), intent(in) :: iout !< unit number for output
    real(DP), dimension(:, :), pointer, contiguous :: dbl2d
    integer(I4B) :: nlay
    integer(I4B) :: nsize1, nsize2
    integer(I4B), dimension(:), allocatable :: array_shape
    integer(I4B), dimension(:), allocatable :: layer_shape
    character(len=LINELENGTH) :: keyword
 
    ! determine the array shape from the input data definition (idt%shape),
    ! which looks like "NCOL, NROW, NLAY"
    call get_shape_from_string(idt%shape, array_shape, mf6_input%mempath)
    nsize1 = array_shape(1)
    nsize2 = array_shape(2)
 
    ! create a new 3d memory managed variable
    call mem_allocate(dbl2d, nsize1, nsize2, idt%mf6varname, mf6_input%mempath)
 
    ! read keyword
    keyword = ''
    call parser%GetStringCaps(keyword)
 
    ! check for "NETCDF" and "LAYERED"
    if (keyword == 'NETCDF') then
      call netcdf_read_array(dbl2d, mshape, idt, mf6_input, nc_vars, &
                             input_fname, iout)
    else if (keyword == 'LAYERED' .and. idt%layered) then
      call get_layered_shape(mshape, nlay, layer_shape)
      call read_dbl2d_layered(parser, dbl2d, idt%mf6varname, nlay, layer_shape)
    else
      call read_dbl2d(parser, dbl2d, idt%mf6varname)
    end if
 
    ! log information on the loaded array to the list file
    call idm_log_var(dbl2d, idt%tagname, mf6_input%mempath, iout)
 
    ! create export file for griddata parameters if optioned
    if (export) then
      if (idt%blockname == 'GRIDDATA') then
        call idm_export(dbl2d, idt%tagname, mf6_input%mempath, idt%shape, iout)
      end if
    end if
  end subroutine load_double2d_type
 
  !> @brief load type 3d double
  !<
  subroutine load_double3d_type(parser, idt, mf6_input, mshape, export, &
                                nc_vars, input_fname, iout)
    use sourcecommonmodule, only: get_shape_from_string, get_layered_shape
    use loadncinputmodule, only: netcdf_read_array
    type(blockparsertype), intent(inout) :: parser !< block parser
    type(inputparamdefinitiontype), intent(in) :: idt !< input data type object describing this record
    type(modflowinputtype), intent(in) :: mf6_input !< description of input
    integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape
    logical(LGP), intent(in) :: export !< export to ascii layer files
    type(ncpackagevarstype), pointer, intent(in) :: nc_vars
    character(len=*), intent(in) :: input_fname !< ascii input file name
    integer(I4B), intent(in) :: iout !< unit number for output
    real(DP), dimension(:, :, :), pointer, contiguous :: dbl3d
    integer(I4B) :: nlay
    integer(I4B) :: nsize1, nsize2, nsize3
    integer(I4B), dimension(:), allocatable :: array_shape
    integer(I4B), dimension(:), allocatable :: layer_shape
    real(DP), dimension(:), pointer, contiguous :: dbl1d_ptr
    character(len=LINELENGTH) :: keyword
 
    ! determine the array shape from the input data definition (idt%shape),
    ! which looks like "NCOL, NROW, NLAY"
    call get_shape_from_string(idt%shape, array_shape, mf6_input%mempath)
    nsize1 = array_shape(1)
    nsize2 = array_shape(2)
    nsize3 = array_shape(3)
 
    ! create a new 3d memory managed variable
    call mem_allocate(dbl3d, nsize1, nsize2, nsize3, idt%mf6varname, &
                      mf6_input%mempath)
 
    ! read keyword
    keyword = ''
    call parser%GetStringCaps(keyword)
 
    ! check for "NETCDF" and "LAYERED"
    if (keyword == 'NETCDF') then
      call netcdf_read_array(dbl3d, mshape, idt, mf6_input, nc_vars, &
                             input_fname, iout)
    else if (keyword == 'LAYERED' .and. idt%layered) then
      call get_layered_shape(mshape, nlay, layer_shape)
      call read_dbl3d_layered(parser, dbl3d, idt%mf6varname, nlay, &
                              layer_shape)
    else
      dbl1d_ptr(1:nsize1 * nsize2 * nsize3) => dbl3d(:, :, :)
      call read_dbl1d(parser, dbl1d_ptr, idt%mf6varname)
    end if
 
    ! log information on the loaded array to the list file
    call idm_log_var(dbl3d, idt%tagname, mf6_input%mempath, iout)
 
    ! create export file for griddata parameters if optioned
    if (export) then
      if (idt%blockname == 'GRIDDATA') then
        call idm_export(dbl3d, idt%tagname, mf6_input%mempath, idt%shape, iout)
      end if
    end if
  end subroutine load_double3d_type
 
  function read_control_record(parser, oc_inunit, iout) result(ibinary)
    use simmodule, only: store_error_unit
    use inputoutputmodule, only: urword
    use inputoutputmodule, only: openfile
    use openspecmodule, only: form, access
    use constantsmodule, only: linelength
    use blockparsermodule, only: blockparsertype
    type(blockparsertype), intent(inout) :: parser
    integer(I4B), intent(inout) :: oc_inunit
    integer(I4B), intent(in) :: iout
    integer(I4B) :: ibinary
    integer(I4B) :: lloc, istart, istop, idum, inunit, itmp, ierr
    integer(I4B) :: nunopn = 99
    character(len=:), allocatable :: line
    character(len=LINELENGTH) :: fname
    logical(LGP) :: exists
    real(dp) :: r
    character(len=*), parameter :: fmtocne = &
      &"('Specified OPEN/CLOSE file ',(A),' does not exist')"
    character(len=*), parameter :: fmtobf = &
      &"(1X,/1X,'OPENING BINARY FILE ON UNIT ',I0,':',/1X,A)"
 
    ! initialize oc_inunit and ibinary
    oc_inunit = 0
    ibinary = 0
    inunit = parser%getunit()
 
    ! Read to the first non-commented line
    lloc = 1
    call parser%line_reader%rdcom(inunit, iout, line, ierr)
    call urword(line, lloc, istart, istop, 1, idum, r, iout, inunit)
 
    if (line(istart:istop) == 'OPEN/CLOSE') then
      ! get filename
      call urword(line, lloc, istart, istop, 0, idum, r, &
                  iout, inunit)
      fname = line(istart:istop)
      ! check to see if file OPEN/CLOSE file exists
      inquire (file=fname, exist=exists)
      if (.not. exists) then
        write (errmsg, fmtocne) line(istart:istop)
        call store_error(errmsg)
        call store_error('Specified OPEN/CLOSE file does not exist')
        call store_error_unit(inunit)
      end if
 
      ! Check for (BINARY) keyword
      call urword(line, lloc, istart, istop, 1, idum, r, &
                  iout, inunit)
 
      if (line(istart:istop) == '(BINARY)') ibinary = 1
      ! Open the file depending on ibinary flag
      if (ibinary == 1) then
        oc_inunit = nunopn
        itmp = iout
        if (iout > 0) then
          itmp = 0
          write (iout, fmtobf) oc_inunit, trim(adjustl(fname))
        end if
        call openfile(oc_inunit, itmp, fname, 'OPEN/CLOSE', &
                      fmtarg_opt=form, accarg_opt=access)
      end if
    end if
 
    if (ibinary == 0) then
      call parser%line_reader%bkspc(parser%getunit())
    end if
  end function read_control_record
 
  !> @brief Save structarray pointer for deferred TS linking
  !!
  !! Saves the current structarray pointer when it contains TS string locs,
  !! then nullifies the pointer so load_block does not destroy it.
  !<
  subroutine save_ts_sa(this)
    class(loadmf6filetype), intent(inout) :: this
    type(structvectortype), pointer :: svect
    type(staticsatype), allocatable :: tmp(:)
    logical(LGP) :: has_ts
    integer(I4B) :: m, n
 
    ! check if any column has deferred TS strlocs
    has_ts = .false.
    do m = 1, this%structarray%count()
      svect => this%structarray%get(m)
      if (svect%idt%timeseries .and. svect%ts_strlocs%count() > 0) then
        has_ts = .true.
        exit
      end if
    end do
 
    if (has_ts) then
      n = size(this%ts_sas)
      allocate (tmp(n + 1))
      tmp(1:n) = this%ts_sas
      tmp(n + 1)%sa => this%structarray
      call move_alloc(tmp, this%ts_sas)
      ! nullify so load_block does not destroy the saved SA
      nullify (this%structarray)
    end if
  end subroutine save_ts_sa
 
  !> @brief Clean up saved static structarrays
  !!
  !! Called from destroy() of dynamic loaders.
  !<
  subroutine cleanup(this)
    use structarraymodule, only: destructstructarray
    class(loadmf6filetype), intent(inout) :: this
    integer(I4B) :: n
    do n = 1, size(this%ts_sas)
      if (associated(this%ts_sas(n)%sa)) then
        call destructstructarray(this%ts_sas(n)%sa)
        nullify (this%ts_sas(n)%sa)
      end if
    end do
    deallocate (this%ts_sas)
  end subroutine cleanup
 
end module loadmf6filemodule