lakmodule Module Reference

Data Types
type	laktabtype
type	lakinputtype
	input context runtime pointers for LAK More...
type	laktype

Functions/Subroutines
subroutine, public	lak_create (packobj, id, ibcnum, inunit, iout, namemodel, pakname, mempath)
	Create a new LAK Package and point bndobj to the new package. More...
subroutine	lak_allocate_scalars (this)
	Allocate scalar members. More...
subroutine	lak_allocate_arrays (this)
	Allocate scalar members. More...
subroutine	lak_source_packagedata (this)
	Source PACKAGEDATA block from input context. More...
subroutine	lak_source_connectiondata (this)
	Source CONNECTIONDATA block from input context. More...
subroutine	lak_source_tables (this)
	Source TABLES block from input context. More...
subroutine	laktables_to_vectors (this, laketables)
	Copy the laketables structure data into flattened vectors that are stored in the memory manager. More...
subroutine	lak_source_table (this, ilak, mempath, laketable)
	Source one lake's table from its LAKTAB6 subpackage. More...
subroutine	lak_source_outlets (this)
	Source OUTLETS block from input context. More...
subroutine	lak_source_dimensions (this)
	Source DIMENSIONS block from input context and load static data. More...
subroutine	lak_read_initial_attr (this)
	Read the initial parameters for this package. More...
subroutine	lak_linear_interpolation (this, n, x, y, z, v)
	Perform linear interpolation of two vectors. More...
subroutine	lak_calculate_sarea (this, ilak, stage, sarea)
	Calculate the surface area of a lake at a given stage. More...
subroutine	lak_calculate_warea (this, ilak, stage, warea, hin)
	Calculate the wetted area of a lake at a given stage. More...
subroutine	lak_calculate_conn_warea (this, ilak, iconn, stage, head, wa)
	Calculate the wetted area of a lake connection at a given stage. More...
subroutine	lak_calculate_vol (this, ilak, stage, volume)
	Calculate the volume of a lake at a given stage. More...
subroutine	lak_calculate_conductance (this, ilak, stage, conductance)
	Calculate the total conductance for a lake at a provided stage. More...
subroutine	lak_calculate_cond_head (this, iconn, stage, head, vv)
	Calculate the controlling lake stage or groundwater head used to calculate the conductance for a lake connection from a provided stage and groundwater head. More...
subroutine	lak_calculate_conn_conductance (this, ilak, iconn, stage, head, cond)
	Calculate the conductance for a lake connection at a provided stage and groundwater head. More...
subroutine	lak_calculate_exchange (this, ilak, stage, totflow)
	Calculate the total groundwater-lake flow at a provided stage. More...
subroutine	lak_calculate_conn_exchange (this, ilak, iconn, stage, head, flow, gwfhcof, gwfrhs)
	Calculate the groundwater-lake flow at a provided stage and groundwater head. More...
subroutine	lak_calculate_conn_exchange_deriv (this, ilak, iconn, stage, head, flow, dqds, dqdh)
	Lakebed seepage and its derivatives for a connection (IMPLICIT) More...
subroutine	lak_estimate_conn_exchange (this, iflag, ilak, iconn, idry, stage, head, flow, source, gwfhcof, gwfrhs)
	Calculate the groundwater-lake flow at a provided stage and groundwater head. More...
subroutine	lak_calculate_storagechange (this, ilak, stage, stage0, delt, dvr)
	Calculate the storage change in a lake based on provided stages and a passed delt. More...
subroutine	lak_calculate_rainfall (this, ilak, stage, ra)
	Calculate the rainfall for a lake. More...
subroutine	lak_calculate_runoff (this, ilak, ro)
	Calculate runoff to a lake. More...
subroutine	lak_calculate_inflow (this, ilak, qin)
	Calculate specified inflow to a lake. More...
subroutine	lak_calculate_external (this, ilak, ex)
	Calculate the external flow terms to a lake. More...
subroutine	lak_calculate_withdrawal (this, ilak, avail, wr)
	Calculate the withdrawal from a lake subject to an available volume. More...
subroutine	lak_calculate_evaporation (this, ilak, stage, avail, ev)
	Calculate the evaporation from a lake at a provided stage subject to an available volume. More...
subroutine	lak_calculate_outlet_inflow (this, ilak, outinf)
	Calculate the outlet inflow to a lake. More...
subroutine	lak_calculate_outlet_outflow (this, ilak, stage, avail, outoutf)
	Calculate the outlet outflow from a lake. More...
subroutine	lak_outlet_outflow_rate (this, ilak, stage, qout)
	Total uncapped outlet outflow rate from a lake at a provided stage. More...
subroutine	lak_get_internal_inlet (this, ilak, outinf)
	Get the outlet inflow to a lake from another lake. More...
subroutine	lak_get_internal_outlet (this, ilak, outoutf)
	Get the outlet from a lake to another lake. More...
subroutine	lak_get_external_outlet (this, ilak, outoutf)
	Get the outlet outflow from a lake to an external boundary. More...
subroutine	lak_get_external_mover (this, ilak, outoutf)
	Get the mover outflow from a lake to an external boundary. More...
subroutine	lak_get_internal_mover (this, ilak, outoutf)
	Get the mover outflow from a lake to another lake. More...
subroutine	lak_get_outlet_tomover (this, ilak, outoutf)
	Get the outlet to mover from a lake. More...
subroutine	lak_vol2stage (this, ilak, vol, stage)
	Determine the stage from a provided volume. More...
integer(i4b) function	lak_check_valid (this, itemno)
	Determine if a valid lake or outlet number has been specified. More...
subroutine	lak_set_nonneg_value (this, n, itemno, setting, value_array, label)
	Apply a PERIOD lake setting and require its value be non-negative. More...
subroutine	lak_set_attribute_error (this, ilak, keyword, msg)
	Issue a parameter error for lake ilak. More...
subroutine	lak_source_options (this)
	Source OPTIONS block from input context. More...
subroutine	lak_input_init (this, mempath, naux, noutlets)
	Bind LakInputType pointers. More...
subroutine	lak_input_destroy (this)
	Nullify all LakInputType pointers. More...
subroutine	lak_ar (this)
	Allocate and Read. More...
subroutine	lak_rp (this)
	Read and prepare (stress period data) More...
subroutine	lak_ad (this)
	Add package connection to matrix. More...
subroutine	lak_cf (this)
	Formulate the HCOF and RHS terms. More...
subroutine	lak_fc (this, rhs, ia, idxglo, matrix_sln)
	Copy rhs and hcof into solution rhs and amat. More...
subroutine	lak_fn (this, rhs, ia, idxglo, matrix_sln)
	Fill newton terms. More...
subroutine	lak_nur (this, neqpak, x, xtemp, dx, inewtonur, dxmax, locmax)
	Apply Newton under-relaxation to the lake stage. More...
subroutine	lak_cc (this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
	Final convergence check for package. More...
subroutine	lak_cq (this, x, flowja, iadv)
	Calculate flows. More...
subroutine	lak_ot_package_flows (this, icbcfl, ibudfl)
	Output LAK package flow terms. More...
subroutine	lak_ot_model_flows (this, icbcfl, ibudfl, icbcun, imap)
	Write flows to binary file and/or print flows to budget. More...
subroutine	lak_ot_dv (this, idvsave, idvprint)
	Save LAK-calculated values to binary file. More...
subroutine	lak_ot_bdsummary (this, kstp, kper, iout, ibudfl)
	Write LAK budget to listing file. More...
subroutine	lak_da (this)
	Deallocate objects. More...
subroutine	define_listlabel (this)
	Define the list heading that is written to iout when PRINT_INPUT option is used. More...
subroutine	lak_set_pointers (this, neq, ibound, xnew, xold, flowja)
	Set pointers to model arrays and variables so that a package has access to these things. More...
subroutine	lak_ac (this, moffset, sparse)
	Add the lake rows and columns to the sparse matrix. More...
subroutine	lak_mc (this, moffset, matrix_sln)
	Find the matrix position of each lake row and connection. More...
logical function	lak_obs_supported (this)
	Procedures related to observations (type-bound) More...
subroutine	lak_df_obs (this)
	Store observation type supported by LAK package. Overrides BndTypebnd_df_obs. More...
subroutine	lak_bd_obs (this)
	Calculate observations this time step and call ObsTypeSaveOneSimval for each LakType observation. More...
subroutine	lak_rp_obs (this)
	Process each observation. More...
subroutine	lak_process_obsid (obsrv, dis, inunitobs, iout)
	This procedure is pointed to by ObsDataTypeProcesssIdPtr. It processes the ID string of an observation definition for LAK package observations. More...
subroutine	lak_accumulate_chterm (this, ilak, rrate, chratin, chratout)
	Accumulate constant head terms for budget. More...
subroutine	lak_bound_update (this)
	Store the lake head and connection conductance in the bound array. More...
subroutine	lak_solve (this, update, only_legacy)
	Solve for lake stage. More...
subroutine	lak_solve_single (this, n, iter, maxiter, ncnv, lupdate)
	Advance one lake stage by a single substitution iteration. More...
subroutine	lak_estimate_seepage_single (this, n, ncnv)
	Estimate the lakebed seepage for a single lake. More...
subroutine	lak_bisection (this, n, ibflg, hlak, temporary_stage, dh, residual)
	@ brief Lake package bisection method More...
subroutine	lak_calculate_available (this, n, hlak, avail, ra, ro, qinf, ex, headp)
	Calculate the available volumetric rate for a lake given a passed stage. More...
subroutine	lak_calculate_residual (this, n, hlak, resid, headp)
	Calculate the residual for a lake given a passed stage. More...
subroutine	lak_setup_budobj (this)
	Set up the budget object that stores all the lake flows. More...
subroutine	lak_fill_budobj (this)
	Copy flow terms into thisbudobj. More...
subroutine	lak_setup_tableobj (this)
	Set up the table object that is used to write the lak stage data. More...
subroutine	lak_activate_density (this)
	Activate addition of density terms. More...
subroutine	lak_activate_viscosity (this)
	Activate viscosity terms. More...
subroutine	lak_calculate_density_exchange (this, iconn, stage, head, cond, botl, flow, gwfhcof, gwfrhs)
	Calculate the groundwater-lake density exchange terms. More...

Variables
character(len=lenftype)	ftype = 'LAK'
character(len=lenpackagename)	text = ' LAK'

Function/Subroutine Documentation

define_listlabel()

subroutine lakmodule::define_listlabel

(

class(laktype), intent(inout)

this

)

Definition at line 4482 of file gwf-lak.f90.

    ! -- 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

lak_ac()

subroutine lakmodule::lak_ac	(	class(laktype), intent(inout)	this,
		integer(i4b), intent(in)	moffset,
		type(sparsematrix), intent(inout)	sparse
	)

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.

Definition at line 4549 of file gwf-lak.f90.

    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

lak_accumulate_chterm()

subroutine lakmodule::lak_accumulate_chterm ( class(laktype)  this, integer(i4b), intent(in)  ilak, real(dp), intent(in)  rrate, real(dp), intent(inout)  chratin, real(dp), intent(inout)  chratout  )

private

Definition at line 5124 of file gwf-lak.f90.

    ! -- 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

lak_activate_density()

subroutine lakmodule::lak_activate_density

(

class(laktype), intent(inout)

this

)

Definition at line 6337 of file gwf-lak.f90.

    ! -- 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))

Here is the call graph for this function:

lak_activate_viscosity()

subroutine lakmodule::lak_activate_viscosity ( class(laktype), intent(inout)  this )

private

Method to activate addition of viscosity terms for a LAK package reach.

Parameters

[in,out]

this

LakType object

Definition at line 6371 of file gwf-lak.f90.

    ! -- 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))

lak_ad()

subroutine lakmodule::lak_ad

(

class(laktype)

this

)

Definition at line 3401 of file gwf-lak.f90.

    ! -- 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

lak_allocate_arrays()

subroutine lakmodule::lak_allocate_arrays ( class(laktype), intent(inout)  this )

private

Definition at line 498 of file gwf-lak.f90.

    ! -- 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)

lak_allocate_scalars()

subroutine lakmodule::lak_allocate_scalars

(

class(laktype), intent(inout)

this

)

Definition at line 436 of file gwf-lak.f90.

    ! -- 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

lak_ar()

subroutine lakmodule::lak_ar ( class(laktype), intent(inout)  this )

private

Create new LAK package and point bndobj to the new package

Definition at line 3128 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_bd_obs()

subroutine lakmodule::lak_bd_obs ( class(laktype)  this )

private

Definition at line 4765 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_bisection()

subroutine lakmodule::lak_bisection ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  n, integer(i4b), intent(inout)  ibflg, real(dp), intent(in)  hlak, real(dp), intent(inout)  temporary_stage, real(dp), intent(inout)  dh, real(dp), intent(inout)  residual  )

private

Use bisection method to find lake stage that reduces the residual

Parameters

[in]	n	lake number
[in,out]	ibflg	bisection flag
[in]	hlak	lake stage
[in,out]	temporary_stage	temporary lake stage
[in,out]	dh	lake stage change
[in,out]	residual	lake residual

Definition at line 5641 of file gwf-lak.f90.

    ! -- 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

lak_bound_update()

subroutine lakmodule::lak_bound_update ( class(laktype), intent(inout)  this )

private

Definition at line 5154 of file gwf-lak.f90.

    ! -- 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

lak_calculate_available()

subroutine lakmodule::lak_calculate_available ( 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  )

private

Definition at line 5689 of file gwf-lak.f90.

    ! -- 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

lak_calculate_cond_head()

subroutine lakmodule::lak_calculate_cond_head ( 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  )

private

Definition at line 2093 of file gwf-lak.f90.

    ! -- 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

lak_calculate_conductance()

subroutine lakmodule::lak_calculate_conductance ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(in)  stage, real(dp), intent(inout)  conductance  )

private

Definition at line 2072 of file gwf-lak.f90.

    ! -- 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

lak_calculate_conn_conductance()

subroutine lakmodule::lak_calculate_conn_conductance ( 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  )

private

Definition at line 2122 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_calculate_conn_exchange()

subroutine lakmodule::lak_calculate_conn_exchange ( 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  )

private

Definition at line 2212 of file gwf-lak.f90.

    ! -- 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

lak_calculate_conn_exchange_deriv()

subroutine lakmodule::lak_calculate_conn_exchange_deriv ( 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  )

private

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.

Definition at line 2281 of file gwf-lak.f90.

    ! -- 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

lak_calculate_conn_warea()

subroutine lakmodule::lak_calculate_conn_warea ( 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  )

private

Definition at line 1968 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_calculate_density_exchange()

subroutine lakmodule::lak_calculate_density_exchange	(	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
	)

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

Definition at line 6411 of file gwf-lak.f90.

    ! -- 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

lak_calculate_evaporation()

subroutine lakmodule::lak_calculate_evaporation ( 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  )

private

Definition at line 2468 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_calculate_exchange()

subroutine lakmodule::lak_calculate_exchange ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(in)  stage, real(dp), intent(inout)  totflow  )

private

Definition at line 2188 of file gwf-lak.f90.

    ! -- 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

lak_calculate_external()

subroutine lakmodule::lak_calculate_external ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  ex  )

private

Definition at line 2430 of file gwf-lak.f90.

    ! -- 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

lak_calculate_inflow()

subroutine lakmodule::lak_calculate_inflow ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  qin  )

private

Definition at line 2418 of file gwf-lak.f90.

    ! -- dummy
    class(LakType), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: qin
    !
    ! -- inflow to lake
    qin = this%inflow(ilak)

lak_calculate_outlet_inflow()

subroutine lakmodule::lak_calculate_outlet_inflow ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  outinf  )

private

Definition at line 2495 of file gwf-lak.f90.

    ! -- 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

lak_calculate_outlet_outflow()

subroutine lakmodule::lak_calculate_outlet_outflow ( 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  )

private

Definition at line 2516 of file gwf-lak.f90.

    ! -- 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

lak_calculate_rainfall()

subroutine lakmodule::lak_calculate_rainfall ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(in)  stage, real(dp), intent(inout)  ra  )

private

Definition at line 2384 of file gwf-lak.f90.

    ! -- 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

lak_calculate_residual()

subroutine lakmodule::lak_calculate_residual	(	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
	)

Definition at line 5754 of file gwf-lak.f90.

    ! -- 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

lak_calculate_runoff()

subroutine lakmodule::lak_calculate_runoff ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  ro  )

private

Definition at line 2406 of file gwf-lak.f90.

    ! -- dummy
    class(LakType), intent(inout) :: this
    integer(I4B), intent(in) :: ilak
    real(DP), intent(inout) :: ro
    !
    ! -- runoff
    ro = this%runoff(ilak)

lak_calculate_sarea()

subroutine lakmodule::lak_calculate_sarea ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(in)  stage, real(dp), intent(inout)  sarea  )

private

Definition at line 1898 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_calculate_storagechange()

subroutine lakmodule::lak_calculate_storagechange ( 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  )

private

Definition at line 2362 of file gwf-lak.f90.

    ! -- 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

lak_calculate_vol()

subroutine lakmodule::lak_calculate_vol ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(in)  stage, real(dp), intent(inout)  volume  )

private

Definition at line 2019 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_calculate_warea()

subroutine lakmodule::lak_calculate_warea ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(in)  stage, real(dp), intent(inout)  warea, real(dp), intent(inout), optional  hin  )

private

Definition at line 1940 of file gwf-lak.f90.

    ! -- 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

lak_calculate_withdrawal()

subroutine lakmodule::lak_calculate_withdrawal ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  avail, real(dp), intent(inout)  wr  )

private

Definition at line 2446 of file gwf-lak.f90.

    ! -- 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

lak_cc()

subroutine lakmodule::lak_cc ( 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  )

private

Definition at line 3721 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_cf()

subroutine lakmodule::lak_cf

(

class(laktype)

this

)

Skip if no lakes, otherwise calculate hcof and rhs

Definition at line 3484 of file gwf-lak.f90.

    ! -- 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()

lak_check_valid()

integer(i4b) function lakmodule::lak_check_valid ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  itemno  )

private

Definition at line 2821 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_cq()

subroutine lakmodule::lak_cq	(	class(laktype), intent(inout)	this,
		real(dp), dimension(:), intent(in)	x,
		real(dp), dimension(:), intent(inout), contiguous	flowja,
		integer(i4b), intent(in), optional	iadv
	)

Definition at line 4033 of file gwf-lak.f90.

    ! -- 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()

Here is the call graph for this function:

lak_create()

subroutine, public lakmodule::lak_create	(	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
	)

Definition at line 394 of file gwf-lak.f90.

    ! -- 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')

Here is the call graph for this function:

Here is the caller graph for this function:

lak_da()

subroutine lakmodule::lak_da

(

class(laktype)

this

)

Definition at line 4301 of file gwf-lak.f90.

    ! -- 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()

lak_df_obs()

subroutine lakmodule::lak_df_obs ( class(laktype)  this )

private

Definition at line 4660 of file gwf-lak.f90.

    ! -- 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

lak_estimate_conn_exchange()

subroutine lakmodule::lak_estimate_conn_exchange ( 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  )

private

Definition at line 2319 of file gwf-lak.f90.

    ! -- 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

lak_estimate_seepage_single()

subroutine lakmodule::lak_estimate_seepage_single	(	class(laktype), intent(inout)	this,
		integer(i4b), intent(in)	n,
		integer(i4b), intent(in)	ncnv
	)

Runs the two-pass connection-seepage estimate for lake n at its current stage and the perturbed stage, accumulating thisseep(n)/thisseep1(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.

Definition at line 5571 of file gwf-lak.f90.

    ! -- 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

lak_fc()

subroutine lakmodule::lak_fc ( 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  )

private

Definition at line 3565 of file gwf-lak.f90.

    ! -- 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

lak_fill_budobj()

subroutine lakmodule::lak_fill_budobj

(

class(laktype)

this

)

Definition at line 6088 of file gwf-lak.f90.

    ! -- 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()

lak_fn()

subroutine lakmodule::lak_fn ( 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  )

private

Definition at line 3606 of file gwf-lak.f90.

    ! -- 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

lak_get_external_mover()

subroutine lakmodule::lak_get_external_mover ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  outoutf  )

private

Definition at line 2678 of file gwf-lak.f90.

    ! -- 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

lak_get_external_outlet()

subroutine lakmodule::lak_get_external_outlet ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  outoutf  )

private

Definition at line 2659 of file gwf-lak.f90.

    ! -- 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

lak_get_internal_inlet()

subroutine lakmodule::lak_get_internal_inlet ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  outinf  )

private

Definition at line 2619 of file gwf-lak.f90.

    ! -- 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

lak_get_internal_mover()

subroutine lakmodule::lak_get_internal_mover ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  outoutf  )

private

Definition at line 2699 of file gwf-lak.f90.

    ! -- 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

lak_get_internal_outlet()

subroutine lakmodule::lak_get_internal_outlet ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  outoutf  )

private

Definition at line 2640 of file gwf-lak.f90.

    ! -- 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

lak_get_outlet_tomover()

subroutine lakmodule::lak_get_outlet_tomover ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(inout)  outoutf  )

private

Definition at line 2720 of file gwf-lak.f90.

    ! -- 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

lak_input_destroy()

subroutine lakmodule::lak_input_destroy

(

class(lakinputtype), intent(inout)

this

)

Definition at line 3107 of file gwf-lak.f90.

    ! -- 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)

lak_input_init()

subroutine lakmodule::lak_input_init ( class(lakinputtype), intent(inout)  this, character(len=*), intent(in)  mempath, integer(i4b), intent(in)  naux, integer(i4b), intent(in)  noutlets  )

private

Called once from lak_ar; pointers remain valid for the simulation lifetime. Input tsmanager advances TS-linked values transparently each timestep.

Definition at line 3078 of file gwf-lak.f90.

    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

lak_linear_interpolation()

subroutine lakmodule::lak_linear_interpolation	(	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
	)

Function assumes x data is sorted in ascending order

Definition at line 1853 of file gwf-lak.f90.

    ! -- 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

lak_mc()

subroutine lakmodule::lak_mc	(	class(laktype), intent(inout)	this,
		integer(i4b), intent(in)	moffset,
		class(matrixbasetype), pointer	matrix_sln
	)

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.

Definition at line 4583 of file gwf-lak.f90.

    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

lak_nur()

subroutine lakmodule::lak_nur ( 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  )

private

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.

Definition at line 3679 of file gwf-lak.f90.

    ! -- 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

lak_obs_supported()

logical function lakmodule::lak_obs_supported

(

class(laktype)

this

)

Return true because LAK package supports observations. Overrides BndTypebnd_obs_supported()

Definition at line 4650 of file gwf-lak.f90.

    ! -- dummy
    class(LakType) :: this
    !
    lak_obs_supported = .true.

lak_ot_bdsummary()

subroutine lakmodule::lak_ot_bdsummary	(	class(laktype)	this,
		integer(i4b), intent(in)	kstp,
		integer(i4b), intent(in)	kper,
		integer(i4b), intent(in)	iout,
		integer(i4b), intent(in)	ibudfl
	)

Parameters

	this	LakType object
[in]	kstp	time step number
[in]	kper	period number
[in]	iout	flag and unit number for the model listing file
[in]	ibudfl	flag indicating budget should be written

Definition at line 4286 of file gwf-lak.f90.

    ! -- 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)

lak_ot_dv()

subroutine lakmodule::lak_ot_dv ( class(laktype)  this, integer(i4b), intent(in)  idvsave, integer(i4b), intent(in)  idvprint  )

private

Definition at line 4215 of file gwf-lak.f90.

    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

Here is the call graph for this function:

lak_ot_model_flows()

subroutine lakmodule::lak_ot_model_flows	(	class(laktype)	this,
		integer(i4b), intent(in)	icbcfl,
		integer(i4b), intent(in)	ibudfl,
		integer(i4b), intent(in)	icbcun,
		integer(i4b), dimension(:), intent(in), optional	imap
	)

Definition at line 4202 of file gwf-lak.f90.

    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)

lak_ot_package_flows()

subroutine lakmodule::lak_ot_package_flows	(	class(laktype)	this,
		integer(i4b), intent(in)	icbcfl,
		integer(i4b), intent(in)	ibudfl
	)

Definition at line 4176 of file gwf-lak.f90.

    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

lak_outlet_outflow_rate()

subroutine lakmodule::lak_outlet_outflow_rate ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(in)  stage, real(dp), intent(inout)  qout  )

private

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.

Definition at line 2580 of file gwf-lak.f90.

    ! -- 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

lak_process_obsid()

subroutine lakmodule::lak_process_obsid	(	type(observetype), intent(inout)	obsrv,
		class(disbasetype), intent(in)	dis,
		integer(i4b), intent(in)	inunitobs,
		integer(i4b), intent(in)	iout
	)

Definition at line 5072 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_read_initial_attr()

subroutine lakmodule::lak_read_initial_attr ( class(laktype), intent(inout)  this )

private

Definition at line 1591 of file gwf-lak.f90.

    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)

Here is the call graph for this function:

lak_rp()

subroutine lakmodule::lak_rp

(

class(laktype), intent(inout)

this

)

Definition at line 3181 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_rp_obs()

subroutine lakmodule::lak_rp_obs ( class(laktype), intent(inout)  this )

private

Only done the first stress period since boundaries are fixed for the simulation

Definition at line 4920 of file gwf-lak.f90.

    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

Here is the call graph for this function:

lak_set_attribute_error()

subroutine lakmodule::lak_set_attribute_error ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, character(len=*), intent(in)  keyword, character(len=*), intent(in)  msg  )

private

Definition at line 2878 of file gwf-lak.f90.

    ! -- 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)

Here is the call graph for this function:

lak_set_nonneg_value()

subroutine lakmodule::lak_set_nonneg_value	(	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
	)

Definition at line 2854 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_set_pointers()

subroutine lakmodule::lak_set_pointers ( 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  )

private

Definition at line 4507 of file gwf-lak.f90.

    ! -- 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

lak_setup_budobj()

subroutine lakmodule::lak_setup_budobj

(

class(laktype)

this

)

Definition at line 5836 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_setup_tableobj()

subroutine lakmodule::lak_setup_tableobj ( class(laktype)  this )

private

The terms listed here must correspond in number and order to the ones written to the stage table in the lak_ot method

Definition at line 6276 of file gwf-lak.f90.

    ! -- 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

Here is the call graph for this function:

lak_solve()

subroutine lakmodule::lak_solve ( class(laktype), intent(inout)  this, logical(lgp), intent(in), optional  update, logical(lgp), intent(in), optional  only_legacy  )

private

Solve the lake stage by substitution. With only_legacy set, solve only the lakes flagged for the legacy solver (thisilegacy /= 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.

Definition at line 5185 of file gwf-lak.f90.

    ! -- 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

lak_solve_single()

subroutine lakmodule::lak_solve_single	(	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
	)

Performs one Newton update (with a bisection backup) of lake n's stage in the default substitution solver, using the per-iteration seepage estimate (thisseep / thisseep1) 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.

Definition at line 5357 of file gwf-lak.f90.

    ! -- 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

lak_source_connectiondata()

subroutine lakmodule::lak_source_connectiondata ( class(laktype), intent(inout)  this )

private

Definition at line 776 of file gwf-lak.f90.

    ! -- 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)

Here is the call graph for this function:

lak_source_dimensions()

subroutine lakmodule::lak_source_dimensions ( class(laktype), intent(inout)  this )

private

Definition at line 1504 of file gwf-lak.f90.

    ! -- 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()

Here is the call graph for this function:

lak_source_options()

subroutine lakmodule::lak_source_options

(

class(laktype), intent(inout)

this

)

Definition at line 2898 of file gwf-lak.f90.

    ! -- 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'

Here is the call graph for this function:

lak_source_outlets()

subroutine lakmodule::lak_source_outlets ( class(laktype), intent(inout)  this )

private

Definition at line 1370 of file gwf-lak.f90.

    ! -- 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)

Here is the call graph for this function:

lak_source_packagedata()

subroutine lakmodule::lak_source_packagedata ( class(laktype), intent(inout)  this )

private

Definition at line 575 of file gwf-lak.f90.

    ! -- 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)

Here is the call graph for this function:

lak_source_table()

subroutine lakmodule::lak_source_table ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, character(len=*), intent(in)  mempath, type(laktabtype), intent(inout)  laketable  )

private

Definition at line 1219 of file gwf-lak.f90.

    ! -- 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)

Here is the call graph for this function:

lak_source_tables()

subroutine lakmodule::lak_source_tables ( class(laktype), intent(inout)  this )

private

Definition at line 1079 of file gwf-lak.f90.

    ! -- 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)

Here is the call graph for this function:

lak_vol2stage()

subroutine lakmodule::lak_vol2stage ( class(laktype), intent(inout)  this, integer(i4b), intent(in)  ilak, real(dp), intent(in)  vol, real(dp), intent(inout)  stage  )

private

Definition at line 2740 of file gwf-lak.f90.

    ! -- 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

laktables_to_vectors()

subroutine lakmodule::laktables_to_vectors ( class(laktype), intent(inout)  this, type(laktabtype), dimension(:), intent(in), contiguous  laketables  )

private

Definition at line 1171 of file gwf-lak.f90.

    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

Variable Documentation

ftype

character(len=lenftype) lakmodule::ftype = 'LAK'

private

Definition at line 46 of file gwf-lak.f90.

  character(len=LENFTYPE) :: ftype = 'LAK'

text

character(len=lenpackagename) lakmodule::text = ' LAK'

private

Definition at line 47 of file gwf-lak.f90.

  character(len=LENPACKAGENAME) :: text = '             LAK'