gwf-lak.f90

Go to the documentation of this file.

module lakmodule
  !
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: linelength, lenboundname, lenmempath, &
                             iwetlake, maxadpit, dzero, dprec, dem30, dem9, &
                             dem6, dem5, dem4, dem2, dem1, dhalf, dp7, dp9, &
                             dp999, done, dtwo, dpi, dthree, deight, dten, &
                             dhundred, dep20, donethird, dtwothirds, &
                             dfivethirds, dgravity, dcd, namedboundflag, &
                             lenftype, lenpackagename, lenpakloc, dnodata, &
                             mnormal, tableft, tabcenter, tabright, &
                             tabstring, tabucstring, tabinteger, tabreal
  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 smoothingmodule, only: squadraticsaturation, sqsaturation, &
                             squadraticsaturationderivative, &
                             sqsaturationderivative
  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_unit, &
                       store_warning, store_error_filename, &
                       deprecation_warning
  use mathutilmodule, only: is_close
  use characterstringmodule, only: characterstringtype
  use openspecmodule, only: access, form
  use simvariablesmodule, only: errmsg, warnmsg
  use matrixbasemodule
  use gwflakinputmodule, only: gwflakparamfoundtype
  !
  implicit none
  !
  private
  public :: laktype
  public :: lak_create
  !
  character(len=LENFTYPE) :: ftype = 'LAK'
  character(len=LENPACKAGENAME) :: text = '             LAK'
  !
  type laktabtype
    real(dp), dimension(:), pointer, contiguous :: tabstage => null()
    real(dp), dimension(:), pointer, contiguous :: tabvolume => null()
    real(dp), dimension(:), pointer, contiguous :: tabsarea => null()
    real(dp), dimension(:), pointer, contiguous :: tabwarea => null()
  end type laktabtype
  !
  !> @brief input context runtime pointers for LAK
  !!
  !! Groups all input context PERIOD and OUTLETS array pointers used at
  !! runtime. All pointers are set once in lak_ar and remain valid for
  !! the simulation lifetime.
  !<
  type :: lakinputtype
    ! -- 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()
    type(characterstringtype), dimension(:), pointer, contiguous :: &
      auxname => null()
    real(dp), dimension(:), pointer, contiguous :: &
      auxval => null()
    ! -- outlets (static, TS-linked)
    integer(I4B), dimension(:), pointer, contiguous :: outletno => null()
    real(dp), dimension(:), pointer, contiguous :: outlets_invert => null()
    real(dp), dimension(:), pointer, contiguous :: outlets_width => null()
    real(dp), dimension(:), pointer, contiguous :: outlets_rough => null()
    real(dp), dimension(:), pointer, contiguous :: outlets_slope => null()
  contains
    procedure :: init => lak_input_init
    procedure :: destroy => lak_input_destroy
  end type lakinputtype
  !
  type, extends(bndexttype) :: laktype
    ! -- scalars
    ! -- characters
    character(len=16), dimension(:), pointer, contiguous :: clakbudget => null()
    character(len=16), dimension(:), pointer, contiguous :: cauxcbc => null()
    ! -- control variables
    ! -- integers
    integer(I4B), pointer :: iprhed => null()
    integer(I4B), pointer :: istageout => null()
    integer(I4B), pointer :: ibudgetout => null()
    integer(I4B), pointer :: ibudcsv => null()
    integer(I4B), pointer :: ipakcsv => null()
    character(len=:), allocatable :: pakcsvfile !< requested PACKAGE_CONVERGENCE file, opened in lak_ar once the IMPLICIT option is known
    integer(I4B), pointer :: cbcauxitems => null()
    integer(I4B), pointer :: nlakes => null()
    integer(I4B), pointer :: noutlets => null()
    integer(I4B), pointer :: ntables => null()
    real(dp), pointer :: convlength => null()
    real(dp), pointer :: convtime => null()
    real(dp), pointer :: outdmax => null()
    integer(I4B), pointer :: igwhcopt => null()
    integer(I4B), pointer :: iconvchk => null()
    integer(I4B), pointer :: maxlakit => null() !< maximum number of iterations in LAK solve
    real(dp), pointer :: surfdep => null()
    real(dp), pointer :: dmaxchg => null()
    real(dp), pointer :: delh => null()
    integer(I4B), pointer :: check_attr => null()
    ! -- implicit formulation: solve the lake stage as an unknown in the
    !    groundwater flow matrix instead of by the legacy substitution iteration
    integer(I4B), pointer :: iimplicit => null() !< flag: solve lake stage in the gwf matrix
    integer(I4B), pointer :: iforceleg => null() !< flag (dev): force every active lake onto the legacy solver
    integer(I4B), pointer :: iforceleglak => null() !< lake (dev) forced onto the legacy solver, 0 if none
    ! -- for budgets
    integer(I4B), pointer :: bditems => null()
    ! -- vectors
    ! -- lake data
    integer(I4B), dimension(:), pointer, contiguous :: nlakeconn => null()
    integer(I4B), dimension(:), pointer, contiguous :: idxlakeconn => null()
    integer(I4B), dimension(:), pointer, contiguous :: ntabrow => null()
    real(dp), dimension(:), pointer, contiguous :: strt => null()
    real(dp), dimension(:), pointer, contiguous :: laketop => null()
    real(dp), dimension(:), pointer, contiguous :: lakebot => null()
    real(dp), dimension(:), pointer, contiguous :: sareamax => null()
    character(len=LENBOUNDNAME), dimension(:), pointer, &
      contiguous :: lakename => null()
    character(len=8), dimension(:), pointer, contiguous :: status => null()
    real(dp), dimension(:), pointer, contiguous :: avail => null()
    real(dp), dimension(:), pointer, contiguous :: lkgwsink => null()
    real(dp), dimension(:), pointer, contiguous :: stage => null()
    real(dp), dimension(:), pointer, contiguous :: rainfall => null()
    real(dp), dimension(:), pointer, contiguous :: evaporation => null()
    real(dp), dimension(:), pointer, contiguous :: runoff => null()
    real(dp), dimension(:), pointer, contiguous :: inflow => null()
    real(dp), dimension(:), pointer, contiguous :: withdrawal => null()
    !
    ! -- table data
    integer(I4B), dimension(:), pointer, contiguous :: ialaktab => null()
    real(dp), dimension(:), pointer, contiguous :: tabstage => null()
    real(dp), dimension(:), pointer, contiguous :: tabvolume => null()
    real(dp), dimension(:), pointer, contiguous :: tabsarea => null()
    real(dp), dimension(:), pointer, contiguous :: tabwarea => null()
    !
    ! -- lake solution data
    integer(I4B), dimension(:), pointer, contiguous :: ncncvr => null()
    ! -- IMPLICIT legacy solve: solve flagged lakes by substitution,
    !    and a per-lake count of consecutive non-converging outer iterations used
    !    to decide when to switch a stalled IMPLICIT lake to the legacy solver
    integer(I4B), dimension(:), pointer, contiguous :: ilegacy => null()
    integer(I4B), dimension(:), pointer, contiguous :: nstuck => null()
    real(dp), dimension(:), pointer, contiguous :: surfin => null()
    real(dp), dimension(:), pointer, contiguous :: surfout => null()
    real(dp), dimension(:), pointer, contiguous :: surfout1 => null()
    real(dp), dimension(:), pointer, contiguous :: precip => null()
    real(dp), dimension(:), pointer, contiguous :: precip1 => null()
    real(dp), dimension(:), pointer, contiguous :: evap => null()
    real(dp), dimension(:), pointer, contiguous :: evap1 => null()
    real(dp), dimension(:), pointer, contiguous :: evapo => null()
    real(dp), dimension(:), pointer, contiguous :: withr => null()
    real(dp), dimension(:), pointer, contiguous :: withr1 => null()
    real(dp), dimension(:), pointer, contiguous :: flwin => null()
    real(dp), dimension(:), pointer, contiguous :: flwiter => null()
    real(dp), dimension(:), pointer, contiguous :: flwiter1 => null()
    real(dp), dimension(:), pointer, contiguous :: seep => null()
    real(dp), dimension(:), pointer, contiguous :: seep1 => null()
    real(dp), dimension(:), pointer, contiguous :: seep0 => null()
    real(dp), dimension(:), pointer, contiguous :: stageiter => null()
    real(dp), dimension(:), pointer, contiguous :: chterm => null()
    !
    ! -- lake convergence
    integer(I4B), dimension(:), pointer, contiguous :: iseepc => null()
    integer(I4B), dimension(:), pointer, contiguous :: idhc => null()
    real(dp), dimension(:), pointer, contiguous :: en1 => null()
    real(dp), dimension(:), pointer, contiguous :: en2 => null()
    real(dp), dimension(:), pointer, contiguous :: r1 => null()
    real(dp), dimension(:), pointer, contiguous :: r2 => null()
    real(dp), dimension(:), pointer, contiguous :: dh0 => null()
    real(dp), dimension(:), pointer, contiguous :: s0 => null()
    real(dp), dimension(:), pointer, contiguous :: qgwf0 => null()
    ! -- groundwater flow matrix bookkeeping for the implicit formulation
    integer(I4B), dimension(:), pointer, contiguous :: idxlocnode => null() !< local index of each lake row in x/rhs
    integer(I4B), dimension(:), pointer, contiguous :: idxdiag => null() !< position of lake-row diagonal (per lake)
    integer(I4B), dimension(:), pointer, contiguous :: idxoffdglo => null() !< position of lake-row -> gwf column (per connection)
    integer(I4B), dimension(:), pointer, contiguous :: idxsymdglo => null() !< position of gwf-row diagonal (per connection)
    integer(I4B), dimension(:), pointer, contiguous :: idxsymoffdglo => null() !< position of gwf-row -> lake column (per connection)
    !
    ! -- lake connection data
    integer(I4B), dimension(:), pointer, contiguous :: imap => null()
    integer(I4B), dimension(:), pointer, contiguous :: gwfconn => null()
    integer(I4B), dimension(:), pointer, contiguous :: nodesontop => null()
    integer(I4B), dimension(:), pointer, contiguous :: ictype => null()
    real(dp), dimension(:), pointer, contiguous :: bedleak => null()
    real(dp), dimension(:), pointer, contiguous :: belev => null()
    real(dp), dimension(:), pointer, contiguous :: telev => null()
    real(dp), dimension(:), pointer, contiguous :: connlength => null()
    real(dp), dimension(:), pointer, contiguous :: connwidth => null()
    real(dp), dimension(:), pointer, contiguous :: sarea => null()
    real(dp), dimension(:), pointer, contiguous :: warea => null()
    real(dp), dimension(:), pointer, contiguous :: satcond => null()
    real(dp), dimension(:), pointer, contiguous :: simcond => null()
    real(dp), dimension(:), pointer, contiguous :: simlakgw => null()
    !
    ! -- lake outlet data
    integer(I4B), dimension(:), pointer, contiguous :: lakein => null()
    integer(I4B), dimension(:), pointer, contiguous :: lakeout => null()
    integer(I4B), dimension(:), pointer, contiguous :: iouttype => null()
    real(dp), dimension(:), pointer, contiguous :: outrate => null()
    real(dp), dimension(:), pointer, contiguous :: outinvert => null()
    real(dp), dimension(:), pointer, contiguous :: outwidth => null()
    real(dp), dimension(:), pointer, contiguous :: outrough => null()
    real(dp), dimension(:), pointer, contiguous :: outslope => null()
    real(dp), dimension(:), pointer, contiguous :: simoutrate => null()
    ! -- true once a PERIOD-block setting has overridden the OUTLETS
    !    block's own (separately time-series-linked) baseline geometry for
    !    that outlet; prevents lak_ad's OUTLETS resync from clobbering it
    logical(LGP), dimension(:), pointer, contiguous :: outlet_invert_set => null()
    logical(LGP), dimension(:), pointer, contiguous :: outlet_width_set => null()
    logical(LGP), dimension(:), pointer, contiguous :: outlet_rough_set => null()
    logical(LGP), dimension(:), pointer, contiguous :: outlet_slope_set => null()
    !
    ! -- lake output data
    real(dp), dimension(:), pointer, contiguous :: qauxcbc => null()
    real(dp), dimension(:), pointer, contiguous :: dbuff => null()
    real(dp), dimension(:), pointer, contiguous :: qleak => null()
    ! -- connected-cell head at the previous outer iteration (per connection),
    !    used by the IMPLICIT legacy-switch detection to compare a lake's stage change
    !    against the change in its connected aquifer heads
    real(dp), dimension(:), pointer, contiguous :: holdconn => null()
    real(dp), dimension(:), pointer, contiguous :: qsto => null()
    !
    ! -- pointer to gwf iss and gwf hk
    integer(I4B), pointer :: gwfiss => null()
    real(dp), dimension(:), pointer, contiguous :: gwfk11 => null()
    real(dp), dimension(:), pointer, contiguous :: gwfk33 => null()
    real(dp), dimension(:), pointer, contiguous :: gwfsat => null()
    integer(I4B), pointer :: gwfik33 => null()
    !
    ! -- package x, xold, and ibound
    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
    !
    ! -- lake budget object
    type(budgetobjecttype), pointer :: budobj => null()
    !
    ! -- input context runtime pointers (PERIOD + OUTLETS arrays)
    type(lakinputtype) :: input
    !
    ! -- lake table objects
    type(tabletype), pointer :: stagetab => null()
    type(tabletype), pointer :: pakcsvtab => null()
    !
    ! -- 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: lak vsc ratio; 2: gwf vsc ratio)
    !
    ! -- type bound procedures
 
  contains
 
    procedure :: lak_allocate_scalars
    procedure :: lak_allocate_arrays
    procedure :: source_options => lak_source_options
    procedure :: source_dimensions => lak_source_dimensions
    procedure :: read_initial_attr => lak_read_initial_attr
    procedure :: set_pointers => lak_set_pointers
    procedure :: bnd_ar => lak_ar
    procedure :: bnd_ac => lak_ac
    procedure :: bnd_mc => lak_mc
    procedure :: bnd_rp => lak_rp
    procedure :: bnd_ad => lak_ad
    procedure :: bnd_cf => lak_cf
    procedure :: bnd_fc => lak_fc
    procedure :: bnd_fn => lak_fn
    procedure :: bnd_nur => lak_nur
    procedure :: bnd_cc => lak_cc
    procedure, private :: lak_set_legacy
    procedure, private :: lak_check_disconnected
    procedure :: bnd_cq => lak_cq
    procedure :: bnd_ot_model_flows => lak_ot_model_flows
    procedure :: bnd_ot_package_flows => lak_ot_package_flows
    procedure :: bnd_ot_dv => lak_ot_dv
    procedure :: bnd_ot_bdsummary => lak_ot_bdsummary
    procedure :: bnd_da => lak_da
    procedure :: define_listlabel
    ! -- methods for observations
    procedure, public :: bnd_obs_supported => lak_obs_supported
    procedure, public :: bnd_df_obs => lak_df_obs
    procedure, public :: bnd_rp_obs => lak_rp_obs
    procedure, public :: bnd_bd_obs => lak_bd_obs
    ! -- private procedures
    procedure, private :: lak_source_packagedata
    procedure, private :: lak_source_connectiondata
    procedure, private :: lak_source_outlets
    procedure, private :: lak_source_tables
    procedure, private :: lak_source_table
    procedure, private :: lak_check_valid
    procedure, private :: lak_set_nonneg_value
    procedure, private :: lak_set_attribute_error
    procedure, private :: lak_bound_update
    procedure, private :: lak_calculate_sarea
    procedure, private :: lak_calculate_warea
    procedure, private :: lak_calculate_conn_warea
    procedure, public :: lak_calculate_vol
    procedure, private :: lak_calculate_conductance
    procedure, private :: lak_calculate_cond_head
    procedure, private :: lak_calculate_conn_conductance
    procedure, private :: lak_calculate_exchange
    procedure, private :: lak_calculate_conn_exchange
    procedure, private :: lak_calculate_conn_exchange_deriv
    procedure, private :: lak_estimate_conn_exchange
    procedure, private :: lak_calculate_storagechange
    procedure, private :: lak_calculate_rainfall
    procedure, private :: lak_calculate_runoff
    procedure, private :: lak_calculate_inflow
    procedure, private :: lak_calculate_external
    procedure, private :: lak_calculate_withdrawal
    procedure, private :: lak_calculate_evaporation
    procedure, private :: lak_calculate_outlet_inflow
    procedure, private :: lak_calculate_outlet_outflow
    procedure, private :: lak_get_internal_inlet
    procedure, private :: lak_get_internal_outlet
    procedure, private :: lak_get_external_outlet
    procedure, private :: lak_get_internal_mover
    procedure, private :: lak_get_external_mover
    procedure, private :: lak_get_outlet_tomover
    procedure, private :: lak_accumulate_chterm
    procedure, private :: lak_vol2stage
    procedure, private :: lak_solve
    procedure, private :: lak_solve_single
    procedure, private :: lak_estimate_seepage_single
    procedure, private :: lak_fc_implicit
    procedure, private :: lak_budget_nogwf
    procedure, private :: lak_outlet_outflow_rate
    procedure, private :: lak_bisection
    procedure, private :: lak_calculate_available
    procedure, private :: lak_calculate_residual
    procedure, private :: lak_linear_interpolation
    procedure, private :: lak_setup_budobj
    procedure, private :: lak_fill_budobj
    procedure, private :: laktables_to_vectors
    ! -- table
    procedure, private :: lak_setup_tableobj
    ! -- density
    procedure :: lak_activate_density
    procedure, private :: lak_calculate_density_exchange
    ! -- viscosity
    procedure :: lak_activate_viscosity
  end type laktype
 
  ! -- implicit-formulation procedures, implemented in the gwf-lak-implicit
  !    submodule (submodules/gwf-lak-implicit.f90)
  interface
    module subroutine lak_budget_nogwf(this, n, stage, b)
      class(laktype), intent(inout) :: this
      integer(I4B), intent(in) :: n
      real(dp), intent(in) :: stage
      real(dp), intent(inout) :: b
    end subroutine
  end interface
 
  interface
    module subroutine lak_fc_implicit(this, rhs, matrix_sln)
      class(laktype) :: this
      real(dp), dimension(:), intent(inout) :: rhs
      class(matrixbasetype), pointer :: matrix_sln
    end subroutine
  end interface
 
  interface
    module subroutine lak_set_legacy(this, kiter, icnvgmod)
      class(laktype), intent(inout) :: this
      integer(I4B), intent(in) :: kiter !< outer (Picard) iteration number
      integer(I4B), intent(in) :: icnvgmod !< 0 if the model has not converged
    end subroutine
  end interface
 
  interface
    module subroutine lak_check_disconnected(this)
      class(laktype), intent(inout) :: this
    end subroutine
  end interface
 
contains
 
  !> @brief Create a new LAK Package and point bndobj to the new package
  !<
  subroutine lak_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
                        mempath)
    ! -- modules
    use memoryhelpermodule, only: create_mem_path
    ! -- 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
    ! -- local
    type(laktype), pointer :: lakobj
    !
    ! -- allocate the object and assign values to object variables
    allocate (lakobj)
    packobj => lakobj
    !
    ! -- create name and memory path
    call packobj%set_names(ibcnum, namemodel, pakname, ftype, mempath)
    packobj%text = text
    !
    ! -- allocate scalars
    call lakobj%lak_allocate_scalars()
    !
    ! -- initialize package
    call packobj%pack_initialize()
    !
    packobj%inunit = inunit
    packobj%iout = iout
    packobj%id = id
    packobj%ibcnum = ibcnum
    packobj%ncolbnd = 3
    packobj%iscloc = 0 ! not supported
    packobj%isadvpak = 1
    packobj%ictMemPath = create_mem_path(namemodel, 'NPF')
  end subroutine lak_create
 
  !> @brief Allocate scalar members
  !<
  subroutine lak_allocate_scalars(this)
    ! -- dummy
    class(laktype), 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%iprhed, 'IPRHED', this%memoryPath)
    call mem_allocate(this%istageout, 'ISTAGEOUT', this%memoryPath)
    call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath)
    call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
    call mem_allocate(this%ipakcsv, 'IPAKCSV', this%memoryPath)
    call mem_allocate(this%nlakes, 'NLAKES', this%memoryPath)
    call mem_allocate(this%noutlets, 'NOUTLETS', this%memoryPath)
    call mem_allocate(this%ntables, 'NTABLES', this%memoryPath)
    call mem_allocate(this%convlength, 'CONVLENGTH', this%memoryPath)
    call mem_allocate(this%convtime, 'CONVTIME', this%memoryPath)
    call mem_allocate(this%outdmax, 'OUTDMAX', this%memoryPath)
    call mem_allocate(this%igwhcopt, 'IGWHCOPT', this%memoryPath)
    call mem_allocate(this%iconvchk, 'ICONVCHK', this%memoryPath)
    call mem_allocate(this%maxlakit, 'MAXLAKIT', this%memoryPath)
    call mem_allocate(this%surfdep, 'SURFDEP', this%memoryPath)
    call mem_allocate(this%dmaxchg, 'DMAXCHG', this%memoryPath)
    call mem_allocate(this%delh, 'DELH', this%memoryPath)
    call mem_allocate(this%check_attr, 'CHECK_ATTR', this%memoryPath)
    call mem_allocate(this%iimplicit, 'IIMPLICIT', this%memoryPath)
    call mem_allocate(this%iforceleg, 'IFORCELEG', this%memoryPath)
    call mem_allocate(this%iforceleglak, 'IFORCELEGLAK', this%memoryPath)
    call mem_allocate(this%bditems, 'BDITEMS', this%memoryPath)
    call mem_allocate(this%cbcauxitems, 'CBCAUXITEMS', this%memoryPath)
    call mem_allocate(this%idense, 'IDENSE', this%memoryPath)
    !
    ! -- Set values
    this%iprhed = 0
    this%istageout = 0
    this%ibudgetout = 0
    this%ibudcsv = 0
    this%ipakcsv = 0
    this%nlakes = 0
    this%noutlets = 0
    this%ntables = 0
    this%convlength = done
    this%convtime = done
    this%outdmax = dzero
    this%igwhcopt = 0
    this%iconvchk = 1
    this%maxlakit = maxadpit
    this%surfdep = dzero
    this%dmaxchg = dem5
    this%delh = dp999 * this%dmaxchg
    this%iimplicit = 0
    this%iforceleg = 0
    this%iforceleglak = 0
    this%bditems = 11
    this%cbcauxitems = 1
    this%idense = 0
    this%ivsc = 0
  end subroutine lak_allocate_scalars
 
  !> @brief Allocate scalar members
  !<
  subroutine lak_allocate_arrays(this)
    ! -- modules
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    integer(I4B) :: i
    !
    ! -- call standard BndExtType allocate arrays
    call this%BndExtType%allocate_arrays()
    !
    ! -- allocate featureauxvar and set pkg_ifno/pkg_auxvar input context pointers
    call this%BndExtType%allocate_featureauxvar()
    !
    ! -- allocate character array for budget text
    allocate (this%clakbudget(this%bditems))
    !
    !-- fill clakbudget
    this%clakbudget(1) = '             GWF'
    this%clakbudget(2) = '        RAINFALL'
    this%clakbudget(3) = '     EVAPORATION'
    this%clakbudget(4) = '          RUNOFF'
    this%clakbudget(5) = '      EXT-INFLOW'
    this%clakbudget(6) = '      WITHDRAWAL'
    this%clakbudget(7) = '     EXT-OUTFLOW'
    this%clakbudget(8) = '         STORAGE'
    this%clakbudget(9) = '        CONSTANT'
    this%clakbudget(10) = '        FROM-MVR'
    this%clakbudget(11) = '          TO-MVR'
    !
    ! -- allocate and initialize dbuff
    if (this%istageout > 0) then
      call mem_allocate(this%dbuff, this%nlakes, 'DBUFF', this%memoryPath)
      do i = 1, this%nlakes
        this%dbuff(i) = dzero
      end do
    else
      call mem_allocate(this%dbuff, 0, 'DBUFF', this%memoryPath)
    end if
    !
    ! -- allocate character array for budget text
    allocate (this%cauxcbc(this%cbcauxitems))
    !
    ! -- allocate and initialize qauxcbc
    call mem_allocate(this%qauxcbc, this%cbcauxitems, 'QAUXCBC', this%memoryPath)
    do i = 1, this%cbcauxitems
      this%qauxcbc(i) = dzero
    end do
    !
    ! -- allocate qleak and qsto
    call mem_allocate(this%qleak, this%maxbound, 'QLEAK', this%memoryPath)
    do i = 1, this%maxbound
      this%qleak(i) = dzero
    end do
    ! -- holdconn is only used by the implicit legacy-switch detector; allocate it at
    !    size 0 otherwise so legacy LAK runs do not pay the maxbound memory cost
    if (this%iimplicit /= 0) then
      call mem_allocate(this%holdconn, this%maxbound, 'HOLDCONN', this%memoryPath)
      do i = 1, this%maxbound
        this%holdconn(i) = dzero
      end do
    else
      call mem_allocate(this%holdconn, 0, 'HOLDCONN', this%memoryPath)
    end if
    call mem_allocate(this%qsto, this%nlakes, 'QSTO', this%memoryPath)
    do i = 1, this%nlakes
      this%qsto(i) = 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 lak_allocate_arrays
 
  !> @brief Source PACKAGEDATA block from input context
  !<
  subroutine lak_source_packagedata(this)
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: ilak, nconn
    character(len=9) :: cno
    character(len=LENBOUNDNAME) :: bndName
    integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
    real(DP), dimension(:), pointer, contiguous :: strt => null()
    integer(I4B), dimension(:), pointer, contiguous :: nlakeconn => null()
    type(characterstringtype), dimension(:), pointer, contiguous :: &
      boundname_cst => null()
    integer(I4B), allocatable :: nboundchk(:)
    !
    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
      ' PACKAGEDATA'
    !
    ! -- allocate lake arrays
    call mem_allocate(this%nlakeconn, this%nlakes, 'NLAKECONN', this%memoryPath)
    call mem_allocate(this%idxlakeconn, this%nlakes + 1, 'IDXLAKECONN', &
                      this%memoryPath)
    call mem_allocate(this%ntabrow, this%nlakes, 'NTABROW', this%memoryPath)
    call mem_allocate(this%strt, this%nlakes, 'STRT', this%memoryPath)
    call mem_allocate(this%laketop, this%nlakes, 'LAKETOP', this%memoryPath)
    call mem_allocate(this%lakebot, this%nlakes, 'LAKEBOT', this%memoryPath)
    call mem_allocate(this%sareamax, this%nlakes, 'SAREAMAX', this%memoryPath)
    ! -- alias into the input context's permanent, feature-indexed arrays
    ! (allocated and DZERO-initialized by the loader)
    call mem_setptr(this%stage, 'STAGE', this%input_mempath)
    call mem_setptr(this%rainfall, 'RAINFALL', this%input_mempath)
    call mem_setptr(this%evaporation, 'EVAPORATION', this%input_mempath)
    call mem_setptr(this%runoff, 'RUNOFF', this%input_mempath)
    call mem_setptr(this%inflow, 'INFLOW', this%input_mempath)
    call mem_setptr(this%withdrawal, 'WITHDRAWAL', this%input_mempath)
    call mem_allocate(this%avail, this%nlakes, 'AVAIL', this%memoryPath)
    call mem_allocate(this%lkgwsink, this%nlakes, 'LKGWSINK', this%memoryPath)
    call mem_allocate(this%ncncvr, this%nlakes, 'NCNCVR', this%memoryPath)
    call mem_allocate(this%ilegacy, this%nlakes, 'ILEGACY', this%memoryPath)
    call mem_allocate(this%nstuck, this%nlakes, 'NSTUCK', this%memoryPath)
    call mem_allocate(this%surfin, this%nlakes, 'SURFIN', this%memoryPath)
    call mem_allocate(this%surfout, this%nlakes, 'SURFOUT', this%memoryPath)
    call mem_allocate(this%surfout1, this%nlakes, 'SURFOUT1', this%memoryPath)
    call mem_allocate(this%precip, this%nlakes, 'PRECIP', this%memoryPath)
    call mem_allocate(this%precip1, this%nlakes, 'PRECIP1', this%memoryPath)
    call mem_allocate(this%evap, this%nlakes, 'EVAP', this%memoryPath)
    call mem_allocate(this%evap1, this%nlakes, 'EVAP1', this%memoryPath)
    call mem_allocate(this%evapo, this%nlakes, 'EVAPO', this%memoryPath)
    call mem_allocate(this%withr, this%nlakes, 'WITHR', this%memoryPath)
    call mem_allocate(this%withr1, this%nlakes, 'WITHR1', this%memoryPath)
    call mem_allocate(this%flwin, this%nlakes, 'FLWIN', this%memoryPath)
    call mem_allocate(this%flwiter, this%nlakes, 'FLWITER', this%memoryPath)
    call mem_allocate(this%flwiter1, this%nlakes, 'FLWITER1', this%memoryPath)
    call mem_allocate(this%seep, this%nlakes, 'SEEP', this%memoryPath)
    call mem_allocate(this%seep1, this%nlakes, 'SEEP1', this%memoryPath)
    call mem_allocate(this%seep0, this%nlakes, 'SEEP0', this%memoryPath)
    call mem_allocate(this%stageiter, this%nlakes, 'STAGEITER', this%memoryPath)
    call mem_allocate(this%chterm, this%nlakes, 'CHTERM', this%memoryPath)
    !
    ! -- lake boundary and stages. In the implicit formulation iboundpak
    !    and xnewpak are aliased to slices of the global ibound/x vectors in
    !    lak_set_pointers (called before this routine), so they are not
    !    allocated here.
    if (this%iimplicit == 0) then
      call mem_allocate(this%iboundpak, this%nlakes, 'IBOUND', this%memoryPath)
      call mem_allocate(this%xnewpak, this%nlakes, 'XNEWPAK', this%memoryPath)
    end if
    call mem_allocate(this%xoldpak, this%nlakes, 'XOLDPAK', this%memoryPath)
    !
    ! -- lake iteration variables
    call mem_allocate(this%iseepc, this%nlakes, 'ISEEPC', this%memoryPath)
    call mem_allocate(this%idhc, this%nlakes, 'IDHC', this%memoryPath)
    call mem_allocate(this%en1, this%nlakes, 'EN1', this%memoryPath)
    call mem_allocate(this%en2, this%nlakes, 'EN2', this%memoryPath)
    call mem_allocate(this%r1, this%nlakes, 'R1', this%memoryPath)
    call mem_allocate(this%r2, this%nlakes, 'R2', this%memoryPath)
    call mem_allocate(this%dh0, this%nlakes, 'DH0', this%memoryPath)
    call mem_allocate(this%s0, this%nlakes, 'S0', this%memoryPath)
    call mem_allocate(this%qgwf0, this%nlakes, 'QGWF0', this%memoryPath)
    !
    ! -- allocate character storage
    allocate (this%lakename(this%nlakes))
    allocate (this%status(this%nlakes))
    !
    ! -- initialize arrays
    do n = 1, this%nlakes
      this%ntabrow(n) = 0
      this%status(n) = 'ACTIVE'
      this%laketop(n) = -dep20
      this%lakebot(n) = dep20
      this%sareamax(n) = dzero
      ! -- in the implicit formulation iboundpak/xnewpak are not yet associated
      !    (they are aliased to the global arrays later in lak_set_pointers and
      !    initialized in lak_read_initial_attr); only touch them here for the
      !    legacy formulation
      if (this%iimplicit == 0) then
        this%iboundpak(n) = 1
        this%xnewpak(n) = dep20
      end if
      this%xoldpak(n) = dep20
      this%rainfall(n) = dzero
      this%evaporation(n) = dzero
      this%runoff(n) = dzero
      this%inflow(n) = dzero
      this%withdrawal(n) = dzero
      this%ilegacy(n) = 0
      this%nstuck(n) = 0
    end do
    !
    ! -- get input context PACKAGEDATA column arrays
    call mem_setptr(ifno, 'PACKAGEDATA_IFNO', this%input_mempath)
    call mem_setptr(strt, 'STRT', this%input_mempath)
    call mem_setptr(nlakeconn, 'NLAKECONN', this%input_mempath)
    if (this%inamedbound /= 0) then
      call mem_setptr(boundname_cst, 'BOUNDNAME', this%input_mempath)
    end if
    !
    ! -- allocate check array
    allocate (nboundchk(this%nlakes))
    nboundchk = 0
    ilak = 0
    nconn = 0
    !
    do n = 1, this%nlakes
      !
      ! -- lake number from input
      ilak = this%validate_ifno(ifno(n), this%nlakes, nboundchk, &
                                'lakeno', 'PACKAGEDATA')
      if (ilak == 0) cycle
      !
      ! -- strt
      this%strt(ilak) = strt(n)
      !
      ! -- nlakeconn
      if (nlakeconn(n) < 0) then
        write (errmsg, '(a,1x,i0)') 'nlakeconn MUST BE >= 0 for lake ', ilak
        call store_error(errmsg)
      end if
      !
      ! -- the IMPLICIT formulation solves the lake stage as a matrix unknown,
      !    which requires at least one groundwater connection to give the lake
      !    row a non-zero diagonal; a lake with no connections is not supported
      if (this%iimplicit /= 0 .and. nlakeconn(n) == 0) then
        write (errmsg, '(a,1x,i0,1x,a)') &
          'lake', ilak, 'has no connections; the IMPLICIT option requires &
          &each lake to have at least one GWF connection.'
        call store_error(errmsg)
      end if
      !
      nconn = nconn + nlakeconn(n)
      this%nlakeconn(ilak) = nlakeconn(n)
      !
      ! -- set default bndName
      write (cno, '(i9.9)') ilak
      bndname = 'Lake'//cno
      !
      ! -- lakename from boundname if provided
      if (this%inamedbound /= 0) then
        bndname = boundname_cst(n)
        if (trim(bndname) == '') then
          write (cno, '(i9.9)') ilak
          bndname = 'Lake'//cno
        end if
      end if
      this%lakename(ilak) = bndname
    end do
    !
    ! -- check for duplicate or missing lakes
    call this%report_ifno_coverage(nboundchk, this%nlakes, 'lake', &
                                   'PACKAGEDATA')
    !
    write (this%iout, '(1x,a)') 'END OF '//trim(adjustl(this%text))// &
      ' PACKAGEDATA'
    !
    ! -- terminate if any errors were detected
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- set MAXBOUND
    this%MAXBOUND = nconn
    write (this%iout, '(//4x,a,i7)') 'MAXBOUND = ', this%maxbound
    !
    ! -- set idxlakeconn
    this%idxlakeconn(1) = 1
    do n = 1, this%nlakes
      this%idxlakeconn(n + 1) = this%idxlakeconn(n) + this%nlakeconn(n)
    end do
    !
    ! -- deallocate check array
    deallocate (nboundchk)
    !
    ! -- release input context packagedata memory; PACKAGEDATA_IFNO excluded
    !    because allocate_featureauxvar retains a pointer to it
    call memorystore_release('STRT', this%input_mempath)
    call memorystore_release('NLAKECONN', this%input_mempath)
    call memorystore_release('BOUNDNAME', this%input_mempath)
  end subroutine lak_source_packagedata
 
  !> @brief Source CONNECTIONDATA block from input context
  !<
  subroutine lak_source_connectiondata(this)
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    integer(I4B) :: n, j, ipos, nn
    integer(I4B) :: ilak, iconn
    integer(I4B) :: ipos0, icellid0, icellid, nn2
    logical(LGP) :: is_lake_bed
    character(len=LINELENGTH) :: cellid_str
    character(len=LINELENGTH) :: ctypenm
    character(len=LINELENGTH) :: cbleak
    real(DP) :: top, bot, rval
    integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
    integer(I4B), dimension(:), pointer, contiguous :: iconns => null()
    integer(I4B), dimension(:, :), pointer, contiguous :: cellid => null()
    type(characterstringtype), dimension(:), pointer, contiguous :: &
      claktype => null()
    type(characterstringtype), dimension(:), pointer, contiguous :: &
      bedleak => null()
    real(DP), dimension(:), pointer, contiguous :: belev => null()
    real(DP), dimension(:), pointer, contiguous :: telev => null()
    real(DP), dimension(:), pointer, contiguous :: connlen => null()
    real(DP), dimension(:), pointer, contiguous :: connwidth => null()
    integer(I4B), allocatable :: nboundchk(:)
    !
    ! -- allocate connection data arrays
    call mem_allocate(this%imap, this%MAXBOUND, 'IMAP', this%memoryPath)
    call mem_allocate(this%gwfconn, this%MAXBOUND, 'GWFCONN', this%memoryPath)
    call mem_allocate(this%nodesontop, this%MAXBOUND, 'NODESONTOP', &
                      this%memoryPath)
    call mem_allocate(this%ictype, this%MAXBOUND, 'ICTYPE', this%memoryPath)
    call mem_allocate(this%bedleak, this%MAXBOUND, 'BEDLEAK', this%memoryPath)
    call mem_allocate(this%belev, this%MAXBOUND, 'BELEV', this%memoryPath)
    call mem_allocate(this%telev, this%MAXBOUND, 'TELEV', this%memoryPath)
    call mem_allocate(this%connlength, this%MAXBOUND, 'CONNLENGTH', &
                      this%memoryPath)
    call mem_allocate(this%connwidth, this%MAXBOUND, 'CONNWIDTH', &
                      this%memoryPath)
    call mem_allocate(this%sarea, this%MAXBOUND, 'SAREA', this%memoryPath)
    call mem_allocate(this%warea, this%MAXBOUND, 'WAREA', 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%simlakgw, this%MAXBOUND, 'SIMLAKGW', this%memoryPath)
    !
    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
      ' CONNECTIONDATA'
    !
    ! -- get input context CONNECTIONDATA column arrays
    call mem_setptr(ifno, 'CONNDATA_IFNO', this%input_mempath)
    call mem_setptr(iconns, 'ICONN', this%input_mempath)
    call mem_setptr(cellid, 'CELLID', this%input_mempath)
    call mem_setptr(claktype, 'CLAKTYPE', this%input_mempath)
    call mem_setptr(bedleak, 'BEDLEAK', this%input_mempath)
    call mem_setptr(belev, 'BELEV', this%input_mempath)
    call mem_setptr(telev, 'TELEV', this%input_mempath)
    call mem_setptr(connlen, 'CONNLEN', this%input_mempath)
    call mem_setptr(connwidth, 'CONNWIDTH', this%input_mempath)
    !
    ! -- allocate check array
    allocate (nboundchk(this%MAXBOUND))
    nboundchk = 0
    !
    do n = 1, this%MAXBOUND
      !
      ilak = ifno(n)
      iconn = iconns(n)
      !
      if (ilak < 1 .or. ilak > this%nlakes) then
        write (errmsg, '(a,1x,i0)') 'lakeno MUST BE > 0 and <= ', this%nlakes
        call store_error(errmsg)
        cycle
      end if
      if (iconn < 1 .or. iconn > this%nlakeconn(ilak)) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
          'iconn FOR LAKE ', ilak, 'MUST BE > 0 and <= ', this%nlakeconn(ilak)
        call store_error(errmsg)
        cycle
      end if
      !
      ipos = this%idxlakeconn(ilak) + iconn - 1
      this%imap(ipos) = ilak
      nboundchk(ipos) = nboundchk(ipos) + 1
      !
      ! -- convert cellid integer array to string and resolve to a 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)
      if (nn < 1) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
          'INVALID cellid FOR LAKE ', ilak, 'connection', iconn
        call store_error(errmsg)
      end if
      this%gwfconn(ipos) = nn
      this%nodesontop(ipos) = nn
      !
      ! -- connection type
      ctypenm = claktype(n)
      select case (trim(ctypenm))
      case ('VERTICAL')
        this%ictype(ipos) = 0
      case ('HORIZONTAL')
        this%ictype(ipos) = 1
      case ('EMBEDDEDH')
        this%ictype(ipos) = 2
      case ('EMBEDDEDV')
        this%ictype(ipos) = 3
      case default
        write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,a,a)') &
          'UNKNOWN ctype FOR LAKE ', ilak, 'connection', iconn, &
          '(', trim(ctypenm), ')'
        call store_error(errmsg)
      end select
      !
      ! -- bed leakance
      cbleak = bedleak(n)
      select case (trim(cbleak))
      case ('NONE')
        is_lake_bed = .false.
        rval = dnodata
        write (warnmsg, '(2(a,1x,i0,1x),a,1pe8.1,a)') &
          'BEDLEAK for connection', iconn, 'in lake', ilak, &
          'is specified to be NONE. Lake connections where the lake-GWF '// &
          'connection conductance is solely a function of aquifer '// &
          'properties in the connected GWF cell should be specified '// &
          'with a DNODATA (', dnodata, ') value.'
        call deprecation_warning('CONNECTIONDATA', 'bedleak=NONE', '6.4.3', &
                                 warnmsg, this%inunit)
      case default
        read (cbleak, *) rval
        if (is_close(rval, dnodata)) then
          is_lake_bed = .false.
        else
          is_lake_bed = .true.
        end if
      end select
      this%bedleak(ipos) = rval
      if (is_lake_bed .and. this%bedleak(ipos) < dzero) then
        write (errmsg, '(a,1x,i0,1x,a)') 'bedleak FOR LAKE ', ilak, &
          'MUST BE >= 0'
        call store_error(errmsg)
      end if
      !
      ! -- belev
      this%belev(ipos) = belev(n)
      !
      ! -- telev
      this%telev(ipos) = telev(n)
      !
      ! -- connection length
      rval = connlen(n)
      if (rval <= dzero) then
        if (this%ictype(ipos) == 1 .or. this%ictype(ipos) == 2 .or. &
            this%ictype(ipos) == 3) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,a,1x,a)') &
            'connection length (connlen) FOR LAKE ', ilak, &
            ', CONNECTION NO.', iconn, ', MUST BE > 0 FOR SPECIFIED ', &
            'connection type (ctype)', trim(ctypenm)
          call store_error(errmsg)
        else
          rval = dzero
        end if
      end if
      this%connlength(ipos) = rval
      !
      ! -- connection width
      rval = connwidth(n)
      if (rval < dzero) then
        if (this%ictype(ipos) == 1) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
            'cell width (connwidth) FOR LAKE ', ilak, &
            ' HORIZONTAL CONNECTION ', iconn, 'MUST BE >= 0'
          call store_error(errmsg)
        else
          rval = dzero
        end if
      end if
      this%connwidth(ipos) = rval
    end do
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' CONNECTIONDATA'
    !
    ! -- terminate if any errors were detected
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- check that embedded lakes have only one connection
    do n = 1, this%nlakes
      j = 0
      do ipos = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        if (this%ictype(ipos) /= 2 .and. this%ictype(ipos) /= 3) cycle
        j = j + 1
        if (j > 1) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
            'nlakeconn FOR LAKE', n, 'EMBEDDED CONNECTION', j, ' EXCEEDS 1.'
          call store_error(errmsg)
        end if
      end do
    end do
    ! -- check that an embedded lake is not in the same cell as a lake
    !   with a vertical connection
    do n = 1, this%nlakes
      ipos0 = this%idxlakeconn(n)
      icellid0 = this%gwfconn(ipos0)
      if (this%ictype(ipos0) /= 2 .and. this%ictype(ipos0) /= 3) cycle
      do nn2 = 1, this%nlakes
        if (nn2 == n) cycle
        j = 0
        do ipos = this%idxlakeconn(nn2), this%idxlakeconn(nn2 + 1) - 1
          j = j + 1
          icellid = this%gwfconn(ipos)
          if (icellid == icellid0) then
            if (this%ictype(ipos) == 0) then
              write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
                'EMBEDDED LAKE', n, &
                'CANNOT COINCIDE WITH VERTICAL CONNECTION', j, &
                'IN LAKE', nn2, '.'
              call store_error(errmsg)
            end if
          end if
        end do
      end do
    end do
    !
    ! -- process the data (post-processing: telev/belev/laketop/lakebot)
    do n = 1, this%nlakes
      j = 0
      do ipos = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        j = j + 1
        nn = this%gwfconn(ipos)
        top = this%dis%top(nn)
        bot = this%dis%bot(nn)
        if (this%ictype(ipos) == 0) then
          this%telev(ipos) = top + this%surfdep
          this%belev(ipos) = top
          this%lakebot(n) = min(this%belev(ipos), this%lakebot(n))
        else if (this%ictype(ipos) == 1) then
          if (this%belev(ipos) == this%telev(ipos)) then
            this%telev(ipos) = top
            this%belev(ipos) = bot
          else
            if (this%belev(ipos) >= this%telev(ipos)) then
              write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
                'telev FOR LAKE ', n, ' HORIZONTAL CONNECTION ', j, &
                'MUST BE >= belev'
              call store_error(errmsg)
            else if (this%belev(ipos) < bot) then
              write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.7,1x,a)') &
                'belev FOR LAKE ', n, ' HORIZONTAL CONNECTION ', j, &
                'MUST BE >= cell bottom (', bot, ')'
              call store_error(errmsg)
            else if (this%telev(ipos) > top) then
              write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.7,1x,a)') &
                'telev FOR LAKE ', n, ' HORIZONTAL CONNECTION ', j, &
                'MUST BE <= cell top (', top, ')'
              call store_error(errmsg)
            end if
          end if
          this%laketop(n) = max(this%telev(ipos), this%laketop(n))
          this%lakebot(n) = min(this%belev(ipos), this%lakebot(n))
        else if (this%ictype(ipos) == 2 .or. this%ictype(ipos) == 3) then
          this%telev(ipos) = top
          this%belev(ipos) = bot
          this%lakebot(n) = bot
        end if
        !
        ! -- check for missing or duplicate connections
        if (nboundchk(ipos) == 0) then
          write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
            'NO DATA SPECIFIED FOR LAKE', 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 LAKE', n, 'CONNECTION', j, &
            'SPECIFIED', nboundchk(ipos), 'TIMES'
          call store_error(errmsg)
        end if
      end do
      if (this%laketop(n) == -dep20) then
        this%laketop(n) = this%lakebot(n) + 100.
      end if
    end do
    !
    deallocate (nboundchk)
    !
    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('ICONN', this%input_mempath)
    call memorystore_release('CELLID', this%input_mempath)
    call memorystore_release('CLAKTYPE', this%input_mempath)
    call memorystore_release('BEDLEAK', this%input_mempath)
    call memorystore_release('BELEV', this%input_mempath)
    call memorystore_release('TELEV', this%input_mempath)
    call memorystore_release('CONNLEN', this%input_mempath)
    call memorystore_release('CONNWIDTH', this%input_mempath)
  end subroutine lak_source_connectiondata
 
  !> @brief Source TABLES block from input context
  !<
  subroutine lak_source_tables(this)
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    integer(I4B) :: n, ntabs
    integer(I4B) :: ilak, iconn
    character(len=LENMEMPATH) :: mempath
    integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
    type(characterstringtype), dimension(:), pointer, contiguous :: &
      laktab6_mempaths => null()
    type(laktabtype), dimension(:), allocatable :: laketables
    integer(I4B), allocatable :: nboundchk(:)
    !
    ! -- skip if no tables
    if (this%ntables < 1) return
    !
    ! -- guard: TABLES block required when ntables > 0
    call get_isize('TABLES_IFNO', this%input_mempath, ntabs)
    if (ntabs < 1) then
      write (errmsg, '(a)') &
        'NTABLES > 0 BUT A TABLES BLOCK WAS NOT FOUND IN THE LAK INPUT.'
      call store_error(errmsg, .true.)
    end if
    !
    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
      ' TABLES'
    !
    ! -- set input context pointers
    call mem_setptr(ifno, 'TABLES_IFNO', this%input_mempath)
    call mem_setptr(laktab6_mempaths, 'LAKTAB6_MEMPATH', this%input_mempath)
    !
    ntabs = size(ifno)
    !
    ! -- allocate
    allocate (laketables(this%nlakes))
    allocate (nboundchk(this%nlakes))
    nboundchk = 0
    !
    do n = 1, ntabs
      ilak = this%validate_ifno(ifno(n), this%nlakes, nboundchk, &
                                'lakeno', 'TABLES')
      if (ilak == 0) cycle
      mempath = laktab6_mempaths(n)
      call this%lak_source_table(ilak, mempath, laketables(ilak))
    end do
    !
    write (this%iout, '(1x,a)') 'END OF '//trim(adjustl(this%text))//' TABLES'
    !
    ! -- check for missing or duplicate table data
    if (ntabs < this%ntables) then
      write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
        'TABLE DATA ARE SPECIFIED', ntabs, &
        'TIMES BUT NTABLES IS SET TO', this%ntables
      call store_error(errmsg)
    end if
    ! -- a table that never loaded any rows is not reported as a duplicate
    do n = 1, this%nlakes
      if (this%ntabrow(n) == 0) nboundchk(n) = 0
    end do
    call this%report_ifno_coverage(nboundchk, this%nlakes, 'lakeno', &
                                   'TABLES', check_missing=.false.)
    !
    deallocate (nboundchk)
    !
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- release input context tables memory
    call memorystore_release('TABLES_IFNO', this%input_mempath)
    !
    ! -- convert laketables to vectors
    call this%laktables_to_vectors(laketables)
    !
    ! -- destroy laketables
    do n = 1, this%nlakes
      if (this%ntabrow(n) > 0) then
        deallocate (laketables(n)%tabstage)
        deallocate (laketables(n)%tabvolume)
        deallocate (laketables(n)%tabsarea)
        iconn = this%idxlakeconn(n)
        if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
          deallocate (laketables(n)%tabwarea)
        end if
      end if
    end do
    deallocate (laketables)
  end subroutine lak_source_tables
 
  !> @brief Copy the laketables structure data into flattened vectors that are
  !! stored in the memory manager
  !<
  subroutine laktables_to_vectors(this, laketables)
    class(laktype), intent(inout) :: this
    type(laktabtype), intent(in), dimension(:), contiguous :: laketables
    integer(I4B) :: n
    integer(I4B) :: ntabrows
    integer(I4B) :: j
    integer(I4B) :: ipos
    integer(I4B) :: iconn
    !
    ! -- allocate index array for lak tables
    call mem_allocate(this%ialaktab, this%nlakes + 1, 'IALAKTAB', this%memoryPath)
    !
    ! -- Move the laktables structure information into flattened arrays
    this%ialaktab(1) = 1
    do n = 1, this%nlakes
      ! -- ialaktab contains a pointer into the flattened lak table data
      this%ialaktab(n + 1) = this%ialaktab(n) + this%ntabrow(n)
    end do
    !
    ! -- Allocate vectors for storing all lake table data
    ntabrows = this%ialaktab(this%nlakes + 1) - 1
    call mem_allocate(this%tabstage, ntabrows, 'TABSTAGE', this%memoryPath)
    call mem_allocate(this%tabvolume, ntabrows, 'TABVOLUME', this%memoryPath)
    call mem_allocate(this%tabsarea, ntabrows, 'TABSAREA', this%memoryPath)
    call mem_allocate(this%tabwarea, ntabrows, 'TABWAREA', this%memoryPath)
    !
    ! -- Copy data from laketables into vectors
    do n = 1, this%nlakes
      j = 1
      do ipos = this%ialaktab(n), this%ialaktab(n + 1) - 1
        this%tabstage(ipos) = laketables(n)%tabstage(j)
        this%tabvolume(ipos) = laketables(n)%tabvolume(j)
        this%tabsarea(ipos) = laketables(n)%tabsarea(j)
        iconn = this%idxlakeconn(n)
        if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
          !
          ! -- tabwarea only filled for ictype 2 and 3
          this%tabwarea(ipos) = laketables(n)%tabwarea(j)
        else
          this%tabwarea(ipos) = dzero
        end if
        j = j + 1
      end do
    end do
  end subroutine laktables_to_vectors
 
  !> @brief Source one lake's table from its LAKTAB6 subpackage
  !<
  subroutine lak_source_table(this, ilak, mempath, laketable)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    character(len=*), intent(in) :: mempath
    type(laktabtype), intent(inout) :: laketable
    ! -- local
    integer(I4B), pointer :: nrow => null()
    integer(I4B), pointer :: ncol => null()
    real(DP), dimension(:), pointer, contiguous :: stage => null()
    real(DP), dimension(:), pointer, contiguous :: volume => null()
    real(DP), dimension(:), pointer, contiguous :: sarea => null()
    real(DP), dimension(:), pointer, contiguous :: barea => null()
    character(len=:), pointer :: tabfname => null()
    integer(I4B) :: jmin, iconn, n, nerr0
    real(DP) :: vol, sa, wa, v, v0
    ! -- formats
    character(len=*), parameter :: fmttaberr = &
      &'(a,1x,i0,1x,a,1x,g15.6,1x,a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.6,1x,a)'
    !
    ! -- this table's own input file, for error attribution below
    call mem_setptr(tabfname, 'INPUT_FNAME', mempath)
    call mem_setptr(nrow, 'NROW', mempath)
    call mem_setptr(ncol, 'NCOL', mempath)
    !
    ! -- errors added below this point belong to this lake's own table
    nerr0 = count_errors()
    !
    if (nrow < 1) then
      write (errmsg, '(a,1x,i0)') &
        'LAKE TABLE NROW MUST BE > 0 FOR LAKE', ilak
      call store_error(errmsg)
    end if
    !
    iconn = this%idxlakeconn(ilak)
    if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
      jmin = 4
    else
      jmin = 3
    end if
    if (ncol < jmin) then
      write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
        'LAKE TABLE NCOL MUST BE >=', jmin, 'FOR LAKE', ilak
      call store_error(errmsg)
    end if
    !
    if (count_errors() > nerr0) then
      call store_error_filename(tabfname, terminate=.false.)
      call memorystore_release('NROW', mempath)
      call memorystore_release('NCOL', mempath)
      call memorystore_release('INPUT_FNAME', mempath)
      return
    end if
    !
    this%ntabrow(ilak) = nrow
    allocate (laketable%tabstage(nrow))
    allocate (laketable%tabvolume(nrow))
    allocate (laketable%tabsarea(nrow))
    call mem_setptr(stage, 'STAGE', mempath)
    call mem_setptr(volume, 'VOLUME', mempath)
    call mem_setptr(sarea, 'SAREA', mempath)
    laketable%tabstage = stage
    laketable%tabvolume = volume
    laketable%tabsarea = sarea
    if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
      allocate (laketable%tabwarea(nrow))
      call mem_setptr(barea, 'BAREA', mempath)
      laketable%tabwarea = barea
    end if
    !
    ! -- set lake bottom based on table if it is an embedded lake
    if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
      do n = 1, this%ntabrow(ilak)
        vol = laketable%tabvolume(n)
        sa = laketable%tabsarea(n)
        wa = laketable%tabwarea(n)
        vol = vol * sa * wa
        ! -- check if all entries are zero
        if (vol > dzero) exit
        ! -- set lake bottom
        this%lakebot(ilak) = laketable%tabstage(n)
        this%belev(ilak) = laketable%tabstage(n)
      end do
      ! -- set maximum surface area for rainfall
      n = this%ntabrow(ilak)
      this%sareamax(ilak) = laketable%tabsarea(n)
    end if
    !
    ! -- verify the table data
    do n = 2, this%ntabrow(ilak)
      v = laketable%tabstage(n)
      v0 = laketable%tabstage(n - 1)
      if (v <= v0) then
        write (errmsg, fmttaberr) &
          'TABLE STAGE ENTRY', n, '(', laketable%tabstage(n), ') FOR LAKE ', &
          ilak, 'MUST BE GREATER THAN THE PREVIOUS STAGE ENTRY', &
          n - 1, '(', laketable%tabstage(n - 1), ')'
        call store_error(errmsg)
      end if
      v = laketable%tabvolume(n)
      v0 = laketable%tabvolume(n - 1)
      if (v <= v0) then
        write (errmsg, fmttaberr) &
          'TABLE VOLUME ENTRY', n, '(', laketable%tabvolume(n), &
          ') FOR LAKE ', &
          ilak, 'MUST BE GREATER THAN THE PREVIOUS VOLUME ENTRY', &
          n - 1, '(', laketable%tabvolume(n - 1), ')'
        call store_error(errmsg)
      end if
      v = laketable%tabsarea(n)
      v0 = laketable%tabsarea(n - 1)
      if (v < v0) then
        write (errmsg, fmttaberr) &
          'TABLE SURFACE AREA ENTRY', n, '(', &
          laketable%tabsarea(n), ') FOR LAKE ', ilak, &
          'MUST BE GREATER THAN OR EQUAL TO THE PREVIOUS SURFACE AREA ENTRY', &
          n - 1, '(', laketable%tabsarea(n - 1), ')'
        call store_error(errmsg)
      end if
      if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
        v = laketable%tabwarea(n)
        v0 = laketable%tabwarea(n - 1)
        if (v < v0) then
          write (errmsg, fmttaberr) &
            'TABLE EXCHANGE AREA ENTRY', n, '(', &
            laketable%tabwarea(n), ') FOR LAKE ', ilak, &
            'MUST BE GREATER THAN OR EQUAL TO THE PREVIOUS EXCHANGE AREA '// &
            'ENTRY', n - 1, '(', laketable%tabwarea(n - 1), ')'
          call store_error(errmsg)
        end if
      end if
    end do
    !
    if (count_errors() > nerr0) then
      call store_error_filename(tabfname, terminate=.false.)
    end if
    !
    ! -- release this table's own input context memory
    call memorystore_release('NROW', mempath)
    call memorystore_release('NCOL', mempath)
    call memorystore_release('STAGE', mempath)
    call memorystore_release('VOLUME', mempath)
    call memorystore_release('SAREA', mempath)
    if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
      call memorystore_release('BAREA', mempath)
    end if
    call memorystore_release('INPUT_FNAME', mempath)
  end subroutine lak_source_table
 
  !> @brief Source OUTLETS block from input context
  !<
  subroutine lak_source_outlets(this)
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    integer(I4B) :: n, isize, ioutlet
    character(len=LINELENGTH) :: keyword
    integer(I4B), dimension(:), pointer, contiguous :: outletno => null()
    integer(I4B), dimension(:), pointer, contiguous :: lakein => null()
    integer(I4B), dimension(:), pointer, contiguous :: lakeout => null()
    type(characterstringtype), dimension(:), pointer, contiguous :: &
      couttype => null()
    real(DP), dimension(:), pointer, contiguous :: outlets_invert => null()
    real(DP), dimension(:), pointer, contiguous :: outlets_width => null()
    real(DP), dimension(:), pointer, contiguous :: outlets_rough => null()
    real(DP), dimension(:), pointer, contiguous :: outlets_slope => null()
    integer(I4B), allocatable :: nboundchk(:)
    !
    ! -- skip if no outlets
    if (this%noutlets < 1) return
    !
    ! -- guard: OUTLETS block required when noutlets > 0
    call get_isize('OUTLETS_OUTLETNO', this%input_mempath, isize)
    if (isize < 1) then
      write (errmsg, '(a)') &
        'NOUTLETS > 0 BUT AN OUTLETS BLOCK WAS NOT FOUND IN THE LAK INPUT.'
      call store_error(errmsg, .true.)
    end if
    !
    ! -- allocate outlet arrays
    call mem_allocate(this%lakein, this%noutlets, 'LAKEIN', this%memoryPath)
    call mem_allocate(this%lakeout, this%noutlets, 'LAKEOUT', this%memoryPath)
    call mem_allocate(this%iouttype, this%noutlets, 'IOUTTYPE', this%memoryPath)
    ! -- alias into the input context's permanent, feature-indexed arrays
    ! (allocated and DZERO-initialized by the loader)
    call mem_setptr(this%outrate, 'RATE', this%input_mempath)
    call mem_setptr(this%outinvert, 'INVERT', this%input_mempath)
    call mem_setptr(this%outwidth, 'WIDTH', this%input_mempath)
    call mem_setptr(this%outrough, 'ROUGH', this%input_mempath)
    call mem_setptr(this%outslope, 'SLOPE', this%input_mempath)
    call mem_allocate(this%simoutrate, this%noutlets, 'SIMOUTRATE', &
                      this%memoryPath)
    call mem_allocate(this%outlet_invert_set, this%noutlets, &
                      'INVERT_SET', this%memoryPath)
    call mem_allocate(this%outlet_width_set, this%noutlets, &
                      'WIDTH_SET', this%memoryPath)
    call mem_allocate(this%outlet_rough_set, this%noutlets, &
                      'ROUGH_SET', this%memoryPath)
    call mem_allocate(this%outlet_slope_set, this%noutlets, &
                      'SLOPE_SET', this%memoryPath)
    this%outlet_invert_set = .false.
    this%outlet_width_set = .false.
    this%outlet_rough_set = .false.
    this%outlet_slope_set = .false.
    !
    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
      ' OUTLETS'
    !
    ! -- set input context pointers
    call mem_setptr(outletno, 'OUTLETS_OUTLETNO', this%input_mempath)
    call mem_setptr(lakein, 'LAKEIN', this%input_mempath)
    call mem_setptr(lakeout, 'LAKEOUT', this%input_mempath)
    call mem_setptr(couttype, 'COUTTYPE', this%input_mempath)
    call mem_setptr(outlets_invert, 'OUTLETS_INVERT', this%input_mempath)
    call mem_setptr(outlets_width, 'OUTLETS_WIDTH', this%input_mempath)
    call mem_setptr(outlets_rough, 'OUTLETS_ROUGH', this%input_mempath)
    call mem_setptr(outlets_slope, 'OUTLETS_SLOPE', this%input_mempath)
    !
    allocate (nboundchk(this%noutlets))
    nboundchk = 0
    !
    do n = 1, size(outletno)
      ioutlet = this%validate_ifno(outletno(n), this%noutlets, &
                                   nboundchk, 'outletno', 'OUTLETS')
      if (ioutlet == 0) cycle
      !
      ! -- lakein
      if (lakein(n) < 1 .or. lakein(n) > this%nlakes) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
          'lakein FOR OUTLET ', ioutlet, 'MUST BE > 0 and <= ', this%nlakes
        call store_error(errmsg)
      end if
      this%lakein(ioutlet) = lakein(n)
      !
      ! -- lakeout
      if (lakeout(n) < 0 .or. lakeout(n) > this%nlakes) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
          'lakeout FOR OUTLET ', ioutlet, 'MUST BE >= 0 and <= ', this%nlakes
        call store_error(errmsg)
      end if
      this%lakeout(ioutlet) = lakeout(n)
      !
      ! -- outlet type
      keyword = couttype(n)
      select case (trim(keyword))
      case ('SPECIFIED')
        this%iouttype(ioutlet) = 0
      case ('MANNING')
        this%iouttype(ioutlet) = 1
      case ('WEIR')
        this%iouttype(ioutlet) = 2
      case default
        write (errmsg, '(a,1x,i0,1x,a,a,a)') &
          'UNKNOWN couttype FOR OUTLET ', ioutlet, '(', trim(keyword), ')'
        call store_error(errmsg)
      end select
      !
      ! -- outlet invert, width, roughness, slope
      this%outinvert(ioutlet) = outlets_invert(n)
      this%outwidth(ioutlet) = outlets_width(n)
      this%outrough(ioutlet) = outlets_rough(n)
      this%outslope(ioutlet) = outlets_slope(n)
    end do
    !
    write (this%iout, '(1x,a)') 'END OF '//trim(adjustl(this%text))//' OUTLETS'
    !
    ! -- check for duplicate or missing outlets
    call this%report_ifno_coverage(nboundchk, this%noutlets, 'outlet', &
                                   'OUTLETS')
    !
    deallocate (nboundchk)
    !
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- release input context outlets memory; OUTLETS_OUTLETNO/INVERT/
    !    WIDTH/ROUGH/SLOPE excluded (bound by lak_input_init for PERIOD use)
    call memorystore_release('LAKEIN', this%input_mempath)
    call memorystore_release('LAKEOUT', this%input_mempath)
    call memorystore_release('COUTTYPE', this%input_mempath)
  end subroutine lak_source_outlets
 
  !> @brief Source DIMENSIONS block from input context and load static data
  !<
  subroutine lak_source_dimensions(this)
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    type(gwflakparamfoundtype) :: found
    !
    ! -- initialize dimensions to -1
    this%nlakes = -1
    this%maxbound = -1
    !
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' DIMENSIONS'
    !
    ! -- NLAKES (mf6internal: nlakes); kept in the input context (not
    !    released) since PERIOD settings resolve it by name via SHAPE
    call mem_set_value(this%nlakes, 'NLAKES', this%input_mempath, &
                       found%nlakes, release=.false.)
    if (found%nlakes) write (this%iout, '(4x,a,i7)') 'NLAKES = ', this%nlakes
    !
    ! -- NOUTLETS (mf6internal: noutlets); kept in the input context (not
    !    released) since PERIOD settings resolve it by name via SHAPE
    call mem_set_value(this%noutlets, 'NOUTLETS', this%input_mempath, &
                       found%noutlets, release=.false.)
    if (found%noutlets) write (this%iout, '(4x,a,i7)') &
      'NOUTLETS = ', this%noutlets
    !
    ! -- NTABLES (mf6internal: ntables)
    call mem_set_value(this%ntables, 'NTABLES', this%input_mempath, found%ntables)
    if (found%ntables) write (this%iout, '(4x,a,i7)') 'NTABLES = ', this%ntables
    !
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%text))//' DIMENSIONS'
    !
    if (this%nlakes < 0) then
      write (errmsg, '(a)') &
        'NLAKES WAS NOT SPECIFIED OR WAS SPECIFIED INCORRECTLY.'
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
      return
    end if
    !
    if (this%iforceleglak /= 0) then
      if (this%iforceleglak < 1 .or. this%iforceleglak > this%nlakes) then
        write (errmsg, '(a,i0,a,i0,a)') &
          'DEV_FORCE_LEGACY_LAKE (', this%iforceleglak, &
          ') MUST BE BETWEEN 1 AND NLAKES (', this%nlakes, ').'
        call store_error(errmsg)
        call store_error_filename(this%input_fname)
        return
      end if
    end if
    !
    ! -- read PACKAGEDATA
    call this%lak_source_packagedata()
    !
    ! -- for the implicit formulation each lake adds one equation (one row and
    !    column) to the groundwater flow matrix
    if (this%iimplicit /= 0) then
      this%npakeq = this%nlakes
      ! -- the implicit lake-aquifer coupling is asymmetric for perched
      !    connections, so flag the coefficient matrix as asymmetric (which
      !    requires the BICGSTAB linear acceleration). Only raise the flag; never
      !    clear an asymmetry already indicated elsewhere.
      this%iasym = 1
    end if
    !
    ! -- read CONNECTIONDATA
    call this%lak_source_connectiondata()
    !
    ! -- read TABLES
    call this%lak_source_tables()
    !
    ! -- read OUTLETS
    call this%lak_source_outlets()
    !
    ! -- Call define_listlabel
    call this%define_listlabel()
    !
    ! -- setup the budget object
    call this%lak_setup_budobj()
    !
    ! -- setup the stage table object
    call this%lak_setup_tableobj()
  end subroutine lak_source_dimensions
 
  !> @brief Read the initial parameters for this package
  !<
  subroutine lak_read_initial_attr(this)
    use constantsmodule, only: linelength
    use memoryhelpermodule, only: create_mem_path
    use simmodule, only: store_error, count_errors
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    integer(I4B) :: j, n
    integer(I4B) :: nn
    integer(I4B) :: idx
    real(DP) :: top
    real(DP) :: bot
    real(DP) :: k
    real(DP) :: area
    real(DP) :: length
    real(DP) :: s
    real(DP) :: dx
    real(DP) :: c
    real(DP) :: sa
    real(DP) :: wa
    real(DP) :: v
    real(DP) :: fact
    real(DP) :: c1
    real(DP) :: c2
    real(DP), allocatable, dimension(:) :: clb, caq
    character(len=14) :: cbedleak
    character(len=14) :: cbedcond
    character(len=10), dimension(0:3) :: ctype
    character(len=15) :: nodestr
    ! -- data
    data ctype(0)/'VERTICAL  '/
    data ctype(1)/'HORIZONTAL'/
    data ctype(2)/'EMBEDDEDH '/
    data ctype(3)/'EMBEDDEDV '/
    !
    ! -- initialize xnewpak and set stage
    do n = 1, this%nlakes
      this%xnewpak(n) = this%strt(n)
      this%stage(n) = this%strt(n)
    end do
    !
    ! -- initialize status (iboundpak) of lakes to active
    do n = 1, this%nlakes
      if (this%status(n) == 'CONSTANT') then
        this%iboundpak(n) = -1
      else if (this%status(n) == 'INACTIVE') then
        this%iboundpak(n) = 0
      else if (this%status(n) == 'ACTIVE ') then
        this%iboundpak(n) = 1
      end if
    end do
    !
    ! -- set boundname for each connection
    if (this%inamedbound /= 0) then
      do n = 1, this%nlakes
        do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
          this%boundname(j) = this%lakename(n)
        end do
      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%gwfk33, 'K33', create_mem_path(this%name_model, 'NPF'))
    call mem_setptr(this%gwfik33, 'IK33', create_mem_path(this%name_model, 'NPF'))
    call mem_setptr(this%gwfsat, 'SAT', create_mem_path(this%name_model, 'NPF'))
    !
    ! -- allocate temporary storage
    allocate (clb(this%MAXBOUND))
    allocate (caq(this%MAXBOUND))
    !
    ! -- calculate saturated conductance for each connection
    do n = 1, this%nlakes
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        nn = this%gwfconn(j)
        top = this%dis%top(nn)
        bot = this%dis%bot(nn)
        ! vertical connection
        if (this%ictype(j) == 0) then
          area = this%dis%area(nn)
          this%sarea(j) = area
          this%warea(j) = area
          this%sareamax(n) = this%sareamax(n) + area
          if (this%gwfik33 == 0) then
            k = this%gwfk11(nn)
          else
            k = this%gwfk33(nn)
          end if
          length = dhalf * (top - bot)
          ! horizontal connection
        else if (this%ictype(j) == 1) then
          area = (this%telev(j) - this%belev(j)) * this%connwidth(j)
          ! -- recalculate area if connected cell is confined and lake
          !    connection top and bot are equal to the cell top and bot
          if (top == this%telev(j) .and. bot == this%belev(j)) then
            if (this%icelltype(nn) == 0) then
              area = this%gwfsat(nn) * (top - bot) * this%connwidth(j)
            end if
          end if
          this%sarea(j) = dzero
          this%warea(j) = area
          this%sareamax(n) = this%sareamax(n) + dzero
          k = this%gwfk11(nn)
          length = this%connlength(j)
          ! embedded horizontal connection
        else if (this%ictype(j) == 2) then
          area = done
          this%sarea(j) = dzero
          this%warea(j) = area
          this%sareamax(n) = this%sareamax(n) + dzero
          k = this%gwfk11(nn)
          length = this%connlength(j)
          ! embedded vertical connection
        else if (this%ictype(j) == 3) then
          area = done
          this%sarea(j) = dzero
          this%warea(j) = area
          this%sareamax(n) = this%sareamax(n) + dzero
          if (this%gwfik33 == 0) then
            k = this%gwfk11(nn)
          else
            k = this%gwfk33(nn)
          end if
          length = this%connlength(j)
        end if
        if (is_close(this%bedleak(j), dnodata)) then
          clb(j) = dnodata
        else if (this%bedleak(j) > dzero) then
          clb(j) = done / this%bedleak(j)
        else
          clb(j) = dzero
        end if
        if (k > dzero) then
          caq(j) = length / k
        else
          caq(j) = dzero
        end if
        if (is_close(this%bedleak(j), dnodata)) then
          this%satcond(j) = area / caq(j)
        else if (clb(j) * caq(j) > dzero) then
          this%satcond(j) = area / (clb(j) + caq(j))
        else
          this%satcond(j) = dzero
        end if
      end do
    end do
    !
    ! -- write a summary of the conductance
    if (this%iprpak > 0) then
      write (this%iout, '(//,29x,a,/)') &
        'INTERFACE CONDUCTANCE BETWEEN LAKE AND AQUIFER CELLS'
      write (this%iout, '(1x,a)') &
     &  '      LAKE CONNECTION                 CONNECTION    LAKEBED'// &
     &  '              C O N D U C T A N C E S        '
      write (this%iout, '(1x,a)') &
     &  '    NUMBER     NUMBER CELLID          DIRECTION    LEAKANCE'// &
     &  '        LAKEBED        AQUIFER       COMBINED'
      write (this%iout, "(1x,108('-'))")
      do n = 1, this%nlakes
        idx = 0
        do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
          idx = idx + 1
          fact = done
          if (this%ictype(j) == 1) then
            fact = this%telev(j) - this%belev(j)
            if (abs(fact) > dzero) then
              fact = done / fact
            end if
          end if
          nn = this%gwfconn(j)
          area = this%warea(j)
          c1 = dzero
          if (is_close(clb(j), dnodata)) then
            cbedleak = '     NONE     '
            cbedcond = '     NONE     '
          else if (clb(j) > dzero) then
            c1 = area * fact / clb(j)
            write (cbedleak, '(g14.5)') this%bedleak(j)
            write (cbedcond, '(g14.5)') c1
          else
            write (cbedleak, '(g14.5)') c1
            write (cbedcond, '(g14.5)') c1
          end if
          c2 = dzero
          if (caq(j) > dzero) then
            c2 = area * fact / caq(j)
          end if
          call this%dis%noder_to_string(nn, nodestr)
          write (this%iout, &
                 '(1x,i10,1x,i10,1x,a15,1x,a10,2(1x,a14),2(1x,g14.5))') &
            n, idx, nodestr, ctype(this%ictype(j)), cbedleak, &
            cbedcond, c2, this%satcond(j) * fact
        end do
      end do
      write (this%iout, "(1x,108('-'))")
      write (this%iout, '(1x,a)') &
        'IF VERTICAL CONNECTION, CONDUCTANCE (L^2/T) IS &
        &BETWEEN AQUIFER CELL AND OVERLYING LAKE CELL.'
      write (this%iout, '(1x,a)') &
        'IF HORIZONTAL CONNECTION, CONDUCTANCES ARE PER &
        &UNIT SATURATED THICKNESS (L/T).'
      write (this%iout, '(1x,a)') &
        'IF EMBEDDED CONNECTION, CONDUCTANCES ARE PER &
        &UNIT EXCHANGE AREA (1/T).'
      !
      !        write(this%iout,*) n, idx, nodestr, this%sarea(j), this%warea(j)
      !
      ! -- calculate stage, surface area, wetted area, volume relation
      do n = 1, this%nlakes
        write (this%iout, '(//1x,a,1x,i10)') 'STAGE/VOLUME RELATION FOR LAKE  ', n
        write (this%iout, '(/1x,5(a14))') '         STAGE', '  SURFACE AREA', &
     &                                    '   WETTED AREA', '   CONDUCTANCE', &
     &                                    '        VOLUME'
        write (this%iout, "(1x,70('-'))")
        dx = (this%laketop(n) - this%lakebot(n)) / 150.
        s = this%lakebot(n)
        do j = 1, 151
          call this%lak_calculate_conductance(n, s, c)
          call this%lak_calculate_sarea(n, s, sa)
          call this%lak_calculate_warea(n, s, wa, s)
          call this%lak_calculate_vol(n, s, v)
          write (this%iout, '(1x,5(E14.5))') s, sa, wa, c, v
          s = s + dx
        end do
        write (this%iout, "(1x,70('-'))")
        !
        write (this%iout, '(//1x,a,1x,i10)') 'STAGE/VOLUME RELATION FOR LAKE  ', n
        write (this%iout, '(/1x,4(a14))') '              ', '              ', &
     &                                    '    CALCULATED', '         STAGE'
        write (this%iout, '(1x,4(a14))') '         STAGE', '        VOLUME', &
    &                                    '         STAGE', '    DIFFERENCE'
        write (this%iout, "(1x,56('-'))")
        s = this%lakebot(n) - dx
        do j = 1, 156
          call this%lak_calculate_vol(n, s, v)
          call this%lak_vol2stage(n, v, c)
          write (this%iout, '(1x,4(E14.5))') s, v, c, s - c
          s = s + dx
        end do
        write (this%iout, "(1x,56('-'))")
      end do
    end if
    !
    ! -- finished with pointer to gwf hydraulic conductivity
    this%gwfk11 => null()
    this%gwfk33 => null()
    this%gwfsat => null()
    this%gwfik33 => null()
    !
    ! -- deallocate temporary storage
    deallocate (clb)
    deallocate (caq)
  end subroutine lak_read_initial_attr
 
  !> @brief Perform linear interpolation of two vectors.
  !!
  !! Function assumes x data is sorted in ascending order
  !<
  subroutine lak_linear_interpolation(this, n, x, y, z, v)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: n
    real(DP), dimension(n), intent(in) :: x
    real(DP), dimension(n), intent(in) :: y
    real(DP), intent(in) :: z
    real(DP), intent(inout) :: v
    ! -- local
    integer(I4B) :: i
    real(DP) :: dx, dydx
    ! code
    v = dzero
    ! below bottom of range - set to lowest value
    if (z <= x(1)) then
      v = y(1)
      ! above highest value
      ! slope calculated from interval between n and n-1
    else if (z > x(n)) then
      dx = x(n) - x(n - 1)
      dydx = dzero
      if (abs(dx) > dzero) then
        dydx = (y(n) - y(n - 1)) / dx
      end if
      dx = (z - x(n))
      v = y(n) + dydx * dx
      ! between lowest and highest value in current interval
    else
      do i = 2, n
        dx = x(i) - x(i - 1)
        dydx = dzero
        if (z >= x(i - 1) .and. z <= x(i)) then
          if (abs(dx) > dzero) then
            dydx = (y(i) - y(i - 1)) / dx
          end if
          dx = (z - x(i - 1))
          v = y(i - 1) + dydx * dx
          exit
        end if
      end do
    end if
  end subroutine lak_linear_interpolation
 
  !> @brief Calculate the surface area of a lake at a given stage
  !<
  subroutine lak_calculate_sarea(this, ilak, stage, sarea)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(inout) :: sarea
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: ifirst
    integer(I4B) :: ilast
    real(DP) :: topl
    real(DP) :: botl
    real(DP) :: sat
    real(DP) :: sa
    !
    sarea = dzero
    i = this%ntabrow(ilak)
    if (i > 0) then
      ifirst = this%ialaktab(ilak)
      ilast = this%ialaktab(ilak + 1) - 1
      if (stage <= this%tabstage(ifirst)) then
        sarea = this%tabsarea(ifirst)
      else if (stage >= this%tabstage(ilast)) then
        sarea = this%tabsarea(ilast)
      else
        call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
                                           this%tabsarea(ifirst:ilast), &
                                           stage, sarea)
      end if
    else
      do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
        topl = this%telev(i)
        botl = this%belev(i)
        sat = squadraticsaturation(topl, botl, stage)
        sa = sat * this%sarea(i)
        sarea = sarea + sa
      end do
    end if
  end subroutine lak_calculate_sarea
 
  !> @brief Calculate the wetted area of a lake at a given stage.
  !<
  subroutine lak_calculate_warea(this, ilak, stage, warea, hin)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(inout) :: warea
    real(DP), optional, intent(inout) :: hin
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: igwfnode
    real(DP) :: head
    real(DP) :: wa
    !
    warea = dzero
    do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
      if (present(hin)) then
        head = hin
      else
        igwfnode = this%gwfconn(i)
        head = this%xnew(igwfnode)
      end if
      call this%lak_calculate_conn_warea(ilak, i, stage, head, wa)
      warea = warea + wa
    end do
  end subroutine lak_calculate_warea
 
  !> @brief Calculate the wetted area of a lake connection at a given stage
  !<
  subroutine lak_calculate_conn_warea(this, ilak, iconn, stage, head, wa)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    integer(I4B), intent(in) :: iconn
    real(DP), intent(in) :: stage
    real(DP), intent(in) :: head
    real(DP), intent(inout) :: wa
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: ifirst
    integer(I4B) :: ilast
    integer(I4B) :: node
    real(DP) :: topl
    real(DP) :: botl
    real(DP) :: vv
    real(DP) :: sat
    !
    wa = dzero
    topl = this%telev(iconn)
    botl = this%belev(iconn)
    call this%lak_calculate_cond_head(iconn, stage, head, vv)
    if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
      if (vv > topl) vv = topl
      i = this%ntabrow(ilak)
      ifirst = this%ialaktab(ilak)
      ilast = this%ialaktab(ilak + 1) - 1
      if (vv <= this%tabstage(ifirst)) then
        wa = this%tabwarea(ifirst)
      else if (vv >= this%tabstage(ilast)) then
        wa = this%tabwarea(ilast)
      else
        call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
                                           this%tabwarea(ifirst:ilast), &
                                           vv, wa)
      end if
    else
      node = this%gwfconn(iconn)
      ! -- confined cell
      if (this%icelltype(node) == 0) then
        sat = done
        ! -- convertible cell
      else
        sat = squadraticsaturation(topl, botl, vv)
      end if
      wa = sat * this%warea(iconn)
    end if
  end subroutine lak_calculate_conn_warea
 
  !> @brief Calculate the volume of a lake at a given stage
  !<
  subroutine lak_calculate_vol(this, ilak, stage, volume)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(inout) :: volume
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: ifirst
    integer(I4B) :: ilast
    real(DP) :: topl
    real(DP) :: botl
    real(DP) :: ds
    real(DP) :: sa
    real(DP) :: v
    real(DP) :: sat
    !
    volume = dzero
    i = this%ntabrow(ilak)
    if (i > 0) then
      ifirst = this%ialaktab(ilak)
      ilast = this%ialaktab(ilak + 1) - 1
      if (stage <= this%tabstage(ifirst)) then
        volume = this%tabvolume(ifirst)
      else if (stage >= this%tabstage(ilast)) then
        ds = stage - this%tabstage(ilast)
        sa = this%tabsarea(ilast)
        volume = this%tabvolume(ilast) + ds * sa
      else
        call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
                                           this%tabvolume(ifirst:ilast), &
                                           stage, volume)
      end if
    else
      do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
        topl = this%telev(i)
        botl = this%belev(i)
        sat = squadraticsaturation(topl, botl, stage)
        sa = sat * this%sarea(i)
        if (stage < botl) then
          v = dzero
        else if (stage > botl .and. stage < topl) then
          v = sa * (stage - botl)
        else
          v = sa * (topl - botl) + sa * (stage - topl)
        end if
        volume = volume + v
      end do
    end if
  end subroutine lak_calculate_vol
 
  !> @brief Calculate the total conductance for a lake at a provided stage
  !<
  subroutine lak_calculate_conductance(this, ilak, stage, conductance)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(inout) :: conductance
    ! -- local
    integer(I4B) :: i
    real(DP) :: c
    !
    conductance = dzero
    do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
      call this%lak_calculate_conn_conductance(ilak, i, stage, stage, c)
      conductance = conductance + c
    end do
  end subroutine lak_calculate_conductance
 
  !> @brief Calculate the controlling lake stage or groundwater head used to
  !! calculate the conductance for a lake connection from a provided stage and
  !! groundwater head
  !<
  subroutine lak_calculate_cond_head(this, iconn, stage, head, vv)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: iconn
    real(DP), intent(in) :: stage
    real(DP), intent(in) :: head
    real(DP), intent(inout) :: vv
    ! -- local
    real(DP) :: ss
    real(DP) :: hh
    real(DP) :: topl
    real(DP) :: botl
    !
    topl = this%telev(iconn)
    botl = this%belev(iconn)
    ss = min(stage, topl)
    hh = min(head, topl)
    if (this%igwhcopt > 0) then
      vv = hh
    else if (this%inewton > 0) then
      vv = max(ss, hh)
    else
      vv = dhalf * (ss + hh)
    end if
  end subroutine lak_calculate_cond_head
 
  !> @brief Calculate the conductance for a lake connection at a provided stage
  !! and groundwater head
  !<
  subroutine lak_calculate_conn_conductance(this, ilak, iconn, stage, head, cond)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    integer(I4B), intent(in) :: iconn
    real(DP), intent(in) :: stage
    real(DP), intent(in) :: head
    real(DP), intent(inout) :: cond
    ! -- local
    integer(I4B) :: node
    !real(DP) :: ss
    !real(DP) :: hh
    real(DP) :: vv
    real(DP) :: topl
    real(DP) :: botl
    real(DP) :: sat
    real(DP) :: wa
    real(DP) :: vscratio
    !
    cond = dzero
    vscratio = done
    topl = this%telev(iconn)
    botl = this%belev(iconn)
    call this%lak_calculate_cond_head(iconn, stage, head, vv)
    sat = squadraticsaturation(topl, botl, vv)
    ! vertical connection
    ! use full saturated conductance if top and bottom of the lake connection
    ! are equal
    if (this%ictype(iconn) == 0) then
      if (abs(topl - botl) < dprec) then
        sat = done
      end if
      ! horizontal connection
      ! use full saturated conductance if the connected cell is not convertible
    else if (this%ictype(iconn) == 1) then
      node = this%gwfconn(iconn)
      if (this%icelltype(node) == 0) then
        sat = done
      end if
      ! embedded connection
    else if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
      node = this%gwfconn(iconn)
      if (this%icelltype(node) == 0) then
        vv = this%telev(iconn)
        call this%lak_calculate_conn_warea(ilak, iconn, vv, vv, wa)
      else
        call this%lak_calculate_conn_warea(ilak, iconn, stage, head, wa)
      end if
      sat = wa
    end if
    !
    ! -- account for viscosity effects (if vsc active)
    if (this%ivsc == 1) then
      ! flow from lake to aquifer
      if (stage > head) then
        vscratio = this%viscratios(1, iconn)
        ! flow from aquifer to lake
      else
        vscratio = this%viscratios(2, iconn)
      end if
    end if
    cond = sat * this%satcond(iconn) * vscratio
  end subroutine lak_calculate_conn_conductance
 
  !> @brief Calculate the total groundwater-lake flow at a provided stage
  !<
  subroutine lak_calculate_exchange(this, ilak, stage, totflow)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(inout) :: totflow
    ! -- local
    integer(I4B) :: j
    integer(I4B) :: igwfnode
    real(DP) :: flow
    real(DP) :: hgwf
    !
    totflow = dzero
    do j = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
      igwfnode = this%gwfconn(j)
      hgwf = this%xnew(igwfnode)
      call this%lak_calculate_conn_exchange(ilak, j, stage, hgwf, flow)
      totflow = totflow + flow
    end do
  end subroutine lak_calculate_exchange
 
  !> @brief Calculate the groundwater-lake flow at a provided stage and
  !! groundwater head
  !<
  subroutine lak_calculate_conn_exchange(this, ilak, iconn, stage, head, flow, &
                                         gwfhcof, gwfrhs)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    integer(I4B), intent(in) :: iconn
    real(DP), intent(in) :: stage
    real(DP), intent(in) :: head
    real(DP), intent(inout) :: flow
    real(DP), intent(inout), optional :: gwfhcof
    real(DP), intent(inout), optional :: gwfrhs
    ! -- local
    real(DP) :: botl
    real(DP) :: cond
    real(DP) :: ss
    real(DP) :: hh
    real(DP) :: gwfhcof0
    real(DP) :: gwfrhs0
    !
    flow = dzero
    call this%lak_calculate_conn_conductance(ilak, iconn, stage, head, cond)
    botl = this%belev(iconn)
    !
    ! -- Set ss to stage or botl
    if (stage >= botl) then
      ss = stage
    else
      ss = botl
    end if
    !
    ! -- set hh to head or botl
    if (head >= botl) then
      hh = head
    else
      hh = botl
    end if
    !
    ! -- calculate flow, positive into lake
    flow = cond * (hh - ss)
    !
    ! -- Calculate gwfhcof and gwfrhs
    if (head >= botl) then
      gwfhcof0 = -cond
      gwfrhs0 = -cond * ss
    else
      gwfhcof0 = dzero
      gwfrhs0 = flow
    end if
    !
    ! Add density contributions, if active
    if (this%idense /= 0) then
      call this%lak_calculate_density_exchange(iconn, stage, head, cond, botl, &
                                               flow, gwfhcof0, gwfrhs0)
    end if
    !
    ! -- If present update gwfhcof and gwfrhs
    if (present(gwfhcof)) gwfhcof = gwfhcof0
    if (present(gwfrhs)) gwfrhs = gwfrhs0
  end subroutine lak_calculate_conn_exchange
 
  !> @brief Lakebed seepage and its derivatives for a connection (IMPLICIT)
  !!
  !! Returns the lakebed seepage for one connection, positive into the lake,
  !! and optionally its stage and head derivatives. The stage and the connected-
  !! cell head are each held at or above the lake bottom, the same wet/dry cutoff
  !! used by the default formulation, so the implicit seepage matches the default
  !! formulation exactly. lak_fc_implicit uses this to assemble the matrix and
  !! lak_cq uses it to report the budget, so the two are consistent.
  !<
  subroutine lak_calculate_conn_exchange_deriv(this, ilak, iconn, stage, &
                                               head, flow, dqds, dqdh)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    integer(I4B), intent(in) :: iconn
    real(DP), intent(in) :: stage
    real(DP), intent(in) :: head
    real(DP), intent(inout) :: flow
    real(DP), intent(inout), optional :: dqds
    real(DP), intent(inout), optional :: dqdh
    ! -- local
    real(DP) :: cond, botl, dps, dph, ss, hh
    !
    call this%lak_calculate_conn_conductance(ilak, iconn, stage, head, cond)
    botl = this%belev(iconn)
    if (stage >= botl) then
      ss = stage
      dps = done
    else
      ss = botl
      dps = dzero
    end if
    if (head >= botl) then
      hh = head
      dph = done
    else
      hh = botl
      dph = dzero
    end if
    flow = cond * (hh - ss)
    if (present(dqds)) dqds = -cond * dps
    if (present(dqdh)) dqdh = cond * dph
  end subroutine lak_calculate_conn_exchange_deriv
 
  !> @brief Calculate the groundwater-lake flow at a provided stage and
  !! groundwater head
  !<
  subroutine lak_estimate_conn_exchange(this, iflag, ilak, iconn, idry, stage, &
                                        head, flow, source, gwfhcof, gwfrhs)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: iflag
    integer(I4B), intent(in) :: ilak
    integer(I4B), intent(in) :: iconn
    integer(I4B), intent(inout) :: idry
    real(DP), intent(in) :: stage
    real(DP), intent(in) :: head
    real(DP), intent(inout) :: flow
    real(DP), intent(inout) :: source
    real(DP), intent(inout), optional :: gwfhcof
    real(DP), intent(inout), optional :: gwfrhs
    ! -- local
    real(DP) :: gwfhcof0, gwfrhs0
    !
    flow = dzero
    idry = 0
    call this%lak_calculate_conn_exchange(ilak, iconn, stage, head, flow, &
                                          gwfhcof0, gwfrhs0)
    if (iflag == 1) then
      if (flow > dzero) then
        source = source + flow
      end if
    else if (iflag == 2) then
      if (-flow > source) then
        flow = -source
        source = dzero
        idry = 1
      else if (flow < dzero) then
        source = source + flow
      end if
    end if
    !
    ! -- Set gwfhcof and gwfrhs if present
    if (present(gwfhcof)) gwfhcof = gwfhcof0
    if (present(gwfrhs)) gwfrhs = gwfrhs0
  end subroutine lak_estimate_conn_exchange
 
  !> @brief Calculate the storage change in a lake based on provided stages
  !! and a passed delt
  !<
  subroutine lak_calculate_storagechange(this, ilak, stage, stage0, delt, dvr)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(in) :: stage0
    real(DP), intent(in) :: delt
    real(DP), intent(inout) :: dvr
    ! -- local
    real(DP) :: v
    real(DP) :: v0
    !
    dvr = dzero
    if (this%gwfiss /= 1) then
      call this%lak_calculate_vol(ilak, stage, v)
      call this%lak_calculate_vol(ilak, stage0, v0)
      dvr = (v0 - v) / delt
    end if
  end subroutine lak_calculate_storagechange
 
  !> @brief Calculate the rainfall for a lake
  !<
  subroutine lak_calculate_rainfall(this, ilak, stage, ra)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(inout) :: ra
    ! -- local
    integer(I4B) :: iconn
    real(DP) :: sa
    !
    ! -- rainfall
    iconn = this%idxlakeconn(ilak)
    if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
      sa = this%sareamax(ilak)
    else
      call this%lak_calculate_sarea(ilak, stage, sa)
    end if
    ra = this%rainfall(ilak) * sa
  end subroutine lak_calculate_rainfall
 
  !> @brief Calculate runoff to a lake
  !<
  subroutine lak_calculate_runoff(this, ilak, ro)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: ro
    !
    ! -- runoff
    ro = this%runoff(ilak)
  end subroutine lak_calculate_runoff
 
  !> @brief Calculate specified inflow to a lake
  !<
  subroutine lak_calculate_inflow(this, ilak, qin)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: qin
    !
    ! -- inflow to lake
    qin = this%inflow(ilak)
  end subroutine lak_calculate_inflow
 
  !> @brief Calculate the external flow terms to a lake
  !<
  subroutine lak_calculate_external(this, ilak, ex)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: ex
    !
    ! -- If mover is active, add receiver water to rhs and
    !    store available water (as positive value)
    ex = dzero
    if (this%imover == 1) then
      ex = this%pakmvrobj%get_qfrommvr(ilak)
    end if
  end subroutine lak_calculate_external
 
  !> @brief Calculate the withdrawal from a lake subject to an available volume
  !<
  subroutine lak_calculate_withdrawal(this, ilak, avail, wr)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: avail
    real(DP), intent(inout) :: wr
    !
    ! -- withdrawals - limit to sum of inflows and available volume
    wr = this%withdrawal(ilak)
    if (wr > avail) then
      wr = -avail
    else
      if (wr > dzero) then
        wr = -wr
      end if
    end if
    avail = avail + wr
  end subroutine lak_calculate_withdrawal
 
  !> @brief Calculate the evaporation from a lake at a provided stage subject
  !! to an available volume
  !<
  subroutine lak_calculate_evaporation(this, ilak, stage, avail, ev)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(inout) :: avail
    real(DP), intent(inout) :: ev
    ! -- local
    real(DP) :: sa
    !
    ! -- evaporation - limit to sum of inflows and available volume
    call this%lak_calculate_sarea(ilak, stage, sa)
    ev = sa * this%evaporation(ilak)
    if (ev > avail) then
      if (is_close(avail, dprec)) then
        ev = dzero
      else
        ev = -avail
      end if
    else
      ev = -ev
    end if
    avail = avail + ev
  end subroutine lak_calculate_evaporation
 
  !> @brief Calculate the outlet inflow to a lake
  !<
  subroutine lak_calculate_outlet_inflow(this, ilak, outinf)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: outinf
    ! -- local
    integer(I4B) :: n
    !
    outinf = dzero
    do n = 1, this%noutlets
      if (this%lakeout(n) == ilak) then
        outinf = outinf - this%simoutrate(n)
        if (this%imover == 1) then
          outinf = outinf - this%pakmvrobj%get_qtomvr(n)
        end if
      end if
    end do
  end subroutine lak_calculate_outlet_inflow
 
  !> @brief Calculate the outlet outflow from a lake
  !<
  subroutine lak_calculate_outlet_outflow(this, ilak, stage, avail, outoutf)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(inout) :: avail
    real(DP), intent(inout) :: outoutf
    ! -- local
    integer(I4B) :: n
    real(DP) :: g
    real(DP) :: d
    real(DP) :: c
    real(DP) :: gsm
    real(DP) :: rate
    !
    outoutf = dzero
    do n = 1, this%noutlets
      if (this%lakein(n) == ilak) then
        rate = dzero
        d = stage - this%outinvert(n)
        if (this%outdmax > dzero) then
          if (d > this%outdmax) d = this%outdmax
        end if
        g = dgravity * this%convlength * this%convtime * this%convtime
        select case (this%iouttype(n))
          ! specified rate
        case (0)
          rate = this%outrate(n)
          if (-rate > avail) then
            rate = -avail
          end if
          ! manning
        case (1)
          if (d > dzero) then
            c = (this%convlength**donethird) * this%convtime
            gsm = dzero
            if (this%outrough(n) > dzero) then
              gsm = done / this%outrough(n)
            end if
            rate = -c * gsm * this%outwidth(n) * (d**dfivethirds) * &
                   sqrt(this%outslope(n))
          end if
          ! weir
        case (2)
          if (d > dzero) then
            rate = -dtwothirds * dcd * this%outwidth(n) * d * &
                   sqrt(dtwo * g * d)
          end if
        end select
        this%simoutrate(n) = rate
        avail = avail + rate
        outoutf = outoutf + rate
      end if
    end do
  end subroutine lak_calculate_outlet_outflow
 
  !> @brief Total uncapped outlet outflow rate from a lake at a provided stage
  !!
  !! Side-effect-free companion to lak_calculate_outlet_outflow used by the
  !! implicit formulation: returns the summed (negative) rating-curve outflow for
  !! all outlets whose source is ilak, without applying the available-water cap or
  !! updating simoutrate. Used to linearize the outlet-outflow sink on the lake
  !! row in stage.
  !<
  subroutine lak_outlet_outflow_rate(this, ilak, stage, qout)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: stage
    real(DP), intent(inout) :: qout
    ! -- local
    integer(I4B) :: n
    real(DP) :: g, d, c, gsm, rate
    !
    qout = dzero
    do n = 1, this%noutlets
      if (this%lakein(n) /= ilak) cycle
      rate = dzero
      d = stage - this%outinvert(n)
      if (this%outdmax > dzero .and. d > this%outdmax) d = this%outdmax
      g = dgravity * this%convlength * this%convtime * this%convtime
      select case (this%iouttype(n))
      case (0) ! specified rate
        rate = this%outrate(n)
      case (1) ! manning
        if (d > dzero) then
          c = (this%convlength**donethird) * this%convtime
          gsm = dzero
          if (this%outrough(n) > dzero) gsm = done / this%outrough(n)
          rate = -c * gsm * this%outwidth(n) * (d**dfivethirds) * &
                 sqrt(this%outslope(n))
        end if
      case (2) ! weir
        if (d > dzero) then
          rate = -dtwothirds * dcd * this%outwidth(n) * d * sqrt(dtwo * g * d)
        end if
      end select
      qout = qout + rate
    end do
  end subroutine lak_outlet_outflow_rate
 
  !> @brief Get the outlet inflow to a lake from another lake
  !<
  subroutine lak_get_internal_inlet(this, ilak, outinf)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: outinf
    ! -- local
    integer(I4B) :: n
    !
    outinf = dzero
    do n = 1, this%noutlets
      if (this%lakeout(n) == ilak) then
        outinf = outinf - this%simoutrate(n)
        if (this%imover == 1) then
          outinf = outinf - this%pakmvrobj%get_qtomvr(n)
        end if
      end if
    end do
  end subroutine lak_get_internal_inlet
 
  !> @brief Get the outlet from a lake to another lake
  !<
  subroutine lak_get_internal_outlet(this, ilak, outoutf)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: outoutf
    ! -- local
    integer(I4B) :: n
    !
    outoutf = dzero
    do n = 1, this%noutlets
      if (this%lakein(n) == ilak) then
        if (this%lakeout(n) < 1) cycle
        outoutf = outoutf + this%simoutrate(n)
      end if
    end do
  end subroutine lak_get_internal_outlet
 
  !> @brief Get the outlet outflow from a lake to an external boundary
  !<
  subroutine lak_get_external_outlet(this, ilak, outoutf)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: outoutf
    ! -- local
    integer(I4B) :: n
    !
    outoutf = dzero
    do n = 1, this%noutlets
      if (this%lakein(n) == ilak) then
        if (this%lakeout(n) > 0) cycle
        outoutf = outoutf + this%simoutrate(n)
      end if
    end do
  end subroutine lak_get_external_outlet
 
  !> @brief Get the mover outflow from a lake to an external boundary
  !<
  subroutine lak_get_external_mover(this, ilak, outoutf)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: outoutf
    ! -- local
    integer(I4B) :: n
    !
    outoutf = dzero
    if (this%imover == 1) then
      do n = 1, this%noutlets
        if (this%lakein(n) == ilak) then
          if (this%lakeout(n) > 0) cycle
          outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
        end if
      end do
    end if
  end subroutine lak_get_external_mover
 
  !> @brief Get the mover outflow from a lake to another lake
  !<
  subroutine lak_get_internal_mover(this, ilak, outoutf)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: outoutf
    ! -- local
    integer(I4B) :: n
    !
    outoutf = dzero
    if (this%imover == 1) then
      do n = 1, this%noutlets
        if (this%lakein(n) == ilak) then
          if (this%lakeout(n) < 1) cycle
          outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
        end if
      end do
    end if
  end subroutine lak_get_internal_mover
 
  !> @brief Get the outlet to mover from a lake
  !<
  subroutine lak_get_outlet_tomover(this, ilak, outoutf)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: outoutf
    ! -- local
    integer(I4B) :: n
    !
    outoutf = dzero
    if (this%imover == 1) then
      do n = 1, this%noutlets
        if (this%lakein(n) == ilak) then
          outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
        end if
      end do
    end if
  end subroutine lak_get_outlet_tomover
 
  !> @brief Determine the stage from a provided volume
  !<
  subroutine lak_vol2stage(this, ilak, vol, stage)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: vol
    real(DP), intent(inout) :: stage
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: ibs
    real(DP) :: s0, s1, sm
    real(DP) :: v0, v1, vm
    real(DP) :: f0, f1, fm
    real(DP) :: sa
    real(DP) :: en0, en1
    real(DP) :: ds, ds0
    real(DP) :: denom
    !
    s0 = this%lakebot(ilak)
    call this%lak_calculate_vol(ilak, s0, v0)
    s1 = this%laketop(ilak)
    call this%lak_calculate_vol(ilak, s1, v1)
    ! -- zero volume
    if (vol <= v0) then
      stage = s0
      ! -- linear relation between stage and volume above top of lake
    else if (vol >= v1) then
      call this%lak_calculate_sarea(ilak, s1, sa)
      stage = s1 + (vol - v1) / sa
      ! -- use combination of secant and bisection
    else
      en0 = s0
      en1 = s1
      ! sm = s1  ! causes divide by zero in 1st line in secantbisection loop
      ! sm = s0  ! causes divide by zero in 1st line in secantbisection loop
      sm = dzero
      f0 = vol - v0
      f1 = vol - v1
      ibs = 0
      secantbisection: do i = 1, 150
        denom = f1 - f0
        if (denom /= dzero) then
          ds = f1 * (s1 - s0) / denom
        else
          ibs = 13
        end if
        if (i == 1) then
          ds0 = ds
        end if
        ! -- use bisection if end points are exceeded
        if (sm < en0 .or. sm > en1) ibs = 13
        ! -- use bisection if secant method stagnates or if
        !    ds exceeds previous ds - bisection would occur
        !    after conditions exceeded in 13 iterations
        if (ds * ds0 < dprec .or. abs(ds) > abs(ds0)) ibs = ibs + 1
        if (ibs > 12) then
          ds = dhalf * (s1 - s0)
          ibs = 0
        end if
        sm = s1 - ds
        if (abs(ds) < dem6) then
          exit secantbisection
        end if
        call this%lak_calculate_vol(ilak, sm, vm)
        fm = vol - vm
        s0 = s1
        f0 = f1
        s1 = sm
        f1 = fm
        ds0 = ds
      end do secantbisection
      stage = sm
      if (abs(ds) >= dem6) then
        write (this%iout, '(1x,a,1x,i0,4(1x,a,1x,g15.6))') &
     &   'LAK_VOL2STAGE failed for lake', ilak, 'volume error =', fm, &
     &   'finding stage (', stage, ') for volume =', vol, &
     &    'final change in stage =', ds
      end if
    end if
  end subroutine lak_vol2stage
 
  !> @brief Determine if a valid lake or outlet number has been specified
  function lak_check_valid(this, itemno) result(ierr)
    ! -- modules
    use simmodule, only: store_error
    ! -- return
    integer(I4B) :: ierr
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: itemno
    ! -- local
    integer(I4B) :: ival
    !
    ierr = 0
    ival = abs(itemno)
    if (itemno > 0) then
      if (ival < 1 .or. ival > this%nlakes) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0,a)') &
          'LAKENO', itemno, 'must be greater than 0 and less than or equal to', &
          this%nlakes, '.'
        call store_error(errmsg)
        ierr = 1
      end if
    else
      if (ival < 1 .or. ival > this%noutlets) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0,a)') &
          'IOUTLET', itemno, 'must be greater than 0 and less than or equal to', &
          this%noutlets, '.'
        call store_error(errmsg)
        ierr = 1
      end if
    end if
  end function lak_check_valid
 
  !> @brief Apply a PERIOD lake setting and require its value be non-negative
  subroutine lak_set_nonneg_value(this, n, itemno, setting, value_array, label)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: itemno
    character(len=*), intent(in) :: setting
    real(DP), dimension(:), intent(inout) :: value_array
    character(len=*), intent(in) :: label
    ! -- local
    integer(I4B) :: istat
    !
    istat = this%lak_check_valid(itemno)
    if (istat /= 0) return
    if (itemno < 1 .or. itemno > this%nlakes) return
    if (value_array(itemno) < dzero) then
      write (errmsg, '(a,i0,a,g0,a)') &
        'Lake ', itemno, ' was assigned a '//trim(label)//' value of ', &
        value_array(itemno), '. Value must be positive.'
      call store_error(errmsg)
    end if
  end subroutine lak_set_nonneg_value
 
  !> @brief Issue a parameter error for lake ilak
  !<
  subroutine lak_set_attribute_error(this, ilak, keyword, msg)
    ! -- modules
    use simmodule, only: store_error
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    character(len=*), intent(in) :: keyword
    character(len=*), intent(in) :: msg
    !
    if (len(msg) == 0) then
      write (errmsg, '(a,1x,a,1x,i0,1x,a)') &
        keyword, ' for LAKE', ilak, 'has already been set.'
    else
      write (errmsg, '(a,1x,a,1x,i0,1x,a)') keyword, ' for LAKE', ilak, msg
    end if
    call store_error(errmsg)
  end subroutine lak_set_attribute_error
 
  !> @brief Source OPTIONS block from input context
  !<
  subroutine lak_source_options(this)
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    type(gwflakparamfoundtype) :: found
    integer(I4B) :: isize
    character(len=LINELENGTH) :: fname
    ! -- formats
    character(len=*), parameter :: fmtlengthconv = &
      &"(4x, 'LENGTH CONVERSION VALUE (',g15.7,') SPECIFIED.')"
    character(len=*), parameter :: fmttimeconv = &
      &"(4x, 'TIME CONVERSION VALUE (',g15.7,') SPECIFIED.')"
    character(len=*), parameter :: fmtlakeopt = &
      &"(4x, 'LAKE ', a, ' VALUE (',g15.7,') SPECIFIED.')"
    character(len=*), parameter :: fmtoutdmax = &
      &"(4x, 'MAXIMUM OUTLET WATER DEPTH (',g15.7,') SPECIFIED.')"
    character(len=*), parameter :: fmtlakbin = &
      "(4x, 'LAK ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', &
      &a, /4x, 'OPENED ON UNIT: ', I0)"
    character(len=*), parameter :: fmtiter = &
      &"(4x, 'MAXIMUM LAK ITERATION VALUE (',i0,') SPECIFIED.')"
    character(len=*), parameter :: fmtdmaxchg = &
      &"(4x, 'MAXIMUM STAGE CHANGE VALUE (',g0,') SPECIFIED.')"
    !
    ! -- call base class to handle AUXILIARY/NAUX, BOUNDNAMES, PRINT_INPUT,
    !    PRINT_FLOWS, SAVE_FLOWS, OBS6_FILENAME
    call this%BndExtType%source_options()
    !
    write (this%iout, '(1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' OPTIONS'
    !
    ! -- PRINT_STAGE
    call mem_set_value(this%iprhed, 'PRINT_STAGE', this%input_mempath, &
                       found%print_stage)
    if (found%print_stage) then
      write (this%iout, '(4x,a)') trim(adjustl(this%text))// &
        ' STAGES WILL BE PRINTED TO LISTING FILE.'
    end if
    !
    ! -- STAGE FILEOUT
    call mem_set_value(fname, 'STAGEFILE', this%input_mempath, found%stagefile)
    if (found%stagefile) then
      this%istageout = getunit()
      call openfile(this%istageout, this%iout, fname, 'DATA(BINARY)', &
                    form, access, 'REPLACE', mode_opt=mnormal)
      write (this%iout, fmtlakbin) 'STAGE', trim(adjustl(fname)), this%istageout
    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, fmtlakbin) '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, fmtlakbin) 'BUDGET CSV', trim(adjustl(fname)), &
        this%ibudcsv
    end if
    !
    ! -- PACKAGE_CONVERGENCE FILEOUT (mf6internal: pkgconvfname)
    ! -- defer opening until lak_ar, where the IMPLICIT option is known. The
    !    file is not opened for the implicit formulation (lak_cc does not
    !    write package convergence in that case).
    call mem_set_value(fname, 'PKGCONVFNAME', this%input_mempath, &
                       found%pkgconvfname)
    if (found%pkgconvfname) then
      this%pakcsvfile = trim(adjustl(fname))
    end if
    !
    ! -- MOVER
    call mem_set_value(this%imover, 'MOVER', this%input_mempath, found%mover)
    if (found%mover) then
      write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
    end if
    !
    ! -- LENGTH_CONVERSION (mf6internal: convlength)
    call mem_set_value(this%convlength, 'CONVLENGTH', this%input_mempath, &
                       found%convlength)
    if (found%convlength) write (this%iout, fmtlengthconv) this%convlength
    !
    ! -- TIME_CONVERSION (mf6internal: convtime)
    call mem_set_value(this%convtime, 'CONVTIME', this%input_mempath, &
                       found%convtime)
    if (found%convtime) write (this%iout, fmttimeconv) this%convtime
    !
    ! -- SURFDEP (mf6internal: surfdep)
    call mem_set_value(this%surfdep, 'SURFDEP', this%input_mempath, found%surfdep)
    if (found%surfdep) then
      if (this%surfdep < dzero) this%surfdep = dzero
      write (this%iout, fmtlakeopt) 'SURFDEP', this%surfdep
    end if
    !
    ! -- MAXIMUM_ITERATIONS (mf6internal: maxlakit)
    call mem_set_value(this%maxlakit, 'MAXLAKIT', this%input_mempath, &
                       found%maxlakit)
    if (found%maxlakit) write (this%iout, fmtiter) this%maxlakit
    !
    ! -- MAXIMUM_STAGE_CHANGE (mf6internal: dmaxchg)
    call mem_set_value(this%dmaxchg, 'DMAXCHG', this%input_mempath, found%dmaxchg)
    if (found%dmaxchg) then
      this%delh = dp999 * this%dmaxchg
      write (this%iout, fmtdmaxchg) this%dmaxchg
    end if
    !
    ! -- DEV_GROUNDWATER_HEAD_CONDUCTANCE (mf6internal: igwhcopt)
    call mem_set_value(this%igwhcopt, 'IGWHCOPT', this%input_mempath, &
                       found%igwhcopt)
    if (found%igwhcopt) then
      write (this%iout, '(4x,a)') &
        'CONDUCTANCE FOR HORIZONTAL CONNECTIONS WILL BE CALCULATED &
        &USING THE GROUNDWATER HEAD'
    end if
    !
    ! -- IMPLICIT (mf6internal: iimplicit)
    call mem_set_value(this%iimplicit, 'IIMPLICIT', this%input_mempath, &
                       found%iimplicit)
    if (found%iimplicit) then
      write (this%iout, '(4x,a)') &
        'LAKE STAGE WILL BE SOLVED AS AN UNKNOWN IN THE GROUNDWATER FLOW '// &
        'MATRIX (IMPLICIT FORMULATION)'
    end if
    !
    ! -- DEV_FORCE_LEGACY (mf6internal: iforceleg)
    call mem_set_value(this%iforceleg, 'IFORCELEG', this%input_mempath, &
                       found%iforceleg)
    !
    ! -- DEV_FORCE_LEGACY_LAKE (mf6internal: iforceleglak); loaded before
    !    either block below so both mutual-exclusion checks see a final value
    call mem_set_value(this%iforceleglak, 'IFORCELEGLAK', this%input_mempath, &
                       found%iforceleglak)
    !
    if (found%iforceleg) then
      write (this%iout, '(4x,a)') &
        'EVERY ACTIVE LAKE WILL BE SOLVED WITH THE LEGACY SUBSTITUTION '// &
        'SOLVER UNDER THE IMPLICIT FORMULATION'
      if (this%iforceleglak /= 0) then
        call store_error('DEV_FORCE_LEGACY and DEV_FORCE_LEGACY_LAKE '// &
                         'cannot both be specified.')
      end if
    end if
    !
    if (found%iforceleglak) then
      write (this%iout, '(4x,a,i0,a)') 'LAKE ', this%iforceleglak, &
        ' WILL BE SOLVED WITH THE LEGACY SUBSTITUTION SOLVER UNDER THE '// &
        'IMPLICIT FORMULATION'
      if (this%iforceleg /= 0) then
        call store_error('DEV_FORCE_LEGACY and DEV_FORCE_LEGACY_LAKE '// &
                         'cannot both be specified.')
      end if
    end if
    !
    ! -- DEV_NO_FINAL_CHECK (mf6internal: iconvchk)
    call get_isize('ICONVCHK', this%input_mempath, isize)
    if (isize > 0) then
      this%iconvchk = 0
      write (this%iout, '(4x,a)') &
        'A FINAL CONVERGENCE CHECK OF THE CHANGE IN LAKE STAGES &
        &WILL NOT BE MADE'
    end if
    !
    ! -- DEV_MAXIMUM_OUTLET_DEPTH (mf6internal: outdmax)
    call mem_set_value(this%outdmax, 'OUTDMAX', this%input_mempath, found%outdmax)
    if (found%outdmax) write (this%iout, fmtoutdmax) this%outdmax
    !
    write (this%iout, '(1x,a)') 'END OF '//trim(adjustl(this%text))//' OPTIONS'
  end subroutine lak_source_options
 
  !> @brief Bind LakInputType pointers.
  !!
  !! Called once from lak_ar; pointers remain valid for the simulation lifetime.
  !! Input tsmanager advances TS-linked values transparently each timestep.
  !<
  subroutine lak_input_init(this, mempath, naux, noutlets)
    use memorymanagermodule, only: mem_setptr
    ! -- dummy
    class(lakinputtype), intent(inout) :: this
    character(len=*), intent(in) :: mempath
    integer(I4B), intent(in) :: naux
    integer(I4B), intent(in) :: noutlets
    !
    ! -- period
    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)
    if (naux > 0) then
      call mem_setptr(this%auxname, 'AUXNAME', mempath)
      call mem_setptr(this%auxval, 'AUXVAL', mempath)
    end if
    ! -- outlets
    if (noutlets > 0) then
      call mem_setptr(this%outletno, 'OUTLETS_OUTLETNO', mempath)
      call mem_setptr(this%outlets_invert, 'OUTLETS_INVERT', mempath)
      call mem_setptr(this%outlets_width, 'OUTLETS_WIDTH', mempath)
      call mem_setptr(this%outlets_rough, 'OUTLETS_ROUGH', mempath)
      call mem_setptr(this%outlets_slope, 'OUTLETS_SLOPE', mempath)
    end if
  end subroutine lak_input_init
 
  !> @brief Nullify all LakInputType pointers.
  !<
  subroutine lak_input_destroy(this)
    ! -- dummy
    class(lakinputtype), intent(inout) :: this
    !
    nullify (this%nbound)
    nullify (this%ifno)
    nullify (this%setting)
    nullify (this%status)
    nullify (this%auxname)
    nullify (this%auxval)
    nullify (this%outletno)
    nullify (this%outlets_invert)
    nullify (this%outlets_width)
    nullify (this%outlets_rough)
    nullify (this%outlets_slope)
  end subroutine lak_input_destroy
 
  !> @brief Allocate and Read
  !!
  !! Create new LAK package and point bndobj to the new package
  !<
  subroutine lak_ar(this)
    ! -- modules
    use constantsmodule, only: mnormal
    use inputoutputmodule, only: getunit, openfile
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- formats
    character(len=*), parameter :: fmtlakbin = &
      "(4x, 'LAK ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', &
      &a, /4x, 'OPENED ON UNIT: ', I0)"
    !
    ! -- open the deferred PACKAGE_CONVERGENCE file now that all options have
    !    been read. lak_cc does not write package convergence for the implicit
    !    formulation, because the lake stage is part of the solver (IMS)
    !    convergence check. With the IMPLICIT option the file is therefore left
    !    unopened (ipakcsv stays 0) and a warning is issued instead.
    if (allocated(this%pakcsvfile)) then
      if (this%iimplicit /= 0) then
        write (warnmsg, '(a)') &
          'PACKAGE_CONVERGENCE output file "'//trim(this%pakcsvfile)// &
          '" is not written when the IMPLICIT option is active; the lake '// &
          'stage is part of the solver (IMS) convergence check.'
        call store_warning(warnmsg)
      else
        this%ipakcsv = getunit()
        call openfile(this%ipakcsv, this%iout, this%pakcsvfile, 'CSV', &
                      filstat_opt='REPLACE', mode_opt=mnormal)
        write (this%iout, fmtlakbin) 'PACKAGE_CONVERGENCE', &
          trim(this%pakcsvfile), this%ipakcsv
      end if
      deallocate (this%pakcsvfile)
    end if
    !
    call this%obs%obs_ar()
    !
    ! -- Allocate arrays in LAK and in package superclass
    call this%lak_allocate_arrays()
    !
    ! -- read optional initial package parameters
    call this%read_initial_attr()
    !
    ! -- bind input context persistent pointers
    call this%input%init(this%input_mempath, this%naux, this%noutlets)
    !
    ! -- setup pakmvrobj
    if (this%imover /= 0) then
      allocate (this%pakmvrobj)
      call this%pakmvrobj%ar(this%noutlets, this%nlakes, this%memoryPath)
    end if
  end subroutine lak_ar
 
  !> @brief Read and prepare (stress period data)
  !<
  subroutine lak_rp(this)
    ! -- modules
    use constantsmodule, only: lenvarname
    use tdismodule, only: kper
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    character(len=LINELENGTH) :: str = ' '
    real(DP) :: auxvaldbl = dzero
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    character(len=LENVARNAME) :: setting
    integer(I4B) :: node, n, j, itemno
    integer(I4B) :: istat
    ! -- formats
    character(len=*), parameter :: fmtlsp = &
      &"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
    !
    ! -- set nbound to maxbound
    this%nbound = this%maxbound
    !
    ! -- check if data has been loaded for this stress period
    if (this%iper == kper) then
      !
      this%check_attr = 1
      !
      ! -- 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, 4, 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 n = 1, 2
          write (text, '(a,1x,i6)') 'VALUE', n
          call this%inputtab%initialize_column(text, 15, alignment=tabcenter)
        end do
      end if
      !
      ! -- process each period data entry
      do n = 1, this%input%nbound
        itemno = this%input%ifno(n)
        !
        ! -- dispatch key for this row
        setting = this%input%setting(n)
        str = ' '
        !
        ! -- STATUS (lake only; not TS-enabled)
        if (trim(setting) == 'STATUS') then
          str = this%input%status(n)
          istat = this%lak_check_valid(itemno)
          if (istat == 0 .and. itemno >= 1 .and. itemno <= this%nlakes) then
            if (trim(str) == 'CONSTANT') then
              this%iboundpak(itemno) = -1
            else if (trim(str) == 'INACTIVE') then
              this%iboundpak(itemno) = 0
            else if (trim(str) == 'ACTIVE') then
              this%iboundpak(itemno) = 1
            else
              write (errmsg, '(a,a)') &
                'Unknown '//trim(this%text)//' STATUS keyword: ', &
                trim(str)//'.'
              call store_error(errmsg)
            end if
            this%status(itemno) = trim(str)
          end if
        end if
        !
        ! -- STAGE (validate itemno; assignment via helper)
        if (trim(setting) == 'STAGE') then
          istat = this%lak_check_valid(itemno)
        end if
        !
        ! -- RAINFALL, EVAPORATION, RUNOFF, INFLOW, WITHDRAWAL (require
        !    non-negative)
        if (trim(setting) == 'RAINFALL') &
          call this%lak_set_nonneg_value(n, itemno, setting, this%rainfall, &
                                         'rainfall')
        if (trim(setting) == 'EVAPORATION') &
          call this%lak_set_nonneg_value(n, itemno, setting, &
                                         this%evaporation, 'evaporation')
        if (trim(setting) == 'RUNOFF') &
          call this%lak_set_nonneg_value(n, itemno, setting, this%runoff, &
                                         'runoff')
        if (trim(setting) == 'INFLOW') &
          call this%lak_set_nonneg_value(n, itemno, setting, this%inflow, &
                                         'inflow')
        if (trim(setting) == 'WITHDRAWAL') &
          call this%lak_set_nonneg_value(n, itemno, setting, &
                                         this%withdrawal, 'withdrawal')
        !
        ! -- RATE, INVERT, WIDTH, ROUGH, SLOPE (outlet; validate with
        !    negative sign convention)
        if (trim(setting) == 'RATE' .or. trim(setting) == 'INVERT' .or. &
            trim(setting) == 'WIDTH' .or. trim(setting) == 'ROUGH' .or. &
            trim(setting) == 'SLOPE') then
          istat = this%lak_check_valid(-itemno)
        end if
        !
        ! -- INVERT, WIDTH, ROUGH, SLOPE latch here so the OUTLETS block's
        !    own (separately time-series-linked) baseline resync in lak_ad
        !    stops overwriting this outlet's PERIOD-block override
        if (istat == 0 .and. itemno >= 1 .and. itemno <= this%noutlets) then
          if (trim(setting) == 'INVERT') this%outlet_invert_set(itemno) = .true.
          if (trim(setting) == 'WIDTH') this%outlet_width_set(itemno) = .true.
          if (trim(setting) == 'ROUGH') this%outlet_rough_set(itemno) = .true.
          if (trim(setting) == 'SLOPE') this%outlet_slope_set(itemno) = .true.
        end if
        !
        ! -- AUXILIARY; featureauxvar set via sync_auxvar()
        if (this%naux > 0) then
          if (trim(setting) == 'PERIOD_AUXILIARY') then
            if (itemno >= 1 .and. itemno <= this%nlakes) then
              str = this%input%auxname(n)
              auxvaldbl = this%input%auxval(n)
            end if
          end if
        end if
        !
        ! -- echo row to period data table
        if (this%iprpak /= 0) then
          call this%inputtab%add_term(itemno)
          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(' ')
          case ('STAGE')
            call this%inputtab%add_term(this%stage(itemno))
            call this%inputtab%add_term(' ')
          case ('RAINFALL')
            call this%inputtab%add_term(this%rainfall(itemno))
            call this%inputtab%add_term(' ')
          case ('EVAPORATION')
            call this%inputtab%add_term(this%evaporation(itemno))
            call this%inputtab%add_term(' ')
          case ('RUNOFF')
            call this%inputtab%add_term(this%runoff(itemno))
            call this%inputtab%add_term(' ')
          case ('INFLOW')
            call this%inputtab%add_term(this%inflow(itemno))
            call this%inputtab%add_term(' ')
          case ('WITHDRAWAL')
            call this%inputtab%add_term(this%withdrawal(itemno))
            call this%inputtab%add_term(' ')
          case ('RATE')
            call this%inputtab%add_term(this%outrate(itemno))
            call this%inputtab%add_term(' ')
          case ('INVERT')
            call this%inputtab%add_term(this%outinvert(itemno))
            call this%inputtab%add_term(' ')
          case ('WIDTH')
            call this%inputtab%add_term(this%outwidth(itemno))
            call this%inputtab%add_term(' ')
          case ('ROUGH')
            call this%inputtab%add_term(this%outrough(itemno))
            call this%inputtab%add_term(' ')
          case ('SLOPE')
            call this%inputtab%add_term(this%outslope(itemno))
            call this%inputtab%add_term(' ')
          case ('PERIOD_AUXILIARY')
            call this%inputtab%add_term(trim(str))
            call this%inputtab%add_term(auxvaldbl)
          case default
            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
      !
      ! -- write summary of lake stress period error messages
      if (count_errors() > 0) then
        call store_error_filename(this%input_fname)
      end if
      !
      ! -- fill bound array with lake stage, conductance, and bottom elevation
      do n = 1, this%nlakes
        do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
          node = this%gwfconn(j)
          this%nodelist(j) = node
          this%bound(1, j) = this%xnewpak(n)
          this%bound(2, j) = this%satcond(j)
          this%bound(3, j) = this%belev(j)
        end do
      end do
      !
      ! -- copy lakein into iprmap so mvr budget contains lake instead of outlet
      if (this%imover == 1) then
        do n = 1, this%noutlets
          this%pakmvrobj%iprmap(n) = this%lakein(n)
        end do
      end if
    else
      write (this%iout, fmtlsp) trim(this%filtyp)
    end if
    !
    ! -- fill nodelist and bound(3) every period; nodelist for non-vertical
    !    connections is not updated by lak_cf so must be set here.
    !    bound(3) (bed elevation) is static and only needs to be set once.
    do n = 1, this%nlakes
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        this%nodelist(j) = this%gwfconn(j)
        this%bound(3, j) = this%belev(j)
      end do
    end do
  end subroutine lak_rp
 
  !> @brief Add package connection to matrix
  !<
  subroutine lak_ad(this)
    ! -- modules
    use simvariablesmodule, only: ifailedstepretry
    ! -- dummy
    class(laktype) :: this
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: j
    integer(I4B) :: iaux
    integer(I4B) :: itemno
    !
    ! -- sync PACKAGEDATA AUX TS and advance observations;
    !    sync_period_auxvar (inside bnd_ad) handles PERIOD AUXILIARY TS
    call this%BndExtType%bnd_ad()
    !
    ! -- re-sync static OUTLETS TS values (invert, width, rough, slope);
    ! -- skip any field a PERIOD-block setting has already overridden for
    !    this outlet -- that override must persist until a later PERIOD
    !    block changes it, not be overwritten by this baseline resync
    if (this%ts_active .and. this%noutlets > 0) then
      do n = 1, size(this%input%outletno)
        itemno = this%input%outletno(n)
        if (itemno < 1 .or. itemno > this%noutlets) cycle
        if (.not. this%outlet_invert_set(itemno)) &
          this%outinvert(itemno) = this%input%outlets_invert(n)
        if (.not. this%outlet_width_set(itemno)) &
          this%outwidth(itemno) = this%input%outlets_width(n)
        if (.not. this%outlet_rough_set(itemno)) &
          this%outrough(itemno) = this%input%outlets_rough(n)
        if (.not. this%outlet_slope_set(itemno)) &
          this%outslope(itemno) = this%input%outlets_slope(n)
      end do
    end if
    !
    ! -- update auxiliary variables
    if (this%naux > 0) then
      do n = 1, this%nlakes
        do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
          do iaux = 1, this%naux
            if (this%noupdateauxvar(iaux) /= 0) cycle
            this%auxvar(iaux, j) = this%featureauxvar(iaux, n)
          end do
        end do
      end do
    end if
    !
    ! -- Update or restore state
    if (ifailedstepretry == 0) then
      !
      ! -- copy xnew into xold and set xnewpak to stage for
      !    constant stage lakes
      do n = 1, this%nlakes
        this%xoldpak(n) = this%xnewpak(n)
        this%stageiter(n) = this%xnewpak(n)
        if (this%iboundpak(n) < 0) then
          this%xnewpak(n) = this%stage(n)
        end if
        this%seep0(n) = dzero
      end do
    else
      !
      ! -- copy xold back into xnew as this is a
      !    retry of this time step
      do n = 1, this%nlakes
        this%xnewpak(n) = this%xoldpak(n)
        this%stageiter(n) = this%xnewpak(n)
        if (this%iboundpak(n) < 0) then
          this%xnewpak(n) = this%stage(n)
        end if
        this%seep0(n) = dzero
      end do
    end if
    !
    ! -- pakmvrobj ad
    if (this%imover == 1) then
      call this%pakmvrobj%ad()
    end if
  end subroutine lak_ad
 
  !> @brief Formulate the HCOF and RHS terms
  !!
  !! Skip if no lakes, otherwise calculate hcof and rhs
  !<
  subroutine lak_cf(this)
    ! -- dummy
    class(laktype) :: this
    ! -- local
    integer(I4B) :: j, n
    integer(I4B) :: igwfnode
    real(DP) :: hlak, bottom_lake
    !
    ! -- save groundwater seepage for lake solution
    do n = 1, this%nlakes
      this%seep0(n) = this%seep(n)
    end do
    !
    ! -- save variables for convergence check
    do n = 1, this%nlakes
      this%s0(n) = this%xnewpak(n)
      call this%lak_calculate_exchange(n, this%s0(n), this%qgwf0(n))
    end do
    !
    ! -- find highest active cell
    do n = 1, this%nlakes
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        ! -- skip horizontal connections
        if (this%ictype(j) /= 0) then
          cycle
        end if
        igwfnode = this%nodesontop(j)
        if (this%ibound(igwfnode) == 0) then
          call this%dis%highest_active(igwfnode, this%ibound)
        end if
        this%nodelist(j) = igwfnode
        this%gwfconn(j) = igwfnode
      end do
    end do
    !
    ! -- reset ibound for cells where lake stage is above the bottom
    !    of the lake in the cell or the lake is inactive - only applied to
    !    vertical connections
    do n = 1, this%nlakes
      !
      hlak = this%xnewpak(n)
      !
      ! -- Go through lake connections
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        !
        ! -- assign gwf node number
        igwfnode = this%gwfconn(j)
        !
        ! -- skip inactive or constant head GWF cells
        if (this%ibound(igwfnode) < 1) then
          cycle
        end if
        !
        ! -- skip horizontal connections
        if (this%ictype(j) /= 0) then
          cycle
        end if
        !
        ! -- skip embedded lakes
        if (this%ictype(j) == 2 .or. this%ictype(j) == 3) then
          cycle
        end if
        !
        ! -- Mark ibound for wet lakes or inactive lakes; reset to 1 otherwise
        bottom_lake = this%belev(j)
        if (hlak > bottom_lake .or. this%iboundpak(n) == 0) then
          this%ibound(igwfnode) = iwetlake
        else
          this%ibound(igwfnode) = 1
        end if
      end do
      !
    end do
    !
    ! -- Store the lake stage and cond in bound array for other
    !    packages, such as the BUY package
    call this%lak_bound_update()
  end subroutine lak_cf
 
  !> @brief Copy rhs and hcof into solution rhs and amat
  !<
  subroutine lak_fc(this, rhs, ia, idxglo, matrix_sln)
    ! -- dummy
    class(laktype) :: 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, n
    integer(I4B) :: igwfnode
    integer(I4B) :: ipossymd
    !
    ! -- pakmvrobj fc
    if (this%imover == 1) then
      call this%pakmvrobj%fc()
    end if
    !
    ! -- implicit formulation: assemble the lake equations directly into the
    !    groundwater flow matrix instead of solving the stage by substitution
    if (this%iimplicit /= 0) then
      call this%lak_fc_implicit(rhs, matrix_sln)
      return
    end if
    !
    ! -- legacy formulation: solve the lake stage by substitution, then add the
    !    resulting lake-aquifer exchange terms to the groundwater flow matrix
    call this%lak_solve()
    do n = 1, this%nlakes
      if (this%iboundpak(n) == 0) cycle
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        igwfnode = this%gwfconn(j)
        if (this%ibound(igwfnode) < 1) cycle
        ipossymd = idxglo(ia(igwfnode))
        call matrix_sln%add_value_pos(ipossymd, this%hcof(j))
        rhs(igwfnode) = rhs(igwfnode) + this%rhs(j)
      end do
    end do
  end subroutine lak_fc
 
  !> @brief Fill newton terms
  !<
  subroutine lak_fn(this, rhs, ia, idxglo, matrix_sln)
    ! -- dummy
    class(laktype) :: 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, n
    integer(I4B) :: ipos
    integer(I4B) :: igwfnode
    integer(I4B) :: idry
    real(DP) :: hlak
    real(DP) :: avail
    real(DP) :: ra
    real(DP) :: ro
    real(DP) :: qinf
    real(DP) :: ex
    real(DP) :: head
    real(DP) :: q
    real(DP) :: q1
    real(DP) :: rterm
    real(DP) :: drterm
    !
    ! -- implicit formulation: the lakebed seepage Jacobian is already
    !    assembled exactly in lak_fc_implicit for constant-conductance
    !    (confined / saturated) connections. Nonlinear conductance Newton terms
    !    are Phase 2.
    if (this%iimplicit /= 0) then
      return
    end if
    !
    do n = 1, this%nlakes
      if (this%iboundpak(n) == 0) cycle
      hlak = this%xnewpak(n)
      call this%lak_calculate_available(n, hlak, avail, &
                                        ra, ro, qinf, ex, this%delh)
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        igwfnode = this%gwfconn(j)
        ipos = ia(igwfnode)
        head = this%xnew(igwfnode)
        if (-this%hcof(j) > dzero) then
          if (this%ibound(igwfnode) > 0) then
            ! -- estimate lake-aquifer exchange with perturbed groundwater head
            !    exchange is relative to the lake
            !avail = DEP20
            call this%lak_estimate_conn_exchange(2, n, j, idry, hlak, &
                                                 head + this%delh, q1, avail)
            q1 = -q1
            ! -- calculate unperturbed lake-aquifer exchange
            q = this%hcof(j) * head - this%rhs(j)
            ! -- calculate rterm
            rterm = this%hcof(j) * head
            ! -- calculate derivative
            drterm = (q1 - q) / this%delh
            ! -- add terms to convert conductance formulation into
            !    newton-raphson formulation
            call matrix_sln%add_value_pos(idxglo(ipos), drterm - this%hcof(j))
            rhs(igwfnode) = rhs(igwfnode) - rterm + drterm * head
          end if
        end if
      end do
    end do
  end subroutine lak_fn
 
  !> @brief Apply Newton under-relaxation to the lake stage
  !!
  !! Implicit formulation only. The model Newton under-relaxation step (gwf_nur)
  !! calls this routine with the lake-stage part of the solution. If an updated
  !! stage falls below the lake bottom, it is relaxed back toward the bottom (90
  !! percent of the way) and its change is zeroed. Active only when the model
  !! uses NEWTON UNDER_RELAXATION.
  !<
  subroutine lak_nur(this, neqpak, x, xtemp, dx, inewtonur, dxmax, locmax)
    ! -- dummy
    class(laktype), 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) :: botl
    real(DP) :: xx
    real(DP) :: dxx
    !
    ! -- only the implicit formulation has the lake stage in the global solution
    if (this%iimplicit == 0) return
    !
    ! -- Newton-Raphson under-relaxation: hold the stage at the lake bottom
    do n = 1, this%nlakes
      if (this%iboundpak(n) < 1) cycle
      botl = this%lakebot(n)
      !
      ! -- only apply under-relaxation if the updated stage is below the bottom
      !    of the lake
      if (x(n) < botl) then
        inewtonur = 1
        xx = xtemp(n) * (done - dp9) + botl * 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
    end do
  end subroutine lak_nur
 
  !> @brief Final convergence check for package
  !<
  subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
    ! -- modules
    use tdismodule, only: totim, kstp, kper, delt
    ! -- dummy
    class(laktype), 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
    character(len=LENPAKLOC) :: cloc
    character(len=LINELENGTH) :: tag
    integer(I4B) :: icheck
    integer(I4B) :: ipakfail
    integer(I4B) :: locdhmax
    integer(I4B) :: locresidmax
    integer(I4B) :: locdgwfmax
    integer(I4B) :: locdqoutmax
    integer(I4B) :: locdqfrommvrmax
    integer(I4B) :: ntabrows
    integer(I4B) :: ntabcols
    integer(I4B) :: n
    real(DP) :: q
    real(DP) :: q0
    real(DP) :: qtolfact
    real(DP) :: area
    real(DP) :: gwf0
    real(DP) :: gwf
    real(DP) :: dh
    real(DP) :: resid
    real(DP) :: dgwf
    real(DP) :: hlak0
    real(DP) :: hlak
    real(DP) :: qout0
    real(DP) :: qout
    real(DP) :: dqout
    real(DP) :: inf
    real(DP) :: ra
    real(DP) :: ro
    real(DP) :: qinf
    real(DP) :: ex
    real(DP) :: dhmax
    real(DP) :: residmax
    real(DP) :: dgwfmax
    real(DP) :: dqoutmax
    real(DP) :: dqfrommvr
    real(DP) :: dqfrommvrmax
    ! -- switch any stalled IMPLICIT lake back to the legacy solver
    call this%lak_set_legacy(kiter, icnvgmod)
    !
    ! -- on the last outer iteration of a solution that has not converged, warn
    !    about a perched (disconnected) lake that has no practical steady state.
    !    Done for both formulations (the implicit path returns just below).
    if (iend /= 0 .and. icnvgmod == 0) then
      call this%lak_check_disconnected()
    end if
    !
    ! -- implicit formulation: lake stage is solved in the global matrix,
    !    so its convergence is governed by the solution's dvclose check on the
    !    stage unknown; no separate package convergence check is performed here.
    if (this%iimplicit /= 0) then
      return
    end if
    !
    ! -- initialize local variables
    icheck = this%iconvchk
    ipakfail = 0
    locdhmax = 0
    locresidmax = 0
    locdgwfmax = 0
    locdqoutmax = 0
    locdqfrommvrmax = 0
    dhmax = dzero
    residmax = dzero
    dgwfmax = dzero
    dqoutmax = dzero
    dqfrommvrmax = dzero
    !
    ! -- if not saving package convergence data on check convergence if
    !    the model is considered converged
    if (this%ipakcsv == 0) then
      if (icnvgmod == 0) then
        icheck = 0
      end if
      !
      ! -- saving package convergence data
    else
      !
      ! -- header for package csv
      if (.not. associated(this%pakcsvtab)) then
        !
        ! -- determine the number of columns and rows
        ntabrows = 1
        ntabcols = 11
        if (this%noutlets > 0) then
          ntabcols = ntabcols + 2
        end if
        if (this%imover == 1) then
          ntabcols = ntabcols + 2
        end if
        !
        ! -- setup table
        call table_cr(this%pakcsvtab, this%packName, '')
        call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, &
                                     lineseparator=.false., separator=',', &
                                     finalize=.false.)
        !
        ! -- add columns to package csv
        tag = 'total_inner_iterations'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'totim'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'kper'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'kstp'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'nouter'
        call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
        tag = 'dvmax'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        tag = 'dvmax_loc'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        tag = 'residmax'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        tag = 'residmax_loc'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        tag = 'dgwfmax'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        tag = 'dgwfmax_loc'
        call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        if (this%noutlets > 0) then
          tag = 'dqoutmax'
          call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
          tag = 'dqoutmax_loc'
          call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
        end if
        if (this%imover == 1) then
          tag = 'dqfrommvrmax'
          call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
          tag = 'dqfrommvrmax_loc'
          call this%pakcsvtab%initialize_column(tag, 16, alignment=tableft)
        end if
      end if
    end if
    !
    ! -- perform package convergence check
    if (icheck /= 0) then
      final_check: do n = 1, this%nlakes
        if (this%iboundpak(n) < 1) cycle
        !
        ! -- set previous and current lake stage
        hlak0 = this%s0(n)
        hlak = this%xnewpak(n)
        !
        ! -- stage difference
        dh = hlak0 - hlak
        !
        ! -- calculate surface area
        call this%lak_calculate_sarea(n, hlak, area)
        !
        ! -- set the Q to length factor
        if (area > dzero) then
          qtolfact = delt / area
        else
          qtolfact = dzero
        end if
        !
        ! -- difference in the residual
        call this%lak_calculate_residual(n, hlak, resid)
        resid = resid * qtolfact
        !
        ! -- change in gwf exchange
        dgwf = dzero
        if (area > dzero) then
          gwf0 = this%qgwf0(n)
          call this%lak_calculate_exchange(n, hlak, gwf)
          dgwf = (gwf0 - gwf) * qtolfact
        end if
        !
        ! -- change in outflows
        dqout = dzero
        if (this%noutlets > 0) then
          if (area > dzero) then
            call this%lak_calculate_available(n, hlak0, inf, ra, ro, qinf, ex)
            call this%lak_calculate_outlet_outflow(n, hlak0, inf, qout0)
            call this%lak_calculate_available(n, hlak, inf, ra, ro, qinf, ex)
            call this%lak_calculate_outlet_outflow(n, hlak, inf, qout)
            dqout = (qout0 - qout) * qtolfact
          end if
        end if
        !
        ! -- q from mvr
        dqfrommvr = dzero
        if (this%imover == 1) then
          q = this%pakmvrobj%get_qfrommvr(n)
          q0 = this%pakmvrobj%get_qfrommvr0(n)
          dqfrommvr = qtolfact * (q0 - q)
        end if
        !
        ! -- evaluate magnitude of differences
        if (n == 1) then
          locdhmax = n
          dhmax = dh
          locdgwfmax = n
          residmax = resid
          locresidmax = n
          dgwfmax = dgwf
          locdqoutmax = n
          dqoutmax = dqout
          dqfrommvrmax = dqfrommvr
          locdqfrommvrmax = n
        else
          if (abs(dh) > abs(dhmax)) then
            locdhmax = n
            dhmax = dh
          end if
          if (abs(resid) > abs(residmax)) then
            locresidmax = n
            residmax = resid
          end if
          if (abs(dgwf) > abs(dgwfmax)) then
            locdgwfmax = n
            dgwfmax = dgwf
          end if
          if (abs(dqout) > abs(dqoutmax)) then
            locdqoutmax = n
            dqoutmax = dqout
          end if
          if (abs(dqfrommvr) > abs(dqfrommvrmax)) then
            dqfrommvrmax = dqfrommvr
            locdqfrommvrmax = n
          end if
        end if
      end do final_check
      !
      ! -- set dpak and cpak
      if (abs(dhmax) > abs(dpak)) then
        ipak = locdhmax
        dpak = dhmax
        write (cloc, "(a,'-',a)") &
          trim(this%packName), 'stage'
        cpak = trim(cloc)
      end if
      if (abs(residmax) > abs(dpak)) then
        ipak = locresidmax
        dpak = residmax
        write (cloc, "(a,'-',a)") &
          trim(this%packName), 'residual'
        cpak = trim(cloc)
      end if
      if (abs(dgwfmax) > abs(dpak)) then
        ipak = locdgwfmax
        dpak = dgwfmax
        write (cloc, "(a,'-',a)") &
          trim(this%packName), 'gwf'
        cpak = trim(cloc)
      end if
      if (this%noutlets > 0) then
        if (abs(dqoutmax) > abs(dpak)) then
          ipak = locdqoutmax
          dpak = dqoutmax
          write (cloc, "(a,'-',a)") &
            trim(this%packName), 'outlet'
          cpak = trim(cloc)
        end if
      end if
      if (this%imover == 1) then
        if (abs(dqfrommvrmax) > abs(dpak)) then
          ipak = locdqfrommvrmax
          dpak = dqfrommvrmax
          write (cloc, "(a,'-',a)") trim(this%packName), 'qfrommvr'
          cpak = trim(cloc)
        end if
      end if
      !
      ! -- write convergence data to package csv
      if (this%ipakcsv /= 0) then
        !
        ! -- write the data
        call this%pakcsvtab%add_term(innertot)
        call this%pakcsvtab%add_term(totim)
        call this%pakcsvtab%add_term(kper)
        call this%pakcsvtab%add_term(kstp)
        call this%pakcsvtab%add_term(kiter)
        call this%pakcsvtab%add_term(dhmax)
        call this%pakcsvtab%add_term(locdhmax)
        call this%pakcsvtab%add_term(residmax)
        call this%pakcsvtab%add_term(locresidmax)
        call this%pakcsvtab%add_term(dgwfmax)
        call this%pakcsvtab%add_term(locdgwfmax)
        if (this%noutlets > 0) then
          call this%pakcsvtab%add_term(dqoutmax)
          call this%pakcsvtab%add_term(locdqoutmax)
        end if
        if (this%imover == 1) then
          call this%pakcsvtab%add_term(dqfrommvrmax)
          call this%pakcsvtab%add_term(locdqfrommvrmax)
        end if
        !
        ! -- finalize the package csv
        if (iend == 1) then
          call this%pakcsvtab%finalize_table()
        end if
      end if
    end if
  end subroutine lak_cc
 
  !> @brief Calculate flows
  !<
  subroutine lak_cq(this, x, flowja, iadv)
    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(laktype), 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) :: chratin, chratout
    ! -- for budget
    integer(I4B) :: j, n, igwfnode
    real(DP) :: hlak, head, flow, dqdh
    real(DP) :: v0, v1, sa, sf
    !
    call this%lak_solve(update=.false.)
    !
    ! -- for the IMPLICIT formulation, report the lake terms with the same
    !    treatment that lak_fc_implicit assembled into the matrix, so the lake
    !    and gwf-cell budgets match the solved flows. lak_solve above set these
    !    from the substitution path (hard-cutoff seepage; availability-limited
    !    losses), which is correct for a lake on the legacy solver but
    !    not for an implicit lake, so overwrite them for the implicit
    !    (non-legacy) lakes only.
    if (this%iimplicit /= 0) then
      do n = 1, this%nlakes
        if (this%iboundpak(n) < 1 .or. this%ilegacy(n) /= 0) cycle
        hlak = this%xnewpak(n)
        !
        ! -- lakebed seepage: use the same exchange the matrix assembled
        do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
          igwfnode = this%gwfconn(j)
          if (this%ibound(igwfnode) < 1) cycle
          head = this%xnew(igwfnode)
          call this%lak_calculate_conn_exchange_deriv(n, j, hlak, head, &
                                                      flow, dqdh=dqdh)
          this%hcof(j) = -dqdh
          this%rhs(j) = -dqdh * head + flow
        end do
        !
        ! -- stage-driven losses: the matrix (lak_budget_nogwf) ramps
        !    evaporation and withdrawal toward zero as the lake approaches its
        !    bottom with the surfdep factor (sf) rather than the substitution
        !    solver's availability limiting, so report them the same way. This
        !    keeps the lake budget closed for a drying lake. (The outlet
        !    outflow is already near zero at the bottom, as its invert is above
        !    the lake bottom, and the storage term reflects the solved stage.)
        call this%lak_calculate_sarea(n, hlak, sa)
        sf = done
        if (this%surfdep > dzero) then
          sf = squadraticsaturation(this%lakebot(n) + this%surfdep, &
                                    this%lakebot(n), hlak)
        end if
        this%evap(n) = -this%evaporation(n) * sa * sf
        this%withr(n) = -this%withdrawal(n) * sf
      end do
    end if
    !
    ! -- call base functionality in bnd_cq.  This will calculate lake-gwf flows
    !    and put them into this%simvals
    call this%BndType%bnd_cq(x, flowja, iadv=1)
    !
    ! -- calculate several budget terms
    chratin = dzero
    chratout = dzero
    do n = 1, this%nlakes
      this%chterm(n) = dzero
      if (this%iboundpak(n) == 0) cycle
      hlak = this%xnewpak(n)
      call this%lak_calculate_vol(n, hlak, v1)
      !
      ! -- add budget terms for active lakes
      if (this%iboundpak(n) /= 0) then
        !
        ! -- rainfall
        rrate = this%precip(n)
        call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
        !
        ! -- evaporation
        rrate = this%evap(n)
        call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
        !
        ! -- runoff
        rrate = this%runoff(n)
        call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
        !
        ! -- inflow
        rrate = this%inflow(n)
        call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
        !
        ! -- withdrawals
        rrate = this%withr(n)
        call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
        !
        ! -- add lake storage changes
        rrate = dzero
        if (this%iboundpak(n) > 0) then
          if (this%gwfiss /= 1) then
            call this%lak_calculate_vol(n, this%xoldpak(n), v0)
            rrate = -(v1 - v0) / delt
            call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
          end if
        end if
        this%qsto(n) = rrate
        !
        ! -- add external outlets
        call this%lak_get_external_outlet(n, rrate)
        call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
        !
        ! -- add mover terms
        if (this%imover == 1) then
          if (this%iboundpak(n) /= 0) then
            rrate = this%pakmvrobj%get_qfrommvr(n)
          else
            rrate = dzero
          end if
          call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
        end if
      end if
    end do
    !
    ! -- gwf flow and constant flow to lake
    do n = 1, this%nlakes
      rrate = dzero
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        ! simvals is from aquifer perspective, and so it is positive
        ! for flow into the aquifer.  Need to switch sign for lake
        ! perspective.
        rrate = -this%simvals(j)
        this%qleak(j) = rrate
        if (this%iboundpak(n) /= 0) then
          call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
        end if
      end do
    end do
    !
    ! -- fill the budget object
    call this%lak_fill_budobj()
  end subroutine lak_cq
 
  !> @brief Output LAK package flow terms
  !<
  subroutine lak_ot_package_flows(this, icbcfl, ibudfl)
    use tdismodule, only: kstp, kper, delt, pertim, totim
    class(laktype) :: 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 lake flows table
    if (ibudfl /= 0 .and. this%iprflow /= 0) then
      call this%budobj%write_flowtable(this%dis, kstp, kper)
    end if
  end subroutine lak_ot_package_flows
 
  !> @brief Write flows to binary file and/or print flows to budget
  !<
  subroutine lak_ot_model_flows(this, icbcfl, ibudfl, icbcun, imap)
    class(laktype) :: 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 lak_ot_model_flows
 
  !> @brief Save LAK-calculated values to binary file
  !<
  subroutine lak_ot_dv(this, idvsave, idvprint)
    use tdismodule, only: kstp, kper, pertim, totim
    use constantsmodule, only: dhnoflo, dhdry
    use inputoutputmodule, only: ulasav
    class(laktype) :: this
    integer(I4B), intent(in) :: idvsave
    integer(I4B), intent(in) :: idvprint
    integer(I4B) :: ibinun
    integer(I4B) :: n
    real(DP) :: v
    real(DP) :: d
    real(DP) :: stage
    real(DP) :: sa
    real(DP) :: wa
    !
    ! -- set unit number for binary dependent variable output
    ibinun = 0
    if (this%istageout /= 0) then
      ibinun = this%istageout
    end if
    if (idvsave == 0) ibinun = 0
    !
    ! -- write lake binary output
    if (ibinun > 0) then
      do n = 1, this%nlakes
        v = this%xnewpak(n)
        d = v - this%lakebot(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, '           STAGE', kstp, kper, pertim, totim, &
                  this%nlakes, 1, 1, ibinun)
    end if
    !
    ! -- Print lake stage table
    if (idvprint /= 0 .and. this%iprhed /= 0) then
      !
      ! -- set table kstp and kper
      call this%stagetab%set_kstpkper(kstp, kper)
      !
      ! -- write data
      do n = 1, this%nlakes
        if (this%iboundpak(n) == 0) then
          stage = dhnoflo
          sa = dhnoflo
          wa = dhnoflo
          v = dhnoflo
        else
          stage = this%xnewpak(n)
          call this%lak_calculate_sarea(n, stage, sa)
          call this%lak_calculate_warea(n, stage, wa)
          call this%lak_calculate_vol(n, stage, v)
        end if
        if (this%inamedbound == 1) then
          call this%stagetab%add_term(this%lakename(n))
        end if
        call this%stagetab%add_term(n)
        call this%stagetab%add_term(stage)
        call this%stagetab%add_term(sa)
        call this%stagetab%add_term(wa)
        call this%stagetab%add_term(v)
      end do
    end if
  end subroutine lak_ot_dv
 
  !> @brief Write LAK budget to listing file
  !<
  subroutine lak_ot_bdsummary(this, kstp, kper, iout, ibudfl)
    ! -- module
    use tdismodule, only: totim, delt
    ! -- dummy
    class(laktype) :: this !< LakType 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 lak_ot_bdsummary
 
  !> @brief Deallocate objects
  !<
  subroutine lak_da(this)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(laktype) :: this
    !
    ! -- nullify input context runtime pointers
    call this%input%destroy()
    !
    ! -- arrays
    deallocate (this%lakename)
    deallocate (this%status)
    deallocate (this%clakbudget)
    call mem_deallocate(this%dbuff)
    deallocate (this%cauxcbc)
    call mem_deallocate(this%qauxcbc)
    call mem_deallocate(this%qleak)
    call mem_deallocate(this%holdconn)
    call mem_deallocate(this%qsto)
    call mem_deallocate(this%denseterms)
    call mem_deallocate(this%viscratios)
    !
    ! -- tables
    if (this%ntables > 0) then
      call mem_deallocate(this%ialaktab)
      call mem_deallocate(this%tabstage)
      call mem_deallocate(this%tabvolume)
      call mem_deallocate(this%tabsarea)
      call mem_deallocate(this%tabwarea)
    end if
    !
    ! -- budobj
    call this%budobj%budgetobject_da()
    deallocate (this%budobj)
    nullify (this%budobj)
    !
    ! -- outlets
    if (this%noutlets > 0) then
      call mem_deallocate(this%lakein)
      call mem_deallocate(this%lakeout)
      call mem_deallocate(this%iouttype)
      ! -- outrate/outinvert/outwidth/outrough/outslope are live aliases
      ! into the input context's permanent arrays, not owned here
      call mem_deallocate(this%simoutrate)
      call mem_deallocate(this%outlet_invert_set)
      call mem_deallocate(this%outlet_width_set)
      call mem_deallocate(this%outlet_rough_set)
      call mem_deallocate(this%outlet_slope_set)
    end if
    !
    ! -- stage table
    if (this%iprhed > 0) then
      call this%stagetab%table_da()
      deallocate (this%stagetab)
      nullify (this%stagetab)
    end if
    !
    ! -- package csv table
    if (this%ipakcsv > 0) then
      if (associated(this%pakcsvtab)) then
        call this%pakcsvtab%table_da()
        deallocate (this%pakcsvtab)
        nullify (this%pakcsvtab)
      end if
    end if
    !
    ! -- scalars
    call mem_deallocate(this%iprhed)
    call mem_deallocate(this%istageout)
    call mem_deallocate(this%ibudgetout)
    call mem_deallocate(this%ibudcsv)
    call mem_deallocate(this%ipakcsv)
    if (allocated(this%pakcsvfile)) deallocate (this%pakcsvfile)
    call mem_deallocate(this%nlakes)
    call mem_deallocate(this%noutlets)
    call mem_deallocate(this%ntables)
    call mem_deallocate(this%convlength)
    call mem_deallocate(this%convtime)
    call mem_deallocate(this%outdmax)
    call mem_deallocate(this%igwhcopt)
    call mem_deallocate(this%iconvchk)
    call mem_deallocate(this%maxlakit)
    call mem_deallocate(this%surfdep)
    call mem_deallocate(this%dmaxchg)
    call mem_deallocate(this%delh)
    call mem_deallocate(this%check_attr)
    call mem_deallocate(this%iimplicit)
    call mem_deallocate(this%iforceleg)
    call mem_deallocate(this%iforceleglak)
    call mem_deallocate(this%bditems)
    call mem_deallocate(this%cbcauxitems)
    call mem_deallocate(this%idense)
    !
    call mem_deallocate(this%nlakeconn)
    call mem_deallocate(this%idxlakeconn)
    call mem_deallocate(this%ntabrow)
    call mem_deallocate(this%strt)
    call mem_deallocate(this%laketop)
    call mem_deallocate(this%lakebot)
    call mem_deallocate(this%sareamax)
    ! -- stage/rainfall/evaporation/runoff/inflow/withdrawal are live aliases
    ! into the input context's permanent arrays, not owned here
    call mem_deallocate(this%avail)
    call mem_deallocate(this%lkgwsink)
    call mem_deallocate(this%ncncvr)
    call mem_deallocate(this%ilegacy)
    call mem_deallocate(this%nstuck)
    call mem_deallocate(this%surfin)
    call mem_deallocate(this%surfout)
    call mem_deallocate(this%surfout1)
    call mem_deallocate(this%precip)
    call mem_deallocate(this%precip1)
    call mem_deallocate(this%evap)
    call mem_deallocate(this%evap1)
    call mem_deallocate(this%evapo)
    call mem_deallocate(this%withr)
    call mem_deallocate(this%withr1)
    call mem_deallocate(this%flwin)
    call mem_deallocate(this%flwiter)
    call mem_deallocate(this%flwiter1)
    call mem_deallocate(this%seep)
    call mem_deallocate(this%seep1)
    call mem_deallocate(this%seep0)
    call mem_deallocate(this%stageiter)
    call mem_deallocate(this%chterm)
    !
    ! -- lake boundary and stages
    if (this%iimplicit == 0) then
      call mem_deallocate(this%iboundpak)
      call mem_deallocate(this%xnewpak)
    else
      ! -- iboundpak aliases the global ibound (not deallocated here); xnewpak
      !    was checked in to the global x vector
      call mem_deallocate(this%xnewpak, 'XNEWPAK', this%memoryPath)
    end if
    call mem_deallocate(this%xoldpak)
    !
    ! -- implicit-formulation matrix-mapping arrays
    call mem_deallocate(this%idxlocnode)
    call mem_deallocate(this%idxdiag)
    call mem_deallocate(this%idxoffdglo)
    call mem_deallocate(this%idxsymdglo)
    call mem_deallocate(this%idxsymoffdglo)
    !
    ! -- lake iteration variables
    call mem_deallocate(this%iseepc)
    call mem_deallocate(this%idhc)
    call mem_deallocate(this%en1)
    call mem_deallocate(this%en2)
    call mem_deallocate(this%r1)
    call mem_deallocate(this%r2)
    call mem_deallocate(this%dh0)
    call mem_deallocate(this%s0)
    call mem_deallocate(this%qgwf0)
    !
    ! -- lake connection variables
    call mem_deallocate(this%imap)
    call mem_deallocate(this%gwfconn)
    call mem_deallocate(this%nodesontop)
    call mem_deallocate(this%ictype)
    call mem_deallocate(this%bedleak)
    call mem_deallocate(this%belev)
    call mem_deallocate(this%telev)
    call mem_deallocate(this%connlength)
    call mem_deallocate(this%connwidth)
    call mem_deallocate(this%sarea)
    call mem_deallocate(this%warea)
    call mem_deallocate(this%satcond)
    call mem_deallocate(this%simcond)
    call mem_deallocate(this%simlakgw)
    !
    ! -- pointers to gwf variables
    nullify (this%gwfiss)
    !
    ! -- Parent object
    call this%BndExtType%bnd_da()
  end subroutine lak_da
 
  !> @brief Define the list heading that is written to iout when PRINT_INPUT
  !! option is used
  !<
  subroutine define_listlabel(this)
    ! -- modules
    class(laktype), 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
  !! access to these things
  !<
  subroutine lak_set_pointers(this, neq, ibound, xnew, xold, flowja)
    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_checkin
    ! -- dummy
    class(laktype) :: 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)
    !
    ! -- for the implicit formulation, point the lake stage and ibound at the
    !    matching slice of the model solution (xnew) and ibound vectors so the
    !    stage is solved directly in the matrix. (The default formulation keeps
    !    these as its own arrays, allocated in lak_allocate_arrays.)
    if (this%iimplicit /= 0) then
      istart = this%dis%nodes + this%ioffset + 1
      iend = istart + this%nlakes - 1
      this%iboundpak => this%ibound(istart:iend)
      this%xnewpak => this%xnew(istart:iend)
      call mem_checkin(this%xnewpak, 'XNEWPAK', this%memoryPath, 'X', &
                       this%memoryPathModel)
      !
      ! -- initialize xnewpak
      do n = 1, this%nlakes
        this%xnewpak(n) = dep20
      end do
    end if
  end subroutine lak_set_pointers
 
  !> @brief Add the lake rows and columns to the sparse matrix
  !!
  !! Implicit formulation only. Each lake adds one equation (row and column) to
  !! the groundwater flow matrix: a diagonal entry, plus a symmetric off-diagonal
  !! pair (lake-to-cell and cell-to-lake) for every lake-cell connection.
  !<
  subroutine lak_ac(this, moffset, sparse)
    use sparsemodule, only: sparsematrix
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: moffset
    type(sparsematrix), intent(inout) :: sparse
    ! -- local
    integer(I4B) :: j, n
    integer(I4B) :: jj
    integer(I4B) :: jglo
    integer(I4B) :: nglo
    !
    ! -- the default formulation adds no rows to the matrix
    if (this%iimplicit == 0) return
    !
    ! -- add a row for each lake and its connections to the cells
    do n = 1, this%nlakes
      nglo = moffset + this%dis%nodes + this%ioffset + n
      call sparse%addconnection(nglo, nglo, 1)
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        jj = this%gwfconn(j)
        jglo = jj + moffset
        call sparse%addconnection(nglo, jglo, 1)
        call sparse%addconnection(jglo, nglo, 1)
      end do
    end do
  end subroutine lak_ac
 
  !> @brief Find the matrix position of each lake row and connection
  !!
  !! Implicit formulation only. Finds the position in the assembled matrix of
  !! each lake-row diagonal, lake-to-cell, cell diagonal, and cell-to-lake entry
  !! and stores them in the idx* arrays so they can be filled during assembly.
  !<
  subroutine lak_mc(this, moffset, matrix_sln)
    use memorymanagermodule, only: mem_allocate
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: moffset
    class(matrixbasetype), pointer :: matrix_sln
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: j
    integer(I4B) :: iglo
    integer(I4B) :: jglo
    integer(I4B) :: ipos
    !
    ! -- the connection-mapping vectors are only used by the implicit
    !    formulation; allocate them at size 0 otherwise so legacy LAK runs do not
    !    pay the maxbound memory cost (lak_da stays symmetric either way)
    if (this%iimplicit == 0) then
      call mem_allocate(this%idxlocnode, 0, 'IDXLOCNODE', this%memoryPath)
      call mem_allocate(this%idxdiag, 0, 'IDXDIAG', this%memoryPath)
      call mem_allocate(this%idxoffdglo, 0, 'IDXOFFDGLO', this%memoryPath)
      call mem_allocate(this%idxsymdglo, 0, 'IDXSYMDGLO', this%memoryPath)
      call mem_allocate(this%idxsymoffdglo, 0, 'IDXSYMOFFDGLO', this%memoryPath)
      return
    end if
    call mem_allocate(this%idxlocnode, this%nlakes, 'IDXLOCNODE', &
                      this%memoryPath)
    call mem_allocate(this%idxdiag, this%nlakes, 'IDXDIAG', 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)
    !
    ! -- lake rows: a per-lake diagonal position and the per-connection
    !    lake->cell off-diagonals. The diagonal is stored per lake (idxdiag) so a
    !    lake with no connections still has a valid diagonal position.
    ipos = 1
    do n = 1, this%nlakes
      iglo = moffset + this%dis%nodes + this%ioffset + n
      this%idxlocnode(n) = this%dis%nodes + this%ioffset + n
      this%idxdiag(n) = matrix_sln%get_position_diag(iglo)
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        jglo = this%gwfconn(j) + moffset
        this%idxoffdglo(ipos) = matrix_sln%get_position(iglo, jglo)
        ipos = ipos + 1
      end do
    end do
    !
    ! -- lake contributions to the gwf portion of the global matrix
    ipos = 1
    do n = 1, this%nlakes
      jglo = moffset + this%dis%nodes + this%ioffset + n
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        iglo = this%gwfconn(j) + moffset
        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 lak_mc
 
  !> @brief Procedures related to observations (type-bound)
  !!
  !! Return true because LAK package supports observations. Overrides
  !! BndType%bnd_obs_supported()
  !<
  logical function lak_obs_supported(this)
    ! -- dummy
    class(laktype) :: this
    !
    lak_obs_supported = .true.
  end function lak_obs_supported
 
  !> @brief Store observation type supported by LAK package. Overrides
  !! BndType%bnd_df_obs
  !<
  subroutine lak_df_obs(this)
    ! -- dummy
    class(laktype) :: this
    ! -- local
    integer(I4B) :: indx
    !
    ! -- Store obs type and assign procedure pointer
    !    for stage observation type.
    call this%obs%StoreObsType('stage', .false., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for ext-inflow observation type.
    call this%obs%StoreObsType('ext-inflow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for outlet-inflow observation type.
    call this%obs%StoreObsType('outlet-inflow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for inflow observation type.
    call this%obs%StoreObsType('inflow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for from-mvr observation type.
    call this%obs%StoreObsType('from-mvr', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for rainfall observation type.
    call this%obs%StoreObsType('rainfall', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for runoff observation type.
    call this%obs%StoreObsType('runoff', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for lak observation type.
    call this%obs%StoreObsType('lak', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for evaporation observation type.
    call this%obs%StoreObsType('evaporation', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for withdrawal observation type.
    call this%obs%StoreObsType('withdrawal', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for ext-outflow observation type.
    call this%obs%StoreObsType('ext-outflow', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for to-mvr observation type.
    call this%obs%StoreObsType('to-mvr', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_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 => lak_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 => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for outlet observation type.
    call this%obs%StoreObsType('outlet', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for volume observation type.
    call this%obs%StoreObsType('volume', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for surface-area observation type.
    call this%obs%StoreObsType('surface-area', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for wetted-area observation type.
    call this%obs%StoreObsType('wetted-area', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
    !
    ! -- Store obs type and assign procedure pointer
    !    for conductance observation type.
    call this%obs%StoreObsType('conductance', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
  end subroutine lak_df_obs
 
  !> @brief Calculate observations this time step and call ObsType%SaveOneSimval
  !! for each LakType observation.
  !<
  subroutine lak_bd_obs(this)
    ! -- dummy
    class(laktype) :: this
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: igwfnode
    integer(I4B) :: j
    integer(I4B) :: jj
    integer(I4B) :: n
    real(DP) :: hgwf
    real(DP) :: hlak
    real(DP) :: v
    real(DP) :: v2
    type(observetype), pointer :: obsrv => null()
    !
    ! Write simulated values for all LAK 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 ('STAGE')
            if (this%iboundpak(jj) /= 0) then
              v = this%xnewpak(jj)
            end if
          case ('EXT-INFLOW')
            if (this%iboundpak(jj) /= 0) then
              call this%lak_calculate_inflow(jj, v)
            end if
          case ('OUTLET-INFLOW')
            if (this%iboundpak(jj) /= 0) then
              call this%lak_calculate_outlet_inflow(jj, v)
            end if
          case ('INFLOW')
            if (this%iboundpak(jj) /= 0) then
              call this%lak_calculate_inflow(jj, v)
              call this%lak_calculate_outlet_inflow(jj, v2)
              v = v + v2
            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 ('RAINFALL')
            if (this%iboundpak(jj) /= 0) then
              v = this%precip(jj)
            end if
          case ('RUNOFF')
            if (this%iboundpak(jj) /= 0) then
              v = this%runoff(jj)
            end if
          case ('LAK')
            n = this%imap(jj)
            if (this%iboundpak(n) /= 0) then
              igwfnode = this%gwfconn(jj)
              hgwf = this%xnew(igwfnode)
              if (this%hcof(jj) /= dzero) then
                v = -(this%hcof(jj) * (this%xnewpak(n) - hgwf))
              else
                v = -this%rhs(jj)
              end if
            end if
          case ('EVAPORATION')
            if (this%iboundpak(jj) /= 0) then
              v = this%evap(jj)
            end if
          case ('WITHDRAWAL')
            if (this%iboundpak(jj) /= 0) then
              v = this%withr(jj)
            end if
          case ('EXT-OUTFLOW')
            n = this%lakein(jj)
            if (this%iboundpak(n) /= 0) then
              if (this%lakeout(jj) == 0) then
                v = this%simoutrate(jj)
                if (v < dzero) then
                  if (this%imover == 1) then
                    v = v + this%pakmvrobj%get_qtomvr(jj)
                  end if
                end if
              end if
            end if
          case ('TO-MVR')
            n = this%lakein(jj)
            if (this%iboundpak(n) /= 0) then
              if (this%imover == 1) then
                v = this%pakmvrobj%get_qtomvr(jj)
                if (v > dzero) then
                  v = -v
                end if
              end if
            end if
          case ('STORAGE')
            if (this%iboundpak(jj) /= 0) then
              v = this%qsto(jj)
            end if
          case ('CONSTANT')
            if (this%iboundpak(jj) /= 0) then
              v = this%chterm(jj)
            end if
          case ('OUTLET')
            n = this%lakein(jj)
            if (this%iboundpak(n) /= 0) then
              v = this%simoutrate(jj)
            end if
          case ('VOLUME')
            if (this%iboundpak(jj) /= 0) then
              call this%lak_calculate_vol(jj, this%xnewpak(jj), v)
            end if
          case ('SURFACE-AREA')
            if (this%iboundpak(jj) /= 0) then
              hlak = this%xnewpak(jj)
              call this%lak_calculate_sarea(jj, hlak, v)
            end if
          case ('WETTED-AREA')
            n = this%imap(jj)
            if (this%iboundpak(n) /= 0) then
              hlak = this%xnewpak(n)
              igwfnode = this%gwfconn(jj)
              hgwf = this%xnew(igwfnode)
              call this%lak_calculate_conn_warea(n, jj, hlak, hgwf, v)
            end if
          case ('CONDUCTANCE')
            n = this%imap(jj)
            if (this%iboundpak(n) /= 0) then
              hlak = this%xnewpak(n)
              igwfnode = this%gwfconn(jj)
              hgwf = this%xnew(igwfnode)
              call this%lak_calculate_conn_conductance(n, jj, hlak, hgwf, v)
            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_unit(this%inunit)
      end if
    end if
  end subroutine lak_bd_obs
 
  !> @brief Process each observation
  !!
  !! Only done the first stress period since boundaries are fixed for the
  !! simulation
  !<
  subroutine lak_rp_obs(this)
    use tdismodule, only: kper
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: nn1
    integer(I4B) :: nn2
    integer(I4B) :: jj
    character(len=LENBOUNDNAME) :: bname
    logical(LGP) :: jfound
    class(observetype), pointer :: obsrv => null()
    ! -- formats
10  format('Boundary "', a, '" for observation "', a, &
           '" is invalid in package "', a, '"')
    !
    ! -- process each package observation
    !    only done the first stress period since boundaries are fixed
    !    for the simulation
    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 lake is based on a boundary name.
            !    Iterate through all lakes to identify and store
            !    corresponding index in bound array.
            jfound = .false.
            if (obsrv%ObsTypeId == 'LAK' .or. &
                obsrv%ObsTypeId == 'CONDUCTANCE' .or. &
                obsrv%ObsTypeId == 'WETTED-AREA') then
              do j = 1, this%nlakes
                do jj = this%idxlakeconn(j), this%idxlakeconn(j + 1) - 1
                  if (this%boundname(jj) == bname) then
                    jfound = .true.
                    call obsrv%AddObsIndex(jj)
                  end if
                end do
              end do
            else if (obsrv%ObsTypeId == 'EXT-OUTFLOW' .or. &
                     obsrv%ObsTypeId == 'TO-MVR' .or. &
                     obsrv%ObsTypeId == 'OUTLET') then
              do j = 1, this%noutlets
                jj = this%lakein(j)
                if (this%lakename(jj) == bname) then
                  jfound = .true.
                  call obsrv%AddObsIndex(j)
                end if
              end do
            else
              do j = 1, this%nlakes
                if (this%lakename(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 == 'LAK' .or. &
                obsrv%ObsTypeId == 'CONDUCTANCE' .or. &
                obsrv%ObsTypeId == 'WETTED-AREA') then
              nn2 = obsrv%NodeNumber2
              j = this%idxlakeconn(nn1) + nn2 - 1
              call obsrv%AddObsIndex(j)
            else
              call obsrv%AddObsIndex(nn1)
            end if
          else
            errmsg = 'Programming error in lak_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 == 'STAGE') 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 lake with a unique boundname.'
            call store_error(errmsg)
          end if
        end if
        !
        ! -- check that index values are valid
        if (obsrv%ObsTypeId == 'TO-MVR' .or. &
            obsrv%ObsTypeId == 'EXT-OUTFLOW' .or. &
            obsrv%ObsTypeId == 'OUTLET') then
          do j = 1, obsrv%indxbnds_count
            nn1 = obsrv%indxbnds(j)
            if (nn1 < 1 .or. nn1 > this%noutlets) then
              write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
                trim(adjustl(obsrv%ObsTypeId)), &
                ' outlet must be > 0 and <=', this%noutlets, &
                '(specified value is ', nn1, ')'
              call store_error(errmsg)
            end if
          end do
        else if (obsrv%ObsTypeId == 'LAK' .or. &
                 obsrv%ObsTypeId == 'CONDUCTANCE' .or. &
                 obsrv%ObsTypeId == 'WETTED-AREA') then
          do j = 1, obsrv%indxbnds_count
            nn1 = obsrv%indxbnds(j)
            if (nn1 < 1 .or. nn1 > this%maxbound) then
              write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
                trim(adjustl(obsrv%ObsTypeId)), &
                'lake connection number must be > 0 and <=', this%maxbound, &
                '(specified value is ', nn1, ')'
              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%nlakes) then
              write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
                trim(adjustl(obsrv%ObsTypeId)), &
                ' lake must be > 0 and <=', this%nlakes, &
                '(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_unit(this%inunit)
      end if
    end if
  end subroutine lak_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 LAK package
  !! observations.
  !<
  subroutine lak_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
    !
    string = obsrv%IDstring
    ! -- Extract lake number from string and store it.
    !    If 1st item is not an integer(I4B), it should be a
    !    lake name--deal with it.
    icol = 1
    ! -- get lake 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 == 'LAK' .or. obsrv%ObsTypeId == 'CONDUCTANCE' .or. &
          obsrv%ObsTypeId == 'WETTED-AREA') 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 (iconn) 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 lake number (NodeNumber)
    obsrv%NodeNumber = nn1
  end subroutine lak_process_obsid
 
  !
  ! -- private LAK methods
  !
 
  !> @brief Accumulate constant head terms for budget
  !<
  subroutine lak_accumulate_chterm(this, ilak, rrate, chratin, chratout)
    ! -- dummy
    class(laktype) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(in) :: rrate
    real(DP), intent(inout) :: chratin
    real(DP), intent(inout) :: chratout
    ! -- locals
    real(DP) :: q
    !
    ! code
    if (this%iboundpak(ilak) < 0) then
      q = -rrate
      this%chterm(ilak) = this%chterm(ilak) + q
      !
      ! -- See if flow is into lake or out of lake.
      if (q < dzero) then
        !
        ! -- Flow is out of lake subtract rate from ratout.
        chratout = chratout - q
      else
        !
        ! -- Flow is into lake; add rate to ratin.
        chratin = chratin + q
      end if
    end if
  end subroutine lak_accumulate_chterm
 
  !> @brief Store the lake head and connection conductance in the bound array
  !<
  subroutine lak_bound_update(this)
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    integer(I4B) :: j, n, node
    real(DP) :: hlak, head, clak
    !
    ! -- Return if no lak lakes
    if (this%nbound == 0) return
    !
    ! -- Calculate hcof and rhs for each lak entry
    do n = 1, this%nlakes
      hlak = this%xnewpak(n)
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        node = this%gwfconn(j)
        head = this%xnew(node)
        call this%lak_calculate_conn_conductance(n, j, hlak, head, clak)
        this%bound(1, j) = hlak
        this%bound(2, j) = clak
      end do
    end do
  end subroutine lak_bound_update
 
  !> @brief Solve for lake stage
  !!
  !! Solve the lake stage by substitution. With only_legacy set, solve only
  !! the lakes flagged for the legacy solver (this%ilegacy /= 0) and
  !! leave the remaining lakes -- whose stage is solved in the global matrix by
  !! the IMPLICIT formulation -- untouched. Without it (the default), every
  !! active lake is solved.
  !<
  subroutine lak_solve(this, update, only_legacy)
    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(laktype), intent(inout) :: this
    logical(LGP), intent(in), optional :: update
    logical(LGP), intent(in), optional :: only_legacy
    ! -- local
    logical(LGP) :: lupdate
    logical(LGP) :: legonly
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: iicnvg
    integer(I4B) :: iter
    integer(I4B) :: maxiter
    integer(I4B) :: ncnv
    real(DP) :: hlak
    real(DP) :: hlak0
    real(DP) :: v0
    real(DP) :: v1
    real(DP) :: ro
    real(DP) :: qinf
    real(DP) :: ex
    real(DP) :: outinf
    real(DP) :: avail
    !
    ! -- set lupdate
    if (present(update)) then
      lupdate = update
    else
      lupdate = .true.
    end if
    !
    ! -- set legonly (solve only the lakes flagged for the legacy solver)
    if (present(only_legacy)) then
      legonly = only_legacy
    else
      legonly = .false.
    end if
    !
    ! -- initialize
    avail = dzero
    !
    ! -- initialize
    do n = 1, this%nlakes
      ! -- a lake not being solved on this call (an IMPLICIT lake when only the
      !    flagged lakes are solved) is treated as already converged and left
      !    untouched
      if (legonly .and. this%ilegacy(n) == 0) then
        this%ncncvr(n) = 1
        cycle
      end if
      this%ncncvr(n) = 0
      this%surfin(n) = dzero
      this%surfout(n) = dzero
      this%surfout1(n) = dzero
      if (this%xnewpak(n) < this%lakebot(n)) then
        this%xnewpak(n) = this%lakebot(n)
      end if
      if (this%gwfiss /= 0) then
        this%xoldpak(n) = this%xnewpak(n)
      end if
      ! -- lake iteration items
      this%iseepc(n) = 0
      this%idhc(n) = 0
      this%en1(n) = this%lakebot(n)
      call this%lak_calculate_residual(n, this%en1(n), this%r1(n))
      this%en2(n) = this%laketop(n)
      call this%lak_calculate_residual(n, this%en2(n), this%r2(n))
    end do
    ! -- a legacy-only solve keeps the outlet rates of the implicit lakes,
    !    which are final and still needed for the mover and downstream lakes
    do n = 1, this%noutlets
      if (legonly .and. this%ilegacy(this%lakein(n)) == 0) cycle
      this%simoutrate(n) = dzero
    end do
    !
    ! -- sum up inflows from mover inflows
    do n = 1, this%nlakes
      call this%lak_calculate_outlet_inflow(n, this%surfin(n))
    end do
    !
    ! -- sum up overland runoff, inflows, and external flows into lake
    !    (includes maximum lake volume)
    do n = 1, this%nlakes
      hlak0 = this%xoldpak(n)
      hlak = this%xnewpak(n)
      call this%lak_calculate_runoff(n, ro)
      call this%lak_calculate_inflow(n, qinf)
      call this%lak_calculate_external(n, ex)
      call this%lak_calculate_vol(n, hlak0, v0)
      call this%lak_calculate_vol(n, hlak, v1)
      this%flwin(n) = this%surfin(n) + ro + qinf + ex + &
                      max(v0, v1) / delt
    end do
    !
    ! -- sum up inflows from upstream outlets
    do n = 1, this%nlakes
      call this%lak_calculate_outlet_inflow(n, outinf)
      this%flwin(n) = this%flwin(n) + outinf
    end do
    !
    iicnvg = 0
    maxiter = this%maxlakit
    !
    ! -- outer loop
    converge: do iter = 1, maxiter
      ncnv = 0
      do n = 1, this%nlakes
        if (this%ncncvr(n) == 0) ncnv = 1
      end do
      if (iter == maxiter) ncnv = 0
      if (ncnv == 0) iicnvg = 1
      !
      ! -- initialize variables
      do n = 1, this%nlakes
        this%evap(n) = dzero
        this%precip(n) = dzero
        this%precip1(n) = dzero
        this%seep(n) = dzero
        this%seep1(n) = dzero
        this%evap(n) = dzero
        this%evap1(n) = dzero
        this%evapo(n) = dzero
        this%withr(n) = dzero
        this%withr1(n) = dzero
        this%flwiter(n) = this%flwin(n)
        this%flwiter1(n) = this%flwin(n)
        if (this%gwfiss /= 0) then
          this%flwiter(n) = dep20
          this%flwiter1(n) = dep20
        end if
        do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
          this%hcof(j) = dzero
          this%rhs(j) = dzero
        end do
      end do
      !
      do n = 1, this%nlakes
        if (legonly .and. this%ilegacy(n) == 0) cycle
        call this%lak_estimate_seepage_single(n, ncnv)
      end do
      !
      laklevel: do n = 1, this%nlakes
        if (legonly .and. this%ilegacy(n) == 0) cycle laklevel
        call this%lak_solve_single(n, iter, maxiter, ncnv, lupdate)
      end do laklevel
      !
      if (iicnvg == 1) exit converge
      !
    end do converge
    !
    ! -- Mover terms: store outflow after diversion loss
    !    as qformvr and reduce outflow (qd)
    !    by how much was actually sent to the mover. A legacy-only solve
    !    skips this because the implicit formulation accumulates the mover
    !    terms for every outlet once the legacy-solved stages are known.
    if (this%imover == 1 .and. .not. legonly) then
      do n = 1, this%noutlets
        call this%pakmvrobj%accumulate_qformvr(n, -this%simoutrate(n))
      end do
    end if
  end subroutine lak_solve
 
  !> @brief Advance one lake stage by a single substitution iteration
  !!
  !! Performs one Newton update (with a bisection backup) of lake n's stage in
  !! the default substitution solver, using the per-iteration seepage estimate
  !! (this%seep / this%seep1) already computed for the lake. Extracted from
  !! lak_solve so the same per-lake update can be reused to solve a single lake
  !! stage on its own -- the per-lake legacy solve for the IMPLICIT formulation.
  !<
  subroutine lak_solve_single(this, n, iter, maxiter, ncnv, lupdate)
    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: iter
    integer(I4B), intent(in) :: maxiter
    integer(I4B), intent(in) :: ncnv
    logical(LGP), intent(in) :: lupdate
    ! -- local
    integer(I4B) :: ibflg
    integer(I4B) :: idhp
    real(DP) :: hlak
    real(DP) :: hlak0
    real(DP) :: ra
    real(DP) :: ro
    real(DP) :: qinf
    real(DP) :: ex
    real(DP) :: ev
    real(DP) :: wr
    real(DP) :: resid
    real(DP) :: resid1
    real(DP) :: residb
    real(DP) :: derv
    real(DP) :: dh
    real(DP) :: ts
    real(DP) :: adh
    real(DP) :: adh0
    real(DP) :: v0
    real(DP) :: v1
    real(DP) :: area
    real(DP) :: qtolfact
    real(DP) :: delh
    !
    delh = this%delh
    !
    ! -- skip inactive lakes
    if (this%iboundpak(n) == 0) then
      this%ncncvr(n) = 1
      return
    end if
    ibflg = 0
    hlak = this%xnewpak(n)
    if (iter < maxiter) then
      this%stageiter(n) = this%xnewpak(n)
    end if
    call this%lak_calculate_rainfall(n, hlak, ra)
    this%precip(n) = ra
    this%flwiter(n) = this%flwiter(n) + ra
    call this%lak_calculate_rainfall(n, hlak + delh, ra)
    this%precip1(n) = ra
    this%flwiter1(n) = this%flwiter1(n) + ra
    !
    ! -- limit withdrawals to lake inflows and lake storage
    call this%lak_calculate_withdrawal(n, this%flwiter(n), wr)
    this%withr(n) = wr
    call this%lak_calculate_withdrawal(n, this%flwiter1(n), wr)
    this%withr1(n) = wr
    !
    ! -- limit evaporation to lake inflows and lake storage
    call this%lak_calculate_evaporation(n, hlak, this%flwiter(n), ev)
    this%evap(n) = ev
    call this%lak_calculate_evaporation(n, hlak + delh, this%flwiter1(n), ev)
    this%evap1(n) = ev
    !
    ! -- no outlet flow if evaporation consumes all water
    call this%lak_calculate_outlet_outflow(n, hlak + delh, &
                                           this%flwiter1(n), &
                                           this%surfout1(n))
    call this%lak_calculate_outlet_outflow(n, hlak, this%flwiter(n), &
                                           this%surfout(n))
    !
    ! -- update the surface inflow values
    call this%lak_calculate_outlet_inflow(n, this%surfin(n))
    !
    !
    if (ncnv == 1) then
      if (this%iboundpak(n) > 0 .and. lupdate .eqv. .true.) then
        !
        ! -- recalculate flwin
        hlak0 = this%xoldpak(n)
        hlak = this%xnewpak(n)
        call this%lak_calculate_vol(n, hlak0, v0)
        call this%lak_calculate_vol(n, hlak, v1)
        call this%lak_calculate_runoff(n, ro)
        call this%lak_calculate_inflow(n, qinf)
        call this%lak_calculate_external(n, ex)
        this%flwin(n) = this%surfin(n) + ro + qinf + ex + &
                        max(v0, v1) / delt
        !
        ! -- compute new lake stage using Newton's method
        resid = this%precip(n) + this%evap(n) + this%withr(n) + ro + &
                qinf + ex + this%surfin(n) + &
                this%surfout(n) + this%seep(n)
        resid1 = this%precip1(n) + this%evap1(n) + this%withr1(n) + ro + &
                 qinf + ex + this%surfin(n) + &
                 this%surfout1(n) + this%seep1(n)
        !
        ! -- add storage changes for transient stress periods
        hlak = this%xnewpak(n)
        if (this%gwfiss /= 1) then
          call this%lak_calculate_vol(n, hlak, v1)
          resid = resid + (v0 - v1) / delt
          call this%lak_calculate_vol(n, hlak + delh, v1)
          resid1 = resid1 + (v0 - v1) / delt
        end if
        !
        ! -- determine the derivative and the stage change
        if (abs(resid1 - resid) > dzero) then
          derv = (resid1 - resid) / delh
          dh = dzero
          if (abs(derv) > dprec) then
            dh = resid / derv
          end if
        else
          if (resid < dzero) then
            resid = dzero
          end if
          call this%lak_vol2stage(n, resid, dh)
          dh = hlak - dh
          this%ncncvr(n) = 1
        end if
        !
        ! -- determine if the updated stage is outside the endpoints
        ts = hlak - dh
        if (iter == 1) this%dh0(n) = dh
        adh = abs(dh)
        adh0 = abs(this%dh0(n))
        if ((ts >= this%en2(n)) .or. (ts < this%en1(n))) then
          ! -- use bisection if dh is increasing or updated stage is below the
          !    bottom of the lake
          if ((adh > adh0) .or. (ts - this%lakebot(n)) < dprec) then
            residb = resid
            call this%lak_bisection(n, ibflg, hlak, ts, dh, residb)
          end if
        end if
        !
        ! -- set seep0 on the first lake iteration
        if (iter == 1) then
          this%seep0(n) = this%seep(n)
        end if
        !
        ! -- check for slow convergence
        if (this%seep(n) * this%seep0(n) < dprec) then
          this%iseepc(n) = this%iseepc(n) + 1
        else
          this%iseepc(n) = 0
        end if
        ! -- determine of convergence is slow and oscillating
        idhp = 0
        if (dh * this%dh0(n) < dprec) idhp = 1
        ! -- determine if stage change is increasing
        adh = abs(dh)
        if (adh > adh0) idhp = 1
        ! -- increment idhc convergence flag
        if (idhp == 1) then
          this%idhc(n) = this%idhc(n) + 1
        end if
        !
        ! -- switch to bisection when the Newton-Raphson method oscillates
        !    or when convergence is slow
        if (ibflg == 1) then
          if (this%iseepc(n) > 7 .or. this%idhc(n) > 12) then
            call this%lak_bisection(n, ibflg, hlak, ts, dh, residb)
          end if
        end if
      else
        dh = dzero
      end if
      !
      ! -- update lake stage
      hlak = hlak - dh
      if (hlak < this%lakebot(n)) then
        hlak = this%lakebot(n)
      end if
      !
      ! -- calculate surface area
      call this%lak_calculate_sarea(n, hlak, area)
      !
      ! -- set the Q to length factor
      if (area > dzero) then
        qtolfact = delt / area
      else
        qtolfact = dzero
      end if
      !
      ! -- recalculate the residual
      call this%lak_calculate_residual(n, hlak, resid)
      !
      ! -- evaluate convergence
      !if (ABS(dh) < delh) then
      if (abs(dh) < delh .and. abs(resid) * qtolfact < this%dmaxchg) then
        this%ncncvr(n) = 1
      end if
      this%xnewpak(n) = hlak
      !
      ! -- save iterates for lake
      this%seep0(n) = this%seep(n)
      this%dh0(n) = dh
    end if
  end subroutine lak_solve_single
 
  !> @brief Estimate the lakebed seepage for a single lake
  !!
  !! Runs the two-pass connection-seepage estimate for lake n at its current
  !! stage and the perturbed stage, accumulating this%seep(n)/this%seep1(n) and,
  !! on the final pass when ncnv == 0, the gwf-cell hcof/rhs contributions. Lakes
  !! do not couple within this estimate (each connection touches only its own
  !! lake's flwiter), so evaluating both passes per lake is equivalent to the
  !! all-lakes-per-pass ordering in lak_solve. Extracted from lak_solve so the
  !! same estimate can drive the single-lake legacy solve for the IMPLICIT
  !! formulation.
  !<
  subroutine lak_estimate_seepage_single(this, n, ncnv)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: n
    integer(I4B), intent(in) :: ncnv
    ! -- local
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: igwfnode
    integer(I4B) :: idry
    integer(I4B) :: idry1
    real(DP) :: hlak
    real(DP) :: head
    real(DP) :: qlakgw
    real(DP) :: qlakgw1
    real(DP) :: gwfhcof
    real(DP) :: gwfrhs
    real(DP) :: delh
    !
    delh = this%delh
    !
    ! -- skip inactive lakes
    if (this%iboundpak(n) == 0) return
    !
    estseep: do i = 1, 2
      ! - set xoldpak to xnewpak if steady-state
      if (this%gwfiss /= 0) then
        this%xoldpak(n) = this%xnewpak(n)
      end if
      hlak = this%xnewpak(n)
      calcconnseep: do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        igwfnode = this%gwfconn(j)
        head = this%xnew(igwfnode)
        if (this%ncncvr(n) /= 2) then
          if (this%ibound(igwfnode) > 0) then
            call this%lak_estimate_conn_exchange(i, n, j, idry, hlak, &
                                                 head, qlakgw, &
                                                 this%flwiter(n), &
                                                 gwfhcof, gwfrhs)
            call this%lak_estimate_conn_exchange(i, n, j, idry1, &
                                                 hlak + delh, head, qlakgw1, &
                                                 this%flwiter1(n))
            !
            ! -- add to gwf matrix
            if (ncnv == 0 .and. i == 2) then
              if (j == this%maxbound) then
                this%ncncvr(n) = 2
              end if
              if (idry /= 1) then
                this%hcof(j) = gwfhcof
                this%rhs(j) = gwfrhs
              else
                this%hcof(j) = dzero
                this%rhs(j) = qlakgw
              end if
            end if
            if (i == 2) then
              this%seep(n) = this%seep(n) + qlakgw
              this%seep1(n) = this%seep1(n) + qlakgw1
            end if
          end if
        end if
      end do calcconnseep
    end do estseep
  end subroutine lak_estimate_seepage_single
 
  !> @ brief Lake package bisection method
  !!
  !!  Use bisection method to find lake stage that reduces the residual
  !<
  subroutine lak_bisection(this, n, ibflg, hlak, temporary_stage, dh, residual)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: n !< lake number
    integer(I4B), intent(inout) :: ibflg !< bisection flag
    real(DP), intent(in) :: hlak !< lake stage
    real(DP), intent(inout) :: temporary_stage !< temporary lake stage
    real(DP), intent(inout) :: dh !< lake stage change
    real(DP), intent(inout) :: residual !< lake residual
    ! -- local
    integer(I4B) :: i
    real(DP) :: temporary_stage0
    real(DP) :: residuala
    real(DP) :: endpoint1
    real(DP) :: endpoint2
    ! -- code
    ibflg = 1
    temporary_stage0 = hlak
    endpoint1 = this%en1(n)
    endpoint2 = this%en2(n)
    call this%lak_calculate_residual(n, temporary_stage, residuala)
    if (hlak > endpoint1 .and. hlak < endpoint2) then
      endpoint2 = hlak
    end if
    do i = 1, this%maxlakit
      temporary_stage = dhalf * (endpoint1 + endpoint2)
      call this%lak_calculate_residual(n, temporary_stage, residual)
      if (abs(residual) == dzero .or. &
          abs(temporary_stage0 - temporary_stage) < this%dmaxchg) then
        exit
      end if
      call this%lak_calculate_residual(n, endpoint1, residuala)
      ! -- change end points
      ! -- root is between temporary_stage and endpoint2
      if (sign(done, residuala) == sign(done, residual)) then
        endpoint1 = temporary_stage
        ! -- root is between endpoint1 and temporary_stage
      else
        endpoint2 = temporary_stage
      end if
      temporary_stage0 = temporary_stage
    end do
    dh = hlak - temporary_stage
  end subroutine lak_bisection
 
  !> @brief Calculate the available volumetric rate for a lake given a passed
  !! stage
  !<
  subroutine lak_calculate_available(this, n, hlak, avail, &
                                     ra, ro, qinf, ex, headp)
    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: n
    real(DP), intent(in) :: hlak
    real(DP), intent(inout) :: avail
    real(DP), intent(inout) :: ra
    real(DP), intent(inout) :: ro
    real(DP), intent(inout) :: qinf
    real(DP), intent(inout) :: ex
    real(DP), intent(in), optional :: headp
    ! -- local
    integer(I4B) :: j
    integer(I4B) :: idry
    integer(I4B) :: igwfnode
    real(DP) :: hp
    real(DP) :: head
    real(DP) :: qlakgw
    real(DP) :: v0
    !
    ! -- set hp
    if (present(headp)) then
      hp = headp
    else
      hp = dzero
    end if
    !
    ! -- initialize
    avail = dzero
    !
    ! -- calculate the aquifer sources to the lake
    do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
      igwfnode = this%gwfconn(j)
      if (this%ibound(igwfnode) == 0) cycle
      head = this%xnew(igwfnode) + hp
      call this%lak_estimate_conn_exchange(1, n, j, idry, hlak, head, qlakgw, &
                                           avail)
    end do
    !
    ! -- add rainfall
    call this%lak_calculate_rainfall(n, hlak, ra)
    avail = avail + ra
    !
    ! -- calculate runoff
    call this%lak_calculate_runoff(n, ro)
    avail = avail + ro
    !
    ! -- calculate inflow
    call this%lak_calculate_inflow(n, qinf)
    avail = avail + qinf
    !
    ! -- calculate external flow terms
    call this%lak_calculate_external(n, ex)
    avail = avail + ex
    !
    ! -- calculate volume available in storage
    call this%lak_calculate_vol(n, this%xoldpak(n), v0)
    avail = avail + v0 / delt
  end subroutine lak_calculate_available
 
  !> @brief Calculate the residual for a lake given a passed stage
  !<
  subroutine lak_calculate_residual(this, n, hlak, resid, headp)
    ! -- modules
    use tdismodule, only: delt
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: n
    real(DP), intent(in) :: hlak
    real(DP), intent(inout) :: resid
    real(DP), intent(in), optional :: headp
    ! -- local
    integer(I4B) :: j
    integer(I4B) :: idry
    integer(I4B) :: igwfnode
    real(DP) :: hp
    real(DP) :: avail
    real(DP) :: head
    real(DP) :: ra
    real(DP) :: ro
    real(DP) :: qinf
    real(DP) :: ex
    real(DP) :: ev
    real(DP) :: wr
    real(DP) :: sout
    real(DP) :: sin
    real(DP) :: qlakgw
    real(DP) :: seep
    real(DP) :: hlak0
    real(DP) :: v0
    real(DP) :: v1
    !
    ! -- set hp
    if (present(headp)) then
      hp = headp
    else
      hp = dzero
    end if
    !
    ! -- initialize
    resid = dzero
    avail = dzero
    seep = dzero
    !
    ! -- calculate the available water
    call this%lak_calculate_available(n, hlak, avail, &
                                      ra, ro, qinf, ex, hp)
    !
    ! -- calculate groundwater seepage
    do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
      igwfnode = this%gwfconn(j)
      if (this%ibound(igwfnode) == 0) cycle
      head = this%xnew(igwfnode) + hp
      call this%lak_estimate_conn_exchange(2, n, j, idry, hlak, head, qlakgw, &
                                           avail)
      seep = seep + qlakgw
    end do
    !
    ! -- limit withdrawals to lake inflows and lake storage
    call this%lak_calculate_withdrawal(n, avail, wr)
    !
    ! -- limit evaporation to lake inflows and lake storage
    call this%lak_calculate_evaporation(n, hlak, avail, ev)
    !
    ! -- no outlet flow if evaporation consumes all water
    call this%lak_calculate_outlet_outflow(n, hlak, avail, sout)
    !
    ! -- update the surface inflow values
    call this%lak_calculate_outlet_inflow(n, sin)
    !
    ! -- calculate residual
    resid = ra + ev + wr + ro + qinf + ex + sin + sout + seep
    !
    ! -- include storage
    if (this%gwfiss /= 1) then
      hlak0 = this%xoldpak(n)
      call this%lak_calculate_vol(n, hlak0, v0)
      call this%lak_calculate_vol(n, hlak, v1)
      resid = resid + (v0 - v1) / delt
    end if
  end subroutine lak_calculate_residual
 
  !> @brief Set up the budget object that stores all the lake flows
  !<
  subroutine lak_setup_budobj(this)
    ! -- modules
    use constantsmodule, only: lenbudtxt
    ! -- dummy
    class(laktype) :: this
    ! -- local
    integer(I4B) :: nbudterm
    integer(I4B) :: nlen
    integer(I4B) :: j, n, n1, n2
    integer(I4B) :: maxlist, naux
    integer(I4B) :: idx
    real(DP) :: q
    character(len=LENBUDTXT) :: text
    character(len=LENBUDTXT), dimension(1) :: auxtxt
    !
    ! -- Determine the number of lake budget terms. These are fixed for
    !    the simulation and cannot change
    nbudterm = 9
    nlen = 0
    do n = 1, this%noutlets
      if (this%lakein(n) > 0 .and. this%lakeout(n) > 0) then
        nlen = nlen + 1
      end if
    end do
    if (nlen > 0) nbudterm = nbudterm + 1
    if (this%imover == 1) nbudterm = nbudterm + 2
    if (this%naux > 0) nbudterm = nbudterm + 1
    !
    ! -- set up budobj
    call budgetobject_cr(this%budobj, this%packName)
    call this%budobj%budgetobject_df(this%nlakes, nbudterm, 0, 0, &
                                     ibudcsv=this%ibudcsv)
    idx = 0
    !
    ! -- Go through and set up each budget term. nlen is the number
    !    of outlets that discharge into another lake
    if (nlen > 0) then
      text = '    FLOW-JA-FACE'
      idx = idx + 1
      maxlist = 2 * nlen
      naux = 0
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux, ordered_id1=.false.)
      !
      ! -- store connectivity
      call this%budobj%budterm(idx)%reset(2 * nlen)
      q = dzero
      do n = 1, this%noutlets
        n1 = this%lakein(n)
        n2 = this%lakeout(n)
        if (n1 > 0 .and. n2 > 0) then
          call this%budobj%budterm(idx)%update_term(n1, n2, q)
          call this%budobj%budterm(idx)%update_term(n2, n1, -q)
        end if
      end do
    end if
    !
    ! --
    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%nlakes
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        n2 = this%gwfconn(j)
        call this%budobj%budterm(idx)%update_term(n, n2, q)
      end do
    end do
    !
    ! --
    text = '        RAINFALL'
    idx = idx + 1
    maxlist = this%nlakes
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '     EVAPORATION'
    idx = idx + 1
    maxlist = this%nlakes
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '          RUNOFF'
    idx = idx + 1
    maxlist = this%nlakes
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '      EXT-INFLOW'
    idx = idx + 1
    maxlist = this%nlakes
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '      WITHDRAWAL'
    idx = idx + 1
    maxlist = this%nlakes
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '     EXT-OUTFLOW'
    idx = idx + 1
    maxlist = this%nlakes
    naux = 0
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux)
    !
    ! --
    text = '         STORAGE'
    idx = idx + 1
    maxlist = this%nlakes
    naux = 1
    auxtxt(1) = '          VOLUME'
    call this%budobj%budterm(idx)%initialize(text, &
                                             this%name_model, &
                                             this%packName, &
                                             this%name_model, &
                                             this%packName, &
                                             maxlist, .false., .false., &
                                             naux, auxtxt)
    !
    ! --
    text = '        CONSTANT'
    idx = idx + 1
    maxlist = this%nlakes
    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%nlakes
      naux = 0
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux)
      !
      ! --
      text = '          TO-MVR'
      idx = idx + 1
      maxlist = this%noutlets
      naux = 0
      call this%budobj%budterm(idx)%initialize(text, &
                                               this%name_model, &
                                               this%packName, &
                                               this%name_model, &
                                               this%packName, &
                                               maxlist, .false., .false., &
                                               naux, ordered_id1=.false.)
      !
      ! -- store to-mvr connection information
      call this%budobj%budterm(idx)%reset(this%noutlets)
      q = dzero
      do n = 1, this%noutlets
        n1 = this%lakein(n)
        call this%budobj%budterm(idx)%update_term(n1, n1, q)
      end do
    end if
    !
    ! --
    naux = this%naux
    if (naux > 0) then
      !
      ! --
      text = '       AUXILIARY'
      idx = idx + 1
      maxlist = this%nlakes
      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 lake 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 lak_setup_budobj
 
  !> @brief Copy flow terms into this%budobj
  !<
  subroutine lak_fill_budobj(this)
    ! -- dummy
    class(laktype) :: this
    ! -- local
    integer(I4B) :: naux
    real(DP), dimension(:), allocatable :: auxvartmp
    !integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: n
    integer(I4B) :: n1
    integer(I4B) :: n2
    integer(I4B) :: ii
    integer(I4B) :: jj
    integer(I4B) :: idx
    integer(I4B) :: nlen
    real(DP) :: v, v1
    real(DP) :: q
    real(DP) :: lkstg, gwhead, wa
    !
    ! -- initialize counter
    idx = 0
 
    ! -- FLOW JA FACE
    nlen = 0
    do n = 1, this%noutlets
      if (this%lakein(n) > 0 .and. this%lakeout(n) > 0) then
        nlen = nlen + 1
      end if
    end do
    if (nlen > 0) then
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(2 * nlen)
      do n = 1, this%noutlets
        n1 = this%lakein(n)
        n2 = this%lakeout(n)
        if (n1 > 0 .and. n2 > 0) then
          q = this%simoutrate(n)
          if (this%imover == 1) then
            q = q + this%pakmvrobj%get_qtomvr(n)
          end if
          call this%budobj%budterm(idx)%update_term(n1, n2, q)
          call this%budobj%budterm(idx)%update_term(n2, n1, -q)
        end if
      end do
    end if
    !
    ! -- GWF (LEAKAGE)
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%maxbound)
    do n = 1, this%nlakes
      do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
        n2 = this%gwfconn(j)
        q = this%qleak(j)
        lkstg = this%xnewpak(n)
        ! -- For the case when the lak stage is exactly equal
        !    to the lake bottom, the wetted area is not returned
        !    equal to 0.0
        gwhead = this%xnew(n2)
        call this%lak_calculate_conn_warea(n, j, lkstg, gwhead, wa)
        ! -- For thermal conduction between a lake and a gw cell,
        !    the shared wetted area should be reset to zero when the lake
        !    stage is below the cell bottom
        if (this%belev(j) > lkstg) wa = dzero
        this%qauxcbc(1) = wa
        call this%budobj%budterm(idx)%update_term(n, n2, q, this%qauxcbc)
      end do
    end do
    !
    ! -- RAIN
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nlakes)
    do n = 1, this%nlakes
      q = this%precip(n)
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- EVAPORATION
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nlakes)
    do n = 1, this%nlakes
      q = this%evap(n)
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- RUNOFF
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nlakes)
    do n = 1, this%nlakes
      q = this%runoff(n)
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- INFLOW
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nlakes)
    do n = 1, this%nlakes
      q = this%inflow(n)
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- WITHDRAWAL
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nlakes)
    do n = 1, this%nlakes
      q = this%withr(n)
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- EXTERNAL OUTFLOW
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nlakes)
    do n = 1, this%nlakes
      call this%lak_get_external_outlet(n, q)
      ! subtract tomover from external outflow
      call this%lak_get_external_mover(n, v)
      q = q + v
      call this%budobj%budterm(idx)%update_term(n, n, q)
    end do
    !
    ! -- STORAGE
    idx = idx + 1
    call this%budobj%budterm(idx)%reset(this%nlakes)
    do n = 1, this%nlakes
      call this%lak_calculate_vol(n, this%xnewpak(n), v1)
      q = this%qsto(n)
      this%qauxcbc(1) = v1
      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%nlakes)
    do n = 1, this%nlakes
      q = this%chterm(n)
      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%nlakes)
      do n = 1, this%nlakes
        q = this%pakmvrobj%get_qfrommvr(n)
        call this%budobj%budterm(idx)%update_term(n, n, q)
      end do
      !
      ! -- TO MOVER
      idx = idx + 1
      call this%budobj%budterm(idx)%reset(this%noutlets)
      do n = 1, this%noutlets
        n1 = this%lakein(n)
        q = this%pakmvrobj%get_qtomvr(n)
        if (q > dzero) then
          q = -q
        end if
        call this%budobj%budterm(idx)%update_term(n1, n1, q)
      end do
      !
    end if
    !
    ! -- AUXILIARY VARIABLES
    naux = this%naux
    if (naux > 0) then
      idx = idx + 1
      allocate (auxvartmp(naux))
      call this%budobj%budterm(idx)%reset(this%nlakes)
      do n = 1, this%nlakes
        q = dzero
        do jj = 1, naux
          ii = n
          auxvartmp(jj) = this%featureauxvar(jj, ii)
        end do
        call this%budobj%budterm(idx)%update_term(n, n, q, auxvartmp)
      end do
      deallocate (auxvartmp)
    end if
    !
    ! --Terms are filled, now accumulate them for this time step
    call this%budobj%accumulate_terms()
  end subroutine lak_fill_budobj
 
  !> @brief Set up the table object that is used to write the lak stage data
  !!
  !! The terms listed here must correspond in number and order to the ones
  !! written to the stage table in the lak_ot method
  !<
  subroutine lak_setup_tableobj(this)
    ! -- modules
    use constantsmodule, only: linelength, lenbudtxt
    ! -- dummy
    class(laktype) :: this
    ! -- local
    integer(I4B) :: nterms
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    !
    ! -- setup stage table
    if (this%iprhed > 0) then
      !
      ! -- Determine the number of lake stage terms. These are fixed for
      !    the simulation and cannot change.  This includes FLOW-JA-FACE
      !    so they can be written to the binary budget files, but these internal
      !    flows are not included as part of the budget table.
      nterms = 5
      if (this%inamedbound == 1) then
        nterms = nterms + 1
      end if
      !
      ! -- set up table title
      title = trim(adjustl(this%text))//' PACKAGE ('// &
              trim(adjustl(this%packName))//') STAGES FOR EACH CONTROL VOLUME'
      !
      ! -- set up stage tableobj
      call table_cr(this%stagetab, this%packName, title)
      call this%stagetab%table_df(this%nlakes, nterms, this%iout, &
                                  transient=.true.)
      !
      ! -- Go through and set up table budget term
      if (this%inamedbound == 1) then
        text = 'NAME'
        call this%stagetab%initialize_column(text, 20, alignment=tableft)
      end if
      !
      ! -- lake number
      text = 'NUMBER'
      call this%stagetab%initialize_column(text, 10, alignment=tabcenter)
      !
      ! -- lake stage
      text = 'STAGE'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
      !
      ! -- lake surface area
      text = 'SURFACE AREA'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
      !
      ! -- lake wetted area
      text = 'WETTED AREA'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
      !
      ! -- lake volume
      text = 'VOLUME'
      call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
    end if
  end subroutine lak_setup_tableobj
 
  !> @brief Activate addition of density terms
  !<
  subroutine lak_activate_density(this)
    ! -- dummy
    class(laktype), intent(inout) :: this
    ! -- local
    integer(I4B) :: i, j
    !
    ! -- the IMPLICIT formulation does not include the density terms yet
    if (this%iimplicit /= 0) then
      write (errmsg, '(a)') &
        'The IMPLICIT option cannot be used with the BUY package until the &
        &implicit formulation includes density terms. Remove the IMPLICIT &
        &option from LAK package '//trim(this%packName)//' to simulate &
        &density.'
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end if
    !
    ! -- 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 LAKE &
      &PACKAGE: '//trim(adjustl(this%packName))
  end subroutine lak_activate_density
 
  !> @brief Activate viscosity terms
  !!
  !! Method to activate addition of viscosity terms for a LAK package reach.
  !<
  subroutine lak_activate_viscosity(this)
    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    ! -- dummy variables
    class(laktype), intent(inout) :: this !< LakType 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 LAK &
      &PACKAGE: '//trim(adjustl(this%packName))
  end subroutine lak_activate_viscosity
 
  !> @brief Calculate the groundwater-lake density exchange terms
  !!
  !! Arguments are as follows:
  !!     iconn       : lak-gwf connection number
  !!     stage       : lake stage
  !!     head        : gwf head
  !!     cond        : conductance
  !!     botl        : bottom elevation of this connection
  !!     flow        : calculated flow, updated here with density terms
  !!     gwfhcof     : gwf head coefficient, updated here with density terms
  !!     gwfrhs      : gwf 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 lake (denselak / denseref)
  !!                     col 2 is relative density of gwf cell (densegwf / denseref)
  !!                     col 3 is elevation of gwf cell
  !<
  subroutine lak_calculate_density_exchange(this, iconn, stage, head, cond, &
                                            botl, flow, gwfhcof, gwfrhs)
    ! -- dummy
    class(laktype), intent(inout) :: this
    integer(I4B), intent(in) :: iconn
    real(DP), intent(in) :: stage
    real(DP), intent(in) :: head
    real(DP), intent(in) :: cond
    real(DP), intent(in) :: botl
    real(DP), intent(inout) :: flow
    real(DP), intent(inout) :: gwfhcof
    real(DP), intent(inout) :: gwfrhs
    ! -- local
    real(DP) :: ss
    real(DP) :: hh
    real(DP) :: havg
    real(DP) :: rdenselak
    real(DP) :: rdensegwf
    real(DP) :: rdenseavg
    real(DP) :: elevlak
    real(DP) :: elevgwf
    real(DP) :: elevavg
    real(DP) :: d1
    real(DP) :: d2
    logical(LGP) :: stage_below_bot
    logical(LGP) :: head_below_bot
    !
    ! -- Set lak density to lak density or gwf density
    if (stage >= botl) then
      ss = stage
      stage_below_bot = .false.
      rdenselak = this%denseterms(1, iconn) ! lak rel density
    else
      ss = botl
      stage_below_bot = .true.
      rdenselak = this%denseterms(2, iconn) ! gwf rel density
    end if
    !
    ! -- set hh to head or botl
    if (head >= botl) then
      hh = head
      head_below_bot = .false.
      rdensegwf = this%denseterms(2, iconn) ! gwf rel density
    else
      hh = botl
      head_below_bot = .true.
      rdensegwf = this%denseterms(1, iconn) ! lak rel density
    end if
    !
    ! -- todo: hack because denseterms not updated in a cf calculation
    if (rdensegwf == dzero) return
    !
    ! -- Update flow
    if (stage_below_bot .and. head_below_bot) then
      !
      ! -- flow is zero, so no terms are updated
      !
    else
      !
      ! -- calculate average relative density
      rdenseavg = dhalf * (rdenselak + rdensegwf)
      !
      ! -- Add contribution of first density term:
      !      cond * (denseavg/denseref - 1) * (hgwf - hlak)
      d1 = cond * (rdenseavg - done)
      gwfhcof = gwfhcof - d1
      gwfrhs = gwfrhs - d1 * ss
      d1 = d1 * (hh - ss)
      flow = flow + d1
      !
      ! -- Add second density term if stage and head not below bottom
      if (.not. stage_below_bot .and. .not. head_below_bot) then
        !
        ! -- Add contribution of second density term:
        !      cond * (havg - elevavg) * (densegwf - denselak) / denseref
        elevgwf = this%denseterms(3, iconn)
        if (this%ictype(iconn) == 0 .or. this%ictype(iconn) == 3) then
          ! -- vertical or embedded vertical connection
          elevlak = botl
        else
          ! -- horizontal or embedded horizontal connection
          elevlak = elevgwf
        end if
        elevavg = dhalf * (elevlak + elevgwf)
        havg = dhalf * (hh + ss)
        d2 = cond * (havg - elevavg) * (rdensegwf - rdenselak)
        gwfrhs = gwfrhs + d2
        flow = flow + d2
      end if
    end if
  end subroutine lak_calculate_density_exchange
 
end module lakmodule