gwf-maw.f90

Go to the documentation of this file.

module mawmodule
  !
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: linelength, lenboundname, &
                             lenbudtxt, dzero, dem9, dem6, dem4, dem2, dquarter, &
                             dhalf, dp7, dp9, done, dtwo, dpi, dtwopi, deight, &
                             dhundred, dep20, namedboundflag, lenpackagename, &
                             lenauxname, lenftype, lenpakloc, dhnoflo, dhdry, &
                             dnodata, &
                             maxcharlen, tableft, tabcenter, tabright, &
                             tabstring, tabucstring, tabinteger, tabreal, mnormal
  use smoothingmodule, only: squadraticsaturation, sqsaturation, &
                             squadraticsaturationderivative, &
                             sqsaturationderivative, &
                             squadratic0sp, &
                             squadratic0spderivative
  use bndmodule, only: bndtype
  use bndextmodule, only: bndexttype
  use budgetobjectmodule, only: budgetobjecttype, budgetobject_cr
  use tablemodule, only: tabletype, table_cr
  use observemodule, only: observetype
  use obsmodule, only: obstype
  use geomutilmodule, only: get_node
  use inputoutputmodule, only: urword, extract_idnum_or_bndname, &
                               getunit, openfile, assign_iounit
  use basedismodule, only: disbasetype
  use simmodule, only: count_errors, store_error, store_error_filename, &
                       store_warning
  use simvariablesmodule, only: errmsg, warnmsg
  use memorymanagermodule, only: mem_allocate, mem_reallocate, mem_setptr, &
                                 mem_deallocate, get_isize
  use memorymanagerextmodule, only: mem_set_value, memorystore_release
  use memoryhelpermodule, only: create_mem_path
  use characterstringmodule, only: characterstringtype
  use openspecmodule, only: access, form
  use sourcecommonmodule, only: filein_fname
  use gwfmawinputmodule, only: gwfmawparamfoundtype
  use matrixbasemodule
  !
  implicit none
 
  public :: mawtype
 
  !
  character(len=LENFTYPE) :: ftype = 'MAW'
  character(len=LENPACKAGENAME) :: text = '             MAW'
 
  private
  public :: maw_create
  public :: maw_damp_weight
  public :: maw_screen_length
  !
  !> @brief input context data pointers for MAW.
  !!
  !! Holds persistent pointers to input context arrays used at
  !! runtime.
  !<
  type :: mawinputtype
    ! -- packagedata (TS-linked; set in maw_source_packagedata)
    real(dp), dimension(:), pointer, contiguous :: strt => null()
    ! -- period
    integer(I4B), pointer :: nbound => null()
    integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: setting => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: status => null()
    integer(I4B), dimension(:), pointer, contiguous :: period_icon => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: connstatus => null()
    real(dp), dimension(:), pointer, contiguous :: period_bottom => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: head_limit => null()
    real(dp), dimension(:), pointer, contiguous :: fwelev => null()
    real(dp), dimension(:), pointer, contiguous :: fwcond => null()
    real(dp), dimension(:), pointer, contiguous :: fwrlen => null()
    real(dp), dimension(:), pointer, contiguous :: minrate => null()
    real(dp), dimension(:), pointer, contiguous :: maxrate => null()
    real(dp), dimension(:), pointer, contiguous :: pump_elevation => null()
    real(dp), dimension(:), pointer, contiguous :: scaling_length => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: auxname => null()
    real(dp), dimension(:), pointer, &
      contiguous :: auxval => null()
  contains
    procedure :: init => maw_input_init
    procedure :: destroy => maw_input_destroy
  end type mawinputtype
  !
  type, extends(bndexttype) :: mawtype
    !
    ! -- scalars
    ! -- characters
    !
    character(len=LENBUDTXT), dimension(:), pointer, &
      contiguous :: cmawbudget => null()
    character(len=LENAUXNAME), dimension(:), pointer, &
      contiguous :: cauxcbc => null()
    !
    ! -- logical
    logical(LGP), pointer :: correct_flow => null()
    !
    ! -- integers
    integer(I4B), pointer :: iprhed => null()
    integer(I4B), pointer :: iheadout => null()
    integer(I4B), pointer :: ibudgetout => null()
    integer(I4B), pointer :: ibudcsv => null()
    integer(I4B), pointer :: cbcauxitems => null()
    integer(I4B), pointer :: iflowingwells => null()
    integer(I4B), pointer :: imawiss => null()
    integer(I4B), pointer :: imawissopt => null()
    integer(I4B), pointer :: nmawwells => null()
    integer(I4B), pointer :: check_attr => null()
    integer(I4B), pointer :: ishutoffcnt => null()
    integer(I4B), pointer :: ieffradopt => null()
    integer(I4B), pointer :: inonvert => null() !< flag indicating that non-vertical (slanted) well connections are simulated
    integer(I4B), pointer :: ioutredflowcsv => null() !< unit number for CSV output file containing MAWs with reduced extraction/injection rates
    real(dp), pointer :: satomega => null()
    !
    ! -- for underrelaxation of estimated well q if using shutoff
    real(dp), pointer :: theta => null()
    real(dp), pointer :: kappa => null()
    !
    ! -- vector data for each well
    character(len=8), dimension(:), pointer, contiguous :: status => null()
    integer(I4B), dimension(:), pointer, contiguous :: ngwfnodes => null()
    integer(I4B), dimension(:), pointer, contiguous :: ieqn => null()
    integer(I4B), dimension(:), pointer, contiguous :: ishutoff => null()
    integer(I4B), dimension(:), pointer, contiguous :: ifwdischarge => null()
    real(dp), dimension(:), pointer, contiguous :: strt => null()
    real(dp), dimension(:), pointer, contiguous :: radius => null()
    real(dp), dimension(:), pointer, contiguous :: area => null()
    real(dp), dimension(:), pointer, contiguous :: pumpelev => null()
    real(dp), dimension(:), pointer, contiguous :: bot => null()
    real(dp), dimension(:), pointer, contiguous :: ratesim => null()
    real(dp), dimension(:), pointer, contiguous :: qsim0 => null() !well rate from the previous outer iteration (rate convergence check)
    real(dp), dimension(:), pointer, contiguous :: reduction_length => null()
    real(dp), dimension(:), pointer, contiguous :: fwelev => null()
    real(dp), dimension(:), pointer, contiguous :: fwcond => null()
    real(dp), dimension(:), pointer, contiguous :: fwrlen => null()
    real(dp), dimension(:), pointer, contiguous :: fwcondsim => null()
    real(dp), dimension(:), pointer, contiguous :: xsto => null()
    real(dp), dimension(:), pointer, contiguous :: xoldsto => null()
    real(dp), dimension(:), pointer, contiguous :: shutoffmin => null()
    real(dp), dimension(:), pointer, contiguous :: shutoffmax => null()
    real(dp), dimension(:), pointer, contiguous :: shutofflevel => null()
    real(dp), dimension(:), pointer, contiguous :: shutoffweight => null()
    real(dp), dimension(:), pointer, contiguous :: shutoffdq => null()
    real(dp), dimension(:), pointer, contiguous :: shutoffqold => null()
    real(dp), dimension(:), pointer, contiguous :: nurdxold => null() !well head change applied in the previous outer iteration
    real(dp), dimension(:), pointer, contiguous :: nurweight => null() !how much each well head change is cut back to stop it from oscillating
    character(len=LENBOUNDNAME), dimension(:), pointer, &
      contiguous :: cmawname => null()
    !
    ! -- time-series aware data
    real(dp), dimension(:), pointer, contiguous :: rate => null()
    real(dp), dimension(:), pointer, contiguous :: well_head => null()
    !
    ! -- .true. once a PERIOD WELL_HEAD entry has been given for a well;
    !    WELL_HEAD then takes precedence over the PACKAGEDATA STRT baseline
    !    for CONSTANT/INACTIVE wells
    logical(LGP), dimension(:), pointer, contiguous :: well_head_set => null()
    !
    ! -- input context data pointers
    type(mawinputtype) :: input
    !
    ! -- ia vector for connections
    integer(I4B), dimension(:), pointer, contiguous :: iaconn => null()
    !
    ! -- vector data for each connections
    integer(I4B), dimension(:), pointer, contiguous :: gwfnodes => null()
    integer(I4B), dimension(:), pointer, contiguous :: iboundconn => null() !< connection status, 0 = inactive; a connection is simulated only if its well is also active
    real(dp), dimension(:), pointer, contiguous :: sradius => null()
    real(dp), dimension(:), pointer, contiguous :: hk => null()
    real(dp), dimension(:), pointer, contiguous :: satcond => null()
    real(dp), dimension(:), pointer, contiguous :: simcond => null()
    real(dp), dimension(:), pointer, contiguous :: topscrn => null()
    real(dp), dimension(:), pointer, contiguous :: botscrn => null()
    real(dp), dimension(:), pointer, contiguous :: angle => null() !< tilt angle of a non-vertical connection, in degrees from vertical (0.0 = vertical)
    real(dp), dimension(:), pointer, contiguous :: connlen => null() !< user-specified in-cell screen length for a non-vertical connection (<= 0.0 = derive from angle and screen elevations)
    real(dp), dimension(:), pointer, contiguous :: usrtopscrn => null() !< user-specified screen top, retained so a non-vertical SPECIFIED connection can honor it
    real(dp), dimension(:), pointer, contiguous :: usrbotscrn => null() !< user-specified screen bottom, retained so a non-vertical SPECIFIED connection can honor it
    !
    ! -- imap vector
    integer(I4B), dimension(:), pointer, contiguous :: imap => null()
    !
    ! -- maw output data
    real(dp), dimension(:), pointer, contiguous :: qauxcbc => null()
    real(dp), dimension(:), pointer, contiguous :: dbuff => null()
    real(dp), dimension(:), pointer, contiguous :: qleak => null()
    real(dp), dimension(:), pointer, contiguous :: qout => null()
    real(dp), dimension(:), pointer, contiguous :: qfw => null()
    real(dp), dimension(:), pointer, contiguous :: qsto => null()
    real(dp), dimension(:), pointer, contiguous :: qconst => null()
    !
    ! -- for budgets
    integer(I4B), pointer :: bditems => null()
    type(budgetobjecttype), pointer :: budobj => null()
    !
    ! -- table objects
    type(tabletype), pointer :: headtab => null()
    !
    ! -- pointer to gwf iss, k11, k22.
    integer(I4B), pointer :: gwfiss => null()
    real(dp), dimension(:), pointer, contiguous :: gwfk11 => null()
    real(dp), dimension(:), pointer, contiguous :: gwfk22 => null()
    integer(I4B), pointer :: gwfik22 => null()
    real(dp), dimension(:), pointer, contiguous :: gwfsat => null()
    !
    ! -- arrays for handling the rows added to the solution matrix
    integer(I4B), dimension(:), pointer, contiguous :: idxlocnode => null() !map position in global rhs and x array of pack entry
    integer(I4B), dimension(:), pointer, contiguous :: idxdglo => null() !map position in global array of package diagonal row entries
    integer(I4B), dimension(:), pointer, contiguous :: idxoffdglo => null() !map position in global array of package off diagonal row entries
    integer(I4B), dimension(:), pointer, contiguous :: idxsymdglo => null() !map position in global array of package diagonal entries to model rows
    integer(I4B), dimension(:), pointer, contiguous :: idxsymoffdglo => null() !map position in global array of package off diagonal entries to model rows
    integer(I4B), dimension(:), pointer, contiguous :: iboundpak => null() !package ibound
    real(dp), dimension(:), pointer, contiguous :: xnewpak => null() !package x vector
    real(dp), dimension(:), pointer, contiguous :: xoldpak => null() !package xold vector
    !
    ! -- density variables
    integer(I4B), pointer :: idense
    real(dp), dimension(:, :), pointer, contiguous :: denseterms => null()
    !
    ! -- viscosity variables
    real(dp), dimension(:, :), pointer, contiguous :: viscratios => null() !< viscosity ratios (1: maw vsc ratio; 2: gwf vsc ratio)
    !
    ! -- type bound procedures
 
  contains
 
    procedure :: maw_allocate_scalars
    procedure :: maw_allocate_well_conn_arrays
    procedure :: maw_allocate_arrays
    procedure :: source_options => maw_source_options
    procedure :: source_dimensions => maw_source_dimensions
    procedure :: read_initial_attr => maw_read_initial_attr
    procedure :: set_pointers => maw_set_pointers
    procedure :: bnd_ac => maw_ac
    procedure :: bnd_mc => maw_mc
    procedure :: bnd_ar => maw_ar
    procedure :: bnd_rp => maw_rp
    procedure :: bnd_ad => maw_ad
    procedure :: bnd_cf => maw_cf
    procedure :: bnd_fc => maw_fc
    procedure :: bnd_fn => maw_fn
    procedure :: bnd_nur => maw_nur
    procedure :: bnd_cc => maw_cc
    procedure :: bnd_cq => maw_cq
    procedure :: bnd_ot_model_flows => maw_ot_model_flows
    procedure :: bnd_ot_package_flows => maw_ot_package_flows
    procedure :: bnd_ot_dv => maw_ot_dv
    procedure :: bnd_ot_bdsummary => maw_ot_bdsummary
    procedure :: bnd_da => maw_da
    procedure :: define_listlabel
    ! -- methods for observations
    procedure, public :: bnd_obs_supported => maw_obs_supported
    procedure, public :: bnd_df_obs => maw_df_obs
    procedure, public :: bnd_rp_obs => maw_rp_obs
    procedure, public :: bnd_bd_obs => maw_bd_obs
    ! -- private procedures
    procedure, private :: maw_source_packagedata
    procedure, private :: maw_source_connectiondata
    procedure, private :: maw_source_angledata
    procedure, private :: maw_cell_extent
    procedure, private :: maw_calc_lcorr
    procedure, private :: maw_check_attributes
    procedure, private :: maw_set_attribute_error
    procedure, private :: maw_set_period_value
    procedure, private :: maw_calculate_saturation
    procedure, private :: maw_calculate_satcond
    procedure, private :: maw_calculate_conn_terms
    procedure, private :: maw_calculate_wellq
    procedure, private :: maw_calculate_qpot
    procedure, private :: maw_cfupdate
    procedure, private :: get_jpos
    procedure, private :: get_gwfnode
    ! -- budget
    procedure, private :: maw_setup_budobj
    procedure, private :: maw_fill_budobj
    ! -- table
    procedure, private :: maw_setup_tableobj
    ! -- density
    procedure :: maw_activate_density
    procedure, private :: maw_calculate_density_exchange
    ! -- MAW reduced flow outputs
    procedure, private :: maw_redflow_csv_init
    procedure, private :: maw_redflow_csv_write
    ! -- viscosity
    procedure :: maw_activate_viscosity
  end type mawtype
 
contains
 
!> @brief Create a New Multi-Aquifer Well (MAW) Package
!!
!! After creating the package object point bndobj to the new package
!<
  subroutine maw_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
                        mempath)
    ! -- dummy
    class(bndtype), pointer :: packobj
    integer(I4B), intent(in) :: id
    integer(I4B), intent(in) :: ibcnum
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    character(len=*), intent(in) :: namemodel
    character(len=*), intent(in) :: pakname
    character(len=*), intent(in) :: mempath
    type(mawtype), pointer :: mawobj
    !
    ! -- allocate the object and assign values to object variables
    allocate (mawobj)
    packobj => mawobj
    !
    ! -- create name and memory path
    call packobj%set_names(ibcnum, namemodel, pakname, ftype, mempath)
    packobj%text = text
    !
    ! -- allocate scalars
    call mawobj%maw_allocate_scalars()
    !
    ! -- initialize package
    call packobj%pack_initialize()
    !
    packobj%inunit = inunit
    packobj%iout = iout
    packobj%id = id
    packobj%ibcnum = ibcnum
    packobj%ncolbnd = 4
    packobj%iscloc = 0 ! not supported
    packobj%isadvpak = 1
    packobj%ictMemPath = create_mem_path(namemodel, 'NPF')
  end subroutine maw_create
 
  !> @brief Allocate scalar members
  !<
  subroutine maw_allocate_scalars(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(mawtype), intent(inout) :: this
    !
    ! -- call standard BndExtType allocate scalars
    call this%BndExtType%allocate_scalars()
    !
    ! -- allocate the object and assign values to object variables
    call mem_allocate(this%correct_flow, 'CORRECT_FLOW', this%memoryPath)
    call mem_allocate(this%iprhed, 'IPRHED', this%memoryPath)
    call mem_allocate(this%iheadout, 'IHEADOUT', this%memoryPath)
    call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath)
    call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
    call mem_allocate(this%iflowingwells, 'IFLOWINGWELLS', this%memoryPath)
    call mem_allocate(this%imawiss, 'IMAWISS', this%memoryPath)
    call mem_allocate(this%imawissopt, 'IMAWISSOPT', this%memoryPath)
    call mem_allocate(this%nmawwells, 'NMAWWELLS', this%memoryPath)
    call mem_allocate(this%check_attr, 'CHECK_ATTR', this%memoryPath)
    call mem_allocate(this%ishutoffcnt, 'ISHUTOFFCNT', this%memoryPath)
    call mem_allocate(this%ieffradopt, 'IEFFRADOPT', this%memoryPath)
    call mem_allocate(this%inonvert, 'INONVERT', this%memoryPath)
    call mem_allocate(this%ioutredflowcsv, 'IOUTREDFLOWCSV', this%memoryPath) !for writing reduced MAW flows to csv file
    call mem_allocate(this%satomega, 'SATOMEGA', this%memoryPath)
    call mem_allocate(this%bditems, 'BDITEMS', this%memoryPath)
    call mem_allocate(this%theta, 'THETA', this%memoryPath)
    call mem_allocate(this%kappa, 'KAPPA', this%memoryPath)
    call mem_allocate(this%cbcauxitems, 'CBCAUXITEMS', this%memoryPath)
    call mem_allocate(this%idense, 'IDENSE', this%memoryPath)
    !
    ! -- Set values
    this%correct_flow = .false.
    this%nmawwells = 0
    this%iprhed = 0
    this%iheadout = 0
    this%ibudgetout = 0
    this%ibudcsv = 0
    this%iflowingwells = 0
    this%imawiss = 0
    this%imawissopt = 0
    this%ieffradopt = 0
    this%inonvert = 0
    this%ioutredflowcsv = 0
    this%satomega = dzero
    this%bditems = 8
    this%theta = dp7
    this%kappa = dem4
    this%cbcauxitems = 1
    this%idense = 0
    this%ivsc = 0
  end subroutine maw_allocate_scalars
 
  !> @brief Allocate well arrays
  !<
  subroutine maw_allocate_well_conn_arrays(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    integer(I4B) :: j
    integer(I4B) :: n
    !
    ! -- allocate character array for budget text
    call mem_allocate(this%cmawbudget, lenbudtxt, this%bditems, 'CMAWBUDGET', &
                      this%memoryPath)
    !
    !-- fill cmawbudget
    this%cmawbudget(1) = '             GWF'
    this%cmawbudget(2) = '            RATE'
    this%cmawbudget(3) = '         STORAGE'
    this%cmawbudget(4) = '        CONSTANT'
    this%cmawbudget(5) = '         FW-RATE'
    this%cmawbudget(6) = '        FROM-MVR'
    this%cmawbudget(7) = '     RATE-TO-MVR'
    this%cmawbudget(8) = '  FW-RATE-TO-MVR'
    !
    ! -- allocate character arrays
    call mem_allocate(this%cmawname, lenboundname, this%nmawwells, 'CMAWNAME', &
                      this%memoryPath)
    call mem_allocate(this%status, 8, this%nmawwells, 'STATUS', this%memoryPath)
    !
    ! -- allocate well data pointers in memory manager
    call mem_allocate(this%ngwfnodes, this%nmawwells, 'NGWFNODES', &
                      this%memoryPath)
    call mem_allocate(this%ieqn, this%nmawwells, 'IEQN', this%memoryPath)
    call mem_allocate(this%ishutoff, this%nmawwells, 'ISHUTOFF', this%memoryPath)
    call mem_allocate(this%ifwdischarge, this%nmawwells, 'IFWDISCHARGE', &
                      this%memoryPath)
    call mem_allocate(this%strt, this%nmawwells, 'STRT', this%memoryPath)
    call mem_allocate(this%radius, this%nmawwells, 'RADIUS', this%memoryPath)
    call mem_allocate(this%area, this%nmawwells, 'AREA', this%memoryPath)
    call mem_allocate(this%pumpelev, this%nmawwells, 'PUMPELEV', this%memoryPath)
    call mem_allocate(this%bot, this%nmawwells, 'BOT', this%memoryPath)
    call mem_allocate(this%ratesim, this%nmawwells, 'RATESIM', this%memoryPath)
    call mem_allocate(this%qsim0, this%nmawwells, 'QSIM0', this%memoryPath)
    call mem_allocate(this%reduction_length, this%nmawwells, 'REDUCTION_LENGTH', &
                      this%memoryPath)
    call mem_allocate(this%fwelev, this%nmawwells, 'FWELEV', this%memoryPath)
    call mem_allocate(this%fwcond, this%nmawwells, 'FWCONDS', this%memoryPath)
    call mem_allocate(this%fwrlen, this%nmawwells, 'FWRLEN', this%memoryPath)
    call mem_allocate(this%fwcondsim, this%nmawwells, 'FWCONDSIM', &
                      this%memoryPath)
    call mem_allocate(this%xsto, this%nmawwells, 'XSTO', this%memoryPath)
    call mem_allocate(this%xoldsto, this%nmawwells, 'XOLDSTO', this%memoryPath)
    call mem_allocate(this%shutoffmin, this%nmawwells, 'SHUTOFFMIN', &
                      this%memoryPath)
    call mem_allocate(this%shutoffmax, this%nmawwells, 'SHUTOFFMAX', &
                      this%memoryPath)
    call mem_allocate(this%shutofflevel, this%nmawwells, 'SHUTOFFLEVEL', &
                      this%memoryPath)
    call mem_allocate(this%shutoffweight, this%nmawwells, 'SHUTOFFWEIGHT', &
                      this%memoryPath)
    call mem_allocate(this%shutoffdq, this%nmawwells, 'SHUTOFFDQ', &
                      this%memoryPath)
    call mem_allocate(this%shutoffqold, this%nmawwells, 'SHUTOFFQOLD', &
                      this%memoryPath)
    call mem_allocate(this%nurdxold, this%nmawwells, 'NURDXOLD', &
                      this%memoryPath)
    call mem_allocate(this%nurweight, this%nmawwells, 'NURWEIGHT', &
                      this%memoryPath)
    !
    ! -- alias into the input context's permanent, feature-indexed arrays
    ! (allocated and DZERO-initialized by the loader)
    call mem_setptr(this%rate, 'RATE', this%input_mempath)
    call mem_setptr(this%well_head, 'WELL_HEAD', this%input_mempath)
    call mem_allocate(this%well_head_set, this%nmawwells, 'WELL_HEAD_SET', &
                      this%memoryPath)
    this%well_head_set = .false.
    !
    ! -- allocate and initialize dbuff
    if (this%iheadout > 0) then
      call mem_allocate(this%dbuff, this%nmawwells, 'DBUFF', this%memoryPath)
    else
      call mem_allocate(this%dbuff, 0, 'DBUFF', this%memoryPath)
    end if
    !
    ! -- allocate iaconn
    call mem_allocate(this%iaconn, this%nmawwells + 1, 'IACONN', this%memoryPath)
    !
    ! -- allocate imap
    call mem_allocate(this%imap, this%MAXBOUND, 'IMAP', this%memoryPath)
    !
    ! -- allocate connection data
    call mem_allocate(this%gwfnodes, this%maxbound, 'GWFNODES', this%memoryPath)
    call mem_allocate(this%iboundconn, this%maxbound, 'IBOUNDCONN', &
                      this%memoryPath)
    call mem_allocate(this%sradius, this%maxbound, 'SRADIUS', this%memoryPath)
    call mem_allocate(this%hk, this%maxbound, 'HK', this%memoryPath)
    call mem_allocate(this%satcond, this%maxbound, 'SATCOND', this%memoryPath)
    call mem_allocate(this%simcond, this%maxbound, 'SIMCOND', this%memoryPath)
    call mem_allocate(this%topscrn, this%maxbound, 'TOPSCRN', this%memoryPath)
    call mem_allocate(this%botscrn, this%maxbound, 'BOTSCRN', this%memoryPath)
    call mem_allocate(this%angle, this%maxbound, 'ANGLE', this%memoryPath)
    call mem_allocate(this%connlen, this%maxbound, 'CONNLEN', this%memoryPath)
    call mem_allocate(this%usrtopscrn, this%maxbound, 'USRTOPSCRN', &
                      this%memoryPath)
    call mem_allocate(this%usrbotscrn, this%maxbound, 'USRBOTSCRN', &
                      this%memoryPath)
    !
    ! -- allocate qleak
    call mem_allocate(this%qleak, this%maxbound, 'QLEAK', this%memoryPath)
    !
    ! -- initialize well data
    do n = 1, this%nmawwells
      this%status(n) = 'ACTIVE'
      this%ngwfnodes(n) = 0
      this%ieqn(n) = 0
      this%ishutoff(n) = 0
      this%ifwdischarge(n) = 0
      this%strt(n) = dep20
      this%radius(n) = dep20
      this%area(n) = dzero
      this%pumpelev(n) = dep20
      this%bot(n) = dep20
      this%ratesim(n) = dzero
      this%qsim0(n) = dzero
      this%reduction_length(n) = dep20
      this%fwelev(n) = dzero
      this%fwcond(n) = dzero
      this%fwrlen(n) = dzero
      this%fwcondsim(n) = dzero
      this%xsto(n) = dzero
      this%xoldsto(n) = dzero
      this%shutoffmin(n) = dzero
      this%shutoffmax(n) = dzero
      this%shutofflevel(n) = dep20
      this%shutoffweight(n) = done
      this%shutoffdq(n) = done
      this%shutoffqold(n) = done
      this%nurdxold(n) = dzero
      this%nurweight(n) = done
      !
      ! -- timeseries aware variables
      this%rate(n) = dzero
      this%well_head(n) = dzero
      !
      ! -- dbuff
      if (this%iheadout > 0) then
        this%dbuff(n) = dzero
      end if
    end do
    !
    ! -- initialize iaconn
    do n = 1, this%nmawwells + 1
      this%iaconn(n) = 0
    end do
    !
    ! -- allocate character array for budget text
    call mem_allocate(this%cauxcbc, lenauxname, this%cbcauxitems, 'CAUXCBC', &
                      this%memoryPath)
    !
    ! -- allocate and initialize qauxcbc
    call mem_allocate(this%qauxcbc, this%cbcauxitems, 'QAUXCBC', this%memoryPath)
    do j = 1, this%cbcauxitems
      this%qauxcbc(j) = dzero
    end do
    !
    ! -- allocate flowing well data
    if (this%iflowingwells /= 0) then
      call mem_allocate(this%qfw, this%nmawwells, 'QFW', this%memoryPath)
    else
      call mem_allocate(this%qfw, 1, 'QFW', this%memoryPath)
    end if
    call mem_allocate(this%qout, this%nmawwells, 'QOUT', this%memoryPath)
    call mem_allocate(this%qsto, this%nmawwells, 'QSTO', this%memoryPath)
    call mem_allocate(this%qconst, this%nmawwells, 'QCONST', this%memoryPath)
    !
    ! -- initialize flowing well, storage, and constant flow terms
    do n = 1, this%nmawwells
      if (this%iflowingwells > 0) then
        this%qfw(n) = dzero
      end if
      this%qsto(n) = dzero
      this%qconst(n) = dzero
    end do
    !
    ! -- initialize connection data
    do j = 1, this%maxbound
      this%imap(j) = 0
      this%gwfnodes(j) = 0
      this%iboundconn(j) = 1
      this%sradius(j) = dzero
      this%hk(j) = dzero
      this%satcond(j) = dzero
      this%simcond(j) = dzero
      this%topscrn(j) = dzero
      this%botscrn(j) = dzero
      this%angle(j) = dzero
      this%connlen(j) = dzero
      this%usrtopscrn(j) = dzero
      this%usrbotscrn(j) = dzero
      this%qleak(j) = dzero
    end do
    !
    ! -- allocate denseterms to size 0
    call mem_allocate(this%denseterms, 3, 0, 'DENSETERMS', this%memoryPath)
    !
    ! -- allocate viscratios to size 0
    call mem_allocate(this%viscratios, 2, 0, 'VISCRATIOS', this%memoryPath)
  end subroutine maw_allocate_well_conn_arrays
 
  !> @brief Allocate arrays
  !<
  subroutine maw_allocate_arrays(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    !
    ! -- allocate BndType connection-level arrays
    call this%BndExtType%allocate_arrays()
    !
    ! -- allocate featureauxvar PACKAGEDATA variables
    call this%BndExtType%allocate_featureauxvar()
  end subroutine maw_allocate_arrays
 
  !> @brief Source OPTIONS block from input context
  !<
  subroutine maw_source_options(this)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    type(gwfmawparamfoundtype) :: found
    integer(I4B) :: isize
    character(len=LINELENGTH) :: fname
    ! -- formats
    character(len=*), parameter :: fmtflowingwells = &
      &"(4x, 'FLOWING WELLS WILL BE SIMULATED.')"
    character(len=*), parameter :: fmtshutdown = &
      &"(4x, 'SHUTDOWN ', a, ' VALUE (',g15.7,') SPECIFIED.')"
    character(len=*), parameter :: fmtnostoragewells = &
      &"(4x, 'WELL STORAGE WILL NOT BE SIMULATED.')"
    character(len=*), parameter :: fmtmawbin = &
      "(4x, 'MAW ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, /4x, &
     &'OPENED ON UNIT: ', I0)"
    !
    ! -- call base class to handle AUXILIARY/NAUX, BOUNDNAMES, PRINT_INPUT,
    !    PRINT_FLOWS, SAVE_FLOWS, OBS6_FILENAME, AUXMULTNAME
    call this%BndExtType%source_options()
    !
    write (this%iout, '(1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' OPTIONS'
    !
    ! -- PRINT_HEAD
    call get_isize('PRINT_HEAD', this%input_mempath, isize)
    if (isize > 0) then
      this%iprhed = 1
      write (this%iout, '(4x,a)') &
        trim(adjustl(this%text))//' heads will be printed to listing file.'
    end if
    !
    ! -- HEAD FILEOUT
    call mem_set_value(fname, 'HEADFILE', this%input_mempath, found%headfile)
    if (found%headfile) then
      call assign_iounit(this%iheadout, this%inunit, "HEAD fileout")
      call openfile(this%iheadout, this%iout, fname, 'DATA(BINARY)', &
                    form, access, 'REPLACE', mode_opt=mnormal)
      write (this%iout, fmtmawbin) 'HEAD', trim(adjustl(fname)), this%iheadout
    end if
    !
    ! -- BUDGET FILEOUT
    call mem_set_value(fname, 'BUDGETFILE', this%input_mempath, found%budgetfile)
    if (found%budgetfile) then
      call assign_iounit(this%ibudgetout, this%inunit, "BUDGET fileout")
      call openfile(this%ibudgetout, this%iout, fname, 'DATA(BINARY)', &
                    form, access, 'REPLACE', mode_opt=mnormal)
      write (this%iout, fmtmawbin) 'BUDGET', trim(adjustl(fname)), this%ibudgetout
    end if
    !
    ! -- BUDGETCSV FILEOUT
    call mem_set_value(fname, 'BUDGETCSVFILE', this%input_mempath, &
                       found%budgetcsvfile)
    if (found%budgetcsvfile) then
      call assign_iounit(this%ibudcsv, this%inunit, "BUDGETCSV fileout")
      call openfile(this%ibudcsv, this%iout, fname, 'CSV', filstat_opt='REPLACE')
      write (this%iout, fmtmawbin) &
        'BUDGET CSV', trim(adjustl(fname)), this%ibudcsv
    end if
    !
    ! -- NO_WELL_STORAGE
    call get_isize('NO_WELL_STORAGE', this%input_mempath, isize)
    if (isize > 0) then
      this%imawissopt = 1
      write (this%iout, fmtnostoragewells)
    end if
    !
    ! -- FLOW_CORRECTION
    call get_isize('FLOW_CORRECTION', this%input_mempath, isize)
    if (isize > 0) then
      this%correct_flow = .true.
      write (this%iout, '(4x,a,/,4x,a)') &
        'MAW-GWF FLOW CORRECTIONS WILL BE APPLIED WHEN MAW HEADS ARE BELOW', &
        'OR GWF HEADS IN CONNECTED CELLS ARE BELOW THE CELL BOTTOM.'
    end if
    !
    ! -- FLOWING_WELLS
    call get_isize('FLOWING_WELLS', this%input_mempath, isize)
    if (isize > 0) then
      this%iflowingwells = 1
      write (this%iout, fmtflowingwells)
    end if
    !
    ! -- SHUTDOWN_THETA
    call mem_set_value(this%theta, 'SHUTDOWN_THETA', this%input_mempath, &
                       found%shutdown_theta)
    if (found%shutdown_theta) then
      write (this%iout, fmtshutdown) 'THETA', this%theta
    end if
    !
    ! -- SHUTDOWN_KAPPA
    call mem_set_value(this%kappa, 'SHUTDOWN_KAPPA', this%input_mempath, &
                       found%shutdown_kappa)
    if (found%shutdown_kappa) then
      write (this%iout, fmtshutdown) 'KAPPA', this%kappa
    end if
    !
    ! -- MAW_FLOW_REDUCE_CSV FILEOUT
    call mem_set_value(fname, 'MFRCSVFILE', this%input_mempath, found%mfrcsvfile)
    if (found%mfrcsvfile) then
      call this%maw_redflow_csv_init(fname)
    end if
    !
    ! -- MOVER
    call get_isize('MOVER', this%input_mempath, isize)
    if (isize > 0) then
      this%imover = 1
      write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
    end if
    !
    ! -- DEV_PEACEMAN_EFFECTIVE_RADIUS
    call get_isize('IEFFRADOPT', this%input_mempath, isize)
    if (isize > 0) then
      this%ieffradopt = 1
      write (this%iout, '(4x,a)') &
        'EFFECTIVE RADIUS FOR STRUCTURED GRIDS WILL BE CALCULATED &
        &USING PEACEMAN 1983'
    end if
    !
    ! -- NON_VERTICAL_WELLS
    call get_isize('INONVERT', this%input_mempath, isize)
    if (isize > 0) then
      this%inonvert = 1
      write (this%iout, '(4x,a)') &
        'NON-VERTICAL (SLANTED) WELL CONNECTIONS ENABLED'
    end if
    !
    write (this%iout, '(1x,a)') 'END OF '//trim(adjustl(this%text))//' OPTIONS'
  end subroutine maw_source_options
 
  !> @brief Source DIMENSIONS block from input context
  !<
  subroutine maw_source_dimensions(this)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    type(gwfmawparamfoundtype) :: found
    !
    ! -- initialize dimensions
    this%nmawwells = -1
    this%maxbound = -1
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' DIMENSIONS'
    !
    ! -- read NMAWWELLS
    call mem_set_value(this%nmawwells, 'MAXBOUND', this%input_mempath, &
                       found%maxbound, release=.false.)
    if (found%maxbound) then
      write (this%iout, '(4x,a,i0)') 'NMAWWELLS = ', this%nmawwells
    end if
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' DIMENSIONS'
    !
    if (this%nmawwells < 0) then
      write (errmsg, '(a)') &
        'NMAWWELLS was not specified or was specified incorrectly.'
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- set npakeq to nmawwells
    this%npakeq = this%nmawwells
    !
    ! -- read PACKAGEDATA
    call this%maw_source_packagedata()
    !
    ! -- read CONNECTIONDATA
    call this%maw_source_connectiondata()
    !
    ! -- read optional ANGLEDATA (non-vertical well connections)
    call this%maw_source_angledata()
    !
    ! -- Call define_listlabel
    call this%define_listlabel()
    !
    ! -- setup the budget object
    call this%maw_setup_budobj()
    !
    ! -- setup the head table object
    call this%maw_setup_tableobj()
  end subroutine maw_source_dimensions
 
  !> @brief Bind MawInputType pointers to input context arrays.
  !!
  !! Values are updated from the input context each stress period
  !! including integrated timeseries for supported arrays.
  !<
  subroutine maw_input_init(this, mempath, naux)
    class(mawinputtype), intent(inout) :: this
    character(len=*), intent(in) :: mempath
    integer(I4B), intent(in) :: naux
    call mem_setptr(this%nbound, 'NBOUND', mempath)
    call mem_setptr(this%ifno, 'IFNO', mempath)
    call mem_setptr(this%setting, 'SETTING', mempath)
    call mem_setptr(this%status, 'STATUS', mempath)
    call mem_setptr(this%period_icon, 'PERIOD_ICON', mempath)
    call mem_setptr(this%connstatus, 'CONNSTATUS', mempath)
    call mem_setptr(this%period_bottom, 'PERIOD_BOTTOM', mempath)
    call mem_setptr(this%head_limit, 'HEAD_LIMIT', mempath)
    call mem_setptr(this%fwelev, 'FWELEV', mempath)
    call mem_setptr(this%fwcond, 'FWCOND', mempath)
    call mem_setptr(this%fwrlen, 'FWRLEN', mempath)
    call mem_setptr(this%minrate, 'MINRATE', mempath)
    call mem_setptr(this%maxrate, 'MAXRATE', mempath)
    call mem_setptr(this%pump_elevation, 'PUMP_ELEVATION', mempath)
    call mem_setptr(this%scaling_length, 'SCALING_LENGTH', mempath)
    if (naux > 0) then
      call mem_setptr(this%auxname, 'AUXNAME', mempath)
      call mem_setptr(this%auxval, 'AUXVAL', mempath)
    end if
  end subroutine maw_input_init
 
  !> @brief Nullify all MawInputType pointers.
  !<
  subroutine maw_input_destroy(this)
    class(mawinputtype), intent(inout) :: this
    nullify (this%strt)
    nullify (this%nbound)
    nullify (this%ifno)
    nullify (this%setting)
    nullify (this%status)
    nullify (this%period_icon)
    nullify (this%connstatus)
    nullify (this%period_bottom)
    nullify (this%head_limit)
    nullify (this%fwelev)
    nullify (this%fwcond)
    nullify (this%fwrlen)
    nullify (this%minrate)
    nullify (this%maxrate)
    nullify (this%pump_elevation)
    nullify (this%scaling_length)
    nullify (this%auxname)
    nullify (this%auxval)
  end subroutine maw_input_destroy
 
  !> @brief Source PACKAGEDATA block from input context
  !<
  subroutine maw_source_packagedata(this)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    integer(I4B) :: n, jj, idx, imaw
    character(len=9) :: cno9
    character(len=LINELENGTH) :: keyword
    real(DP) :: rval
    integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
    real(DP), dimension(:), pointer, contiguous :: radius => null()
    real(DP), dimension(:), pointer, contiguous :: bottom => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: condeqn => null()
    integer(I4B), dimension(:), pointer, contiguous :: ngwfnodes => null()
    type(characterstringtype), dimension(:), pointer, &
      contiguous :: boundname => null()
    integer(I4B), allocatable :: nboundchk(:)
    ! -- formats
    character(len=*), parameter :: fmthdbot = &
      "('well head (', G0, ') must be greater than or equal to the &
      &BOTTOM_ELEVATION (', G0, ').')"
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' PACKAGEDATA'
    !
    ! -- get input context PACKAGEDATA column arrays
    call mem_setptr(ifno, 'PACKAGEDATA_IFNO', this%input_mempath)
    call mem_setptr(radius, 'RADIUS', this%input_mempath)
    call mem_setptr(bottom, 'BOTTOM', this%input_mempath)
    call mem_setptr(this%input%strt, 'STRT', this%input_mempath)
    call mem_setptr(condeqn, 'CONDEQN', this%input_mempath)
    call mem_setptr(ngwfnodes, 'NGWFNODES', this%input_mempath)
    if (this%inamedbound /= 0) then
      call mem_setptr(boundname, 'BOUNDNAME', this%input_mempath)
    end if
    !
    ! -- total GWF connections is the sum of per-well node counts;
    ! -- allocate well and connection arrays before filling them
    this%maxbound = sum(ngwfnodes)
    write (this%iout, '(//4x,a,i7)') 'MAXBOUND = ', this%maxbound
    call this%maw_allocate_well_conn_arrays()
    !
    ! -- validate and fill allocated arrays in feature-number order
    allocate (nboundchk(this%nmawwells))
    nboundchk = 0
    do n = 1, this%nmawwells
      imaw = this%validate_ifno(ifno(n), this%nmawwells, nboundchk, &
                                'IMAW', 'PACKAGEDATA')
      if (imaw == 0) cycle
      !
      ! -- radius
      rval = radius(n)
      if (rval <= dzero) then
        write (errmsg, '(a,1x,i0,1x,a)') &
          'Radius for well', imaw, 'must be greater than zero.'
        call store_error(errmsg)
      end if
      this%radius(imaw) = rval
      this%area(imaw) = dpi * rval**dtwo
      !
      ! -- well bottom
      this%bot(imaw) = bottom(n)
      !
      ! -- conductance equation
      keyword = condeqn(n)
      select case (trim(keyword))
      case ('SPECIFIED')
        this%ieqn(imaw) = 0
      case ('THIEM')
        this%ieqn(imaw) = 1
      case ('THEIM') ! # codespell:ignore
        this%ieqn(imaw) = 1
        write (warnmsg, '(a,a,a,a,a,a)') &
          "CONDEQN in '", trim(this%packName), "' should be ", &
          "corrected from '", trim(keyword), "' to 'THIEM'."
        call store_warning(warnmsg)
      case ('SKIN')
        this%ieqn(imaw) = 2
      case ('CUMULATIVE')
        this%ieqn(imaw) = 3
      case ('MEAN')
        this%ieqn(imaw) = 4
      case default
        write (errmsg, '(a,1x,i0,1x,a)') &
          'CONDEQN for well', imaw, &
          "must be 'CUMULATIVE', 'THIEM', 'MEAN', or 'SKIN'."
        call store_error(errmsg)
        this%ieqn(imaw) = 0
      end select
      !
      ! -- ngwfnodes
      if (ngwfnodes(n) < 1) then
        write (errmsg, '(a,1x,i0,1x,a)') &
          'NGWFNODES for well', imaw, 'must be greater than zero.'
        call store_error(errmsg)
      end if
      this%ngwfnodes(imaw) = ngwfnodes(n)
      !
      ! -- read the TS integrated numeric value directly
      this%well_head(imaw) = this%input%strt(n)
      this%strt(imaw) = this%input%strt(n)
      if (this%strt(imaw) < this%bot(imaw)) then
        write (keyword, fmthdbot) this%strt(imaw), this%bot(imaw)
        call this%maw_set_attribute_error(imaw, 'STRT', trim(keyword))
      end if
      !
      ! -- well name / boundname
      write (cno9, '(i9.9)') imaw
      this%cmawname(imaw) = 'MAWWELL'//cno9
      if (this%inamedbound /= 0) then
        text = boundname(n)
        if (trim(text) /= '') then
          this%cmawname(imaw) = trim(text)
        end if
      end if
    end do
    !
    ! -- check for duplicate or missing wells
    call this%report_ifno_coverage(nboundchk, this%nmawwells, 'well', &
                                   'PACKAGEDATA')
    deallocate (nboundchk)
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' PACKAGEDATA'
    !
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- set iaconn and imap for each connection
    idx = 0
    this%iaconn(1) = 1
    do n = 1, this%nmawwells
      do jj = 1, this%ngwfnodes(n)
        idx = idx + 1
        this%imap(idx) = n
      end do
      this%iaconn(n + 1) = idx + 1
    end do
    !
    ! -- release input context packagedata memory;
    !    exclude STRT and PACKAGEDATA_IFNO
    call memorystore_release('RADIUS', this%input_mempath)
    call memorystore_release('BOTTOM', this%input_mempath)
    call memorystore_release('CONDEQN', this%input_mempath)
    call memorystore_release('NGWFNODES', this%input_mempath)
    call memorystore_release('BOUNDNAME', this%input_mempath)
  end subroutine maw_source_packagedata
 
  !> @brief Source CONNECTIONDATA block from input context
  !<
  subroutine maw_source_connectiondata(this)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    integer(I4B) :: n, nconn, imaw, iconn, ipos, jpos, j, jj
    integer(I4B) :: nn, nn2
    integer(I4B) :: ireset_scrntop, ireset_scrnbot, ireset_wellbot
    real(DP) :: rval, topnn, botnn, botw
    character(len=LINELENGTH) :: cellid_str
    character(len=30) :: nodestr
    integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
    integer(I4B), dimension(:), pointer, contiguous :: icon => null()
    integer(I4B), dimension(:, :), pointer, contiguous :: cellid => null()
    real(DP), dimension(:), pointer, contiguous :: scrn_top => null()
    real(DP), dimension(:), pointer, contiguous :: scrn_bot => null()
    real(DP), dimension(:), pointer, contiguous :: hk_skin => null()
    real(DP), dimension(:), pointer, contiguous :: radius_skin => null()
    integer(I4B), allocatable :: iachk(:)
    integer(I4B), allocatable :: nboundchk(:)
    !
    ireset_scrntop = 0
    ireset_scrnbot = 0
    ireset_wellbot = 0
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' CONNECTIONDATA'
    !
    ! -- get input context CONNECTIONDATA columns
    call mem_setptr(ifno, 'CONNDATA_IFNO', this%input_mempath)
    call mem_setptr(icon, 'ICON', this%input_mempath)
    call mem_setptr(cellid, 'CELLID', this%input_mempath)
    call mem_setptr(scrn_top, 'SCRN_TOP', this%input_mempath)
    call mem_setptr(scrn_bot, 'SCRN_BOT', this%input_mempath)
    call mem_setptr(hk_skin, 'HK_SKIN', this%input_mempath)
    call mem_setptr(radius_skin, 'RADIUS_SKIN', this%input_mempath)
    !
    nconn = size(ifno)
    !
    ! -- allocate check arrays
    allocate (iachk(this%nmawwells + 1))
    iachk(1) = 1
    do n = 1, this%nmawwells
      iachk(n + 1) = iachk(n) + this%ngwfnodes(n)
    end do
    allocate (nboundchk(this%maxbound))
    nboundchk = 0
    !
    do n = 1, nconn
      imaw = ifno(n)
      iconn = icon(n)
      !
      ! -- check for error condition
      if (imaw < 1 .or. imaw > this%nmawwells) then
        write (errmsg, '(a,1x,i0,a)') &
          'IMAW must be greater than 0 and less than or equal to ', &
          this%nmawwells, '.'
        call store_error(errmsg)
        cycle
      end if
      if (iconn < 1 .or. iconn > this%ngwfnodes(imaw)) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0,a)') &
          'JCONN for well ', imaw, &
          'must be greater than 1 and less than or equal to ', &
          this%ngwfnodes(imaw), '.'
        call store_error(errmsg)
        cycle
      end if
      !
      ipos = iachk(imaw) + iconn - 1
      nboundchk(ipos) = nboundchk(ipos) + 1
      j = iconn
      jpos = this%get_jpos(imaw, iconn)
      !
      ! -- convert cellid integer array to string and get node number
      write (cellid_str, '(10(i0,1x))') cellid(1:this%dis%ndim, n)
      nn = this%dis%noder_from_cellid(cellid_str, this%inunit, this%iout)
      topnn = this%dis%top(nn)
      botnn = this%dis%bot(nn)
      botw = this%bot(imaw)
      !
      ! -- set gwf node number for connection
      this%gwfnodes(jpos) = nn
      !
      ! -- top of screen; retain user value before clipping for non-vertical SPECIFIED connections
      this%usrtopscrn(jpos) = scrn_top(n)
      rval = scrn_top(n)
      if (this%ieqn(imaw) /= 4) then
        rval = topnn
      else
        if (rval > topnn) then
          ireset_scrntop = ireset_scrntop + 1
          rval = topnn
        end if
      end if
      this%topscrn(jpos) = rval
      !
      ! -- bottom of screen; retain user value before clipping for non-vertical SPECIFIED connections
      this%usrbotscrn(jpos) = scrn_bot(n)
      rval = scrn_bot(n)
      if (this%ieqn(imaw) /= 4) then
        rval = botnn
      else
        if (rval < botnn) then
          ireset_scrnbot = ireset_scrnbot + 1
          rval = botnn
        end if
      end if
      this%botscrn(jpos) = rval
      !
      ! -- adjust the bottom of the well for all conductance approaches
      !    except for "mean"
      if (rval < botw) then
        if (this%ieqn(imaw) /= 4) then
          ireset_wellbot = ireset_wellbot + 1
          botw = rval
          this%bot(imaw) = rval
        else
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,g0,a,g0,a)') &
            'Screen bottom for maw well', imaw, 'connection', j, '(', &
            this%botscrn(jpos), ') is less than the well bottom (', &
            this%bot(imaw), ').'
          call store_error(errmsg)
        end if
      end if
      !
      ! -- hydraulic conductivity or conductance
      rval = hk_skin(n)
      if (this%ieqn(imaw) == 0) then
        this%satcond(jpos) = rval
      else if (this%ieqn(imaw) == 2 .OR. this%ieqn(imaw) == 3 .OR. &
               this%ieqn(imaw) == 4) then
        this%hk(jpos) = rval
      end if
      !
      ! -- skin radius
      rval = radius_skin(n)
      if (this%ieqn(imaw) == 2 .OR. this%ieqn(imaw) == 3 .OR. &
          this%ieqn(imaw) == 4) then
        this%sradius(jpos) = rval
        if (this%sradius(jpos) <= this%radius(imaw)) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,g0,a,g0,a)') &
            'Screen radius for maw well', imaw, 'connection', j, '(', &
            this%sradius(jpos), &
            ') is less than or equal to the well radius (', &
            this%radius(imaw), ').'
          call store_error(errmsg)
        end if
      end if
    end do
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' CONNECTIONDATA'
    !
    ! -- check for missing or duplicate connections
    ipos = 0
    do n = 1, this%nmawwells
      do j = 1, this%ngwfnodes(n)
        ipos = ipos + 1
        if (nboundchk(ipos) == 0) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,a)') &
            'No data specified for maw well', n, 'connection', j, '.'
          call store_error(errmsg)
        else if (nboundchk(ipos) > 1) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
            'Data for maw well', n, 'connection', j, &
            'specified', nboundchk(ipos), 'times.'
          call store_error(errmsg)
        end if
      end do
    end do
    !
    ! -- make sure that more than one connection per cell is only specified
    !    for wells using the mean conductance type
    do n = 1, this%nmawwells
      if (this%ieqn(n) /= 4) then
        do j = 1, this%ngwfnodes(n)
          nn = this%get_gwfnode(n, j)
          do jj = 1, this%ngwfnodes(n)
            if (jj == j) cycle
            nn2 = this%get_gwfnode(n, jj)
            if (nn2 == nn) then
              call this%dis%noder_to_string(nn, nodestr)
              write (errmsg, '(a,1x,i0,1x,a,1x,i0,3(1x,a))') &
                'Only one connection can be specified for maw well', &
                n, 'connection', j, 'to gwf cell', trim(adjustl(nodestr)), &
                'unless the mean condeqn is specified.'
              call store_error(errmsg)
            end if
          end do
        end do
      end if
    end do
    !
    ! -- deallocate local storage
    deallocate (iachk)
    deallocate (nboundchk)
    !
    ! -- add warning messages
    if (ireset_scrntop > 0) then
      write (warnmsg, '(a,1x,a,1x,a,1x,i0,1x,a)') &
        'The screen tops in multi-aquifer well package', trim(this%packName), &
        'were reset to the top of the connected cell', ireset_scrntop, 'times.'
      call store_warning(warnmsg)
    end if
    if (ireset_scrnbot > 0) then
      write (warnmsg, '(a,1x,a,1x,a,1x,i0,1x,a)') &
        'The screen bottoms in multi-aquifer well package', trim(this%packName), &
        'were reset to the bottom of the connected cell', ireset_scrnbot, &
        'times.'
      call store_warning(warnmsg)
    end if
    if (ireset_wellbot > 0) then
      write (warnmsg, '(a,1x,a,1x,a,1x,i0,1x,a)') &
        'The well bottoms in multi-aquifer well package', trim(this%packName), &
        'were reset to the bottom of the connected cell', ireset_wellbot, &
        'times.'
      call store_warning(warnmsg)
    end if
    !
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- release input context connectiondata memory
    call memorystore_release('CONNDATA_IFNO', this%input_mempath)
    call memorystore_release('ICON', this%input_mempath)
    call memorystore_release('CELLID', this%input_mempath)
    call memorystore_release('SCRN_TOP', this%input_mempath)
    call memorystore_release('SCRN_BOT', this%input_mempath)
    call memorystore_release('HK_SKIN', this%input_mempath)
    call memorystore_release('RADIUS_SKIN', this%input_mempath)
  end subroutine maw_source_connectiondata
 
  !> @brief Source ANGLEDATA block from input context for non-vertical MAW connections
  !<
  subroutine maw_source_angledata(this)
    use constantsmodule, only: dzero, dem6, dtwo, dpio180, dnodata
    use mathutilmodule, only: is_close
    use memorymanagerextmodule, only: memorystore_release
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    integer(I4B) :: n, imaw, iconn, jpos, ichk, nangle, node
    real(DP) :: angle, dz, omega, lw, topexp
    real(DP) :: hlen
    real(DP) :: extent
    character(len=LINELENGTH) :: cndmsg
    character(len=LINELENGTH) :: extmsg
    logical(LGP) :: estimated
    integer(I4B), dimension(:), pointer, contiguous :: angledata_ifno => null()
    integer(I4B), dimension(:), pointer, contiguous :: angledata_icon => null()
    real(DP), dimension(:), pointer, contiguous :: angle_arr => null()
    real(DP), dimension(:), pointer, contiguous :: connlen_arr => null()
    integer(I4B), allocatable :: iachk(:)
    integer(I4B), allocatable :: nboundchk(:)
    real(DP), parameter :: coszero = dem6
    real(DP), parameter :: dninety = 9.0d1
    !
    call get_isize('ANGLEDATA_IFNO', this%input_mempath, nangle)
    !
    if (nangle > 0) then
      !
      if (this%inonvert == 0) then
        call store_error('An ANGLEDATA block was specified but the '// &
                         'NON_VERTICAL_WELLS option was not specified in the '// &
                         'OPTIONS block.')
        call store_error_filename(this%input_fname)
      end if
      !
      call mem_setptr(angledata_ifno, 'ANGLEDATA_IFNO', this%input_mempath)
      call mem_setptr(angledata_icon, 'ANGLEDATA_ICON', this%input_mempath)
      call mem_setptr(angle_arr, 'ANGLE', this%input_mempath)
      call mem_setptr(connlen_arr, 'CONNLEN', this%input_mempath)
      !
      allocate (iachk(this%nmawwells + 1))
      iachk(1) = 1
      do n = 1, this%nmawwells
        iachk(n + 1) = iachk(n) + this%ngwfnodes(n)
      end do
      allocate (nboundchk(this%maxbound))
      nboundchk = 0
      !
      write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
        ' ANGLEDATA'
      !
      do n = 1, nangle
        imaw = angledata_ifno(n)
        if (imaw < 1 .or. imaw > this%nmawwells) then
          write (errmsg, '(a,1x,i0,a)') &
            'IFNO must be greater than 0 and less than or equal to ', &
            this%nmawwells, '.'
          call store_error(errmsg)
          cycle
        end if
        !
        iconn = angledata_icon(n)
        if (iconn < 1 .or. iconn > this%ngwfnodes(imaw)) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,a)') &
            'ICON for well ', imaw, &
            'must be greater than 0 and less than or equal to ', &
            this%ngwfnodes(imaw), '.'
          call store_error(errmsg)
          cycle
        end if
        jpos = this%get_jpos(imaw, iconn)
        !
        ichk = iachk(imaw) + iconn - 1
        nboundchk(ichk) = nboundchk(ichk) + 1
        if (nboundchk(ichk) > 1) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
            'ANGLEDATA for maw well', imaw, 'connection', iconn, &
            'is specified more than once.'
          call store_error(errmsg)
        end if
        !
        angle = angle_arr(n)
        this%angle(jpos) = angle
        if (connlen_arr(n) /= dnodata) then
          this%connlen(jpos) = max(connlen_arr(n), dzero)
        end if
        !
        if (angle < dzero .or. angle > dninety) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,g0,a)') &
            'ANGLE for maw well', imaw, 'connection', iconn, '(', angle, &
            ') must be greater than or equal to 0.0 and less than or '// &
            'equal to 90.0 degrees.'
          call store_error(errmsg)
          cycle
        end if
        !
        node = this%get_gwfnode(imaw, iconn)
        !
        ! -- a SPECIFIED connection provides the saturated conductance directly,
        !    so the length correction is not applied; the user-specified screen
        !    elevations (which are otherwise reset to the cell top and bottom)
        !    are restored here, clamped to the cell, so the connection
        !    saturation is calculated over the correct interval. It is the
        !    user's responsibility to calculate the correct conductance.
        if (this%ieqn(imaw) == 0) then
          this%topscrn(jpos) = min(this%usrtopscrn(jpos), this%dis%top(node))
          this%botscrn(jpos) = max(this%usrbotscrn(jpos), this%dis%bot(node))
        end if
        !
        dz = this%topscrn(jpos) - this%botscrn(jpos)
        if (dz <= dzero) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
            'The screen top must be greater than the screen bottom for maw '// &
            'well', imaw, 'connection', iconn, 'listed in the ANGLEDATA block.'
          call store_error(errmsg)
          cycle
        end if
        !
        omega = angle * dpio180
        !
        ! -- the in-cell screen length is specified directly or derived from
        !    the screen elevations, well radius, and tilt angle; it cannot be
        !    derived for a (near) horizontal connection
        if (this%connlen(jpos) > dzero) then
          lw = this%connlen(jpos)
        else if (cos(omega) > coszero) then
          lw = maw_screen_length(dz, this%radius(imaw), omega)
        else
          lw = dzero
        end if
        !
        ! -- the in-cell screen length is used only by the conductance
        !    equations calculated by the program (it is not used by the
        !    SPECIFIED equation). For those equations a (near) horizontal
        !    connection requires a connection length, because the length cannot
        !    be derived from the screen elevations, and the calculated length
        !    must be positive.
        if (this%ieqn(imaw) /= 0) then
          if (cos(omega) <= coszero .and. this%connlen(jpos) <= dzero) then
            write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
              'A connection length must be specified for the (near) '// &
              'horizontal maw well', imaw, 'connection', iconn, &
              'listed in the ANGLEDATA block.'
            call store_error(errmsg)
            cycle
          end if
          if (lw <= dzero) then
            write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,g0,a)') &
              'The calculated in-cell screen length for maw well', imaw, &
              'connection', iconn, '(', lw, &
              ') is not greater than zero. Specify a connection length in '// &
              'the ANGLEDATA block.'
            call store_error(errmsg)
            cycle
          end if
        end if
        !
        ! -- a connection spans a horizontal distance of lw * sin(omega), which
        !    cannot exceed the horizontal extent of the connected cell
        if (lw > dzero) then
          hlen = lw * sin(omega)
          extent = this%maw_cell_extent(node, estimated)
          if (hlen > extent) then
            !
            ! -- the extent is estimated from the cell area where the cell
            !    vertices are not defined
            if (estimated) then
              extmsg = 'estimated maximum horizontal extent'
            else
              extmsg = 'maximum horizontal extent'
            end if
            !
            ! -- the length correction is not applied to a SPECIFIED
            !    connection, so the conductance of one is not the program's
            !    to describe as too large
            if (this%ieqn(imaw) == 0) then
              cndmsg = 'The specified saturated conductance is applied to a '// &
                       'screen that is longer than the cell.'
            else
              cndmsg = 'The calculated saturated conductance is '// &
                       'correspondingly too large.'
            end if
            write (warnmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,g0,a,g0,a)') &
              'The horizontal distance spanned by maw well', imaw, &
              'connection', iconn, '(', hlen, &
              ') is greater than the '//trim(extmsg)//' of the '// &
              'connected cell (', extent, &
              '), so the screen extends beyond the cell it is connected '// &
              'to. '//trim(cndmsg)//' Reduce ANGLE or CONN_LENGTH, or '// &
              'specify a separate connection to each cell the connection '// &
              'passes through.'
            call store_warning(warnmsg)
          end if
        end if
        !
        ! -- the radial (THIEM, SKIN, and CUMULATIVE) conductance equations use
        !    the horizontal-plane flow geometry, so a (near) horizontal
        !    connection is an approximation; recommend the MEAN equation
        if (cos(omega) <= coszero .and. &
            this%ieqn(imaw) /= 4 .and. this%ieqn(imaw) /= 0) then
          write (warnmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
            'The (near) horizontal maw well', imaw, 'connection', iconn, &
            'uses a radial conductance equation (THIEM, SKIN, or '// &
            'CUMULATIVE). The calculated conductance is an approximation '// &
            'for horizontal connections; the MEAN conductance equation is '// &
            'recommended for horizontal connections.'
          call store_warning(warnmsg)
        end if
        !
        if (cos(omega) <= coszero .and. &
            (this%ieqn(imaw) == 4 .or. this%ieqn(imaw) == 0)) then
          topexp = this%botscrn(jpos) + dtwo * this%radius(imaw)
          if (is_close(this%topscrn(jpos), topexp)) then
            this%topscrn(jpos) = topexp
          else
            write (warnmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,g0,a)') &
              'The vertical screen extent (SCRN_TOP - SCRN_BOT) for the '// &
              '(near) horizontal maw well', imaw, 'connection', iconn, &
              'is not equal to the well diameter (2 * RADIUS = ', &
              dtwo * this%radius(imaw), &
              '). The vertical screen extent is used to determine the '// &
              'saturation of a connection and should equal the well '// &
              'diameter for a horizontal connection.'
            call store_warning(warnmsg)
          end if
        end if
      end do
      !
      write (this%iout, '(1x,a)') &
        'END OF '//trim(adjustl(this%text))//' ANGLEDATA'
      !
      deallocate (iachk)
      deallocate (nboundchk)
      !
      call memorystore_release('ANGLEDATA_IFNO', this%input_mempath)
      call memorystore_release('ANGLEDATA_ICON', this%input_mempath)
      call memorystore_release('ANGLE', this%input_mempath)
      call memorystore_release('CONNLEN', this%input_mempath)
    else
      if (this%inonvert /= 0) then
        write (warnmsg, '(a)') &
          'The NON_VERTICAL_WELLS option was specified but an ANGLEDATA '// &
          'block was not found. All multi-aquifer well connections will be '// &
          'treated as vertical.'
        call store_warning(warnmsg)
      end if
    end if
    !
    if (count_errors() > 0) call store_error_filename(this%input_fname)
  end subroutine maw_source_angledata
 
  !> @brief In-cell screen length of a non-vertical well connection
  !!
  !! The length is derived from the vertical screen extent, the well radius,
  !! and the tilt angle, and grows without bound as the connection approaches
  !! horizontal.
  !<
  pure function maw_screen_length(dz, radius, omega) result(lw)
    ! -- dummy
    real(dp), intent(in) :: dz !< vertical screen extent
    real(dp), intent(in) :: radius !< well radius
    real(dp), intent(in) :: omega !< tilt angle from vertical, in radians
    ! -- return
    real(dp) :: lw
    !
    lw = (dz - dtwo * radius * sin(omega)) / cos(omega)
  end function maw_screen_length
 
  !> @brief Maximum horizontal extent of a cell
  !!
  !! The extent is the largest distance between two cell vertices where they
  !! are defined, and is estimated as the diagonal of a square with the same
  !! area where they are not.
  !<
  function maw_cell_extent(this, node, estimated) result(extent)
    use constantsmodule, only: dis, disv, disu, dtwo
    use geomutilmodule, only: polygon_extent
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: node !< reduced node number of the connected cell
    logical(LGP), intent(out) :: estimated !< extent is estimated from the cell area
    ! -- return
    real(dp) :: extent
    ! -- local
    real(dp), allocatable, dimension(:, :) :: polyverts
    integer(I4B) :: nverts
    !
    ! -- the vertices are optional for a disu grid and are not defined for
    !    grids that have no cell polygons
    nverts = 0
    select case (this%dis%get_dis_enum())
    case (dis, disv, disu)
      nverts = this%dis%get_npolyverts(node)
    end select
    !
    estimated = nverts == 0
    if (estimated) then
      extent = sqrt(dtwo * this%dis%area(node))
    else
      call this%dis%get_polyverts(node, polyverts)
      extent = polygon_extent(polyverts(1, :), polyverts(2, :))
      deallocate (polyverts)
    end if
  end function maw_cell_extent
 
  !> @brief Calculate the length correction factor for a multi-aquifer well
  !! connection
  !!
  !! The saturated conductance of a non-vertical (slanted) connection is scaled
  !! by the ratio of the in-cell screen length to the vertical screen thickness.
  !! The factor is 1.0 for vertical connections (angle = 0 and no connection
  !! length specified). The in-cell screen length is either specified directly
  !! (CONNLEN > 0) or derived from the screen elevations and the tilt angle,
  !! accounting for the vertical band occupied by the finite-radius borehole.
  !<
  function maw_calc_lcorr(this, i, jpos) result(lcorr)
    use constantsmodule, only: dzero, done, dtwo, dpio180
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: i !< well number
    integer(I4B), intent(in) :: jpos !< connection position
    ! -- return
    real(dp) :: lcorr
    ! -- local
    real(dp) :: dz
    real(dp) :: omega
    real(dp) :: lw
    !
    ! -- default to no correction (vertical connection)
    lcorr = done
    if (this%angle(jpos) == dzero .and. this%connlen(jpos) <= dzero) then
      return
    end if
    !
    ! -- vertical screen thickness
    dz = this%topscrn(jpos) - this%botscrn(jpos)
    if (dz <= dzero) then
      return
    end if
    !
    ! -- in-cell screen length: specified directly or derived from the screen
    !    elevations and the tilt angle
    if (this%connlen(jpos) > dzero) then
      lw = this%connlen(jpos)
    else
      omega = this%angle(jpos) * dpio180
      lw = maw_screen_length(dz, this%radius(i), omega)
    end if
    !
    lcorr = lw / dz
  end function maw_calc_lcorr
 
  !> @brief Read the initial parameters for this package
  !<
  subroutine maw_read_initial_attr(this)
    ! -- modules
    use constantsmodule, only: linelength
    use memorymanagermodule, only: mem_setptr
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    integer(I4B) :: ntabcols
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: nn
    integer(I4B) :: jpos
    integer(I4B) :: inode
    integer(I4B) :: idx
    real(DP) :: k11
    real(DP) :: k22
    character(len=10), dimension(0:4) :: ccond
    character(len=30) :: nodestr
    ! -- data
    data ccond(0)/'SPECIFIED '/
    data ccond(1)/'THIEM     '/
    data ccond(2)/'SKIN      '/
    data ccond(3)/'CUMULATIVE'/
    data ccond(4)/'MEAN      '/
    ! -- format
    character(len=*), parameter :: fmtwelln = &
      "(1X,//43X,'MULTI-AQUIFER WELL DATA'&
      &/1X,109('-'),&
      &/1X,7(A10,1X),A16)"
    character(len=*), parameter :: fmtwelld = &
      &"(1X,I10,1X,4(G10.3,1X),I10,1X,A10,1X,A16)"
    character(len=*), parameter :: fmtline = &
      &"(1X,119('-'),//)"
    character(len=*), parameter :: fmtwellcn = &
      "(1X,//37X,'MULTI-AQUIFER WELL CONNECTION DATA'&
      &/1X,119('-'),&
      &/1X,2(A10,1X),A20,7(A10,1X))"
    character(len=*), parameter :: fmtwellcd = &
      &"(1X,2(I10,1X),A20,1X,2(G10.3,1X),2(A10,1X),3(G10.3,1X))"
    !
    ! -- initialize xnewpak
    do n = 1, this%nmawwells
      this%xnewpak(n) = this%strt(n)
      this%xsto(n) = this%strt(n)
    end do
    !
    ! -- initialize status (iboundpak) of maw wells to active
    do n = 1, this%nmawwells
      select case (this%status(n))
      case ('CONSTANT')
        this%iboundpak(n) = -1
      case ('INACTIVE')
        this%iboundpak(n) = 0
      case ('ACTIVE')
        this%iboundpak(n) = 1
      end select
    end do
    !
    ! -- set imap and boundname for each connection
    if (this%inamedbound /= 0) then
      idx = 0
      do n = 1, this%nmawwells
        do j = 1, this%ngwfnodes(n)
          idx = idx + 1
          this%boundname(idx) = this%cmawname(n)
          this%imap(idx) = n
        end do
      end do
    else
      do n = 1, this%nmawwells
        this%cmawname(n) = ''
      end do
    end if
    !
    ! -- copy boundname into boundname_cst
    call this%copy_boundname()
    !
    ! -- set pointer to gwf iss and gwf hk
    call mem_setptr(this%gwfiss, 'ISS', create_mem_path(this%name_model))
    call mem_setptr(this%gwfk11, 'K11', create_mem_path(this%name_model, 'NPF'))
    call mem_setptr(this%gwfk22, 'K22', create_mem_path(this%name_model, 'NPF'))
    call mem_setptr(this%gwfik22, 'IK22', create_mem_path(this%name_model, 'NPF'))
    call mem_setptr(this%gwfsat, 'SAT', create_mem_path(this%name_model, 'NPF'))
    !
    ! -- qa data
    call this%maw_check_attributes()
    !
    ! -- Calculate the saturated conductance
    do n = 1, this%nmawwells
      !
      ! -- calculate saturated conductance only if CONDUCTANCE was not
      !    specified for each maw-gwf connection (CONDUCTANCE keyword).
      do j = 1, this%ngwfnodes(n)
        if (this%ieqn(n) /= 0) then
          inode = this%get_gwfnode(n, j)
          call this%maw_calculate_satcond(n, j, inode)
        end if
      end do
    end do
    !
    ! -- write summary of static well data
    ! -- write well data
    if (this%iprpak /= 0) then
      ntabcols = 7
      if (this%inamedbound /= 0) then
        ntabcols = ntabcols + 1
      end if
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') STATIC WELL DATA'
      call table_cr(this%inputtab, this%packName, title)
      call this%inputtab%table_df(this%nmawwells, ntabcols, this%iout)
      text = 'NUMBER'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'RADIUS'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'AREA'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'WELL BOTTOM'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'STARTING HEAD'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'NUMBER OF GWF NODES'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'CONDUCT. EQUATION'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      if (this%inamedbound /= 0) then
        text = 'NAME'
        call this%inputtab%initialize_column(text, 20, alignment=tableft)
      end if
      do n = 1, this%nmawwells
        call this%inputtab%add_term(n)
        call this%inputtab%add_term(this%radius(n))
        call this%inputtab%add_term(this%area(n))
        call this%inputtab%add_term(this%bot(n))
        call this%inputtab%add_term(this%strt(n))
        call this%inputtab%add_term(this%ngwfnodes(n))
        call this%inputtab%add_term(ccond(this%ieqn(n)))
        if (this%inamedbound /= 0) then
          call this%inputtab%add_term(this%cmawname(n))
        end if
      end do
    end if
    !
    ! -- write well connection data
    if (this%iprpak /= 0) then
      ntabcols = 10
      if (this%inonvert /= 0) then
        ntabcols = ntabcols + 1
      end if
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') STATIC WELL CONNECTION DATA'
      call table_cr(this%inputtab, this%packName, title)
      call this%inputtab%table_df(this%maxbound, ntabcols, this%iout)
      text = 'NUMBER'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'WELL CONNECTION'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'CELLID'
      call this%inputtab%initialize_column(text, 20, alignment=tableft)
      text = 'TOP OF SCREEN'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'BOTTOM OF SCREEN'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'SKIN RADIUS'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'SKIN K'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'K11'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'K22'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'SATURATED WELL CONDUCT.'
      call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      if (this%inonvert /= 0) then
        text = 'ANGLE (DEG)'
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
      end if
      !
      ! -- write the data to the table
      do n = 1, this%nmawwells
        do j = 1, this%ngwfnodes(n)
          call this%inputtab%add_term(n)
          call this%inputtab%add_term(j)
          jpos = this%get_jpos(n, j)
          nn = this%get_gwfnode(n, j)
          call this%dis%noder_to_string(nn, nodestr)
          call this%inputtab%add_term(nodestr)
          call this%inputtab%add_term(this%topscrn(jpos))
          call this%inputtab%add_term(this%botscrn(jpos))
          if (this%ieqn(n) == 2 .or. &
              this%ieqn(n) == 3 .or. &
              this%ieqn(n) == 4) then
            call this%inputtab%add_term(this%sradius(jpos))
            call this%inputtab%add_term(this%hk(jpos))
          else
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
          end if
          if (this%ieqn(n) == 1 .or. &
              this%ieqn(n) == 2 .or. &
              this%ieqn(n) == 3) then
            k11 = this%gwfk11(nn)
            if (this%gwfik22 == 0) then
              k22 = this%gwfk11(nn)
            else
              k22 = this%gwfk22(nn)
            end if
            call this%inputtab%add_term(k11)
            call this%inputtab%add_term(k22)
          else
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
          end if
          call this%inputtab%add_term(this%satcond(jpos))
          if (this%inonvert /= 0) then
            call this%inputtab%add_term(this%angle(jpos))
          end if
        end do
      end do
    end if
    !
    ! -- finished with pointer to gwf hydraulic conductivity
    this%gwfk11 => null()
    this%gwfk22 => null()
    this%gwfik22 => null()
    this%gwfsat => null()
    !
    ! -- check for any error conditions
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
  end subroutine maw_read_initial_attr
 
  !> @brief Issue a parameter error for mawweslls(imaw)
  !<
  subroutine maw_set_attribute_error(this, imaw, keyword, msg)
    use simmodule, only: store_error
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: imaw
    character(len=*), intent(in) :: keyword
    character(len=*), intent(in) :: msg
    ! -- local
    ! -- formats
    !
    if (len(msg) == 0) then
      write (errmsg, '(a,1x,a,1x,i0,1x,a)') &
        keyword, ' for MAW well', imaw, 'has already been set.'
    else
      write (errmsg, '(a,1x,a,1x,i0,1x,a)') &
        keyword, ' for MAW well', imaw, msg
    end if
    call store_error(errmsg)
  end subroutine maw_set_attribute_error
 
  !> @brief Apply a PERIOD setting value to well imaw
  !!
  !! Shared by maw_rp (once per period) and maw_ad (per-timestep resync of
  !! TS-linked settings); validation stays with the caller.
  !<
  subroutine maw_set_period_value(this, n, imaw, setting)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: n !< row index into input context arrays
    integer(I4B), intent(in) :: imaw !< well number
    character(len=*), intent(in) :: setting !< uppercase tagname dispatch key
    !
    select case (trim(setting))
    case ('WELL_HEAD')
      this%well_head_set(imaw) = .true.
      this%xnewpak(imaw) = this%well_head(imaw)
    end select
  end subroutine maw_set_period_value
 
  !> @brief Issue parameter errors for mawwells(imaw)
  !<
  subroutine maw_check_attributes(this)
    use simmodule, only: store_error
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    character(len=LINELENGTH) :: cgwfnode
    integer(I4B) :: idx
    integer(I4B) :: n
    integer(I4B) :: j
    integer(I4B) :: jpos
    ! -- formats
    !
    idx = 1
    do n = 1, this%nmawwells
      if (this%ngwfnodes(n) < 1) then
        call this%maw_set_attribute_error(n, 'NGWFNODES', 'must be greater '// &
                                          'than 0.')
      end if
      if (this%radius(n) == dep20) then
        call this%maw_set_attribute_error(n, 'RADIUS', 'has not been specified.')
      end if
      if (this%shutoffmin(n) > dzero) then
        if (this%shutoffmin(n) >= this%shutoffmax(n)) then
          call this%maw_set_attribute_error(n, 'SHUT_OFF', 'shutoffmax must '// &
                                            'be greater than shutoffmin.')
        end if
      end if
      do j = 1, this%ngwfnodes(n)
        !
        ! -- calculate jpos
        jpos = this%get_jpos(n, j)
        !
        ! -- write gwfnode number
        write (cgwfnode, '(a,i0,a)') 'gwfnode(', j, ')'
        !
        ! -- connection screen data
        if (this%botscrn(jpos) >= this%topscrn(jpos)) then
          call this%maw_set_attribute_error(n, 'SCREEN_TOP', 'screen bottom '// &
                                            'must be less than screen top. '// &
                                            trim(cgwfnode))
        end if
        !
        ! -- connection skin hydraulic conductivity
        if (this%ieqn(n) == 2 .OR. this%ieqn(n) == 3 .OR. &
            this%ieqn(n) == 4) then
          if (this%hk(jpos) <= dzero) then
            call this%maw_set_attribute_error(n, 'HK_SKIN', 'skin hyraulic '// &
                                              'conductivity must be greater '// &
                                              'than zero. '//trim(cgwfnode))
          end if
        else if (this%ieqn(n) == 0) then
          !
          ! -- saturated conductance
          if (this%satcond(jpos) < dzero) then
            call this%maw_set_attribute_error(n, 'HK_SKIN', &
                                              'skin hyraulic conductivity '// &
                                              'must be greater than or '// &
                                              'equal to zero when using '// &
                                              'SPECIFIED condeqn. '// &
                                              trim(cgwfnode))
          end if
        end if
        idx = idx + 1
      end do
    end do
    ! -- reset check_attr
    this%check_attr = 0
  end subroutine maw_check_attributes
 
  !> @brief Add package connection to matrix
  !<
  subroutine maw_ac(this, moffset, sparse)
    use sparsemodule, only: sparsematrix
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: moffset
    type(sparsematrix), intent(inout) :: sparse
    ! -- local
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: jj
    integer(I4B) :: jglo
    integer(I4B) :: nglo
    ! -- format
    !
    ! -- Add package rows to sparse
    do n = 1, this%nmawwells
      nglo = moffset + this%dis%nodes + this%ioffset + n
      call sparse%addconnection(nglo, nglo, 1)
      do j = 1, this%ngwfnodes(n)
        jj = this%get_gwfnode(n, j)
        jglo = jj + moffset
        call sparse%addconnection(nglo, jglo, 1)
        call sparse%addconnection(jglo, nglo, 1)
      end do
 
    end do
  end subroutine maw_ac
 
  !> @brief Map package connection to matrix
  !<
  subroutine maw_mc(this, moffset, matrix_sln)
    use sparsemodule, only: sparsematrix
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: moffset
    class(matrixbasetype), pointer :: matrix_sln
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: j
    integer(I4B) :: ii
    integer(I4B) :: iglo
    integer(I4B) :: jglo
    integer(I4B) :: ipos
    ! -- format
    !
    ! -- allocate connection mapping vectors
    call mem_allocate(this%idxlocnode, this%nmawwells, 'IDXLOCNODE', &
                      this%memoryPath)
    call mem_allocate(this%idxdglo, this%maxbound, 'IDXDGLO', this%memoryPath)
    call mem_allocate(this%idxoffdglo, this%maxbound, 'IDXOFFDGLO', &
                      this%memoryPath)
    call mem_allocate(this%idxsymdglo, this%maxbound, 'IDXSYMDGLO', &
                      this%memoryPath)
    call mem_allocate(this%idxsymoffdglo, this%maxbound, 'IDXSYMOFFDGLO', &
                      this%memoryPath)
    !
    ! -- Find the position of each connection in the global ia, ja structure
    !    and store them in idxglo.  idxglo allows this model to insert or
    !    retrieve values into or from the global A matrix
    ! -- maw rows
    ipos = 1
    do n = 1, this%nmawwells
      iglo = moffset + this%dis%nodes + this%ioffset + n
      this%idxlocnode(n) = this%dis%nodes + this%ioffset + n
      do ii = 1, this%ngwfnodes(n)
        j = this%get_gwfnode(n, ii)
        jglo = j + moffset
        this%idxdglo(ipos) = matrix_sln%get_position_diag(iglo)
        this%idxoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
        ipos = ipos + 1
      end do
    end do
    ! -- maw contributions gwf portion of global matrix
    ipos = 1
    do n = 1, this%nmawwells
      do ii = 1, this%ngwfnodes(n)
        iglo = this%get_gwfnode(n, ii) + moffset
        jglo = moffset + this%dis%nodes + this%ioffset + n
        this%idxsymdglo(ipos) = matrix_sln%get_position_diag(iglo)
        this%idxsymoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
        ipos = ipos + 1
      end do
    end do
  end subroutine maw_mc
 
  !> @brief Allocate and Read
  !!
  !! Create new MAW package and point bndobj to the new package
  !<
  subroutine maw_ar(this)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    ! -- format
    !
    call this%obs%obs_ar()
    !
    ! -- set omega value used for saturation calculations
    if (this%inewton > 0) then
      this%satomega = dem6
    end if
    !
    ! -- Allocate connection arrays in MAW and in package superclass
    call this%maw_allocate_arrays()
    !
    ! -- read optional initial package parameters
    call this%read_initial_attr()
    !
    ! -- bind persistent pointers to input context data
    call this%input%init(this%input_mempath, this%naux)
    !
    ! -- setup pakmvrobj
    if (this%imover /= 0) then
      allocate (this%pakmvrobj)
      call this%pakmvrobj%ar(this%nmawwells, this%nmawwells, this%memoryPath)
    end if
  end subroutine maw_ar
 
  !> @brief Read and Prepare
  !!
  !! Read itmp and new boundaries if itmp > 0
  !<
  subroutine maw_rp(this)
    use constantsmodule, only: lenvarname
    use tdismodule, only: kper
    use memorymanagermodule, only: mem_setptr
    use characterstringmodule, only: characterstringtype
    use simmodule, only: store_error_filename
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    character(len=LINELENGTH) :: csteady
    character(len=LINELENGTH) :: str
    real(DP) :: auxvaldbl
    character(len=LENVARNAME) :: setting
    character(len=LINELENGTH) :: cstr
    character(len=LINELENGTH) :: errmsgr
    integer(I4B) :: node
    integer(I4B) :: n
    integer(I4B) :: ntabcols
    integer(I4B) :: ntabrows
    integer(I4B) :: imaw
    integer(I4B) :: ibnd
    integer(I4B) :: j
    integer(I4B) :: jpos
    integer(I4B) :: jj
    integer(I4B) :: nactive
    integer(I4B) :: icon
    integer(I4B), dimension(:), allocatable :: ibotset
    integer(I4B), dimension(:), allocatable :: iheadset
    integer(I4B), dimension(:), allocatable :: ilimset
    integer(I4B) :: istat
    integer(I4B) :: iheadlimit_warning
    ! -- formats
    character(len=*), parameter :: fmtlsp = &
      &"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
    character(len=*), parameter :: fmthdbot = &
      &"('well head (',G0,') must be >= BOTTOM_ELEVATION (',G0, ').')"
    !
    ! -- initialize counters
    iheadlimit_warning = 0
    !
    ! -- set steady-state flag based on gwfiss
    this%imawiss = this%gwfiss
    !
    ! -- reset maw steady flag if 'STEADY-STATE' specified in the OPTIONS block
    if (this%imawissopt == 1) then
      this%imawiss = 1
    end if
    !
    ! -- set nbound to maxbound
    this%nbound = this%maxbound
    !
    ! -- check if input context has been updated for this stress period
    if (this%iper == kper) then
      !
      ! -- setup table to echo period data
      if (this%iprpak /= 0) then
        title = trim(adjustl(this%text))//' PACKAGE ('// &
                trim(adjustl(this%packName))//') DATA FOR PERIOD'
        write (title, '(a,1x,i6)') trim(adjustl(title)), kper
        call table_cr(this%inputtab, this%packName, title)
        call this%inputtab%table_df(1, 5, this%iout, finalize=.false.)
        text = 'NUMBER'
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
        text = 'KEYWORD'
        call this%inputtab%initialize_column(text, 20, alignment=tableft)
        do jj = 1, 3
          write (text, '(a,1x,i6)') 'VALUE', jj
          call this%inputtab%initialize_column(text, 15, alignment=tabcenter)
        end do
      end if
      !
      ! -- set flag to check attributes
      this%check_attr = 1
      ! -- wells with BOTTOM, WELL_HEAD, or HEAD_LIMIT set this period, which
      !    are checked after all of the period's settings are applied
      allocate (ibotset(this%nmawwells))
      allocate (iheadset(this%nmawwells))
      allocate (ilimset(this%nmawwells))
      do imaw = 1, this%nmawwells
        ibotset(imaw) = 0
        iheadset(imaw) = 0
        ilimset(imaw) = 0
      end do
      !
      do n = 1, this%input%nbound
        imaw = this%input%ifno(n)
        if (imaw < 1 .or. imaw > this%nmawwells) then
          write (errmsg, '(2(a,1x),i0,a)') &
            'IMAW must be greater than 0 and', &
            'less than or equal to ', this%nmawwells, '.'
          call store_error(errmsg)
          cycle
        end if
        !
        ! -- dispatch key for this row
        setting = this%input%setting(n)
        !
        ! -- STATUS
        if (trim(setting) == 'STATUS') then
          str = this%input%status(n)
          this%status(imaw) = str(1:8)
          select case (trim(str))
          case ('CONSTANT')
            this%iboundpak(imaw) = -1
          case ('INACTIVE')
            this%iboundpak(imaw) = 0
          case ('ACTIVE')
            this%iboundpak(imaw) = 1
          case default
            write (errmsg, '(2a)') &
              'Unknown '//trim(this%text)//" maw status keyword: '", &
              trim(str)//"'."
            call store_error(errmsg)
          end select
        end if
        !
        ! -- CONNECTION_STATUS (compound group)
        if (trim(setting) == 'CONN_STATUS') then
          icon = this%input%period_icon(n)
          str = this%input%connstatus(n)
          if (icon < 1 .or. icon > this%ngwfnodes(imaw)) then
            write (errmsg, '(2(a,1x),i0,1x,a,1x,i0,a)') &
              'ICON must be greater than 0 and', &
              'less than or equal to ', this%ngwfnodes(imaw), &
              'for maw well', imaw, '.'
            call store_error(errmsg)
          else
            jpos = this%get_jpos(imaw, icon)
            select case (trim(str))
            case ('INACTIVE')
              this%iboundconn(jpos) = 0
            case ('ACTIVE')
              this%iboundconn(jpos) = 1
            case default
              write (errmsg, '(2a)') &
                'Unknown '//trim(this%text)//" maw connection status "// &
                "keyword: '", trim(str)//"'."
              call store_error(errmsg)
            end select
          end if
        end if
        !
        ! -- BOTTOM
        if (trim(setting) == 'PERIOD_BOTTOM') then
          this%bot(imaw) = this%input%period_bottom(n)
          ibotset(imaw) = 1
        end if
        !
        ! -- RATE / WELL_HEAD
        call this%maw_set_period_value(n, imaw, setting)
        !
        ! -- WELL_HEAD
        if (trim(setting) == 'WELL_HEAD') then
          iheadset(imaw) = 1
        end if
        !
        ! -- HEAD_LIMIT
        if (trim(setting) == 'HEAD_LIMIT') then
          str = this%input%head_limit(n)
          if (trim(str) == 'OFF') then
            this%shutofflevel(imaw) = dep20
          else
            read (str, *, iostat=istat, iomsg=errmsgr) &
              this%shutofflevel(imaw)
            if (istat /= 0) then
              errmsg = 'Could not read HEAD_LIMIT value. '//trim(errmsgr)
              call store_error(errmsg)
            end if
          end if
          ilimset(imaw) = 1
        end if
        !
        ! -- FLOWING_WELL (compound group)
        if (trim(setting) == 'FLOWING_WELL') then
          this%fwelev(imaw) = this%input%fwelev(n)
          this%fwcond(imaw) = this%input%fwcond(n)
          this%fwrlen(imaw) = this%input%fwrlen(n)
          if (this%iflowingwells == 0) then
            this%iflowingwells = -1
            text = 'Flowing well data is specified in the '// &
                   trim(this%packName)//' package but FLOWING_WELL was'// &
                   ' not specified in the OPTIONS block.'
            call store_warning(text)
          end if
        end if
        !
        ! -- SHUT_OFF (compound group)
        if (trim(setting) == 'SHUT_OFF') then
          this%shutoffmin(imaw) = this%input%minrate(n)
          this%shutoffmax(imaw) = this%input%maxrate(n)
        end if
        !
        ! -- RATE_SCALING (compound group)
        if (trim(setting) == 'RATE_SCALING') then
          this%pumpelev(imaw) = this%input%pump_elevation(n)
          this%reduction_length(imaw) = this%input%scaling_length(n)
          if (this%reduction_length(imaw) < dzero) then
            call this%maw_set_attribute_error(imaw, 'RATE_SCALING', &
              &'must be greater than or equal to 0.')
          end if
        end if
        !
        ! -- AUXILIARY (compound group); the current value is already set --
        ! -- only extract auxname here for echoing this period's row below
        if (this%naux > 0) then
          if (trim(setting) == 'PERIOD_AUXILIARY') then
            str = this%input%auxname(n)
            auxvaldbl = this%input%auxval(n)
          end if
        end if
        !
        ! -- echo row to period data table
        if (this%iprpak /= 0) then
          call this%inputtab%add_term(imaw)
          if (trim(setting) == 'PERIOD_AUXILIARY') then
            call this%inputtab%add_term('AUXILIARY')
          else
            call this%inputtab%add_term(trim(setting))
          end if
          select case (trim(setting))
          case ('STATUS')
            call this%inputtab%add_term(trim(str))
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
          case ('CONN_STATUS')
            call this%inputtab%add_term(icon)
            call this%inputtab%add_term(trim(str))
            call this%inputtab%add_term(' ')
          case ('PERIOD_BOTTOM')
            call this%inputtab%add_term(this%bot(imaw))
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
          case ('RATE')
            call this%inputtab%add_term(this%rate(imaw))
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
          case ('WELL_HEAD')
            call this%inputtab%add_term(this%well_head(imaw))
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
          case ('HEAD_LIMIT')
            call this%inputtab%add_term(trim(str))
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
          case ('FLOWING_WELL')
            call this%inputtab%add_term(this%fwelev(imaw))
            call this%inputtab%add_term(this%fwcond(imaw))
            call this%inputtab%add_term(this%fwrlen(imaw))
          case ('SHUT_OFF')
            call this%inputtab%add_term(this%shutoffmin(imaw))
            call this%inputtab%add_term(this%shutoffmax(imaw))
            call this%inputtab%add_term(' ')
          case ('RATE_SCALING')
            call this%inputtab%add_term(this%pumpelev(imaw))
            call this%inputtab%add_term(this%reduction_length(imaw))
            call this%inputtab%add_term(' ')
          case ('PERIOD_AUXILIARY')
            call this%inputtab%add_term(trim(str))
            call this%inputtab%add_term(auxvaldbl)
            call this%inputtab%add_term(' ')
          case default
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
            call this%inputtab%add_term(' ')
          end select
        end if
      end do
      if (this%iprpak /= 0) then
        call this%inputtab%finalize_table()
      end if
      !
      ! -- the well bottom, the datum for well storage, may not be above the
      !    head or the screen bottom of an active connection; a CONSTANT well
      !    is checked against its well head
      do imaw = 1, this%nmawwells
        if (ibotset(imaw) /= 0 .and. this%iboundpak(imaw) >= 0) then
          if (this%bot(imaw) > this%xnewpak(imaw)) then
            write (cstr, fmthdbot) this%xnewpak(imaw), this%bot(imaw)
            call this%maw_set_attribute_error(imaw, 'BOTTOM', trim(cstr))
          end if
        end if
        do jj = 1, this%ngwfnodes(imaw)
          jpos = this%get_jpos(imaw, jj)
          if (this%iboundconn(jpos) == 0) cycle
          if (this%bot(imaw) > this%botscrn(jpos)) then
            write (errmsg, '(a,g0,a,1x,i0,1x,a,1x,i0,1x,a,g0,a)') &
              'BOTTOM (', this%bot(imaw), ') for maw well', imaw, &
              'is above the screen bottom of active connection', jj, '(', &
              this%botscrn(jpos), ').'
            call store_error(errmsg)
            exit
          end if
        end do
        if (iheadset(imaw) /= 0 .or. this%iboundpak(imaw) < 0) then
          if (this%well_head(imaw) < this%bot(imaw)) then
            write (cstr, fmthdbot) this%well_head(imaw), this%bot(imaw)
            call this%maw_set_attribute_error(imaw, 'WELL HEAD', trim(cstr))
          end if
        end if
        if (ibotset(imaw) /= 0 .or. ilimset(imaw) /= 0) then
          if (this%shutofflevel(imaw) <= this%bot(imaw)) then
            iheadlimit_warning = iheadlimit_warning + 1
          end if
        end if
      end do
      deallocate (ibotset)
      deallocate (iheadset)
      deallocate (ilimset)
      !
      ! -- a well must keep at least one active connection; STATUS INACTIVE,
      !    not CONNECTION_STATUS, is used to deactivate every connection
      do imaw = 1, this%nmawwells
        nactive = 0
        do jj = 1, this%ngwfnodes(imaw)
          jpos = this%get_jpos(imaw, jj)
          if (this%iboundconn(jpos) /= 0) nactive = nactive + 1
        end do
        if (nactive == 0) then
          write (errmsg, '(a,1x,i0,1x,a)') &
            'Every connection of maw well', imaw, 'is inactive. Use '// &
            'STATUS INACTIVE to deactivate all connections of a well.'
          call store_error(errmsg)
        end if
      end do
      !
      ! -- using data from the last stress period
    else
      write (this%iout, fmtlsp) trim(this%filtyp)
    end if
    !
    ! -- issue warning messages
    if (iheadlimit_warning > 0) then
      write (warnmsg, '(a,a,a,1x,a,1x,a)') &
        "HEAD_LIMIT in '", trim(this%packName), "' was below the well bottom", &
        "for one or more multi-aquifer well(s). This may result in", &
        "convergence failures for some models."
      call store_warning(warnmsg, substring=warnmsg(:50))
    end if
    !
    ! -- write summary of maw well stress period error messages
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- qa data if necessary
    if (this%check_attr /= 0) then
      call this%maw_check_attributes()
 
      ! -- write summary of stress period data for MAW
      if (this%iprpak == 1) then
        if (this%imawiss /= 0) then
          csteady = 'STEADY-STATE    '
        else
          csteady = 'TRANSIENT       '
        end if
        !
        ! -- reset the input table object for rate data
        title = trim(adjustl(this%text))//' PACKAGE ('// &
                trim(adjustl(this%packName))//') '//trim(adjustl(csteady))// &
                ' RATE DATA FOR PERIOD'
        write (title, '(a,1x,i6)') trim(adjustl(title)), kper
        ntabcols = 6
        call table_cr(this%inputtab, this%packName, title)
        call this%inputtab%table_df(this%nmawwells, ntabcols, this%iout)
        text = 'NUMBER'
        call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
        text = 'STATUS'
        call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
        text = 'RATE'
        call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
        text = 'SPECIFIED HEAD'
        call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
        text = 'PUMP ELEVATION'
        call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
        text = 'REDUCTION LENGTH'
        call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
        do n = 1, this%nmawwells
          call this%inputtab%add_term(n)
          call this%inputtab%add_term(this%status(n))
          call this%inputtab%add_term(this%rate(n))
          if (this%iboundpak(n) < 0) then
            call this%inputtab%add_term(this%well_head(n))
          else
            call this%inputtab%add_term(' ')
          end if
          call this%inputtab%add_term(this%pumpelev(n))
          if (this%reduction_length(n) /= dep20) then
            call this%inputtab%add_term(this%reduction_length(n))
          else
            call this%inputtab%add_term(' ')
          end if
        end do
        !
        ! -- flowing wells
        if (this%iflowingwells > 0) then
          !
          ! -- reset the input table object for flowing well data
          title = trim(adjustl(this%text))//' PACKAGE ('// &
                  trim(adjustl(this%packName))//') '//trim(adjustl(csteady))// &
                  ' FLOWING WELL DATA FOR PERIOD'
          write (title, '(a,1x,i6)') trim(adjustl(title)), kper
          ntabcols = 4
          ntabrows = 0
          do n = 1, this%nmawwells
            if (this%fwcond(n) > dzero) then
              ntabrows = ntabrows + 1
            end if
          end do
          if (ntabrows > 0) then
            call table_cr(this%inputtab, this%packName, title)
            call this%inputtab%table_df(ntabrows, ntabcols, this%iout)
            text = 'NUMBER'
            call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
            text = 'ELEVATION'
            call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
            text = 'CONDUCT.'
            call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
            text = 'REDUCTION LENGTH'
            call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
            do n = 1, this%nmawwells
              if (this%fwcond(n) > dzero) then
                call this%inputtab%add_term(n)
                call this%inputtab%add_term(this%fwelev(n))
                call this%inputtab%add_term(this%fwcond(n))
                call this%inputtab%add_term(this%fwrlen(n))
              end if
            end do
          end if
        end if
        !
        ! -- reset the input table object for shutoff data
        title = trim(adjustl(this%text))//' PACKAGE ('// &
                trim(adjustl(this%packName))//') '//trim(adjustl(csteady))// &
                ' WELL SHUTOFF DATA FOR PERIOD'
        write (title, '(a,1x,i6)') trim(adjustl(title)), kper
        ntabcols = 4
        ntabrows = 0
        do n = 1, this%nmawwells
          if (this%shutofflevel(n) /= dep20) then
            ntabrows = ntabrows + 1
          end if
        end do
        if (ntabrows > 0) then
          call table_cr(this%inputtab, this%packName, title)
          call this%inputtab%table_df(ntabrows, ntabcols, this%iout)
          text = 'NUMBER'
          call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
          text = 'ELEVATION'
          call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
          text = 'MINIMUM. Q'
          call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
          text = 'MAXIMUM Q'
          call this%inputtab%initialize_column(text, 12, alignment=tabcenter)
          do n = 1, this%nmawwells
            if (this%shutofflevel(n) /= dep20) then
              call this%inputtab%add_term(n)
              call this%inputtab%add_term(this%shutofflevel(n))
              call this%inputtab%add_term(this%shutoffmin(n))
              call this%inputtab%add_term(this%shutoffmax(n))
            end if
          end do
        end if
      end if
    end if
    !
    ! -- fill arrays
    ibnd = 1
    do n = 1, this%nmawwells
      do j = 1, this%ngwfnodes(n)
        jpos = this%get_jpos(n, j)
        node = this%get_gwfnode(n, j)
        this%nodelist(ibnd) = node
        this%bound(1, ibnd) = this%xnewpak(n)
        this%bound(2, ibnd) = this%satcond(jpos)
        this%bound(3, ibnd) = this%botscrn(jpos)
        if (this%iboundpak(n) > 0) then
          this%bound(4, ibnd) = this%rate(n)
        else
          this%bound(4, ibnd) = dzero
        end if
        ibnd = ibnd + 1
      end do
    end do
  end subroutine maw_rp
 
  !> @brief Add package connection to matrix
  !<
  subroutine maw_ad(this)
    use constantsmodule, only: lenvarname
    use tdismodule, only: kper, kstp
    ! -- dummy
    class(mawtype) :: this
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: n
    integer(I4B) :: j
    integer(I4B) :: jj
    integer(I4B) :: ibnd
    integer(I4B) :: imaw
    character(len=LENVARNAME) :: setting
    !
    ! -- sync PACKAGEDATA AUX TS and advance observations;
    ! -- step auxvar period block overrides
    call this%BndExtType%bnd_ad()
    !
    ! -- sync PERIOD TS from input context
    if (this%ts_active .and. this%input%nbound > 0) then
      do n = 1, this%input%nbound
        imaw = this%input%ifno(n)
        if (imaw < 1 .or. imaw > this%nmawwells) cycle
        setting = this%input%setting(n)
        call this%maw_set_period_value(n, imaw, setting)
      end do
    end if
    !
    ! -- update auxiliary variables
    if (this%naux > 0) then
      ibnd = 1
      do n = 1, this%nmawwells
        do j = 1, this%ngwfnodes(n)
          do jj = 1, this%naux
            if (this%noupdateauxvar(jj) /= 0) cycle
            this%auxvar(jj, ibnd) = this%featureauxvar(jj, n)
          end do
          ibnd = ibnd + 1
        end do
      end do
    end if
    !
    ! -- copy xnew into xold; CONSTANT well heads track the PACKAGEDATA STRT
    !    baseline until a PERIOD WELL_HEAD entry overrides it
    do i = 1, this%nmawwells
      n = this%pkg_ifno(i)
      this%xoldpak(n) = this%xnewpak(n)
      this%xoldsto(n) = this%xsto(n)
      if (this%iboundpak(n) < 0) then
        if (.not. this%well_head_set(n)) then
          this%well_head(n) = this%input%strt(i)
        end if
        this%xnewpak(n) = this%well_head(n)
      end if
      !
      ! -- start each time step with damping turned off. Clearing the weight
      !    and the previous head change stops a leftover value from the last
      !    time step from triggering damping on the first iteration.
      this%nurdxold(n) = dzero
      this%nurweight(n) = done
    end do
    !
    !--use the appropriate xoldsto if initial heads are above the
    !  specified flowing well discharge elevation
    if (kper == 1 .and. kstp == 1) then
      do n = 1, this%nmawwells
        if (this%fwcond(n) > dzero) then
          if (this%xoldsto(n) > this%fwelev(n)) then
            this%xoldsto(n) = this%fwelev(n)
          end if
        end if
      end do
    end if
    !
    ! -- reset ishutoffcnt (equivalent to kiter) to zero
    this%ishutoffcnt = 0
    !
    ! -- pakmvrobj ad
    if (this%imover == 1) then
      call this%pakmvrobj%ad()
    end if
  end subroutine maw_ad
 
  !> @brief Formulate the HCOF and RHS terms
  !!
  !! Skip if no multi-aquifer wells, otherwise, calculate hcof and rhs
  !<
  subroutine maw_cf(this)
    ! -- dummy
    class(mawtype) :: this
    ! -- local
    !
    ! -- Calculate maw conductance and update package RHS and HCOF
    call this%maw_cfupdate()
  end subroutine maw_cf
 
  !> @brief Copy rhs and hcof into solution rhs and amat
  !<
  subroutine maw_fc(this, rhs, ia, idxglo, matrix_sln)
    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(mawtype) :: this
    real(DP), dimension(:), intent(inout) :: rhs
    integer(I4B), dimension(:), intent(in) :: ia
    integer(I4B), dimension(:), intent(in) :: idxglo
    class(matrixbasetype), pointer :: matrix_sln
    ! -- local
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: idx
    integer(I4B) :: iloc
    integer(I4B) :: isymloc
    integer(I4B) :: igwfnode
    integer(I4B) :: iposd
    integer(I4B) :: iposoffd
    integer(I4B) :: isymnode
    integer(I4B) :: ipossymd
    integer(I4B) :: ipossymoffd
    integer(I4B) :: jpos
    integer(I4B) :: icflow
    real(DP) :: hmaw
    real(DP) :: hgwf
    real(DP) :: cfw
    real(DP) :: cmaw
    real(DP) :: cterm
    real(DP) :: term
    real(DP) :: scale
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: rate
    real(DP) :: ratefw
    real(DP) :: flow
    real(DP) :: tled
    real(DP) :: sd
    real(DP) :: ss
    real(DP) :: ssold
    !
    ! -- pakmvrobj fc
    if (this%imover == 1) then
      call this%pakmvrobj%fc()
    end if
    !
    ! -- Copy package rhs and hcof into solution rhs and amat
    idx = 1
    do n = 1, this%nmawwells
      iloc = this%idxlocnode(n)
      !
      ! -- save the well rate from the previous outer iteration before it is
      !    recalculated below, for the rate convergence check in maw_cc
      this%qsim0(n) = this%ratesim(n)
      !
      ! -- update head value for constant head maw wells
      if (this%iboundpak(n) < 0) then
        this%xnewpak(n) = this%well_head(n)
      end if
      hmaw = this%xnewpak(n)
      !
      ! -- add pumping rate to active or constant maw well
      if (this%iboundpak(n) == 0) then
        this%ratesim(n) = dzero
      else
        call this%maw_calculate_wellq(n, hmaw, rate)
        this%ratesim(n) = rate
        rhs(iloc) = rhs(iloc) - rate
        !
        ! -- location of diagonal for maw row
        iposd = this%idxdglo(idx)
        !
        ! -- add flowing well
        this%xsto(n) = hmaw
        ratefw = dzero
        if (this%iflowingwells > 0) then
          if (this%fwcond(n) > dzero) then
            bt = this%fwelev(n)
            tp = bt + this%fwrlen(n)
            scale = sqsaturation(tp, bt, hmaw)
            cfw = scale * this%fwcond(n)
            this%ifwdischarge(n) = 0
            if (cfw > dzero) then
              this%ifwdischarge(n) = 1
              this%xsto(n) = bt
            end if
            this%fwcondsim(n) = cfw
            call matrix_sln%add_value_pos(iposd, -cfw)
            rhs(iloc) = rhs(iloc) - cfw * bt
            ratefw = cfw * (bt - hmaw)
          end if
        end if
        !
        ! -- add maw storage changes
        if (this%imawiss /= 1) then
          if (this%ifwdischarge(n) /= 1) then
            ! -- well storage. The storage water level is not allowed to drop
            !    below the bottom of the well: ss is a smooth version of
            !    max(hmaw, well bottom). The storage release then fades out as
            !    the well empties instead of drawing water from below the well
            !    bottom. While the head is above the bottom this is the original
            !    term (sd = 1, ss = hmaw, giving -area/delt).
            tled = this%area(n) / delt
            sd = squadratic0spderivative(hmaw, this%bot(n), this%satomega)
            ss = squadratic0sp(hmaw, this%bot(n), this%satomega)
            ssold = squadratic0sp(this%xoldsto(n), this%bot(n), this%satomega)
            call matrix_sln%add_value_pos(iposd, -tled * sd)
            rhs(iloc) = rhs(iloc) - tled * (sd * hmaw - ss + ssold)
          else
            cterm = this%xoldsto(n) - this%fwelev(n)
            rhs(iloc) = rhs(iloc) - (this%area(n) * cterm / delt)
          end if
        end if
        !
        ! -- If mover is active, add receiver water to rhs and
        !    store available water (as positive value)
        if (this%imover == 1) then
          rhs(iloc) = rhs(iloc) - this%pakmvrobj%get_qfrommvr(n)
          !
          ! -- add pumping rate to mover if not injection
          if (rate < 0) then
            call this%pakmvrobj%accumulate_qformvr(n, -rate) !pumped water
          end if
          !
          ! -- add flowing well flow to mover
          call this%pakmvrobj%accumulate_qformvr(n, -ratefw) !flowing water
        end if
        !
      end if
      !
      ! -- process each maw/gwf connection
      do j = 1, this%ngwfnodes(n)
        jpos = this%get_jpos(n, j)
        if (this%iboundpak(n) /= 0 .and. this%iboundconn(jpos) /= 0) then
          igwfnode = this%get_gwfnode(n, j)
          hgwf = this%xnew(igwfnode)
          !
          ! -- calculate connection terms
          call this%maw_calculate_conn_terms(n, j, icflow, cmaw, cterm, term, &
                                             flow)
          this%simcond(jpos) = cmaw
          !
          ! -- add to maw row
          iposd = this%idxdglo(idx)
          iposoffd = this%idxoffdglo(idx)
          call matrix_sln%add_value_pos(iposd, -term)
          call matrix_sln%set_value_pos(iposoffd, term)
          !
          ! -- add correction term
          rhs(iloc) = rhs(iloc) - cterm
          !
          ! -- add to gwf row for maw connection
          isymnode = this%get_gwfnode(n, j)
          isymloc = ia(isymnode)
          ipossymd = this%idxsymdglo(idx)
          ipossymoffd = this%idxsymoffdglo(idx)
          call matrix_sln%add_value_pos(ipossymd, -term)
          call matrix_sln%set_value_pos(ipossymoffd, term)
          !
          ! -- add correction term to gwf row
          rhs(isymnode) = rhs(isymnode) + cterm
        end if
        !
        ! -- increment maw connection counter
        idx = idx + 1
      end do
    end do
  end subroutine maw_fc
 
  !> @brief Fill newton terms
  !<
  subroutine maw_fn(this, rhs, ia, idxglo, matrix_sln)
    ! -- dummy
    class(mawtype) :: this
    real(DP), dimension(:), intent(inout) :: rhs
    integer(I4B), dimension(:), intent(in) :: ia
    integer(I4B), dimension(:), intent(in) :: idxglo
    class(matrixbasetype), pointer :: matrix_sln
    ! -- local
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: idx
    integer(I4B) :: iloc
    integer(I4B) :: isymloc
    integer(I4B) :: igwfnode
    integer(I4B) :: iposd
    integer(I4B) :: iposoffd
    integer(I4B) :: isymnode
    integer(I4B) :: ipossymd
    integer(I4B) :: ipossymoffd
    integer(I4B) :: jpos
    integer(I4B) :: icflow
    real(DP) :: hmaw
    real(DP) :: hgwf
    real(DP) :: scale
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: cfw
    real(DP) :: rate
    real(DP) :: rate2
    real(DP) :: rterm
    real(DP) :: derv
    real(DP) :: drterm
    real(DP) :: cmaw
    real(DP) :: cterm
    real(DP) :: term
    real(DP) :: flow
    real(DP) :: term2
    real(DP) :: rhsterm
    !
    ! -- Calculate Newton-Raphson corrections
    idx = 1
    do n = 1, this%nmawwells
      iloc = this%idxlocnode(n)
      hmaw = this%xnewpak(n)
      !
      ! -- add pumping rate to active or constant maw well
      if (this%iboundpak(n) /= 0) then
        iposd = this%idxdglo(idx)
        scale = done
        drterm = dzero
        rate = this%ratesim(n)
        !
        !-- calculate final derivative for pumping rate
        call this%maw_calculate_wellq(n, hmaw + dem4, rate2)
        drterm = (rate2 - rate) / dem4
        !
        !-- fill amat and rhs with newton-raphson terms
        call matrix_sln%add_value_pos(iposd, drterm)
        rhs(iloc) = rhs(iloc) + drterm * hmaw
        !
        ! -- add flowing well
        if (this%iflowingwells > 0) then
          if (this%fwcond(n) > dzero) then
            bt = this%fwelev(n)
            tp = bt + this%fwrlen(n)
            scale = sqsaturation(tp, bt, hmaw)
            cfw = scale * this%fwcond(n)
            this%ifwdischarge(n) = 0
            if (cfw > dzero) then
              this%ifwdischarge(n) = 1
            end if
            this%fwcondsim(n) = cfw
            rate = cfw * (bt - hmaw)
            rterm = -cfw * hmaw
            !
            ! --calculate derivative for flowing well
            if (hmaw < tp) then
              derv = sqsaturationderivative(tp, bt, hmaw)
              drterm = -(cfw + this%fwcond(n) * derv * (hmaw - bt))
              !
              ! -- fill amat and rhs with newton-raphson terms
              call matrix_sln%add_value_pos(iposd, &
                                            -this%fwcond(n) * derv * (hmaw - bt))
              rhs(iloc) = rhs(iloc) - rterm + drterm * hmaw
            end if
          end if
        end if
      end if
      !
      ! -- process each maw/gwf connection
      do j = 1, this%ngwfnodes(n)
        jpos = this%get_jpos(n, j)
        if (this%iboundpak(n) /= 0 .and. this%iboundconn(jpos) /= 0) then
          igwfnode = this%get_gwfnode(n, j)
          hgwf = this%xnew(igwfnode)
          !
          ! -- add to maw row
          iposd = this%idxdglo(idx)
          iposoffd = this%idxoffdglo(idx)
          !
          ! -- add to gwf row for maw connection
          isymnode = this%get_gwfnode(n, j)
          isymloc = ia(isymnode)
          ipossymd = this%idxsymdglo(idx)
          ipossymoffd = this%idxsymoffdglo(idx)
          !
          ! -- calculate newton terms
          call this%maw_calculate_conn_terms(n, j, icflow, cmaw, cterm, term, &
                                             flow, term2)
          !
          ! -- maw is upstream
          if (hmaw > hgwf) then
            if (icflow /= 0) then
              rhsterm = term2 * hgwf + term * hmaw
              rhs(iloc) = rhs(iloc) + rhsterm
              rhs(isymnode) = rhs(isymnode) - rhsterm
              if (this%iboundpak(n) > 0) then
                call matrix_sln%add_value_pos(iposd, term)
                call matrix_sln%add_value_pos(iposoffd, term2)
              end if
              call matrix_sln%add_value_pos(ipossymd, -term2)
              call matrix_sln%add_value_pos(ipossymoffd, -term)
            else
              rhs(iloc) = rhs(iloc) + term * hmaw
              rhs(isymnode) = rhs(isymnode) - term * hmaw
              call matrix_sln%add_value_pos(iposd, term)
              if (this%ibound(igwfnode) > 0) then
                call matrix_sln%add_value_pos(ipossymoffd, -term)
              end if
            end if
            !
            ! -- gwf is upstream
          else
            if (icflow /= 0) then
              rhsterm = term2 * hmaw + term * hgwf
              rhs(iloc) = rhs(iloc) + rhsterm
              rhs(isymnode) = rhs(isymnode) - rhsterm
              if (this%iboundpak(n) > 0) then
                call matrix_sln%add_value_pos(iposd, term2)
                call matrix_sln%add_value_pos(iposoffd, term)
              end if
              call matrix_sln%add_value_pos(ipossymd, -term)
              call matrix_sln%add_value_pos(ipossymoffd, -term2)
            else
              rhs(iloc) = rhs(iloc) + term * hgwf
              rhs(isymnode) = rhs(isymnode) - term * hgwf
              if (this%iboundpak(n) > 0) then
                call matrix_sln%add_value_pos(iposoffd, term)
              end if
              call matrix_sln%add_value_pos(ipossymd, -term)
            end if
          end if
        end if
        !
        ! -- increment maw connection counter
        idx = idx + 1
      end do
    end do
  end subroutine maw_fn
 
  !> @brief Apply Newton under-relaxation to the MAW Package well heads.
  !!
  !! Two corrections are applied to each well head, in order:
  !!   1. Oscillation damping. If the well head change reverses direction and is
  !!      not getting smaller, it is cut back by a per-well weight. The weight
  !!      shrinks while the head oscillates and grows back when the head moves
  !!      steadily in one direction. This settles poorly connected wells (for
  !!      example a low-K cell or a single connection) whose full Newton step is
  !!      too large.
  !!   2. Bottom limit. If the head is still below the well bottom, it is moved
  !!      back up toward the well bottom (the original behavior).
  !<
  subroutine maw_nur(this, neqpak, x, xtemp, dx, inewtonur, dxmax, locmax)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: neqpak
    real(DP), dimension(neqpak), intent(inout) :: x
    real(DP), dimension(neqpak), intent(in) :: xtemp
    real(DP), dimension(neqpak), intent(inout) :: dx
    integer(I4B), intent(inout) :: inewtonur
    real(DP), intent(inout) :: dxmax
    integer(I4B), intent(inout) :: locmax
    ! -- local
    integer(I4B) :: n
    real(DP) :: botw
    real(DP) :: xx
    real(DP) :: dxx
    real(DP) :: dxprop
    real(DP) :: weight
    ! -- parameters
    real(DP), parameter :: damptheta = dp7 !factor the weight is cut by when the head oscillates
    real(DP), parameter :: damptol = dhalf !only damp if the change is at least this fraction of the last one
    real(DP), parameter :: weightmin = dem2 !smallest allowed weight
    real(DP), parameter :: recover = 0.2_dp !amount the weight grows back each iteration
    !
    ! -- Newton-Raphson under-relaxation
    do n = 1, this%nmawwells
      if (this%iboundpak(n) < 1) cycle
      botw = this%bot(n)
      !
      ! -- full Newton head change proposed for this outer iteration
      dxprop = x(n) - xtemp(n)
      !
      ! -- oscillation damping. Cut the weight only when the head change flips
      !    direction and is not already shrinking on its own (it is still at
      !    least damptol times the previous change). Otherwise let the weight
      !    grow back toward one, so a normal, converging well is not slowed
      !    down. This catches oscillations above the well bottom, which the
      !    bottom limit below cannot.
      weight = maw_damp_weight(dxprop, this%nurdxold(n), this%nurweight(n), &
                               damptheta, damptol, weightmin, recover)
      this%nurweight(n) = weight
      !
      ! -- apply the reduced step when damping is turned on
      if (weight < done) then
        inewtonur = 1
        xx = xtemp(n) + weight * dxprop
        dxx = x(n) - xx
        if (abs(dxx) > abs(dxmax)) then
          locmax = n
          dxmax = dxx
        end if
        x(n) = xx
        dx(n) = weight * dxprop
      end if
      !
      ! -- bottom limit: if the (reduced) head is still below the bottom of the
      !    well, move it back up toward the well bottom
      if (x(n) < botw) then
        inewtonur = 1
        xx = xtemp(n) * (done - dp9) + botw * dp9
        dxx = x(n) - xx
        if (abs(dxx) > abs(dxmax)) then
          locmax = n
          dxmax = dxx
        end if
        x(n) = xx
        dx(n) = dzero
      end if
      !
      ! -- save the head change actually applied this iteration (after damping
      !    and the bottom limit). Saving it after the bottom limit lets the
      !    next iteration see the downward push from the limit, so the upward
      !    swing that follows can be damped.
      this%nurdxold(n) = x(n) - xtemp(n)
    end do
  end subroutine maw_nur
 
  !> @brief Update the oscillation-damping weight for a single MAW well head.
  !!
  !! This is the decision used by maw_nur. The weight is cut back (multiplied by
  !! damptheta, but never below weightmin) only when the proposed head change
  !! reversed direction and is not already getting smaller (its size is at least
  !! damptol times the previous change). Otherwise the weight grows back toward
  !! one by recover. It is a pure function so the logic can be unit tested.
  !<
  pure function maw_damp_weight(dxprop, dxold, weight, damptheta, damptol, &
                                weightmin, recover) result(new_weight)
    ! -- dummy
    real(dp), intent(in) :: dxprop !proposed well head change this outer iteration
    real(dp), intent(in) :: dxold !well head change applied last outer iteration
    real(dp), intent(in) :: weight !current damping weight
    real(dp), intent(in) :: damptheta !factor the weight is cut by when oscillating
    real(dp), intent(in) :: damptol !smallest |dxprop|/|dxold| ratio that is damped
    real(dp), intent(in) :: weightmin !smallest allowed weight
    real(dp), intent(in) :: recover !amount the weight grows back each iteration
    ! -- return
    real(dp) :: new_weight
    !
    if (dxprop * dxold < dzero .and. abs(dxprop) > damptol * abs(dxold)) then
      new_weight = max(damptheta * weight, weightmin)
    else
      new_weight = min(weight + recover, done)
    end if
  end function maw_damp_weight
 
  !> @brief Extra convergence check for the MAW package.
  !!
  !! This does two things:
  !!   1. While a well is being damped by Newton under-relaxation, it keeps the
  !!      model from converging until the well rate has also stopped changing
  !!      (reported through dpak). Wells that are not being damped, and models
  !!      that do not use Newton under-relaxation, are not affected.
  !!   2. When the model fails to converge, it warns about any pumping well that
  !!      is asking for more water than the aquifer can supply, a likely reason
  !!      for the failure.
  !<
  subroutine maw_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: innertot
    integer(I4B), intent(in) :: kiter
    integer(I4B), intent(in) :: iend
    integer(I4B), intent(in) :: icnvgmod
    character(len=LENPAKLOC), intent(inout) :: cpak
    integer(I4B), intent(inout) :: ipak
    real(DP), intent(inout) :: dpak
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: j
    integer(I4B) :: jpos
    integer(I4B) :: igwfnode
    integer(I4B) :: locdpak
    real(DP) :: botw
    real(DP) :: bmaw
    real(DP) :: hgwf
    real(DP) :: hv
    real(DP) :: sat
    real(DP) :: cmaw
    real(DP) :: qmax
    real(DP) :: qreq
    real(DP) :: qtolfact
    real(DP) :: dq
    real(DP) :: dpakmax
    character(len=LENPAKLOC) :: cloc
    ! -- parameters
    real(DP), parameter :: qtol = 1.001_dp !the request must exceed supply by this factor to warn
    ! -- formats
    character(len=*), parameter :: fmtwarn = &
      "('MAW well ', a, ' requests an extraction rate (', g0.5, &
      &') larger than the maximum rate (', g0.5, ') the aquifer can supply &
      &with the well head at its bottom. This may be preventing convergence. &
      &Consider reducing the requested rate, applying or widening RATE_SCALING, &
      &or reviewing the connection conductance (e.g. aquifer K).')"
    !
    ! -- rate convergence check. While a well is being damped (nurweight < 1)
    !    its head barely moves, so also require the well rate to stop changing
    !    before the model is allowed to converge. The rate change is turned
    !    into an equivalent head change (using the well area and time step) so
    !    it can be compared with the solver tolerance, and is passed back in
    !    dpak. Wells that are not being damped are left unchanged.
    dpakmax = dzero
    locdpak = 0
    do n = 1, this%nmawwells
      if (this%iboundpak(n) < 1) cycle
      if (this%nurweight(n) >= done) cycle
      if (this%area(n) > dzero) then
        qtolfact = delt / this%area(n)
      else
        qtolfact = dzero
      end if
      ! -- qsim0 and ratesim both come from maw_fc (before the solve), so dq is
      !    the rate change from the previous outer iteration -- one step behind
      !    the newest heads. That is fine here: a lagged rate can only delay
      !    convergence, never accept a bad answer, because the solver already
      !    checks the heads and steady heads mean a steady rate. The up-to-date
      !    rate is not used because recomputing it would change the well's
      !    shutoff and RATE_SCALING state.
      dq = (this%qsim0(n) - this%ratesim(n)) * qtolfact
      if (abs(dq) > abs(dpakmax)) then
        dpakmax = dq
        locdpak = n
      end if
    end do
    if (locdpak > 0 .and. abs(dpakmax) > abs(dpak)) then
      ipak = locdpak
      dpak = dpakmax
      write (cloc, "(a,'-',a)") trim(this%packName), 'rate'
      cpak = trim(cloc)
    end if
    !
    ! -- the over-demand warning is only written on the last outer iteration of
    !    a model that did not converge
    if (iend == 0) return
    if (icnvgmod /= 0) return
    !
    ! -- look at each active pumping (extraction) well
    do n = 1, this%nmawwells
      if (this%iboundpak(n) < 1) cycle
      if (this%rate(n) >= dzero) cycle
      botw = this%bot(n)
      !
      ! -- largest rate the aquifer could supply if the well head dropped all
      !    the way to the bottom of the well
      qmax = dzero
      do j = 1, this%ngwfnodes(n)
        jpos = this%get_jpos(n, j)
        if (this%iboundconn(jpos) == 0) cycle
        igwfnode = this%get_gwfnode(n, j)
        hgwf = this%xnew(igwfnode)
        bmaw = this%botscrn(jpos)
        hv = max(botw, bmaw)
        !
        ! -- connection saturation evaluated with the well head at the bottom of
        !    the well (hv), consistent with this maximum-supply estimate, rather
        !    than at the current well head
        call this%maw_calculate_saturation(n, j, igwfnode, sat, hv)
        cmaw = this%satcond(jpos) * sat
        !
        ! -- only count connections that can supply water. A connection whose
        !    aquifer head is below the well bottom would take water in, so skip
        !    it instead of letting it lower the maximum supply.
        qmax = qmax + cmaw * max(hgwf - hv, dzero)
      end do
      !
      ! -- warn if the well is asking for more than the aquifer can supply
      qreq = -this%rate(n)
      if (qmax >= dzero .and. qreq > qtol * qmax) then
        write (cloc, '(a, a, a, a, i0, a, i0, a)') trim(this%name_model), '-(', &
          trim(this%filtyp), '_', this%ibcnum, '-', n, ')'
        write (warnmsg, fmtwarn) trim(cloc), qreq, qmax
        call store_warning(warnmsg)
      end if
    end do
  end subroutine maw_cc
 
  !> @brief Calculate flows
  !<
  subroutine maw_cq(this, x, flowja, iadv)
    ! -- modules
    use tdismodule, only: delt
    use constantsmodule, only: lenboundname
    use budgetmodule, only: budgettype
    ! -- dummy
    class(mawtype), intent(inout) :: this
    real(DP), dimension(:), intent(in) :: x
    real(DP), dimension(:), contiguous, intent(inout) :: flowja
    integer(I4B), optional, intent(in) :: iadv
    ! -- local
    real(DP) :: rrate
    real(DP) :: ss
    real(DP) :: ssold
    ! -- for budget
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: ibnd
    real(DP) :: hmaw
    real(DP) :: cfw
    ! -- for observations
    ! -- formats
    !
    ! -- recalculate package HCOF and RHS terms with latest groundwater and
    !    maw heads prior to calling base budget functionality
    call this%maw_cfupdate()
    !
    ! -- call base functionality in bnd_cq.  This will calculate maw-gwf flows
    !    and put them into this%simvals
    call this%BndType%bnd_cq(x, flowja, iadv=1)
    !
    ! -- calculate maw budget flow and storage terms
    do n = 1, this%nmawwells
      this%qout(n) = dzero
      this%qsto(n) = dzero
      if (this%iflowingwells > 0) then
        this%qfw(n) = dzero
      end if
      if (this%iboundpak(n) == 0) then
        cycle
      end if
      !
      ! -- set hmaw and xsto
      hmaw = this%xnewpak(n)
      this%xsto(n) = hmaw
      !
      ! -- add pumping rate to active maw well
      rrate = this%ratesim(n)
      !
      ! -- If flow is out of maw set qout to rrate.
      if (rrate < dzero) then
        this%qout(n) = rrate
      end if
      !
      ! -- add flowing well
      if (this%iflowingwells > 0) then
        if (this%fwcond(n) > dzero) then
          cfw = this%fwcondsim(n)
          ! -- only raise the storage level to the flowing-well elevation when
          !    the well is actually discharging (ifwdischarge == 1), matching
          !    the storage term assembled in maw_fc
          if (this%ifwdischarge(n) == 1) then
            this%xsto(n) = this%fwelev(n)
          end if
          rrate = cfw * (this%fwelev(n) - hmaw)
          this%qfw(n) = rrate
          !
          ! -- Subtract flowing well rrate from qout.
          this%qout(n) = this%qout(n) + rrate
        end if
      end if
      !
      ! -- Calculate qsto so it matches the storage term built in maw_fc. A
      !    flowing well that is actively discharging (ifwdischarge == 1) uses
      !    the flowing-well elevation (xsto = fwelev). Otherwise the storage
      !    water level is not allowed to drop below the bottom of the well.
      !    Testing ifwdischarge (not fwcond) keeps the budget in step with the
      !    matrix when a flowing well is at or below its discharge elevation.
      if (this%imawiss /= 1) then
        if (this%iflowingwells > 0 .and. this%ifwdischarge(n) == 1) then
          rrate = -this%area(n) * (this%xsto(n) - this%xoldsto(n)) / delt
        else
          ss = squadratic0sp(hmaw, this%bot(n), this%satomega)
          ssold = squadratic0sp(this%xoldsto(n), this%bot(n), this%satomega)
          rrate = -this%area(n) * (ss - ssold) / delt
        end if
        this%qsto(n) = rrate
      end if
    end do
    !
    ! -- gwf and constant flow
    ibnd = 1
    do n = 1, this%nmawwells
      hmaw = this%xnewpak(n)
      this%qconst(n) = dzero
      do j = 1, this%ngwfnodes(n)
        rrate = -this%simvals(ibnd)
        this%qleak(ibnd) = rrate
        if (this%iboundpak(n) < 0) then
          this%qconst(n) = this%qconst(n) - rrate
          !
          ! -- If flow is out increment qout by -rrate.
          if (-rrate < dzero) then
            this%qout(n) = this%qout(n) - rrate
          end if
        end if
        !
        ! -- increment ibnd counter
        ibnd = ibnd + 1
      end do
      !
      ! -- add additional flow terms to constant head term
      if (this%iboundpak(n) < 0) then
        !
        ! -- add well pumping rate
        this%qconst(n) = this%qconst(n) - this%ratesim(n)
        !
        ! -- add flowing well rate
        if (this%iflowingwells > 0) then
          this%qconst(n) = this%qconst(n) - this%qfw(n)
        end if
        !
        ! -- add storage term
        if (this%imawiss /= 1) then
          this%qconst(n) = this%qconst(n) - this%qsto(n)
        end if
      end if
    end do
    !
    ! -- fill the budget object
    call this%maw_fill_budobj()
  end subroutine maw_cq
 
  !> @brief Write flows to binary file and/or print flows to budget
  !<
  subroutine maw_ot_model_flows(this, icbcfl, ibudfl, icbcun, imap)
    ! -- dummy
    class(mawtype) :: this
    integer(I4B), intent(in) :: icbcfl
    integer(I4B), intent(in) :: ibudfl
    integer(I4B), intent(in) :: icbcun
    integer(I4B), dimension(:), optional, intent(in) :: imap
    !
    ! -- write the flows from the budobj
    call this%BndType%bnd_ot_model_flows(icbcfl, ibudfl, icbcun, this%imap)
  end subroutine maw_ot_model_flows
 
  !> @brief Output MAW package flow terms.
  !<
  subroutine maw_ot_package_flows(this, icbcfl, ibudfl)
    use tdismodule, only: kstp, kper, delt, pertim, totim
    class(mawtype) :: this
    integer(I4B), intent(in) :: icbcfl
    integer(I4B), intent(in) :: ibudfl
    integer(I4B) :: ibinun
    !
    ! -- write the flows from the budobj
    ibinun = 0
    if (this%ibudgetout /= 0) then
      ibinun = this%ibudgetout
    end if
    if (icbcfl == 0) ibinun = 0
    if (ibinun > 0) then
      call this%budobj%save_flows(this%dis, ibinun, kstp, kper, delt, &
                                  pertim, totim, this%iout)
    end if
    !
    ! -- Print maw flows table
    if (ibudfl /= 0 .and. this%iprflow /= 0) then
      call this%budobj%write_flowtable(this%dis, kstp, kper)
    end if
  end subroutine maw_ot_package_flows
 
  !> @brief Save maw-calculated values to binary file
  !<
  subroutine maw_ot_dv(this, idvsave, idvprint)
    use tdismodule, only: kstp, kper, pertim, totim
    use constantsmodule, only: dhnoflo, dhdry
    use inputoutputmodule, only: ulasav
    class(mawtype) :: this
    integer(I4B), intent(in) :: idvsave
    integer(I4B), intent(in) :: idvprint
    integer(I4B) :: ibinun
    integer(I4B) :: n
    real(DP) :: v
    real(DP) :: d
    !
    ! -- set unit number for binary dependent variable output
    ibinun = 0
    if (this%iheadout /= 0) then
      ibinun = this%iheadout
    end if
    if (idvsave == 0) ibinun = 0
    !
    ! -- write maw binary output
    if (ibinun > 0) then
      do n = 1, this%nmawwells
        v = this%xnewpak(n)
        d = v - this%bot(n)
        if (this%iboundpak(n) == 0) then
          v = dhnoflo
        else if (d <= dzero) then
          v = dhdry
        end if
        this%dbuff(n) = v
      end do
      call ulasav(this%dbuff, '            HEAD', &
                  kstp, kper, pertim, totim, &
                  this%nmawwells, 1, 1, ibinun)
    end if
    !
    ! -- write maw head table
    if (idvprint /= 0 .and. this%iprhed /= 0) then
      !
      ! -- set table kstp and kper
      call this%headtab%set_kstpkper(kstp, kper)
      !
      ! -- fill stage data
      do n = 1, this%nmawwells
        if (this%inamedbound == 1) then
          call this%headtab%add_term(this%cmawname(n))
        end if
        call this%headtab%add_term(n)
        call this%headtab%add_term(this%xnewpak(n))
      end do
    end if
  end subroutine maw_ot_dv
 
  !> @brief Write MAW budget to listing file
  !<
  subroutine maw_ot_bdsummary(this, kstp, kper, iout, ibudfl)
    ! -- module
    use tdismodule, only: totim, delt
    ! -- dummy
    class(mawtype) :: this !< MawType object
    integer(I4B), intent(in) :: kstp !< time step number
    integer(I4B), intent(in) :: kper !< period number
    integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
    integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
    !
    call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim, delt)
  end subroutine maw_ot_bdsummary
 
  !> @brief Deallocate memory
  !<
  subroutine maw_da(this)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(mawtype) :: this
    ! -- local
    !
    ! -- budobj
    call this%budobj%budgetobject_da()
    deallocate (this%budobj)
    nullify (this%budobj)
    !
    ! -- head table
    if (this%iprhed > 0) then
      call this%headtab%table_da()
      deallocate (this%headtab)
      nullify (this%headtab)
    end if
    !
    ! -- character arrays
    call mem_deallocate(this%cmawbudget, 'CMAWBUDGET', this%memoryPath)
    call mem_deallocate(this%cmawname, 'CMAWNAME', this%memoryPath)
    call mem_deallocate(this%status, 'STATUS', this%memoryPath)
    !
    ! -- deallocate well data pointers in memory manager
    call mem_deallocate(this%ngwfnodes)
    call mem_deallocate(this%ieqn)
    call mem_deallocate(this%ishutoff)
    call mem_deallocate(this%ifwdischarge)
    call mem_deallocate(this%strt)
    call mem_deallocate(this%radius)
    call mem_deallocate(this%area)
    call mem_deallocate(this%pumpelev)
    call mem_deallocate(this%bot)
    call mem_deallocate(this%ratesim)
    call mem_deallocate(this%qsim0)
    call mem_deallocate(this%reduction_length)
    call mem_deallocate(this%fwelev)
    call mem_deallocate(this%fwcond)
    call mem_deallocate(this%fwrlen)
    call mem_deallocate(this%fwcondsim)
    call mem_deallocate(this%xsto)
    call mem_deallocate(this%xoldsto)
    call mem_deallocate(this%shutoffmin)
    call mem_deallocate(this%shutoffmax)
    call mem_deallocate(this%shutofflevel)
    call mem_deallocate(this%shutoffweight)
    call mem_deallocate(this%shutoffdq)
    call mem_deallocate(this%shutoffqold)
    call mem_deallocate(this%nurdxold)
    call mem_deallocate(this%nurweight)
    !
    ! -- input-context-owned aliases, not package-allocated
    nullify (this%rate)
    nullify (this%well_head)
    call mem_deallocate(this%well_head_set)
    !
    ! -- connection data
    call mem_deallocate(this%iaconn)
    call mem_deallocate(this%gwfnodes)
    call mem_deallocate(this%sradius)
    call mem_deallocate(this%hk)
    call mem_deallocate(this%iboundconn)
    call mem_deallocate(this%satcond)
    call mem_deallocate(this%simcond)
    call mem_deallocate(this%topscrn)
    call mem_deallocate(this%botscrn)
    call mem_deallocate(this%angle)
    call mem_deallocate(this%connlen)
    call mem_deallocate(this%usrtopscrn)
    call mem_deallocate(this%usrbotscrn)
    !
    ! -- imap vector
    call mem_deallocate(this%imap)
    call mem_deallocate(this%dbuff)
    call mem_deallocate(this%cauxcbc, 'CAUXCBC', this%memoryPath)
    call mem_deallocate(this%qauxcbc)
    call mem_deallocate(this%qleak)
    call mem_deallocate(this%qfw)
    call mem_deallocate(this%qout)
    call mem_deallocate(this%qsto)
    call mem_deallocate(this%qconst)
    call mem_deallocate(this%denseterms)
    call mem_deallocate(this%viscratios)
    call mem_deallocate(this%idxlocnode)
    call mem_deallocate(this%idxdglo)
    call mem_deallocate(this%idxoffdglo)
    call mem_deallocate(this%idxsymdglo)
    call mem_deallocate(this%idxsymoffdglo)
    call mem_deallocate(this%xoldpak)
    !
    ! -- nullify pointers
    call mem_deallocate(this%xnewpak, 'HEAD', this%memoryPath)
    !
    ! -- scalars
    call mem_deallocate(this%correct_flow)
    call mem_deallocate(this%iprhed)
    call mem_deallocate(this%iheadout)
    call mem_deallocate(this%ibudgetout)
    call mem_deallocate(this%ibudcsv)
    call mem_deallocate(this%iflowingwells)
    call mem_deallocate(this%imawiss)
    call mem_deallocate(this%imawissopt)
    call mem_deallocate(this%nmawwells)
    call mem_deallocate(this%check_attr)
    call mem_deallocate(this%ishutoffcnt)
    call mem_deallocate(this%ieffradopt)
    call mem_deallocate(this%inonvert)
    call mem_deallocate(this%ioutredflowcsv)
    call mem_deallocate(this%satomega)
    call mem_deallocate(this%bditems)
    call mem_deallocate(this%theta)
    call mem_deallocate(this%kappa)
    call mem_deallocate(this%cbcauxitems)
    call mem_deallocate(this%idense)
    !
    ! -- pointers to gwf variables
    nullify (this%gwfiss)
    !
    ! -- nullify MAW input context pointers
    call this%input%destroy()
    !
    ! -- deallocate base
    call this%BndExtType%bnd_da()
  end subroutine maw_da
 
  !> @brief Define the list heading that is written to iout when PRINT_INPUT
  !! option is used.
  !<
  subroutine define_listlabel(this)
    class(mawtype), intent(inout) :: this
    !
    ! -- 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
 
  !> @brief Set pointers to model arrays and variables so that a package has
  !! has access to these things.
  !<
  subroutine maw_set_pointers(this, neq, ibound, xnew, xold, flowja)
    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_checkin
    ! -- dummy
    class(mawtype) :: this
    integer(I4B), pointer :: neq
    integer(I4B), dimension(:), pointer, contiguous :: ibound
    real(DP), dimension(:), pointer, contiguous :: xnew
    real(DP), dimension(:), pointer, contiguous :: xold
    real(DP), dimension(:), pointer, contiguous :: flowja
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: istart, iend
    !
    ! -- call base BndType set_pointers
    call this%BndType%set_pointers(neq, ibound, xnew, xold, flowja)
    !
    ! -- Set the MAW pointers
    !
    ! -- set package pointers
    istart = this%dis%nodes + this%ioffset + 1
    iend = istart + this%nmawwells - 1
    this%iboundpak => this%ibound(istart:iend)
    this%xnewpak => this%xnew(istart:iend)
    call mem_checkin(this%xnewpak, 'HEAD', this%memoryPath, 'X', &
                     this%memoryPathModel)
    call mem_allocate(this%xoldpak, this%nmawwells, 'XOLDPAK', this%memoryPath)
    !
    ! -- initialize xnewpak
    do n = 1, this%nmawwells
      this%xnewpak(n) = dep20
    end do
  end subroutine maw_set_pointers
 
  ! -- Procedures related to observations (type-bound)
 
  !> @brief Return true because MAW package supports observations
  !!
  !! Overrides BndType%bnd_obs_supported()
  !<
  logical function maw_obs_supported(this)
    class(mawtype) :: this
    !
    maw_obs_supported = .true.
  end function maw_obs_supported
 
  !> @brief Store observation type supported by MAW package
  !!
  !! Overrides BndType%bnd_df_obs
  !<
  subroutine maw_df_obs(this)
    ! -- dummy
    class(mawtype) :: this
    ! -- local
    integer(I4B) :: indx
    !
    ! -- Store obs type and assign procedure pointer
    !    for head observation type.
    call this%obs%StoreObsType('head', .false., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for frommvr observation type.
    call this%obs%StoreObsType('from-mvr', .false., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for conn-rate observation type.
    call this%obs%StoreObsType('maw', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for rate observation type.
    call this%obs%StoreObsType('rate', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for rate-to-mvr observation type.
    call this%obs%StoreObsType('rate-to-mvr', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for fw-rate observation type.
    call this%obs%StoreObsType('fw-rate', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for rate-to-mvr observation type.
    call this%obs%StoreObsType('fw-to-mvr', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for storage observation type.
    call this%obs%StoreObsType('storage', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for constant observation type.
    call this%obs%StoreObsType('constant', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for cond observation type.
    call this%obs%StoreObsType('conductance', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for fw-conductance observation type.
    call this%obs%StoreObsType('fw-conductance', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => maw_process_obsid
  end subroutine maw_df_obs
 
  !> @brief Calculate observations this time step and call
  !! ObsType%SaveOneSimval for each MawType observation.
  !<
  subroutine maw_bd_obs(this)
    ! -- dummy
    class(mawtype) :: this
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: jj
    integer(I4B) :: n
    integer(I4B) :: nn
    integer(I4B) :: jpos
    real(DP) :: cmaw
    real(DP) :: hmaw
    real(DP) :: v
    real(DP) :: qfact
    type(observetype), pointer :: obsrv => null()
    !
    ! Calculate, save, and write simulated values for all MAW observations
    if (this%obs%npakobs > 0) then
      call this%obs%obs_bd_clear()
      do i = 1, this%obs%npakobs
        obsrv => this%obs%pakobs(i)%obsrv
        do j = 1, obsrv%indxbnds_count
          v = dnodata
          jj = obsrv%indxbnds(j)
          select case (obsrv%ObsTypeId)
          case ('HEAD')
            if (this%iboundpak(jj) /= 0) then
              v = this%xnewpak(jj)
            end if
          case ('FROM-MVR')
            if (this%iboundpak(jj) /= 0) then
              if (this%imover == 1) then
                v = this%pakmvrobj%get_qfrommvr(jj)
              end if
            end if
          case ('MAW')
            ! -- an inactive connection reports DNODATA, as an inactive well
            n = this%imap(jj)
            nn = jj - this%iaconn(n) + 1
            jpos = this%get_jpos(n, nn)
            if (this%iboundpak(n) /= 0 .and. this%iboundconn(jpos) /= 0) then
              v = this%qleak(jj)
            end if
          case ('RATE')
            if (this%iboundpak(jj) /= 0) then
              v = this%ratesim(jj)
              if (v < dzero .and. this%qout(jj) < dzero) then
                qfact = v / this%qout(jj)
                if (this%imover == 1) then
                  v = v + this%pakmvrobj%get_qtomvr(jj) * qfact
                end if
              end if
            end if
          case ('RATE-TO-MVR')
            if (this%iboundpak(jj) /= 0) then
              if (this%imover == 1) then
                v = this%ratesim(jj)
                qfact = dzero
                if (v < dzero .and. this%qout(jj) < dzero) then
                  qfact = v / this%qout(jj)
                end if
                v = this%pakmvrobj%get_qtomvr(jj) * qfact
                if (v > dzero) then
                  v = -v
                end if
              end if
            end if
          case ('FW-RATE')
            if (this%iboundpak(jj) /= 0 .and. this%iflowingwells > 0) then
              hmaw = this%xnewpak(jj)
              cmaw = this%fwcondsim(jj)
              v = cmaw * (this%fwelev(jj) - hmaw)
              if (v < dzero .and. this%qout(jj) < dzero) then
                qfact = v / this%qout(jj)
                if (this%imover == 1) then
                  v = v + this%pakmvrobj%get_qtomvr(jj) * qfact
                end if
              end if
            end if
          case ('FW-TO-MVR')
            if (this%iboundpak(jj) /= 0 .and. this%iflowingwells > 0) then
              if (this%imover == 1) then
                hmaw = this%xnewpak(jj)
                cmaw = this%fwcondsim(jj)
                v = cmaw * (this%fwelev(jj) - hmaw)
                qfact = dzero
                if (v < dzero .and. this%qout(jj) < dzero) then
                  qfact = v / this%qout(jj)
                end if
                v = this%pakmvrobj%get_qtomvr(jj) * qfact
                if (v > dzero) then
                  v = -v
                end if
              end if
            end if
          case ('STORAGE')
            if (this%iboundpak(jj) /= 0 .and. this%imawissopt /= 1) then
              v = this%qsto(jj)
            end if
          case ('CONSTANT')
            if (this%iboundpak(jj) /= 0) then
              v = this%qconst(jj)
            end if
          case ('CONDUCTANCE')
            n = this%imap(jj)
            nn = jj - this%iaconn(n) + 1
            jpos = this%get_jpos(n, nn)
            if (this%iboundpak(n) /= 0 .and. this%iboundconn(jpos) /= 0) then
              v = this%simcond(jpos)
            end if
          case ('FW-CONDUCTANCE')
            if (this%iboundpak(jj) /= 0) then
              v = this%fwcondsim(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
      end do
      !
      ! -- write summary of error messages
      if (count_errors() > 0) then
        call store_error_filename(this%input_fname)
      end if
    end if
    !
    ! -- Write the MAW reduced flows to csv file entries for this step
    if (this%ioutredflowcsv > 0) then
      call this%maw_redflow_csv_write()
    end if
  end subroutine maw_bd_obs
 
  !> @brief Process each observation
  !!
  !! Only done the first stress period since boundaries are fixed for the
  !! simulation
  !<
  subroutine maw_rp_obs(this)
    use tdismodule, only: kper
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: nn1
    integer(I4B) :: nn2
    integer(I4B) :: jj
    character(len=LENBOUNDNAME) :: bname
    logical :: jfound
    class(observetype), pointer :: obsrv => null()
    ! -- formats
10  format('Boundary "', a, '" for observation "', a, &
           '" is invalid in package "', a, '"')
    !
    if (kper == 1) then
      do i = 1, this%obs%npakobs
        obsrv => this%obs%pakobs(i)%obsrv
        !
        ! -- get node number 1
        nn1 = obsrv%NodeNumber
        if (nn1 == namedboundflag) then
          bname = obsrv%FeatureName
          if (bname /= '') then
            ! -- Observation maw is based on a boundary name.
            !    Iterate through all multi-aquifer wells to identify and store
            !    corresponding index in bound array.
            jfound = .false.
            if (obsrv%ObsTypeId == 'MAW' .or. &
                obsrv%ObsTypeId == 'CONDUCTANCE') then
              do j = 1, this%nmawwells
                do jj = this%iaconn(j), this%iaconn(j + 1) - 1
                  if (this%boundname(jj) == bname) then
                    jfound = .true.
                    call obsrv%AddObsIndex(jj)
                  end if
                end do
              end do
            else
              do j = 1, this%nmawwells
                if (this%cmawname(j) == bname) then
                  jfound = .true.
                  call obsrv%AddObsIndex(j)
                end if
              end do
            end if
            if (.not. jfound) then
              write (errmsg, 10) &
                trim(bname), trim(obsrv%Name), trim(this%packName)
              call store_error(errmsg)
            end if
          end if
        else
          if (obsrv%indxbnds_count == 0) then
            if (obsrv%ObsTypeId == 'MAW' .or. &
                obsrv%ObsTypeId == 'CONDUCTANCE') then
              nn2 = obsrv%NodeNumber2
              j = this%iaconn(nn1) + nn2 - 1
              call obsrv%AddObsIndex(j)
            else
              call obsrv%AddObsIndex(nn1)
            end if
          else
            errmsg = 'Programming error in maw_rp_obs'
            call store_error(errmsg)
          end if
        end if
        !
        ! -- catch non-cumulative observation assigned to observation defined
        !    by a boundname that is assigned to more than one element
        if (obsrv%ObsTypeId == 'HEAD') then
          if (obsrv%indxbnds_count > 1) then
            write (errmsg, '(a,3(1x,a))') &
              trim(adjustl(obsrv%ObsTypeId)), &
              'for observation', trim(adjustl(obsrv%Name)), &
              'must be assigned to a multi-aquifer well with a unique boundname.'
            call store_error(errmsg)
          end if
        end if
        !
        ! -- check that index values are valid
        if (obsrv%ObsTypeId == 'MAW' .or. &
            obsrv%ObsTypeId == 'CONDUCTANCE') then
          do j = 1, obsrv%indxbnds_count
            nn1 = obsrv%indxbnds(j)
            n = this%imap(nn1)
            nn2 = nn1 - this%iaconn(n) + 1
            jj = this%iaconn(n + 1) - this%iaconn(n)
            if (nn1 < 1 .or. nn1 > this%maxbound) then
              write (errmsg, '(3(a,1x),i0,1x,a,i0,a)') &
                trim(adjustl(obsrv%ObsTypeId)), &
                'multi-aquifer well connection number must be greater than 0', &
                'and less than', jj, '(specified value is ', nn2, ').'
              call store_error(errmsg)
            end if
          end do
        else
          do j = 1, obsrv%indxbnds_count
            nn1 = obsrv%indxbnds(j)
            if (nn1 < 1 .or. nn1 > this%nmawwells) then
              write (errmsg, '(3(a,1x),i0,1x,a,i0,a)') &
                trim(adjustl(obsrv%ObsTypeId)), &
                'multi-aquifer well must be greater than 0 ', &
                'and less than or equal to', this%nmawwells, &
                '(specified value is ', nn1, ').'
              call store_error(errmsg)
            end if
          end do
        end if
      end do
      !
      ! -- evaluate if there are any observation errors
      if (count_errors() > 0) then
        call store_error_filename(this%input_fname)
      end if
    end if
  end subroutine maw_rp_obs
 
  !
  ! -- Procedures related to observations (NOT type-bound)
 
  !> @brief This procedure is pointed to by ObsDataType%ProcesssIdPtr. It
  !! processes the ID string of an observation definition for MAW package
  !! observations.
  !<
  subroutine maw_process_obsid(obsrv, dis, inunitobs, iout)
    ! -- dummy
    type(observetype), intent(inout) :: obsrv
    class(disbasetype), intent(in) :: dis
    integer(I4B), intent(in) :: inunitobs
    integer(I4B), intent(in) :: iout
    ! -- local
    integer(I4B) :: nn1, nn2
    integer(I4B) :: icol, istart, istop
    character(len=LINELENGTH) :: string
    character(len=LENBOUNDNAME) :: bndname
    ! formats
    !
    string = obsrv%IDstring
    ! -- Extract multi-aquifer well number from string and store it.
    !    If 1st item is not an integer(I4B), it should be a
    !    maw name--deal with it.
    icol = 1
    ! -- get multi-aquifer well number or boundary name
    call extract_idnum_or_bndname(string, icol, istart, istop, nn1, bndname)
    if (nn1 == namedboundflag) then
      obsrv%FeatureName = bndname
    else
      if (obsrv%ObsTypeId == 'MAW' .or. &
          obsrv%ObsTypeId == 'CONDUCTANCE') then
        call extract_idnum_or_bndname(string, icol, istart, istop, nn2, bndname)
        if (len_trim(bndname) < 1 .and. nn2 < 0) then
          write (errmsg, '(a,1x,a,a,1x,a,1x,a)') &
            'For observation type', trim(adjustl(obsrv%ObsTypeId)), &
            ', ID given as an integer and not as boundname,', &
            'but ID2 (icon) is missing.  Either change ID to valid', &
            'boundname or supply valid entry for ID2.'
          call store_error(errmsg)
        end if
        if (nn2 == namedboundflag) then
          obsrv%FeatureName = bndname
          ! -- reset nn1
          nn1 = nn2
        else
          obsrv%NodeNumber2 = nn2
        end if
      end if
    end if
    ! -- store multi-aquifer well number (NodeNumber)
    obsrv%NodeNumber = nn1
  end subroutine maw_process_obsid
 
  !
  ! -- private MAW methods
 
  !> @brief Initialize the auto flow reduce csv output file
  !<
  subroutine maw_redflow_csv_init(this, fname)
    ! -- dummy variables
    class(mawtype), intent(inout) :: this !< MawType object
    character(len=*), intent(in) :: fname
    ! -- format
    character(len=*), parameter :: fmtredflowcsv = &
      "(4x, 'MAW REDUCED FLOW INFORMATION WILL BE SAVED TO FILE: ', a, /4x, &
    &'OPENED ON UNIT: ', I0)"
 
    this%ioutredflowcsv = getunit()
    call openfile(this%ioutredflowcsv, this%iout, fname, 'CSV', &
                  filstat_opt='REPLACE')
    write (this%iout, fmtredflowcsv) trim(adjustl(fname)), &
      this%ioutredflowcsv
    write (this%ioutredflowcsv, '(a)') &
      'time,period,step,MAWnumber,rate-requested,rate-actual,maw-reduction'
  end subroutine maw_redflow_csv_init
 
  !> @brief MAW reduced flows only when & where they occur
  !<
  subroutine maw_redflow_csv_write(this)
    ! -- modules
    use tdismodule, only: totim, kstp, kper
    ! -- dummy variables
    class(mawtype), intent(inout) :: this !< MawType object
    ! -- local
    integer(I4B) :: n
    !integer(I4B) :: nodereduced
    !integer(I4B) :: nodeuser
    real(DP) :: v
    ! -- format
    do n = 1, this%nmawwells
      !
      ! -- test if node is constant or inactive
      if (this%status(n) .ne. 'ACTIVE') then
        cycle
      end if
      v = this%rate(n) - this%ratesim(n) !reductions in extraction will be negative and reductions in injection will be positive; follows convention of WEL AUTO_FLOW_REDUCE_CSV
      if (abs(v) > dem9) then !need to check absolute value of difference for both extraction and injection; using 1e-9 as epsilon value but could be tweaked
        write (this%ioutredflowcsv, '(*(G0,:,","))') &
          totim, kper, kstp, n, this%rate(n), this%ratesim(n), v
      end if
    end do
  end subroutine maw_redflow_csv_write
 
  !> @brief Calculate the appropriate saturated conductance to use based on
  !! aquifer and multi-aquifer well characteristics
  !<
  subroutine maw_calculate_satcond(this, i, j, node)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: i
    integer(I4B), intent(in) :: j
    integer(I4B), intent(in) :: node
    ! -- local
    integer(I4B) :: iTcontrastErr
    integer(I4B) :: jpos
    real(DP) :: c
    real(DP) :: k11
    real(DP) :: k22
    real(DP) :: sqrtk11k22
    real(DP) :: hks
    real(DP) :: area
    real(DP) :: eradius
    real(DP) :: topw
    real(DP) :: botw
    real(DP) :: tthkw
    real(DP) :: tthka
    real(DP) :: Tcontrast
    real(DP) :: skin
    real(DP) :: ravg
    real(DP) :: slen
    real(DP) :: pavg
    real(DP) :: gwfsat
    real(DP) :: gwftop
    real(DP) :: gwfbot
    real(DP) :: lc1
    real(DP) :: lc2
    real(DP) :: dx
    real(DP) :: dy
    real(DP) :: Txx
    real(DP) :: Tyy
    real(DP) :: T2pi
    real(DP) :: yx4
    real(DP) :: xy4
    ! -- formats
    !
    ! -- initialize conductance variables
    itcontrasterr = 0
    lc1 = dzero
    lc2 = dzero
    !
    ! -- calculate connection position
    jpos = this%get_jpos(i, j)
    !
    ! -- set K11 and K22
    k11 = this%gwfk11(node)
    if (this%gwfik22 == 0) then
      k22 = this%gwfk11(node)
    else
      k22 = this%gwfk22(node)
    end if
    sqrtk11k22 = sqrt(k11 * k22)
    !
    ! -- set gwftop, gwfbot, and gwfsat
    gwftop = this%dis%top(node)
    gwfbot = this%dis%bot(node)
    tthka = gwftop - gwfbot
    gwfsat = this%gwfsat(node)
    !
    ! -- set top and bottom of well screen
    c = dzero
    topw = this%topscrn(jpos)
    botw = this%botscrn(jpos)
    tthkw = topw - botw
    !
    ! -- scale screen thickness using gwfsat (for NPF Package THICKSTRT)
    if (gwftop == topw .and. gwfbot == botw) then
      if (this%icelltype(node) == 0) then
        tthkw = tthkw * gwfsat
        tthka = tthka * gwfsat
      end if
    end if
    !
    ! -- calculate the aquifer transmissivity (T2pi)
    t2pi = dtwopi * tthka * sqrtk11k22
    !
    ! -- calculate effective radius
    if (this%dis%ndim == 3 .and. this%ieffradopt /= 0) then
      txx = k11 * tthka
      tyy = k22 * tthka
      dx = sqrt(this%dis%area(node))
      dy = dx
      yx4 = (tyy / txx)**dquarter
      xy4 = (txx / tyy)**dquarter
      eradius = 0.28_dp * ((yx4 * dx)**dtwo + &
                           (xy4 * dy)**dtwo)**dhalf / (yx4 + xy4)
    else
      area = this%dis%area(node)
      eradius = sqrt(area / (deight * dpi))
    end if
    !
    ! -- conductance calculations
    ! -- Thiem equation (1) and cumulative Thiem and skin equations (3)
    if (this%ieqn(i) == 1 .or. this%ieqn(i) == 3) then
      lc1 = log(eradius / this%radius(i)) / t2pi
    end if
    !
    ! -- skin equation (2) and cumulative Thiem and skin equations (3)
    if (this%ieqn(i) == 2 .or. this%ieqn(i) == 3) then
      hks = this%hk(jpos)
      if (tthkw * hks > dzero) then
        tcontrast = (sqrtk11k22 * tthka) / (hks * tthkw)
        skin = (tcontrast - done) * log(this%sradius(jpos) / this%radius(i))
        !
        ! -- trap invalid transmissvity contrast if using skin equation (2).
        !    Not trapped for cumulative Thiem and skin equations (3)
        !    because the MNW2 package allowed this condition (for
        !    backward compatibility with the MNW2 package for
        !    MODFLOW-2005, MODFLOW-NWT, and MODFLOW-USG).
        if (tcontrast <= 1 .and. this%ieqn(i) == 2) then
          itcontrasterr = 1
          write (errmsg, '(a,g0,a,1x,i0,1x,a,1x,i0,a,4(1x,a))') &
            'Invalid calculated transmissivity contrast (', tcontrast, &
            ') for maw well', i, 'connection', j, '.', 'This happens when the', &
            'skin transmissivity equals or exceeds the aquifer transmissivity.', &
            'Consider decreasing HK_SKIN for the connection or using the', &
            'CUMULATIVE or MEAN conductance equations.'
          call store_error(errmsg)
        else
          lc2 = skin / t2pi
        end if
      end if
    end if
    ! -- conductance using screen elevations, hk, well radius,
    !    and screen radius
    if (this%ieqn(i) == 4) then
      hks = this%hk(jpos)
      ravg = dhalf * (this%radius(i) + this%sradius(jpos))
      slen = this%sradius(jpos) - this%radius(i)
      pavg = dtwopi * ravg
      c = hks * pavg * tthkw / slen
    end if
    !
    ! -- calculate final conductance for Thiem (1), Skin (2), and
    ! and cumulative Thiem and skin equations (3)
    if (this%ieqn(i) < 4) then
      if (lc1 + lc2 /= dzero) then
        c = done / (lc1 + lc2)
      else
        c = -dnodata
      end if
    end if
    !
    ! -- ensure that the conductance is not negative. Only write error message
    !    if error condition has not occurred for skin calculations (LC2)
    if (c < dzero .and. itcontrasterr == 0) then
      write (errmsg, '(a,g0,a,1x,i0,1x,a,1x,i0,a,4(1x,a))') &
        'Invalid calculated negative conductance (', c, &
        ') for maw well', i, 'connection', j, '.', 'this happens when the', &
        'skin transmissivity equals or exceeds the aquifer transmissivity.', &
        'consider decreasing hk_skin for the connection or using the', &
        'mean conductance equation.'
      call store_error(errmsg)
    end if
    !
    ! -- scale the saturated conductance by the length correction factor for
    !    non-vertical (slanted) connections. The factor is 1.0 for vertical
    !    connections, so this has no effect unless an ANGLEDATA block was read.
    if (this%inonvert /= 0) then
      c = c * this%maw_calc_lcorr(i, jpos)
    end if
    !
    ! -- set saturated conductance
    this%satcond(jpos) = c
  end subroutine maw_calculate_satcond
 
  !> @brief Calculate the saturation between the aquifer maw well_head
  !<
  subroutine maw_calculate_saturation(this, n, j, node, sat, hwell_in)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: j
    integer(I4B), intent(in) :: node
    real(DP), intent(inout) :: sat
    real(DP), intent(in), optional :: hwell_in !well head to use instead of the current well head
    ! -- local
    integer(I4B) :: jpos
    real(DP) :: h_temp
    real(DP) :: hwell
    real(DP) :: topw
    real(DP) :: botw
    ! -- formats
    !
    ! -- initialize saturation
    sat = dzero
    !
    ! -- calculate current saturation for convertible cells
    if (this%icelltype(node) /= 0) then
      !
      ! -- set hwell (use the caller-supplied head if provided)
      hwell = this%xnewpak(n)
      if (present(hwell_in)) then
        hwell = hwell_in
      end if
      !
      ! -- set connection position
      jpos = this%get_jpos(n, j)
      !
      ! -- set top and bottom of the well connection
      topw = this%topscrn(jpos)
      botw = this%botscrn(jpos)
      !
      ! -- calculate appropriate saturation
      if (this%inewton /= 1) then
        h_temp = this%xnew(node)
        if (h_temp < botw) then
          h_temp = botw
        end if
        if (hwell < botw) then
          hwell = botw
        end if
        h_temp = dhalf * (h_temp + hwell)
      else
        h_temp = this%xnew(node)
        if (hwell > h_temp) then
          h_temp = hwell
        end if
        if (h_temp < botw) then
          h_temp = botw
        end if
      end if
      ! -- calculate saturation
      sat = squadraticsaturation(topw, botw, h_temp, this%satomega)
    else
      sat = done
    end if
  end subroutine maw_calculate_saturation
 
  !> @brief Calculate matrix terms for a multi-aquifer well connection. Terms
  !! for fc and fn methods are calculated based on whether term2 is passed
  !! Arguments are as follows:
  !!     n       : maw well number
  !!     j       : connection number for well n
  !!     icflow  : flag indicating that flow should be corrected
  !!     cmaw    : maw-gwf conducance
  !!     cterm   : correction term for flow to dry cell
  !!     term    : xxx
  !!     flow    : calculated flow for this connection, positive into well
  !!     term2   : xxx
  !<
  subroutine maw_calculate_conn_terms(this, n, j, icflow, cmaw, cterm, term, &
                                      flow, term2)
    ! -- dummy
    class(mawtype) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: j
    integer(I4B), intent(inout) :: icflow
    real(DP), intent(inout) :: cmaw
    real(DP), intent(inout) :: cterm
    real(DP), intent(inout) :: term
    real(DP), intent(inout) :: flow
    real(DP), intent(inout), optional :: term2
    ! -- local
    logical(LGP) :: correct_flow
    integer(I4B) :: inewton
    integer(I4B) :: jpos
    integer(I4B) :: igwfnode
    real(DP) :: hmaw
    real(DP) :: hgwf
    real(DP) :: hups
    real(DP) :: hdowns
    real(DP) :: sat
    real(DP) :: tmaw
    real(DP) :: bmaw
    real(DP) :: en
    real(DP) :: hbar
    real(DP) :: drterm
    real(DP) :: dhbarterm
    real(DP) :: vscratio
    !
    ! -- initialize terms
    cterm = dzero
    vscratio = done
    icflow = 0
    if (present(term2)) then
      inewton = 1
    else
      inewton = 0
    end if
    !
    ! -- set common terms
    jpos = this%get_jpos(n, j)
    igwfnode = this%get_gwfnode(n, j)
    hgwf = this%xnew(igwfnode)
    hmaw = this%xnewpak(n)
    tmaw = this%topscrn(jpos)
    bmaw = this%botscrn(jpos)
    !
    ! -- if vsc active, select appropriate viscosity ratio
    if (this%ivsc == 1) then
      ! flow out of well (flow is negative)
      if (flow < 0) then
        vscratio = this%viscratios(1, n)
      else
        vscratio = this%viscratios(2, n)
      end if
    end if
    !
    ! -- calculate saturation
    call this%maw_calculate_saturation(n, j, igwfnode, sat)
    cmaw = this%satcond(jpos) * vscratio * sat
    !
    ! -- set upstream head, term, and term2 if returning newton terms
    if (inewton == 1) then
      term = dzero
      term2 = dzero
      hups = hmaw
      if (hgwf > hups) then
        hups = hgwf
      end if
      !
      ! -- slope of the saturation with head. In a confined (non-convertible)
      !    cell the saturation is always one and does not change with head, so
      !    this slope must be zero. Otherwise the matrix terms would not match
      !    the flow that is actually calculated for the well.
      if (this%icelltype(igwfnode) /= 0) then
        drterm = squadraticsaturationderivative(tmaw, bmaw, hups, this%satomega)
      else
        drterm = dzero
      end if
    else
      term = cmaw
    end if
    !
    ! -- calculate correction term if flow_correction option specified
    if (this%correct_flow) then
      !
      ! -- set bmaw, determine en, and set correct_flow flag
      en = max(bmaw, this%dis%bot(igwfnode))
      correct_flow = .false.
      if (hmaw < en) then
        correct_flow = .true.
      end if
      if (hgwf < en .and. this%icelltype(igwfnode) /= 0) then
        correct_flow = .true.
      end if
      !
      ! -- if flow should be corrected because hgwf or hmaw is below bottom
      !    then calculate correction term (cterm)
      if (correct_flow) then
        icflow = 1
        hdowns = min(hmaw, hgwf)
        hbar = squadratic0sp(hdowns, en, this%satomega)
        if (hgwf > hmaw) then
          cterm = cmaw * (hmaw - hbar)
        else
          cterm = cmaw * (hbar - hgwf)
        end if
      end if
      !
      ! -- if newton formulation then calculate newton terms
      if (inewton /= 0) then
        !
        ! -- maw is upstream
        if (hmaw > hgwf) then
          hbar = squadratic0sp(hgwf, en, this%satomega)
          term = drterm * this%satcond(jpos) * vscratio * (hbar - hmaw)
          dhbarterm = squadratic0spderivative(hgwf, en, this%satomega)
          term2 = cmaw * (dhbarterm - done)
          !
          ! -- gwf is upstream
        else
          hbar = squadratic0sp(hmaw, en, this%satomega)
          term = -drterm * this%satcond(jpos) * vscratio * (hgwf - hbar)
          dhbarterm = squadratic0spderivative(hmaw, en, this%satomega)
          term2 = cmaw * (done - dhbarterm)
        end if
      end if
    else
      !
      ! -- flow is not corrected, so calculate term for newton formulation
      if (inewton /= 0) then
        term = drterm * this%satcond(jpos) * vscratio * (hgwf - hmaw)
      end if
    end if
    !
    ! -- calculate flow relative to maw for fc and bd
    flow = dzero
    if (inewton == 0) then
      flow = term * (hgwf - hmaw) + cterm
    end if
    !
    ! -- add density part here
    if (this%idense /= 0 .and. inewton == 0) then
      call this%maw_calculate_density_exchange(jpos, hmaw, hgwf, cmaw, &
                                               bmaw, flow, term, cterm)
    end if
  end subroutine maw_calculate_conn_terms
 
  !> @brief Calculate well pumping rate based on constraints
  !<
  subroutine maw_calculate_wellq(this, n, hmaw, q)
    ! -- dummy
    class(mawtype) :: this
    integer(I4B), intent(in) :: n
    real(DP), intent(in) :: hmaw
    real(DP), intent(inout) :: q
    ! -- local
    real(DP) :: scale
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: rate
    real(DP) :: weight
    real(DP) :: dq
    !
    ! -- Initialize q
    q = dzero
    !
    ! -- Assign rate as the user-provided base pumping rate
    rate = this%rate(n)
    !
    ! -- Assign q differently depending on whether this is an extraction well
    !    (rate < 0) or an injection well (rate > 0).
    if (rate < dzero) then
      !
      ! -- If well shut off is activated, then turn off well if necessary,
      !    or if shut off is not activated then check to see if rate scaling
      !    is on.
      if (this%shutofflevel(n) /= dep20) then
        call this%maw_calculate_qpot(n, q)
        if (q < dzero) q = dzero
        if (q > -rate) q = -rate
 
        if (this%ishutoffcnt == 1) then
          this%shutoffweight(n) = done
          this%shutoffdq(n) = dzero
          this%shutoffqold(n) = q
        end if
 
        dq = q - this%shutoffqold(n)
        weight = this%shutoffweight(n)
        !
        ! -- for oscillating condition, decrease factor
        if (this%shutoffdq(n) * dq < dzero) then
          weight = this%theta * this%shutoffweight(n)
          !
          ! -- when change is of same sign, increase factor
        else
          weight = this%shutoffweight(n) + this%kappa
        end if
        if (weight > done) weight = done
 
        q = this%shutoffqold(n) + weight * dq
 
        this%shutoffqold(n) = q
        this%shutoffdq(n) = dq
        this%shutoffweight(n) = weight
        !
        ! -- If shutoffmin and shutoffmax are specified then apply
        !    additional checks for when to shut off the well.
        if (this%shutoffmin(n) > dzero) then
          if (hmaw < this%shutofflevel(n)) then
            !
            ! -- calculate adjusted well rate subject to constraints
            ! -- well is shutoff
            if (this%ishutoff(n) /= 0) then
              q = dzero
              !
              ! --- well is not shut off
            else
              ! -- turn off well if q is less than the minimum rate and
              !    reset the ishutoff flag if at least on iteration 3
              if (q < this%shutoffmin(n)) then
                if (this%ishutoffcnt > 2) then
                  this%ishutoff(n) = 1
                end if
                q = dzero
                !
                ! -- leave well on and use the specified rate
                !    or the potential rate
              end if
            end if
            !
            ! -- try to use the specified rate or the potential rate
          else
            if (q > this%shutoffmax(n)) then
              if (this%ishutoffcnt <= 2) then
                this%ishutoff(n) = 0
              end if
            end if
            if (this%ishutoff(n) /= 0) then
              q = dzero
            end if
          end if
        end if
 
        if (q /= dzero) q = -q
 
      else
        scale = done
        !
        ! -- Apply rate scaling by reducing pumpage when hmaw is less than the
        !    sum of maw pump elevation (pumpelev) and the specified reduction
        !    length.  The rate will go to zero as hmaw drops to the pump
        !    elevation.
        if (this%reduction_length(n) /= dep20) then
          bt = this%pumpelev(n)
          tp = bt + this%reduction_length(n)
          scale = sqsaturation(tp, bt, hmaw)
        end if
        q = scale * rate
      end if
      !
    else
      !
      ! -- Handle the injection case (rate > 0) differently than extraction.
      q = rate
      if (this%shutofflevel(n) /= dep20) then
        call this%maw_calculate_qpot(n, q)
        q = -q
        if (q < dzero) q = dzero
        if (q > rate) q = rate
 
        if (this%ishutoffcnt == 1) then
          this%shutoffweight(n) = done
          this%shutoffdq(n) = dzero
          this%shutoffqold(n) = q
        end if
 
        dq = q - this%shutoffqold(n)
        weight = this%shutoffweight(n)
        !
        ! -- for oscillating condition, decrease factor
        if (this%shutoffdq(n) * dq < dzero) then
          weight = this%theta * this%shutoffweight(n)
          !
          ! -- when change is of same sign, increase factor
        else
          weight = this%shutoffweight(n) + this%kappa
        end if
        if (weight > done) weight = done
 
        q = this%shutoffqold(n) + weight * dq
 
        this%shutoffqold(n) = q
        this%shutoffdq(n) = dq
        this%shutoffweight(n) = weight
 
      else
        scale = done
        !
        ! -- Apply rate scaling for an injection well by reducing the
        !    injection rate as hmaw rises above the pump elevation.  The rate
        !    will approach zero as hmaw approaches pumpelev + reduction_length.
        if (this%reduction_length(n) /= dep20) then
          bt = this%pumpelev(n)
          tp = bt + this%reduction_length(n)
          scale = done - sqsaturation(tp, bt, hmaw)
        end if
        q = scale * rate
      end if
    end if
  end subroutine maw_calculate_wellq
 
  !> @brief Calculate groundwater inflow to a maw well
  !<
  subroutine maw_calculate_qpot(this, n, qnet)
    use tdismodule, only: delt
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: n
    real(DP), intent(inout) :: qnet
    ! -- local
    integer(I4B) :: j
    integer(I4B) :: jpos
    integer(I4B) :: igwfnode
    real(DP) :: bt
    real(DP) :: tp
    real(DP) :: scale
    real(DP) :: cfw
    real(DP) :: hdterm
    real(DP) :: sat
    real(DP) :: cmaw
    real(DP) :: hgwf
    real(DP) :: bmaw
    real(DP) :: h_temp
    real(DP) :: hv
    real(DP) :: vscratio
    ! -- format
    !
    ! -- initialize qnet and h_temp
    qnet = dzero
    vscratio = done
    h_temp = this%shutofflevel(n)
    !
    ! -- if vsc active, select appropriate viscosity ratio
    if (this%ivsc == 1) then
      ! flow out of well (flow is negative)
      if (qnet < 0) then
        vscratio = this%viscratios(1, n)
      else
        vscratio = this%viscratios(2, n)
      end if
    end if
    !
    ! -- calculate discharge to flowing wells
    if (this%iflowingwells > 0) then
      if (this%fwcond(n) > dzero) then
        bt = this%fwelev(n)
        tp = bt + this%fwrlen(n)
        scale = sqsaturation(tp, bt, h_temp)
        cfw = scale * this%fwcond(n) * this%viscratios(2, n)
        this%ifwdischarge(n) = 0
        if (cfw > dzero) then
          this%ifwdischarge(n) = 1
          this%xsto(n) = bt
        end if
        qnet = qnet + cfw * (bt - h_temp)
      end if
    end if
    !
    ! -- calculate maw storage changes
    if (this%imawiss /= 1) then
      if (this%ifwdischarge(n) /= 1) then
        hdterm = this%xoldsto(n) - h_temp
      else
        hdterm = this%xoldsto(n) - this%fwelev(n)
      end if
      qnet = qnet - (this%area(n) * hdterm / delt)
    end if
    !
    ! -- calculate inflow from aquifer
    do j = 1, this%ngwfnodes(n)
      jpos = this%get_jpos(n, j)
      if (this%iboundconn(jpos) == 0) cycle
      igwfnode = this%get_gwfnode(n, j)
      call this%maw_calculate_saturation(n, j, igwfnode, sat)
      cmaw = this%satcond(jpos) * vscratio * sat
      hgwf = this%xnew(igwfnode)
      bmaw = this%botscrn(jpos)
      hv = h_temp
      if (hv < bmaw) then
        hv = bmaw
      end if
      if (hgwf < bmaw) then
        hgwf = bmaw
      end if
      qnet = qnet + cmaw * (hgwf - hv)
    end do
  end subroutine maw_calculate_qpot
 
  !> @brief Update MAW satcond and package rhs and hcof
  !<
  subroutine maw_cfupdate(this)
    class(mawtype) :: this
    ! -- dummy
    ! -- local
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: jpos
    integer(I4B) :: icflow
    integer(I4B) :: ibnd
    real(DP) :: flow
    real(DP) :: cmaw
    real(DP) :: hmaw
    real(DP) :: cterm
    real(DP) :: term
    !
    ! -- Return if no maw wells
    if (this%nbound .eq. 0) return
    !
    ! -- Update shutoff count
    this%ishutoffcnt = this%ishutoffcnt + 1
    !
    ! -- Calculate hcof and rhs for each maw entry
    ibnd = 1
    do n = 1, this%nmawwells
      hmaw = this%xnewpak(n)
      do j = 1, this%ngwfnodes(n)
        jpos = this%get_jpos(n, j)
        this%hcof(ibnd) = dzero
        this%rhs(ibnd) = dzero
        !
        ! -- set bound, hcof, and rhs components
        !
        ! -- use connection method so the gwf-maw budget flows
        !    are consistent with the maw-gwf budget flows
        if (this%iboundpak(n) == 0 .or. this%iboundconn(jpos) == 0) then
          cmaw = dzero
          term = dzero
          cterm = dzero
        else
          call this%maw_calculate_conn_terms(n, j, icflow, cmaw, cterm, &
                                             term, flow)
        end if
        this%simcond(jpos) = cmaw
        this%bound(2, ibnd) = cmaw
        this%hcof(ibnd) = -term
        this%rhs(ibnd) = -term * hmaw + cterm
        !
        ! -- increment boundary number
        ibnd = ibnd + 1
      end do
    end do
  end subroutine maw_cfupdate
 
  !> @brief Set up the budget object that stores all the maw flows
  !! The terms listed here must correspond in number and order to the ones
  !! listed in the maw_fill_budobj routine.
  !<
  subroutine maw_setup_budobj(this)
    ! -- modules
    use constantsmodule, only: lenbudtxt
    ! -- dummy
    class(mawtype) :: this
    ! -- local
    integer(I4B) :: nbudterm
    integer(I4B) :: n, j, n2
    real(DP) :: q
    integer(I4B) :: maxlist, naux
    integer(I4B) :: idx
    character(len=LENBUDTXT) :: text
    character(len=LENBUDTXT), dimension(1) :: auxtxt
    !
    ! -- Determine the number of maw budget terms. These are fixed for
    !    the simulation and cannot change.
    ! gwf rate [flowing_well] storage constant_flow [frommvr tomvr tomvrcf [tomvrfw]] [aux]
    nbudterm = 4
    if (this%iflowingwells > 0) then
      nbudterm = nbudterm + 1
    end if
    if (this%imover == 1) then
      nbudterm = nbudterm + 3
      if (this%iflowingwells > 0) then
        nbudterm = nbudterm + 1
      end if
    end if
    if (this%naux > 0) nbudterm = nbudterm + 1
    !
    ! -- set up budobj
    call budgetobject_cr(this%budobj, this%packName)
    call this%budobj%budgetobject_df(this%nmawwells, nbudterm, 0, 0, &
                                     ibudcsv=this%ibudcsv)
    idx = 0
    !
    ! -- Go through and set up each budget term
    !
    ! --
    text = '             GWF'
    idx = idx + 1
    maxlist = this%maxbound
    naux = 1
    auxtxt(1) = '       FLOW-AREA'
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%name_model, &
                                             maxlist, .false., .true., &
                                             naux, auxtxt)
    call this%budobj%budterm(idx)%reset(this%maxbound)
    q = dzero
    do n = 1, this%nmawwells
      do j = 1, this%ngwfnodes(n)
        n2 = this%get_gwfnode(n, j)
        call this%budobj%budterm(idx)%update_term(n, n2, q)
      end do
    end do
    !
    ! --
    text = '            RATE'
    idx = idx + 1
    maxlist = this%nmawwells
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    if (this%iflowingwells > 0) then
      text = '         FW-RATE'
      idx = idx + 1
      maxlist = this%nmawwells
      naux = 0
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux)
    end if
    !
    ! --
    text = '         STORAGE'
    idx = idx + 1
    maxlist = this%nmawwells
    naux = 1
    auxtxt(1) = '          VOLUME'
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%name_model, &
                                             maxlist, .false., .true., &
                                             naux, auxtxt)
    !
    ! --
    text = '        CONSTANT'
    idx = idx + 1
    maxlist = this%nmawwells
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    if (this%imover == 1) then
      !
      ! --
      text = '        FROM-MVR'
      idx = idx + 1
      maxlist = this%nmawwells
      naux = 0
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux)
      !
      ! --
      text = '     RATE-TO-MVR'
      idx = idx + 1
      maxlist = this%nmawwells
      naux = 0
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux)
      !
      ! -- constant-head flow to mover
      text = ' CONSTANT-TO-MVR'
      idx = idx + 1
      maxlist = this%nmawwells
      naux = 0
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux)
      !
      ! -- flowing-well flow to mover
      if (this%iflowingwells > 0) then
        !
        ! --
        text = '  FW-RATE-TO-MVR'
        idx = idx + 1
        maxlist = this%nmawwells
        naux = 0
        call this%budobj%budterm(idx)%initialize(text, &
                                                 this%name_model, &
                                                 this%packName, &
                                                 this%name_model, &
                                                 this%packName, &
                                                 maxlist, .false., .false., &
                                                 naux)
      end if
    end if
    !
    ! -- auxiliary variable
    naux = this%naux
    if (naux > 0) then
      !
      ! --
      text = '       AUXILIARY'
      idx = idx + 1
      maxlist = this%maxbound
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux, this%auxname)
    end if
    !
    ! -- if maw flow for each reach are written to the listing file
    if (this%iprflow /= 0) then
      call this%budobj%flowtable_df(this%iout)
    end if
  end subroutine maw_setup_budobj
 
  !> @brief Copy flow terms into this%budobj
  !!
  !! terms include a combination of the following:
  !! gwf rate [flowing_well] [storage] constant_flow [frommvr tomvr tomvrcf [tomvrfw]] [aux]
  !<
  subroutine maw_fill_budobj(this)
    ! -- modules
    ! -- dummy
    class(mawtype) :: this
    ! -- local
    integer(I4B) :: naux
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: n2
    integer(I4B) :: jpos
    integer(I4B) :: idx
    integer(I4B) :: ibnd
    real(DP) :: q
    real(DP) :: tmaw
    real(DP) :: bmaw
    real(DP) :: sat
    real(DP) :: qfact
    real(DP) :: q2
    real(DP) :: b
    real(DP) :: v
    ! -- formats
    !
    ! -- initialize counter
    idx = 0
    !
    ! -- GWF (LEAKAGE) and connection surface area (aux)
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%maxbound)
    ibnd = 1
    do n = 1, this%nmawwells
      do j = 1, this%ngwfnodes(n)
        jpos = this%get_jpos(n, j)
        n2 = this%get_gwfnode(n, j)
        !
        ! -- a connection that is not simulated has no wetted area, so it
        !    exchanges no heat by conduction with the aquifer (MWE)
        if (this%iboundpak(n) == 0 .or. this%iboundconn(jpos) == 0) then
          this%qauxcbc(1) = dzero
        else
          tmaw = this%topscrn(jpos)
          bmaw = this%botscrn(jpos)
          call this%maw_calculate_saturation(n, j, n2, sat)
          this%qauxcbc(1) = dtwo * dpi * this%radius(n) * sat * (tmaw - bmaw)
        end if
        q = this%qleak(ibnd)
        call this%budobj%budterm(idx)%update_term(n, n2, q, this%qauxcbc)
        ibnd = ibnd + 1
      end do
    end do
    !
    ! -- RATE (WITHDRAWAL RATE)
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nmawwells)
    do n = 1, this%nmawwells
      q = this%ratesim(n)
      !
      ! -- adjust if well rate is an outflow
      if (this%imover == 1 .and. q < dzero) then
        qfact = done
        if (this%qout(n) < dzero) then
          qfact = q / this%qout(n)
        end if
        q = q + qfact * this%pakmvrobj%get_qtomvr(n)
      end if
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- FLOWING WELL
    if (this%iflowingwells > 0) then
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(this%nmawwells)
      do n = 1, this%nmawwells
        q = this%qfw(n)
        if (this%imover == 1) then
          qfact = done
          !
          ! -- adjust if well rate is an outflow
          if (this%qout(n) < dzero) then
            qfact = q / this%qout(n)
          end if
          q = q + qfact * this%pakmvrobj%get_qtomvr(n)
        end if
        call this%budobj%budterm(idx)%update_term(n, n, q)
      end do
    end if
    !
    ! -- STORAGE (AND VOLUME AS AUX)
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nmawwells)
    do n = 1, this%nmawwells
      b = this%xsto(n) - this%bot(n)
      if (b < dzero) then
        b = dzero
      end if
      v = this%area(n) * b
      if (this%imawissopt /= 1) then
        q = this%qsto(n)
      else
        q = dzero
      end if
      this%qauxcbc(1) = v
      call this%budobj%budterm(idx)%update_term(n, n, q, this%qauxcbc)
    end do
    !
    ! -- CONSTANT FLOW
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nmawwells)
    do n = 1, this%nmawwells
      q = this%qconst(n)
      !
      ! -- adjust if constant-flow rate is an outflow
      if (this%imover == 1 .and. q < dzero) then
        qfact = done
        if (this%qout(n) < dzero) then
          qfact = q / this%qout(n)
        end if
        q = q + qfact * this%pakmvrobj%get_qtomvr(n)
      end if
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- MOVER
    if (this%imover == 1) then
      !
      ! -- FROM MOVER
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(this%nmawwells)
      do n = 1, this%nmawwells
        if (this%iboundpak(n) == 0) then
          q = dzero
        else
          q = this%pakmvrobj%get_qfrommvr(n)
        end if
        call this%budobj%budterm(idx)%update_term(n, n, q)
      end do
      !
      ! -- RATE TO MOVER
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(this%nmawwells)
      do n = 1, this%nmawwells
        q = this%pakmvrobj%get_qtomvr(n)
        if (q > dzero) then
          q = -q
          q2 = this%ratesim(n)
          !
          ! -- adjust TO MOVER if well rate is outflow
          if (q2 < dzero) then
            qfact = q2 / this%qout(n)
            q = q * qfact
          else
            q = dzero
          end if
        end if
        call this%budobj%budterm(idx)%update_term(n, n, q)
      end do
      !
      ! -- CONSTANT TO MOVER
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(this%nmawwells)
      do n = 1, this%nmawwells
        q = this%pakmvrobj%get_qtomvr(n)
        if (q > dzero) then
          q = -q
          q2 = this%qconst(n)
          ! -- adjust TO MOVER if well rate is outflow
          if (q2 < dzero) then
            qfact = q2 / this%qout(n)
            q = q * qfact
          else
            q = dzero
          end if
        end if
        call this%budobj%budterm(idx)%update_term(n, n, q)
      end do
      !
      ! -- FLOWING WELL TO MOVER
      if (this%iflowingwells > 0) then
        idx = idx + 1
        call this%budobj%budterm(idx)%reset(this%nmawwells)
        do n = 1, this%nmawwells
          q = this%pakmvrobj%get_qtomvr(n)
          if (q > dzero) then
            q = -q
            q2 = this%ratesim(n)
            !
            ! -- adjust TO MOVER if well rate is outflow
            qfact = done
            if (this%qout(n) < dzero) then
              qfact = this%qfw(n) / this%qout(n)
            end if
            q = q * qfact
          end if
          call this%budobj%budterm(idx)%update_term(n, n, q)
        end do
      end if
 
    end if
    !
    ! -- AUXILIARY VARIABLES
    naux = this%naux
    if (naux > 0) then
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(this%nmawwells)
      do n = 1, this%nmawwells
        q = dzero
        call this%budobj%budterm(idx)%update_term(n, n, q, this%auxvar(:, n))
      end do
    end if
    !
    ! --Terms are filled, now accumulate them for this time step
    call this%budobj%accumulate_terms()
  end subroutine maw_fill_budobj
 
  !> @brief Set up the table object that is used to write the maw head data
  !!
  !! The terms listed here must correspond in number and order to the ones
  !! written to the head table in the maw_ot method.
  !<
  subroutine maw_setup_tableobj(this)
    ! -- modules
    use constantsmodule, only: linelength, lenbudtxt
    ! -- dummy
    class(mawtype) :: this
    ! -- local
    integer(I4B) :: nterms
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    !
    ! -- setup well head table
    if (this%iprhed > 0) then
      !
      ! -- Determine the number of head table columns
      nterms = 2
      if (this%inamedbound == 1) nterms = nterms + 1
      !
      ! -- set up table title
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') HEADS FOR EACH CONTROL VOLUME'
      !
      ! -- set up head tableobj
      call table_cr(this%headtab, this%packName, title)
      call this%headtab%table_df(this%nmawwells, nterms, this%iout, &
                                 transient=.true.)
      !
      ! -- Go through and set up table budget term
      if (this%inamedbound == 1) then
        text = 'NAME'
        call this%headtab%initialize_column(text, 20, alignment=tableft)
      end if
      !
      ! -- reach number
      text = 'NUMBER'
      call this%headtab%initialize_column(text, 10, alignment=tabcenter)
      !
      ! -- reach stage
      text = 'HEAD'
      call this%headtab%initialize_column(text, 12, alignment=tabcenter)
    end if
  end subroutine maw_setup_tableobj
 
  !> @brief Get position of value in connection data
  !<
  function get_jpos(this, n, j) result(jpos)
    ! -- return variable
    integer(I4B) :: jpos
    ! -- dummy
    class(mawtype) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: j
    ! -- local
    !
    ! -- set jpos
    jpos = this%iaconn(n) + j - 1
  end function get_jpos
 
  !> @brief Get the gwfnode for connection
  !<
  function get_gwfnode(this, n, j) result(igwfnode)
    ! -- return variable
    integer(I4B) :: igwfnode
    ! -- dummy
    class(mawtype) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: j
    ! -- local
    integer(I4B) :: jpos
    !
    ! -- set jpos
    jpos = this%get_jpos(n, j)
    igwfnode = this%gwfnodes(jpos)
  end function get_gwfnode
 
  !> @brief Activate density terms
  !<
  subroutine maw_activate_density(this)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    ! -- local
    integer(I4B) :: i, j
    ! -- formats
    !
    ! -- Set idense and reallocate denseterms to be of size MAXBOUND
    this%idense = 1
    call mem_reallocate(this%denseterms, 3, this%MAXBOUND, 'DENSETERMS', &
                        this%memoryPath)
    do i = 1, this%maxbound
      do j = 1, 3
        this%denseterms(j, i) = dzero
      end do
    end do
    write (this%iout, '(/1x,a)') 'DENSITY TERMS HAVE BEEN ACTIVATED FOR MAW &
      &PACKAGE: '//trim(adjustl(this%packName))
  end subroutine maw_activate_density
 
  !> @brief Activate viscosity terms
  !!
  !! Method to activate addition of viscosity terms for a MAW package reach.
  !<
  subroutine maw_activate_viscosity(this)
    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    ! -- dummy variables
    class(mawtype), intent(inout) :: this !< MawType object
    ! -- local variables
    integer(I4B) :: i
    integer(I4B) :: j
    !
    ! -- Set ivsc and reallocate viscratios to be of size MAXBOUND
    this%ivsc = 1
    call mem_reallocate(this%viscratios, 2, this%MAXBOUND, 'VISCRATIOS', &
                        this%memoryPath)
    do i = 1, this%maxbound
      do j = 1, 2
        this%viscratios(j, i) = done
      end do
    end do
    write (this%iout, '(/1x,a)') 'VISCOSITY HAS BEEN ACTIVATED FOR MAW &
      &PACKAGE: '//trim(adjustl(this%packName))
  end subroutine maw_activate_viscosity
 
  !> @brief Calculate the groundwater-maw density exchange terms
  !!
  !! Arguments are as follows:
  !!     iconn       : maw-gwf connection number
  !!     hmaw        : maw head
  !!     hgwf        : gwf head
  !!     cond        : conductance
  !!     bmaw        : bottom elevation of this connection
  !!     flow        : calculated flow, updated here with density terms, + into maw
  !!     hcofterm    : head coefficient term
  !!     rhsterm     : right-hand-side value, updated here with density terms
  !!
  !! Member variable used here
  !!     denseterms  : shape (3, MAXBOUND), filled by buoyancy package
  !!                     col 1 is relative density of maw (densemaw / denseref)
  !!                     col 2 is relative density of gwf cell (densegwf / denseref)
  !!                     col 3 is elevation of gwf cell
  !!
  !! Upon return, amat and rhs for maw row should be updated as:
  !!    amat(idiag) = amat(idiag) - hcofterm
  !!    rhs(n) = rhs(n) + rhsterm
  !<
  subroutine maw_calculate_density_exchange(this, iconn, hmaw, hgwf, cond, &
                                            bmaw, flow, hcofterm, rhsterm)
    ! -- dummy
    class(mawtype), intent(inout) :: this
    integer(I4B), intent(in) :: iconn
    real(DP), intent(in) :: hmaw
    real(DP), intent(in) :: hgwf
    real(DP), intent(in) :: cond
    real(DP), intent(in) :: bmaw
    real(DP), intent(inout) :: flow
    real(DP), intent(inout) :: hcofterm
    real(DP), intent(inout) :: rhsterm
    ! -- local
    real(DP) :: t
    real(DP) :: havg
    real(DP) :: rdensemaw
    real(DP) :: rdensegwf
    real(DP) :: rdenseavg
    real(DP) :: elevavg
    ! -- formats
    !
    ! -- assign relative density terms, return if zero which means not avail yet
    rdensemaw = this%denseterms(1, iconn)
    rdensegwf = this%denseterms(2, iconn)
    if (rdensegwf == dzero) return
    !
    ! -- update rhsterm with density contribution
    if (hmaw > bmaw .and. hgwf > bmaw) then
      !
      ! -- hmaw and hgwf both above bmaw
      rdenseavg = dhalf * (rdensemaw + rdensegwf)
      !
      ! -- update rhsterm with first density term
      t = cond * (rdenseavg - done) * (hgwf - hmaw)
      rhsterm = rhsterm + t
      flow = flow + t
      !
      ! -- update rhterm with second density term
      havg = dhalf * (hgwf + hmaw)
      elevavg = this%denseterms(3, iconn)
      t = cond * (havg - elevavg) * (rdensegwf - rdensemaw)
      rhsterm = rhsterm + t
      flow = flow + t
    else if (hmaw > bmaw) then
      !
      ! -- if only hmaw is above bmaw, then increase correction term by density
      t = (rdensemaw - done) * rhsterm
      rhsterm = rhsterm + t
      !
    else if (hgwf > bmaw) then
      !
      ! -- if only hgwf is above bmaw, then increase correction term by density
      t = (rdensegwf - done) * rhsterm
      rhsterm = rhsterm + t
      !
    else
      !
      ! -- Flow should be zero so do nothing
    end if
  end subroutine maw_calculate_density_exchange
 
end module mawmodule