gwf-wel.f90

Go to the documentation of this file.

!> @brief This module contains the WEL package methods
!!
!! This module contains the overridden methods for the standard WEL package.
!! Several methods need to be overridden because of the AUTO_FLOW_REDUCE
!! option. Overridden methods include:
!!   - bnd_cf (AUTO_FLOW_REDUCE)
!!   - bnd_fc (AUTO_FLOW_REDUCE)
!!   - bnd_fn (AUTO_FLOW_REDUCE Newton-Raphson terms)
!!   - bnd_ot_package_flows (write AUTO_FLOW_REDUCE terms to csv file)
!!   - bnd_da (deallocate AUTO_FLOW_REDUCE variables)
!!   - bnd_bd_obs (wel-reduction observation added)
!!
!<
 
module welmodule
  ! -- modules used by WelModule methods
  use kindmodule, only: dp, i4b
  use constantsmodule, only: dzero, dem1, done, lenftype, dnodata, linelength, &
                             lenauxname
  use simvariablesmodule, only: errmsg, warnmsg
  use simmodule, only: count_errors, store_error, store_error_filename, &
                       store_warning
  use memoryhelpermodule, only: create_mem_path
  use bndmodule, only: bndtype
  use bndextmodule, only: bndexttype
  use obsmodule, only: defaultobsidprocessor
  use smoothingmodule, only: sqsaturation, sqsaturationderivative
  use observemodule, only: observetype
  use blockparsermodule, only: blockparsertype
  use inputoutputmodule, only: getunit, openfile
  use matrixbasemodule
  use flowreductioninterfacemodule, only: flowreductiontype
  !
  implicit none
  !
  private
  public :: wel_create
  public :: weltype
  !
  character(len=LENFTYPE) :: ftype = 'WEL' !< package ftype
  character(len=16) :: text = '             WEL' !< package flow text string
  !
  type, extends(bndexttype) :: weltype
    real(dp), dimension(:), pointer, contiguous :: q => null() !< volumetric well rate
    integer(I4B), pointer :: iflowred => null() !< flag indicating if the AUTO_FLOW_REDUCE option is active
    real(dp), pointer :: flowred => null() !< AUTO_FLOW_REDUCE variable
    integer(I4B), pointer :: ioutafrcsv => null() !< unit number for CSV output file containing wells with reduced puping rates
    integer(I4B), pointer :: iflowredlen => null() !< flag indicating flowred variable is a length value
    integer(I4B), pointer :: iafrauxcol => null() !< column index in auxvar for AUTO_FLOW_REDUCE_AUXNAME (0 if not active)
    class(flowreductiontype), pointer :: flowred_provider => null() !< supplies the flow-reduction factor; null when no formulation does
  contains
    procedure :: allocate_scalars => wel_allocate_scalars
    procedure :: allocate_arrays => wel_allocate_arrays
    procedure :: source_options => wel_options
    procedure :: log_wel_options
    procedure :: bnd_ck => wel_ck
    procedure :: bnd_cf => wel_cf
    procedure, private :: reduction_factor
    procedure :: bnd_fc => wel_fc
    procedure :: bnd_fn => wel_fn
    procedure :: bnd_da => wel_da
    procedure :: define_listlabel
    procedure :: bound_value => wel_bound_value
    procedure :: q_mult
    ! -- methods for observations
    procedure, public :: bnd_obs_supported => wel_obs_supported
    procedure, public :: bnd_df_obs => wel_df_obs
    procedure, public :: bnd_bd_obs => wel_bd_obs
    ! -- afr
    procedure, private :: wel_afr_csv_init
    procedure, private :: wel_afr_csv_write
  end type weltype
 
contains
 
  !> @ brief Create a new package object
  !!
  !!  Create a new WEL Package object
  !!
  !<
  subroutine wel_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
                        mempath)
    ! -- dummy variables
    class(bndtype), pointer :: packobj !< pointer to default package type
    integer(I4B), intent(in) :: id !< package id
    integer(I4B), intent(in) :: ibcnum !< boundary condition number
    integer(I4B), intent(in) :: inunit !< unit number of WEL package input file
    integer(I4B), intent(in) :: iout !< unit number of model listing file
    character(len=*), intent(in) :: namemodel !< model name
    character(len=*), intent(in) :: pakname !< package name
    character(len=*), intent(in) :: mempath !< input mempath
    ! -- local variables
    type(weltype), pointer :: welobj
    !
    ! -- allocate the object and assign values to object variables
    allocate (welobj)
    packobj => welobj
    !
    ! -- create name and memory path
    call packobj%set_names(ibcnum, namemodel, pakname, ftype, mempath)
    packobj%text = text
    !
    ! -- allocate scalars
    call welobj%allocate_scalars()
    !
    ! -- initialize package
    call packobj%pack_initialize()
 
    packobj%inunit = inunit
    packobj%iout = iout
    packobj%id = id
    packobj%ibcnum = ibcnum
    packobj%ncolbnd = 1
    packobj%ictMemPath = create_mem_path(namemodel, 'NPF')
  end subroutine wel_create
 
  !> @ brief Deallocate package memory
    !!
    !!  Deallocate WEL package scalars and arrays.
    !!
  !<
  subroutine wel_da(this)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy variables
    class(weltype) :: this !< WelType object
    !
    ! -- Deallocate parent package
    call this%BndExtType%bnd_da()
    !
    ! -- scalars
    call mem_deallocate(this%iflowred)
    call mem_deallocate(this%flowred)
    call mem_deallocate(this%ioutafrcsv)
    call mem_deallocate(this%iflowredlen)
    call mem_deallocate(this%iafrauxcol)
    call mem_deallocate(this%q, 'Q', this%memoryPath)
    nullify (this%flowred_provider)
  end subroutine wel_da
 
  !> @ brief Allocate scalars
    !!
    !! Allocate and initialize scalars for the WEL package. The base model
    !! allocate scalars method is also called.
    !!
  !<
  subroutine wel_allocate_scalars(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy variables
    class(weltype) :: this !< WelType object
    !
    ! -- call base type allocate scalars
    call this%BndExtType%allocate_scalars()
    !
    ! -- allocate the object and assign values to object variables
    call mem_allocate(this%iflowred, 'IFLOWRED', this%memoryPath)
    call mem_allocate(this%flowred, 'FLOWRED', this%memoryPath)
    call mem_allocate(this%ioutafrcsv, 'IOUTAFRCSV', this%memoryPath)
    call mem_allocate(this%iflowredlen, 'IFLOWREDLEN', this%memoryPath)
    call mem_allocate(this%iafrauxcol, 'IAFRAUXCOL', this%memoryPath)
    !
    ! -- Set values
    this%iflowred = 0
    this%ioutafrcsv = 0
    this%flowred = dzero
    this%iflowredlen = 0
    this%iafrauxcol = 0
  end subroutine wel_allocate_scalars
 
  !> @ brief Allocate arrays
    !!
    !! Allocate and initialize arrays for the WEL package
    !!
  !<
  subroutine wel_allocate_arrays(this, nodelist, auxvar)
    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr, mem_checkin
    ! -- dummy
    class(weltype) :: this
    integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
    real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
    ! -- local
    !
    ! -- call BndExtType allocate scalars
    call this%BndExtType%allocate_arrays(nodelist, auxvar)
    !
    ! -- set constant head array input context pointer
    call mem_setptr(this%q, 'Q', this%input_mempath)
    !
    ! -- checkin constant head array input context pointer
    call mem_checkin(this%q, 'Q', this%memoryPath, &
                     'Q', this%input_mempath)
  end subroutine wel_allocate_arrays
 
  !> @ brief Source additional options for package
    !!
    !!  Source additional options for WEL package.
    !!
  !<
  subroutine wel_options(this)
    ! -- modules
    use inputoutputmodule, only: urword
    use memorymanagerextmodule, only: mem_set_value
    use gwfwelinputmodule, only: gwfwelparamfoundtype
    ! -- dummy variables
    class(weltype), intent(inout) :: this !< WelType object
    ! -- local variables
    character(len=LINELENGTH) :: fname
    character(len=LENAUXNAME) :: afrauxname
    type(gwfwelparamfoundtype) :: found
    integer(I4B) :: n
    ! -- formats
    character(len=*), parameter :: fmtflowred = &
      &"(4x, 'AUTOMATIC FLOW REDUCTION OF WELLS IMPLEMENTED.')"
    character(len=*), parameter :: fmtflowredv = &
      &"(4x, 'AUTOMATIC FLOW REDUCTION FRACTION (',g15.7,').')"
    !
    ! -- source base BndExtType options
    call this%BndExtType%source_options()
    !
    ! -- source well options from input context
    call mem_set_value(this%flowred, 'FLOWRED', this%input_mempath, found%flowred)
    call mem_set_value(fname, 'AFRCSVFILE', this%input_mempath, found%afrcsvfile)
    call mem_set_value(this%imover, 'MOVER', this%input_mempath, found%mover)
    call mem_set_value(this%iflowredlen, 'IFLOWREDLEN', this%input_mempath, &
                       found%iflowredlen)
    call mem_set_value(afrauxname, 'AFRAUXNAME', this%input_mempath, &
                       found%afrauxname)
 
    if (found%iflowredlen) then
      if (found%flowred .eqv. .false.) then
        write (warnmsg, '(a)') &
          'FLOW_REDUCTION_LENGTH option specified but a AUTO_FLOW_REDUCTION value &
          &is not specified. The FLOW_REDUCTION_LENGTH option will be ignored.'
        call store_warning(warnmsg)
      else
        this%iflowredlen = 1
      end if
    end if
 
    if (found%flowred) then
      this%iflowred = 1
      if (this%flowred <= dzero) then
        if (found%iflowredlen) then
          write (errmsg, '(a)') &
            'An AUTO_FLOW_REDUCTION value less than or equal to zero cannot be &
            &specified if the FLOW_REDUCTION_LENGTH option is specified.'
          call store_error(errmsg)
        else
          this%flowred = dem1
        end if
      else if (this%flowred > done .and. this%iflowredlen == 0) then
        this%flowred = done
      end if
    end if
 
    if (found%afrcsvfile) then
      call this%wel_afr_csv_init(fname)
    end if
 
    if (found%mover) then
      this%imover = 1
    end if
 
    if (found%afrauxname) then
      if (.not. found%flowred) then
        write (warnmsg, '(a)') &
          'AUTO_FLOW_REDUCE_AUXNAME is specified but AUTO_FLOW_REDUCE is not &
          &specified. The AUTO_FLOW_REDUCE_AUXNAME option will be ignored.'
        call store_warning(warnmsg)
      end if
      if (.not. found%iflowredlen) then
        write (warnmsg, '(a)') &
          'AUTO_FLOW_REDUCE_AUXNAME is specified but FLOW_REDUCTION_LENGTH is &
          &not specified. The AUTO_FLOW_REDUCE_AUXNAME value will be &
          &interpreted as a fraction of the cell thickness.'
        call store_warning(warnmsg)
      end if
      if (found%flowred) then
        if (this%naux == 0) then
          write (errmsg, '(a,2(1x,a))') &
            'AUTO_FLOW_REDUCE_AUXNAME was specified as', &
            trim(adjustl(afrauxname)), 'but no AUX variables specified.'
          call store_error(errmsg)
        end if
        this%iafrauxcol = 0
        do n = 1, this%naux
          if (afrauxname == this%auxname(n)) then
            this%iafrauxcol = n
            exit
          end if
        end do
        if (this%iafrauxcol == 0) then
          write (errmsg, '(a,2(1x,a))') &
            'AUTO_FLOW_REDUCE_AUXNAME was specified as', &
            trim(adjustl(afrauxname)), &
            'but no AUX variable found with this name.'
          call store_error(errmsg)
        end if
      end if
    end if
 
    ! -- log WEL specific options
    call this%log_wel_options(found)
    !
    ! -- terminate if errors were detected
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
  end subroutine wel_options
 
  !> @ brief Log WEL specific package options
  !<
  subroutine log_wel_options(this, found)
    ! -- modules
    use gwfwelinputmodule, only: gwfwelparamfoundtype
    ! -- dummy variables
    class(weltype), intent(inout) :: this
    type(gwfwelparamfoundtype), intent(in) :: found
    ! -- local variables
    ! -- format
    character(len=*), parameter :: fmtflowred = &
      &"(4x, 'AUTOMATIC FLOW REDUCTION OF WELLS IMPLEMENTED.')"
    character(len=*), parameter :: fmtflowredv = &
      &"(4x, 'AUTOMATIC FLOW REDUCTION FRACTION (',g15.7,').')"
    character(len=*), parameter :: fmtflowredl = &
      &"(4x, 'AUTOMATIC FLOW REDUCTION LENGTH (',g15.7,').')"
    !
    ! -- log found options
    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text)) &
      //' OPTIONS'
 
    if (found%iflowredlen) then
      write (this%iout, fmtflowred)
      write (this%iout, '(4x,A)') &
        'AUTOMATIC FLOW REDUCTION FRACTION INTERPRETED AS A LENGTH'
    end if
 
    if (found%flowred) then
      if (this%iflowredlen == 0) then
        write (this%iout, fmtflowredv) this%flowred
      else
        write (this%iout, fmtflowredl) this%flowred
      end if
    end if
    !
    if (found%afrcsvfile) then
      ! -- currently no-op
    end if
 
    if (found%mover) then
      write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
    end if
 
    if (found%afrauxname) then
      write (this%iout, '(4x,A)') &
        'AUTO_FLOW_REDUCE_AUXNAME OPTION ENABLED FOR PER-WELL FLOW REDUCTION'
    end if
    !
    ! -- close logging block
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' OPTIONS'
  end subroutine log_wel_options
 
  !> @ brief Check WEL period data.
    !!
    !!  Verify that the per-well AUTO_FLOW_REDUCE_AUXNAME auxiliary values are
    !!  within valid bounds. When FLOW_REDUCTION_LENGTH is not specified the
    !!  value is interpreted as a fraction of the cell thickness and must be
    !!  between 0 and 1. When FLOW_REDUCTION_LENGTH is specified the value is
    !!  interpreted as a length above the cell bottom and must be between 0 and
    !!  the cell thickness. The check runs each stress period because the
    !!  auxiliary values may change between periods.
    !!
  !<
  subroutine wel_ck(this)
    ! -- modules
    use simmodule, only: count_errors, store_error_unit
    ! -- dummy variables
    class(weltype), intent(inout) :: this !< WelType object
    ! -- local variables
    character(len=LINELENGTH) :: errmsg
    integer(I4B) :: i
    integer(I4B) :: node
    real(DP) :: afraux
    real(DP) :: thick
    ! -- formats
    character(len=*), parameter :: fmtfracerr = &
      "('WELL (',i0,') AUTO_FLOW_REDUCE_AUXNAME value (',g0,') must be greater &
      &than 0 and less than or equal to 1 when FLOW_REDUCTION_LENGTH is not &
      &specified (it is interpreted as a fraction of the cell thickness).')"
    character(len=*), parameter :: fmtlenerr = &
      "('WELL (',i0,') AUTO_FLOW_REDUCE_AUXNAME value (',g0,') must be greater &
      &than 0 and less than or equal to the cell thickness (',g0,') when &
      &FLOW_REDUCTION_LENGTH is specified.')"
    !
    ! -- nothing to check unless AUTO_FLOW_REDUCE_AUXNAME is active
    if (this%iflowred == 0 .or. this%iafrauxcol == 0) return
    !
    ! -- check the per-well auxiliary flow reduction values. The value must be
    !    strictly greater than zero: a value of zero places the turnoff
    !    threshold at the cell bottom, giving a zero-width smoothing interval
    !    (a divide-by-zero in sQSaturation), consistent with the global
    !    AUTO_FLOW_REDUCE length-mode check.
    do i = 1, this%nbound
      node = this%nodelist(i)
      if (node == 0) cycle
      ! -- without a provider the reduction applies only to convertible
      !    cells; with one, any cell can be reduced, so the per-well value
      !    must be valid there too
      if (this%icelltype(node) == 0 .and. &
          .not. associated(this%flowred_provider)) cycle
      afraux = this%auxvar(this%iafrauxcol, i)
      if (this%iflowredlen == 0) then
        ! -- value is a fraction of the cell thickness: valid range (0, 1]
        if (afraux <= dzero .or. afraux > done) then
          write (errmsg, fmt=fmtfracerr) i, afraux
          call store_error(errmsg)
        end if
      else
        ! -- value is a length above the cell bottom: valid range (0, thick]
        thick = this%dis%top(node) - this%dis%bot(node)
        if (afraux <= dzero .or. afraux > thick) then
          write (errmsg, fmt=fmtlenerr) i, afraux, thick
          call store_error(errmsg)
        end if
      end if
    end do
    !
    ! -- terminate if any errors were detected
    if (count_errors() > 0) then
      call store_error_unit(this%inunit)
    end if
  end subroutine wel_ck
 
  !> @ brief Formulate the package hcof and rhs terms.
    !!
    !!  Formulate the hcof and rhs terms for the WEL package that will be
    !!  added to the coefficient matrix and right-hand side vector.
    !!
  !<
  subroutine wel_cf(this)
    ! -- dummy variables
    class(weltype) :: this !< WelType object
    ! -- local variables
    integer(I4B) :: i, node
    real(DP) :: f
    real(DP) :: dfdh
    real(DP) :: q
    !
    ! -- Return if no wells
    if (this%nbound == 0) return
    !
    ! -- Calculate hcof and rhs for each well entry
    do i = 1, this%nbound
      node = this%nodelist(i)
      this%hcof(i) = dzero
      if (this%ibound(node) <= 0) then
        this%rhs(i) = dzero
        cycle
      end if
      q = this%q_mult(i)
      if (this%iflowred /= 0 .and. q < dzero) then
        call this%reduction_factor(i, node, f, dfdh)
        q = q * f
      end if
      this%rhs(i) = -q
    end do
  end subroutine wel_cf
 
  !> @brief Factor by which a withdrawal is reduced, and its derivative
  !!
  !! The taper length is what the user asked for: the per-well auxiliary
  !! value or the package value, a length as given or a fraction of the cell
  !! thickness. Without a provider the factor tapers up from the cell bottom
  !! over that length and applies to convertible cells only, as it always
  !! has. A provider decides for itself how much of the fluid being withdrawn
  !! the cell holds and returns the factor directly, so with one a cell can
  !! be reduced whatever its cell type.
  !<
  subroutine reduction_factor(this, i, node, f, dfdh)
    class(weltype), intent(inout) :: this !< this WEL instance
    integer(I4B), intent(in) :: i !< well number
    integer(I4B), intent(in) :: node !< reduced node number
    real(DP), intent(out) :: f !< factor to apply to the requested rate
    real(DP), intent(out) :: dfdh !< derivative of the factor with respect to head
    ! -- local variables
    real(DP) :: taper, bt, tp
 
    f = done
    dfdh = dzero
 
    ! -- the length to taper over
    if (this%iafrauxcol > 0) then
      taper = this%auxvar(this%iafrauxcol, i)
    else
      taper = this%flowred
    end if
    if (this%iflowredlen == 0) then
      taper = taper * (this%dis%top(node) - this%dis%bot(node))
    end if
 
    ! -- a formulation that knows what the cell holds answers for it
    if (associated(this%flowred_provider)) then
      call this%flowred_provider%reduce(node, this%xnew(node), taper, f, dfdh)
      return
    end if
 
    ! -- otherwise the cell itself: taper up from the bottom, convertible only
    if (this%icelltype(node) == 0) return
    bt = this%dis%bot(node)
    tp = bt + taper
    if (tp <= bt) return
    f = sqsaturation(tp, bt, this%xnew(node))
    dfdh = sqsaturationderivative(tp, bt, this%xnew(node))
 
  end subroutine reduction_factor
 
  !> @ brief Copy hcof and rhs terms into solution.
    !!
    !!  Add the hcof and rhs terms for the WEL package to the
    !!  coefficient matrix and right-hand side vector.
    !!
  !<
  subroutine wel_fc(this, rhs, ia, idxglo, matrix_sln)
    ! -- dummy variables
    class(weltype) :: this !< WelType object
    real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector for model
    integer(I4B), dimension(:), intent(in) :: ia !< solution CRS row pointers
    integer(I4B), dimension(:), intent(in) :: idxglo !< mapping vector for model (local) to solution (global)
    class(matrixbasetype), pointer :: matrix_sln !< solution coefficient matrix
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: n
    integer(I4B) :: ipos
    !
    ! -- pakmvrobj fc
    if (this%imover == 1) then
      call this%pakmvrobj%fc()
    end if
    !
    ! -- Copy package rhs and hcof into solution rhs and amat
    do i = 1, this%nbound
      n = this%nodelist(i)
      rhs(n) = rhs(n) + this%rhs(i)
      ipos = ia(n)
      call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i))
      !
      ! -- If mover is active and this well is discharging,
      !    store available water (as positive value).
      if (this%imover == 1 .and. this%rhs(i) > dzero) then
        call this%pakmvrobj%accumulate_qformvr(i, this%rhs(i))
      end if
    end do
  end subroutine wel_fc
 
  !> @ brief Add Newton-Raphson terms for package into solution.
    !!
    !!  Calculate and add the Newton-Raphson terms for the WEL package to the
    !!  coefficient matrix and right-hand side vector.
    !!
  !<
  subroutine wel_fn(this, rhs, ia, idxglo, matrix_sln)
    ! -- dummy variables
    class(weltype) :: this !< WelType object
    real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector for model
    integer(I4B), dimension(:), intent(in) :: ia !< solution CRS row pointers
    integer(I4B), dimension(:), intent(in) :: idxglo !< mapping vector for model (local) to solution (global)
    class(matrixbasetype), pointer :: matrix_sln !< solution coefficient matrix
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: node
    integer(I4B) :: ipos
    real(DP) :: f
    real(DP) :: dfdh
    real(DP) :: drterm
    real(DP) :: q
    !
    ! -- Copy package rhs and hcof into solution rhs and amat
    do i = 1, this%nbound
      node = this%nodelist(i)
      !
      ! -- test if node is constant or inactive
      if (this%ibound(node) <= 0) then
        cycle
      end if
      !
      ! -- well rate is possibly head dependent
      if (this%iflowred /= 0) then
        q = -this%rhs(i)
        if (q < dzero) then
          ! -- calculate derivative for well
          call this%reduction_factor(i, node, f, dfdh)
          drterm = dfdh * this%q_mult(i)
          if (drterm == dzero) cycle
          !--fill amat and rhs with newton-raphson terms
          ipos = ia(node)
          call matrix_sln%add_value_pos(idxglo(ipos), drterm)
          rhs(node) = rhs(node) + drterm * this%xnew(node)
        end if
      end if
    end do
  end subroutine wel_fn
 
  !> @brief Initialize the auto flow reduce csv output file
  subroutine wel_afr_csv_init(this, fname)
    ! -- dummy variables
    class(weltype), intent(inout) :: this !< WelType object
    character(len=*), intent(in) :: fname
    ! -- format
    character(len=*), parameter :: fmtafrcsv = &
      "(4x, 'AUTO FLOW REDUCE INFORMATION WILL BE SAVED TO FILE: ', a, /4x, &
      &'OPENED ON UNIT: ', I0)"
 
    this%ioutafrcsv = getunit()
    call openfile(this%ioutafrcsv, this%iout, fname, 'CSV', &
                  filstat_opt='REPLACE')
    write (this%iout, fmtafrcsv) trim(adjustl(fname)), &
      this%ioutafrcsv
    write (this%ioutafrcsv, '(a)') &
      'time,period,step,boundnumber,cellnumber,rate-requested,&
      &rate-actual,wel-reduction'
  end subroutine wel_afr_csv_init
 
  !> @brief Write out auto flow reductions only when & where they occur
  subroutine wel_afr_csv_write(this)
    ! -- modules
    use tdismodule, only: totim, kstp, kper
    ! -- dummy variables
    class(weltype), intent(inout) :: this !< WelType object
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: nodereduced
    integer(I4B) :: nodeuser
    real(DP) :: v
    ! -- format
    do i = 1, this%nbound
      nodereduced = this%nodelist(i)
      !
      ! -- test if node is constant or inactive
      if (this%ibound(nodereduced) <= 0) then
        cycle
      end if
      v = this%q_mult(i) + this%rhs(i)
      if (v < dzero) then
        nodeuser = this%dis%get_nodeuser(nodereduced)
        write (this%ioutafrcsv, '(*(G0,:,","))') &
          totim, kper, kstp, i, nodeuser, this%q_mult(i), this%simvals(i), v
      end if
    end do
  end subroutine wel_afr_csv_write
 
  !> @ brief Define the list label for the package
    !!
    !!  Method defined the list label for the WEL package. The list label is
    !!  the heading that is written to iout when PRINT_INPUT option is used.
    !!
  !<
  subroutine define_listlabel(this)
    ! -- dummy variables
    class(weltype), intent(inout) :: this !< WelType object
    !
    ! -- create the header list label
    this%listlabel = trim(this%filtyp)//' NO.'
    if (this%dis%ndim == 3) then
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW'
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL'
    elseif (this%dis%ndim == 2) then
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D'
    else
      write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE'
    end if
    write (this%listlabel, '(a, a16)') trim(this%listlabel), 'STRESS RATE'
    if (this%inamedbound == 1) then
      write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
    end if
  end subroutine define_listlabel
 
  ! -- Procedures related to observations
 
  !> @brief Determine if observations are supported.
    !!
    !! Function to determine if observations are supported by the WEL package.
    !! Observations are supported by the WEL package.
    !!
    !! @return  wel_obs_supported       boolean indicating if observations are supported
    !!
  !<
  logical function wel_obs_supported(this)
    ! -- dummy variables
    class(weltype) :: this !< WelType object
    !
    ! -- set boolean
    wel_obs_supported = .true.
  end function wel_obs_supported
 
  !> @brief Define the observation types available in the package
    !!
    !! Method to define the observation types available in the WEL package.
    !!
  !<
  subroutine wel_df_obs(this)
    ! -- dummy variables
    class(weltype) :: this !< WelType object
    ! -- local variables
    integer(I4B) :: indx
    !
    ! -- initialize observations
    call this%obs%StoreObsType('wel', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
    !
    ! -- Store obs type and assign procedure pointer
    !    for to-mvr observation type.
    call this%obs%StoreObsType('to-mvr', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
    !
    ! -- Store obs type and assign procedure pointer
    !    for wel-reduction observation type.
    call this%obs%StoreObsType('wel-reduction', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
  end subroutine wel_df_obs
 
  !> @brief Save observations for the package
    !!
    !! Method to save simulated values for the WEL package.
    !!
  !<
  subroutine wel_bd_obs(this)
    ! -- dummy variables
    class(weltype) :: this !< WelType object
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: n
    integer(I4B) :: jj
    real(DP) :: v
    type(observetype), pointer :: obsrv => null()
    !
    ! -- clear the observations
    call this%obs%obs_bd_clear()
    !
    ! -- Save simulated values for all of package's observations.
    do i = 1, this%obs%npakobs
      obsrv => this%obs%pakobs(i)%obsrv
      if (obsrv%BndFound) then
        do n = 1, obsrv%indxbnds_count
          v = dnodata
          jj = obsrv%indxbnds(n)
          select case (obsrv%ObsTypeId)
          case ('TO-MVR')
            if (this%imover == 1) then
              v = this%pakmvrobj%get_qtomvr(jj)
              if (v > dzero) then
                v = -v
              end if
            end if
          case ('WEL')
            v = this%simvals(jj)
          case ('WEL-REDUCTION')
            if (this%iflowred > 0) then
              v = this%q_mult(jj) + this%rhs(jj)
            end if
          case default
            errmsg = 'Unrecognized observation type: '//trim(obsrv%ObsTypeId)
            call store_error(errmsg)
          end select
          call this%obs%SaveOneSimval(obsrv, v)
        end do
      else
        call this%obs%SaveOneSimval(obsrv, dnodata)
      end if
    end do
    !
    ! -- Write the auto flow reduce csv file entries for this step
    if (this%ioutafrcsv > 0) then
      call this%wel_afr_csv_write()
    end if
  end subroutine wel_bd_obs
 
  function q_mult(this, row) result(q)
    ! -- modules
    use constantsmodule, only: dzero
    ! -- dummy variables
    class(weltype), intent(inout) :: this
    integer(I4B), intent(in) :: row
    ! -- result
    real(dp) :: q
    !
    if (this%iauxmultcol > 0) then
      q = this%q(row) * this%auxvar(this%iauxmultcol, row)
    else
      q = this%q(row)
    end if
  end function q_mult
 
  !> @ brief Return a bound value
    !!
    !!  Return a bound value associated with an ncolbnd index
    !!  and row.
    !!
  !<
  function wel_bound_value(this, col, row) result(bndval)
    ! -- modules
    use constantsmodule, only: dzero
    ! -- dummy variables
    class(weltype), intent(inout) :: this
    integer(I4B), intent(in) :: col
    integer(I4B), intent(in) :: row
    ! -- result
    real(dp) :: bndval
    !
    select case (col)
    case (1)
      bndval = this%q_mult(row)
    case default
      errmsg = 'Programming error. WEL bound value requested column '&
               &'outside range of ncolbnd (1).'
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end select
  end function wel_bound_value
 
end module welmodule