Mf6FileGridArray.f90

Go to the documentation of this file.

!> @brief This module contains the GridArrayLoadModule
!!
!! This module contains the routines for reading period block
!! grid array-based input for stress packages that use the
!! READARRAYGRID option (CHD, WEL, DRN, RIV, GHB).
!!
!<
module gridarrayloadmodule
 
  use kindmodule, only: i4b, dp, lgp
  use constantsmodule, only: linelength, lenvarname
  use simvariablesmodule, only: errmsg
  use simmodule, only: store_error, store_error_filename
  use inputdefinitionmodule, only: inputparamdefinitiontype
  use memorymanagermodule, only: mem_allocate, mem_setptr
  use blockparsermodule, only: blockparsertype
  use modflowinputmodule, only: modflowinputtype
  use loadcontextmodule, only: loadcontexttype, readstatevartype
  use asciiinputloadtypemodule, only: asciidynamicpkgloadbasetype
 
  implicit none
  private
  public :: gridarrayloadtype
 
  !> @brief Ascii grid based dynamic loader type
  !<
  type, extends(asciidynamicpkgloadbasetype) :: gridarrayloadtype
    type(readstatevartype), dimension(:), allocatable :: param_reads !< read states for current load
    type(loadcontexttype) :: ctx
    integer(I4B), dimension(:), pointer, contiguous :: nodeulist
  contains
    procedure :: ainit
    procedure :: df
    procedure :: ts_advance
    procedure :: rp
    procedure :: destroy
    procedure :: reset
    procedure :: params_alloc
    procedure :: param_load
  end type gridarrayloadtype
 
contains
 
  subroutine ainit(this, mf6_input, component_name, &
                   component_input_name, input_name, &
                   iperblock, parser, iout)
    use memorymanagermodule, only: get_isize
    use blockparsermodule, only: blockparsertype
    use loadmf6filemodule, only: loadmf6filetype
    class(gridarrayloadtype), intent(inout) :: this
    type(modflowinputtype), intent(in) :: mf6_input
    character(len=*), intent(in) :: component_name
    character(len=*), intent(in) :: component_input_name
    character(len=*), intent(in) :: input_name
    integer(I4B), intent(in) :: iperblock
    type(blockparsertype), pointer, intent(inout) :: parser
    integer(I4B), intent(in) :: iout
    type(loadmf6filetype) :: loader
    integer(I4B) :: isize
    integer(I4B), pointer :: maxbound
 
    ! initialize base type
    call this%DynamicPkgLoadType%init(mf6_input, component_name, &
                                      component_input_name, &
                                      input_name, iperblock, iout)
    ! initialize
    this%iout = iout
 
    ! load static input
    call loader%load(parser, mf6_input, this%nc_vars, this%input_name, iout)
 
    ! maxbound is optional
    call get_isize('MAXBOUND', mf6_input%mempath, isize)
    if (isize < 0) then
      ! set maxbound to grid nodes
      call mem_allocate(maxbound, 'MAXBOUND', mf6_input%mempath)
      maxbound = product(loader%mshape)
    end if
 
    ! initialize input context memory
    call this%ctx%init(mf6_input)
 
    ! allocate user nodelist
    call mem_allocate(this%nodeulist, this%ctx%maxbound, &
                      'NODEULIST', mf6_input%mempath)
 
    ! allocate dfn params
    call this%params_alloc()
  end subroutine ainit
 
  subroutine df(this)
    class(gridarrayloadtype), intent(inout) :: this
  end subroutine df
 
  subroutine ts_advance(this)
    class(gridarrayloadtype), intent(inout) :: this
  end subroutine ts_advance
 
  subroutine rp(this, parser)
    use blockparsermodule, only: blockparsertype
    use inputdefinitionmodule, only: inputparamdefinitiontype
    use definitionselectmodule, only: get_param_definition_type
    use arrayhandlersmodule, only: ifind
    use sourcecommonmodule, only: ifind_charstr
    use idmloggermodule, only: idm_log_header, idm_log_close
    class(gridarrayloadtype), intent(inout) :: this
    type(blockparsertype), pointer, intent(inout) :: parser
    logical(LGP) :: endOfBlock, netcdf, layered
    character(len=LINELENGTH) :: keyword, param_tag
    type(inputparamdefinitiontype), pointer :: idt
    integer(I4B) :: iaux
 
    ! reset for this period
    call this%reset()
 
    ! log lst file header
    call idm_log_header(this%mf6_input%component_name, &
                        this%mf6_input%subcomponent_name, this%iout)
 
    ! read array block
    do
      ! initialize
      iaux = 0
      netcdf = .false.
      layered = .false.
 
      ! read next line
      call parser%GetNextLine(endofblock)
      if (endofblock) exit
      ! read param_tag
      call parser%GetStringCaps(param_tag)
 
      ! is param tag an auxvar?
      iaux = ifind_charstr(this%ctx%auxname_cst, param_tag)
 
      ! any auxvar corresponds to the definition tag 'AUX'
      if (iaux > 0) param_tag = 'AUX'
 
      ! set input definition
      idt => get_param_definition_type(this%mf6_input%param_dfns, &
                                       this%mf6_input%component_type, &
                                       this%mf6_input%subcomponent_type, &
                                       'PERIOD', param_tag, this%input_name)
      ! look for Layered and NetCDF keywords
      call parser%GetStringCaps(keyword)
      if (keyword == 'LAYERED' .and. idt%layered) then
        layered = .true.
      else if (keyword == 'NETCDF') then
        netcdf = .true.
      end if
 
      ! read and load the parameter
      call this%param_load(parser, idt, this%mf6_input%mempath, layered, &
                           netcdf, iaux)
    end do
 
    ! log lst file header
    call idm_log_close(this%mf6_input%component_name, &
                       this%mf6_input%subcomponent_name, this%iout)
  end subroutine rp
 
  subroutine destroy(this)
    use memorymanagermodule, only: mem_deallocate
    class(gridarrayloadtype), intent(inout) :: this
    call mem_deallocate(this%nodeulist)
  end subroutine destroy
 
  subroutine reset(this)
    class(gridarrayloadtype), intent(inout) :: this
    integer(I4B) :: n
 
    this%ctx%nbound = 0
 
    do n = 1, this%nparam
      ! reset read state
      this%param_reads(n)%invar = 0
    end do
  end subroutine reset
 
  subroutine params_alloc(this)
    class(gridarrayloadtype), intent(inout) :: this
    character(len=LENVARNAME) :: rs_varname
    integer(I4B), pointer :: intvar
    integer(I4B) :: iparam
 
    ! set in scope param names
    call this%ctx%tags(this%param_names, this%nparam, this%input_name, &
                       create=.true.)
    call this%ctx%allocate_arrays()
 
    ! allocate and set param_reads pointer array
    allocate (this%param_reads(this%nparam))
 
    ! store read state variable pointers
    do iparam = 1, this%nparam
      ! allocate and store name of read state variable
      rs_varname = this%ctx%rsv_alloc(this%param_names(iparam))
      call mem_setptr(intvar, rs_varname, this%mf6_input%mempath)
      this%param_reads(iparam)%invar => intvar
      this%param_reads(iparam)%invar = 0
    end do
  end subroutine params_alloc
 
  subroutine param_load(this, parser, idt, mempath, layered, netcdf, iaux)
    use tdismodule, only: kper
    use constantsmodule, only: dnodata
    use arrayhandlersmodule, only: ifind
    use inputdefinitionmodule, only: inputparamdefinitiontype
    use double1dreadermodule, only: read_dbl1d
    use layeredarrayreadermodule, only: read_dbl1d_layered
    use loadncinputmodule, only: netcdf_read_array
    use sourcecommonmodule, only: get_layered_shape
    use idmloggermodule, only: idm_log_var
    class(gridarrayloadtype), intent(inout) :: this
    type(blockparsertype), intent(in) :: parser
    type(inputparamdefinitiontype), intent(in) :: idt
    character(len=*), intent(in) :: mempath
    logical(LGP), intent(in) :: layered
    logical(LGP), intent(in) :: netcdf
    integer(I4B), intent(in) :: iaux
    real(DP), dimension(:), pointer, contiguous :: dbl1d, nodes
    real(DP), dimension(:, :), pointer, contiguous :: dbl2d
    integer(I4B), dimension(:), allocatable :: layer_shape
    integer(I4B) :: iparam, n, nlay, nnode
 
    nnode = 0
 
    select case (idt%datatype)
    case ('DOUBLE1D')
      call mem_setptr(dbl1d, idt%mf6varname, mempath)
      allocate (nodes(this%ctx%nodes))
      if (netcdf) then
        call netcdf_read_array(nodes, this%ctx%mshape, idt, &
                               this%mf6_input, this%nc_vars, this%input_name, &
                               this%iout, kper)
      else if (layered) then
        call get_layered_shape(this%ctx%mshape, nlay, layer_shape)
        call read_dbl1d_layered(parser, nodes, idt%mf6varname, nlay, layer_shape)
      else
        call read_dbl1d(parser, nodes, idt%mf6varname)
      end if
 
      call idm_log_var(nodes, idt%tagname, mempath, this%iout)
 
      if (this%ctx%nbound > 0) then
        ! nodeulist already established: extract values at known positions
        do n = 1, this%ctx%nbound
          dbl1d(n) = nodes(this%nodeulist(n))
        end do
      else
        ! first array: filter by DNODATA to establish nodeulist and nbound
        do n = 1, this%ctx%nodes
          if (nodes(n) /= dnodata) then
            nnode = nnode + 1
            dbl1d(nnode) = nodes(n)
            this%nodeulist(nnode) = n
          end if
        end do
        this%ctx%nbound = nnode
      end if
      deallocate (nodes)
    case ('DOUBLE2D')
      call mem_setptr(dbl2d, idt%mf6varname, mempath)
      allocate (nodes(this%ctx%nodes))
 
      if (netcdf) then
        call netcdf_read_array(nodes, this%ctx%mshape, idt, &
                               this%mf6_input, this%nc_vars, this%input_name, &
                               this%iout, kper, iaux)
      else if (layered) then
        call get_layered_shape(this%ctx%mshape, nlay, layer_shape)
        call read_dbl1d_layered(parser, nodes, idt%mf6varname, nlay, layer_shape)
      else
        call read_dbl1d(parser, nodes, idt%mf6varname)
      end if
 
      call idm_log_var(nodes, idt%tagname, mempath, this%iout)
 
      if (this%ctx%nbound > 0) then
        ! nodeulist already established: extract values at known positions
        do n = 1, this%ctx%nbound
          dbl2d(iaux, n) = nodes(this%nodeulist(n))
        end do
      else
        ! first array: filter by DNODATA to establish nodeulist and nbound
        do n = 1, this%ctx%nodes
          if (nodes(n) /= dnodata) then
            nnode = nnode + 1
            dbl2d(iaux, nnode) = nodes(n)
            this%nodeulist(nnode) = n
          end if
        end do
        this%ctx%nbound = nnode
      end if
      deallocate (nodes)
    case default
      errmsg = 'IDM unimplemented. GridArrayLoad::param_load &
               &datatype='//trim(idt%datatype)
      call store_error(errmsg)
      call store_error_filename(this%input_name)
    end select
 
    ! if param is tracked set read state
    iparam = ifind(this%param_names, idt%tagname)
    if (iparam > 0) then
      this%param_reads(iparam)%invar = 1
    end if
  end subroutine param_load
 
end module gridarrayloadmodule