MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
gwf-lak.f90
Go to the documentation of this file.
1 module lakmodule
2  !
3  use kindmodule, only: dp, i4b, lgp
4  use constantsmodule, only: linelength, lenboundname, lentimeseriesname, &
5  iwetlake, maxadpit, &
6  dzero, dprec, dem30, dem9, dem6, dem5, &
7  dem4, dem2, dem1, dhalf, dp7, dp999, done, &
8  dtwo, dpi, dthree, deight, dten, dhundred, dep20, &
9  donethird, dtwothirds, dfivethirds, &
10  dgravity, dcd, &
11  namedboundflag, lenftype, lenpackagename, &
12  lenpakloc, dnodata, &
13  tableft, tabcenter, tabright, &
14  tabstring, tabucstring, tabinteger, tabreal
15  use memorymanagermodule, only: mem_allocate, mem_reallocate, mem_setptr, &
16  mem_deallocate
17  use memoryhelpermodule, only: create_mem_path
18  use smoothingmodule, only: squadraticsaturation, sqsaturation, &
19  squadraticsaturationderivative, &
20  sqsaturationderivative
21  use bndmodule, only: bndtype
22  use budgetobjectmodule, only: budgetobjecttype, budgetobject_cr
23  use tablemodule, only: tabletype, table_cr
24  use observemodule, only: observetype
25  use obsmodule, only: obstype
26  use geomutilmodule, only: get_node
27  use inputoutputmodule, only: urword, extract_idnum_or_bndname
28  use basedismodule, only: disbasetype
29  use simmodule, only: count_errors, store_error, store_error_unit, &
30  deprecation_warning
31  use mathutilmodule, only: is_close
32  use blockparsermodule, only: blockparsertype
33  use basedismodule, only: disbasetype
34  use simvariablesmodule, only: errmsg, warnmsg
35  use matrixbasemodule
36  !
37  implicit none
38  !
39  private
40  public :: laktype
41  public :: lak_create
42  !
43  character(len=LENFTYPE) :: ftype = 'LAK'
44  character(len=LENPACKAGENAME) :: text = ' LAK'
45  !
46  type laktabtype
47  real(dp), dimension(:), pointer, contiguous :: tabstage => null()
48  real(dp), dimension(:), pointer, contiguous :: tabvolume => null()
49  real(dp), dimension(:), pointer, contiguous :: tabsarea => null()
50  real(dp), dimension(:), pointer, contiguous :: tabwarea => null()
51  end type laktabtype
52  !
53  type, extends(bndtype) :: laktype
54  ! -- scalars
55  ! -- characters
56  character(len=16), dimension(:), pointer, contiguous :: clakbudget => null()
57  character(len=16), dimension(:), pointer, contiguous :: cauxcbc => null()
58  ! -- control variables
59  ! -- integers
60  integer(I4B), pointer :: iprhed => null()
61  integer(I4B), pointer :: istageout => null()
62  integer(I4B), pointer :: ibudgetout => null()
63  integer(I4B), pointer :: ibudcsv => null()
64  integer(I4B), pointer :: ipakcsv => null()
65  integer(I4B), pointer :: cbcauxitems => null()
66  integer(I4B), pointer :: nlakes => null()
67  integer(I4B), pointer :: noutlets => null()
68  integer(I4B), pointer :: ntables => null()
69  real(dp), pointer :: convlength => null()
70  real(dp), pointer :: convtime => null()
71  real(dp), pointer :: outdmax => null()
72  integer(I4B), pointer :: igwhcopt => null()
73  integer(I4B), pointer :: iconvchk => null()
74  integer(I4B), pointer :: iconvresidchk => null()
75  integer(I4B), pointer :: maxlakit => null() !< maximum number of iterations in LAK solve
76  real(dp), pointer :: surfdep => null()
77  real(dp), pointer :: dmaxchg => null()
78  real(dp), pointer :: delh => null()
79  real(dp), pointer :: pdmax => null()
80  integer(I4B), pointer :: check_attr => null()
81  ! -- for budgets
82  integer(I4B), pointer :: bditems => null()
83  ! -- vectors
84  ! -- lake data
85  integer(I4B), dimension(:), pointer, contiguous :: nlakeconn => null()
86  integer(I4B), dimension(:), pointer, contiguous :: idxlakeconn => null()
87  integer(I4B), dimension(:), pointer, contiguous :: ntabrow => null()
88  real(dp), dimension(:), pointer, contiguous :: strt => null()
89  real(dp), dimension(:), pointer, contiguous :: laketop => null()
90  real(dp), dimension(:), pointer, contiguous :: lakebot => null()
91  real(dp), dimension(:), pointer, contiguous :: sareamax => null()
92  character(len=LENBOUNDNAME), dimension(:), pointer, &
93  contiguous :: lakename => null()
94  character(len=8), dimension(:), pointer, contiguous :: status => null()
95  real(dp), dimension(:), pointer, contiguous :: avail => null()
96  real(dp), dimension(:), pointer, contiguous :: lkgwsink => null()
97  real(dp), dimension(:), pointer, contiguous :: stage => null()
98  real(dp), dimension(:), pointer, contiguous :: rainfall => null()
99  real(dp), dimension(:), pointer, contiguous :: evaporation => null()
100  real(dp), dimension(:), pointer, contiguous :: runoff => null()
101  real(dp), dimension(:), pointer, contiguous :: inflow => null()
102  real(dp), dimension(:), pointer, contiguous :: withdrawal => null()
103  real(dp), dimension(:, :), pointer, contiguous :: lauxvar => null()
104  !
105  ! -- table data
106  integer(I4B), dimension(:), pointer, contiguous :: ialaktab => null()
107  real(dp), dimension(:), pointer, contiguous :: tabstage => null()
108  real(dp), dimension(:), pointer, contiguous :: tabvolume => null()
109  real(dp), dimension(:), pointer, contiguous :: tabsarea => null()
110  real(dp), dimension(:), pointer, contiguous :: tabwarea => null()
111  !
112  ! -- lake solution data
113  integer(I4B), dimension(:), pointer, contiguous :: ncncvr => null()
114  real(dp), dimension(:), pointer, contiguous :: surfin => null()
115  real(dp), dimension(:), pointer, contiguous :: surfout => null()
116  real(dp), dimension(:), pointer, contiguous :: surfout1 => null()
117  real(dp), dimension(:), pointer, contiguous :: precip => null()
118  real(dp), dimension(:), pointer, contiguous :: precip1 => null()
119  real(dp), dimension(:), pointer, contiguous :: evap => null()
120  real(dp), dimension(:), pointer, contiguous :: evap1 => null()
121  real(dp), dimension(:), pointer, contiguous :: evapo => null()
122  real(dp), dimension(:), pointer, contiguous :: withr => null()
123  real(dp), dimension(:), pointer, contiguous :: withr1 => null()
124  real(dp), dimension(:), pointer, contiguous :: flwin => null()
125  real(dp), dimension(:), pointer, contiguous :: flwiter => null()
126  real(dp), dimension(:), pointer, contiguous :: flwiter1 => null()
127  real(dp), dimension(:), pointer, contiguous :: seep => null()
128  real(dp), dimension(:), pointer, contiguous :: seep1 => null()
129  real(dp), dimension(:), pointer, contiguous :: seep0 => null()
130  real(dp), dimension(:), pointer, contiguous :: stageiter => null()
131  real(dp), dimension(:), pointer, contiguous :: chterm => null()
132  !
133  ! -- lake convergence
134  integer(I4B), dimension(:), pointer, contiguous :: iseepc => null()
135  integer(I4B), dimension(:), pointer, contiguous :: idhc => null()
136  real(dp), dimension(:), pointer, contiguous :: en1 => null()
137  real(dp), dimension(:), pointer, contiguous :: en2 => null()
138  real(dp), dimension(:), pointer, contiguous :: r1 => null()
139  real(dp), dimension(:), pointer, contiguous :: r2 => null()
140  real(dp), dimension(:), pointer, contiguous :: dh0 => null()
141  real(dp), dimension(:), pointer, contiguous :: s0 => null()
142  real(dp), dimension(:), pointer, contiguous :: qgwf0 => null()
143  !
144  ! -- lake connection data
145  integer(I4B), dimension(:), pointer, contiguous :: imap => null()
146  integer(I4B), dimension(:), pointer, contiguous :: cellid => null()
147  integer(I4B), dimension(:), pointer, contiguous :: nodesontop => null()
148  integer(I4B), dimension(:), pointer, contiguous :: ictype => null()
149  real(dp), dimension(:), pointer, contiguous :: bedleak => null()
150  real(dp), dimension(:), pointer, contiguous :: belev => null()
151  real(dp), dimension(:), pointer, contiguous :: telev => null()
152  real(dp), dimension(:), pointer, contiguous :: connlength => null()
153  real(dp), dimension(:), pointer, contiguous :: connwidth => null()
154  real(dp), dimension(:), pointer, contiguous :: sarea => null()
155  real(dp), dimension(:), pointer, contiguous :: warea => null()
156  real(dp), dimension(:), pointer, contiguous :: satcond => null()
157  real(dp), dimension(:), pointer, contiguous :: simcond => null()
158  real(dp), dimension(:), pointer, contiguous :: simlakgw => null()
159  !
160  ! -- lake outlet data
161  integer(I4B), dimension(:), pointer, contiguous :: lakein => null()
162  integer(I4B), dimension(:), pointer, contiguous :: lakeout => null()
163  integer(I4B), dimension(:), pointer, contiguous :: iouttype => null()
164  real(dp), dimension(:), pointer, contiguous :: outrate => null()
165  real(dp), dimension(:), pointer, contiguous :: outinvert => null()
166  real(dp), dimension(:), pointer, contiguous :: outwidth => null()
167  real(dp), dimension(:), pointer, contiguous :: outrough => null()
168  real(dp), dimension(:), pointer, contiguous :: outslope => null()
169  real(dp), dimension(:), pointer, contiguous :: simoutrate => null()
170  !
171  ! -- lake output data
172  real(dp), dimension(:), pointer, contiguous :: qauxcbc => null()
173  real(dp), dimension(:), pointer, contiguous :: dbuff => null()
174  real(dp), dimension(:), pointer, contiguous :: qleak => null()
175  real(dp), dimension(:), pointer, contiguous :: qsto => null()
176  !
177  ! -- pointer to gwf iss and gwf hk
178  integer(I4B), pointer :: gwfiss => null()
179  real(dp), dimension(:), pointer, contiguous :: gwfk11 => null()
180  real(dp), dimension(:), pointer, contiguous :: gwfk33 => null()
181  real(dp), dimension(:), pointer, contiguous :: gwfsat => null()
182  integer(I4B), pointer :: gwfik33 => null()
183  !
184  ! -- package x, xold, and ibound
185  integer(I4B), dimension(:), pointer, contiguous :: iboundpak => null() !package ibound
186  real(dp), dimension(:), pointer, contiguous :: xnewpak => null() !package x vector
187  real(dp), dimension(:), pointer, contiguous :: xoldpak => null() !package xold vector
188  !
189  ! -- lake budget object
190  type(budgetobjecttype), pointer :: budobj => null()
191  !
192  ! -- lake table objects
193  type(tabletype), pointer :: stagetab => null()
194  type(tabletype), pointer :: pakcsvtab => null()
195  !
196  ! -- density variables
197  integer(I4B), pointer :: idense
198  real(dp), dimension(:, :), pointer, contiguous :: denseterms => null()
199  !
200  ! -- viscosity variables
201  real(dp), dimension(:, :), pointer, contiguous :: viscratios => null() !< viscosity ratios (1: lak vsc ratio; 2: gwf vsc ratio)
202  !
203  ! -- type bound procedures
204 
205  contains
206 
207  procedure :: lak_allocate_scalars
208  procedure :: lak_allocate_arrays
209  procedure :: bnd_options => lak_options
210  procedure :: read_dimensions => lak_read_dimensions
211  procedure :: read_initial_attr => lak_read_initial_attr
212  procedure :: set_pointers => lak_set_pointers
213  procedure :: bnd_ar => lak_ar
214  procedure :: bnd_rp => lak_rp
215  procedure :: bnd_ad => lak_ad
216  procedure :: bnd_cf => lak_cf
217  procedure :: bnd_fc => lak_fc
218  procedure :: bnd_fn => lak_fn
219  procedure :: bnd_cc => lak_cc
220  procedure :: bnd_cq => lak_cq
221  procedure :: bnd_ot_model_flows => lak_ot_model_flows
222  procedure :: bnd_ot_package_flows => lak_ot_package_flows
223  procedure :: bnd_ot_dv => lak_ot_dv
224  procedure :: bnd_ot_bdsummary => lak_ot_bdsummary
225  procedure :: bnd_da => lak_da
226  procedure :: define_listlabel
227  ! -- methods for observations
228  procedure, public :: bnd_obs_supported => lak_obs_supported
229  procedure, public :: bnd_df_obs => lak_df_obs
230  procedure, public :: bnd_rp_obs => lak_rp_obs
231  procedure, public :: bnd_bd_obs => lak_bd_obs
232  ! -- private procedures
233  procedure, private :: lak_read_lakes
234  procedure, private :: lak_read_lake_connections
235  procedure, private :: lak_read_outlets
236  procedure, private :: lak_read_tables
237  procedure, private :: lak_read_table
238  procedure, private :: lak_check_valid
239  procedure, private :: lak_set_stressperiod
240  procedure, private :: lak_set_attribute_error
241  procedure, private :: lak_bound_update
242  procedure, private :: lak_calculate_sarea
243  procedure, private :: lak_calculate_warea
244  procedure, private :: lak_calculate_conn_warea
245  procedure, public :: lak_calculate_vol
246  procedure, private :: lak_calculate_conductance
247  procedure, private :: lak_calculate_cond_head
248  procedure, private :: lak_calculate_conn_conductance
249  procedure, private :: lak_calculate_exchange
250  procedure, private :: lak_calculate_conn_exchange
251  procedure, private :: lak_estimate_conn_exchange
252  procedure, private :: lak_calculate_storagechange
253  procedure, private :: lak_calculate_rainfall
254  procedure, private :: lak_calculate_runoff
255  procedure, private :: lak_calculate_inflow
256  procedure, private :: lak_calculate_external
257  procedure, private :: lak_calculate_withdrawal
258  procedure, private :: lak_calculate_evaporation
259  procedure, private :: lak_calculate_outlet_inflow
260  procedure, private :: lak_calculate_outlet_outflow
261  procedure, private :: lak_get_internal_inlet
262  procedure, private :: lak_get_internal_outlet
263  procedure, private :: lak_get_external_outlet
264  procedure, private :: lak_get_internal_mover
265  procedure, private :: lak_get_external_mover
266  procedure, private :: lak_get_outlet_tomover
267  procedure, private :: lak_accumulate_chterm
268  procedure, private :: lak_vol2stage
269  procedure, private :: lak_solve
270  procedure, private :: lak_bisection
271  procedure, private :: lak_calculate_available
272  procedure, private :: lak_calculate_residual
273  procedure, private :: lak_linear_interpolation
274  procedure, private :: lak_setup_budobj
275  procedure, private :: lak_fill_budobj
276  procedure, private :: laktables_to_vectors
277  ! -- table
278  procedure, private :: lak_setup_tableobj
279  ! -- density
280  procedure :: lak_activate_density
281  procedure, private :: lak_calculate_density_exchange
282  ! -- viscosity
283  procedure :: lak_activate_viscosity
284  end type laktype
285 
286 contains
287 
288  !> @brief Create a new LAK Package and point bndobj to the new package
289  !<
290  subroutine lak_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
291  ! -- dummy
292  class(bndtype), pointer :: packobj
293  integer(I4B), intent(in) :: id
294  integer(I4B), intent(in) :: ibcnum
295  integer(I4B), intent(in) :: inunit
296  integer(I4B), intent(in) :: iout
297  character(len=*), intent(in) :: namemodel
298  character(len=*), intent(in) :: pakname
299  ! -- local
300  type(laktype), pointer :: lakobj
301  !
302  ! -- allocate the object and assign values to object variables
303  allocate (lakobj)
304  packobj => lakobj
305  !
306  ! -- create name and memory path
307  call packobj%set_names(ibcnum, namemodel, pakname, ftype)
308  packobj%text = text
309  !
310  ! -- allocate scalars
311  call lakobj%lak_allocate_scalars()
312  !
313  ! -- initialize package
314  call packobj%pack_initialize()
315  !
316  packobj%inunit = inunit
317  packobj%iout = iout
318  packobj%id = id
319  packobj%ibcnum = ibcnum
320  packobj%ncolbnd = 3
321  packobj%iscloc = 0 ! not supported
322  packobj%isadvpak = 1
323  packobj%ictMemPath = create_mem_path(namemodel, 'NPF')
324  end subroutine lak_create
325 
326  !> @brief Allocate scalar members
327  !<
328  subroutine lak_allocate_scalars(this)
329  ! -- dummy
330  class(laktype), intent(inout) :: this
331  !
332  ! -- call standard BndType allocate scalars
333  call this%BndType%allocate_scalars()
334  !
335  ! -- allocate the object and assign values to object variables
336  call mem_allocate(this%iprhed, 'IPRHED', this%memoryPath)
337  call mem_allocate(this%istageout, 'ISTAGEOUT', this%memoryPath)
338  call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath)
339  call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
340  call mem_allocate(this%ipakcsv, 'IPAKCSV', this%memoryPath)
341  call mem_allocate(this%nlakes, 'NLAKES', this%memoryPath)
342  call mem_allocate(this%noutlets, 'NOUTLETS', this%memoryPath)
343  call mem_allocate(this%ntables, 'NTABLES', this%memoryPath)
344  call mem_allocate(this%convlength, 'CONVLENGTH', this%memoryPath)
345  call mem_allocate(this%convtime, 'CONVTIME', this%memoryPath)
346  call mem_allocate(this%outdmax, 'OUTDMAX', this%memoryPath)
347  call mem_allocate(this%igwhcopt, 'IGWHCOPT', this%memoryPath)
348  call mem_allocate(this%iconvchk, 'ICONVCHK', this%memoryPath)
349  call mem_allocate(this%iconvresidchk, 'ICONVRESIDCHK', this%memoryPath)
350  call mem_allocate(this%maxlakit, 'MAXLAKIT', this%memoryPath)
351  call mem_allocate(this%surfdep, 'SURFDEP', this%memoryPath)
352  call mem_allocate(this%dmaxchg, 'DMAXCHG', this%memoryPath)
353  call mem_allocate(this%delh, 'DELH', this%memoryPath)
354  call mem_allocate(this%pdmax, 'PDMAX', this%memoryPath)
355  call mem_allocate(this%check_attr, 'CHECK_ATTR', this%memoryPath)
356  call mem_allocate(this%bditems, 'BDITEMS', this%memoryPath)
357  call mem_allocate(this%cbcauxitems, 'CBCAUXITEMS', this%memoryPath)
358  call mem_allocate(this%idense, 'IDENSE', this%memoryPath)
359  !
360  ! -- Set values
361  this%iprhed = 0
362  this%istageout = 0
363  this%ibudgetout = 0
364  this%ibudcsv = 0
365  this%ipakcsv = 0
366  this%nlakes = 0
367  this%noutlets = 0
368  this%ntables = 0
369  this%convlength = done
370  this%convtime = done
371  this%outdmax = dzero
372  this%igwhcopt = 0
373  this%iconvchk = 1
374  this%iconvresidchk = 1
375  this%maxlakit = maxadpit
376  this%surfdep = dzero
377  this%dmaxchg = dem5
378  this%delh = dp999 * this%dmaxchg
379  this%pdmax = dem1
380  this%bditems = 11
381  this%cbcauxitems = 1
382  this%idense = 0
383  this%ivsc = 0
384  end subroutine lak_allocate_scalars
385 
386  !> @brief Allocate scalar members
387  !<
388  subroutine lak_allocate_arrays(this)
389  ! -- modules
390  ! -- dummy
391  class(laktype), intent(inout) :: this
392  ! -- local
393  integer(I4B) :: i
394  !
395  ! -- call standard BndType allocate scalars
396  call this%BndType%allocate_arrays()
397  !
398  ! -- allocate character array for budget text
399  allocate (this%clakbudget(this%bditems))
400  !
401  !-- fill clakbudget
402  this%clakbudget(1) = ' GWF'
403  this%clakbudget(2) = ' RAINFALL'
404  this%clakbudget(3) = ' EVAPORATION'
405  this%clakbudget(4) = ' RUNOFF'
406  this%clakbudget(5) = ' EXT-INFLOW'
407  this%clakbudget(6) = ' WITHDRAWAL'
408  this%clakbudget(7) = ' EXT-OUTFLOW'
409  this%clakbudget(8) = ' STORAGE'
410  this%clakbudget(9) = ' CONSTANT'
411  this%clakbudget(10) = ' FROM-MVR'
412  this%clakbudget(11) = ' TO-MVR'
413  !
414  ! -- allocate and initialize dbuff
415  if (this%istageout > 0) then
416  call mem_allocate(this%dbuff, this%nlakes, 'DBUFF', this%memoryPath)
417  do i = 1, this%nlakes
418  this%dbuff(i) = dzero
419  end do
420  else
421  call mem_allocate(this%dbuff, 0, 'DBUFF', this%memoryPath)
422  end if
423  !
424  ! -- allocate character array for budget text
425  allocate (this%cauxcbc(this%cbcauxitems))
426  !
427  ! -- allocate and initialize qauxcbc
428  call mem_allocate(this%qauxcbc, this%cbcauxitems, 'QAUXCBC', this%memoryPath)
429  do i = 1, this%cbcauxitems
430  this%qauxcbc(i) = dzero
431  end do
432  !
433  ! -- allocate qleak and qsto
434  call mem_allocate(this%qleak, this%maxbound, 'QLEAK', this%memoryPath)
435  do i = 1, this%maxbound
436  this%qleak(i) = dzero
437  end do
438  call mem_allocate(this%qsto, this%nlakes, 'QSTO', this%memoryPath)
439  do i = 1, this%nlakes
440  this%qsto(i) = dzero
441  end do
442  !
443  ! -- allocate denseterms to size 0
444  call mem_allocate(this%denseterms, 3, 0, 'DENSETERMS', this%memoryPath)
445  !
446  ! -- allocate viscratios to size 0
447  call mem_allocate(this%viscratios, 2, 0, 'VISCRATIOS', this%memoryPath)
448  end subroutine lak_allocate_arrays
449 
450  !> @brief Read the dimensions for this package
451  !<
452  subroutine lak_read_lakes(this)
453  ! -- modules
454  use constantsmodule, only: linelength
455  use simmodule, only: store_error, count_errors, store_error_unit
456  use timeseriesmanagermodule, only: read_value_or_time_series_adv
457  ! -- dummy
458  class(laktype), intent(inout) :: this
459  ! -- local
460  character(len=LINELENGTH) :: text
461  character(len=LENBOUNDNAME) :: bndName, bndNameTemp
462  character(len=9) :: cno
463  character(len=50), dimension(:), allocatable :: caux
464  integer(I4B) :: ierr, ival
465  logical(LGP) :: isfound, endOfBlock
466  integer(I4B) :: n
467  integer(I4B) :: ii, jj
468  integer(I4B) :: iaux
469  integer(I4B) :: itmp
470  integer(I4B) :: nlak
471  integer(I4B) :: nconn
472  integer(I4B), dimension(:), pointer, contiguous :: nboundchk
473  real(DP), pointer :: bndElem => null()
474  !
475  ! -- initialize itmp
476  itmp = 0
477  !
478  ! -- allocate lake data
479  call mem_allocate(this%nlakeconn, this%nlakes, 'NLAKECONN', this%memoryPath)
480  call mem_allocate(this%idxlakeconn, this%nlakes + 1, 'IDXLAKECONN', &
481  this%memoryPath)
482  call mem_allocate(this%ntabrow, this%nlakes, 'NTABROW', this%memoryPath)
483  call mem_allocate(this%strt, this%nlakes, 'STRT', this%memoryPath)
484  call mem_allocate(this%laketop, this%nlakes, 'LAKETOP', this%memoryPath)
485  call mem_allocate(this%lakebot, this%nlakes, 'LAKEBOT', this%memoryPath)
486  call mem_allocate(this%sareamax, this%nlakes, 'SAREAMAX', this%memoryPath)
487  call mem_allocate(this%stage, this%nlakes, 'STAGE', this%memoryPath)
488  call mem_allocate(this%rainfall, this%nlakes, 'RAINFALL', this%memoryPath)
489  call mem_allocate(this%evaporation, this%nlakes, 'EVAPORATION', &
490  this%memoryPath)
491  call mem_allocate(this%runoff, this%nlakes, 'RUNOFF', this%memoryPath)
492  call mem_allocate(this%inflow, this%nlakes, 'INFLOW', this%memoryPath)
493  call mem_allocate(this%withdrawal, this%nlakes, 'WITHDRAWAL', this%memoryPath)
494  call mem_allocate(this%lauxvar, this%naux, this%nlakes, 'LAUXVAR', &
495  this%memoryPath)
496  call mem_allocate(this%avail, this%nlakes, 'AVAIL', this%memoryPath)
497  call mem_allocate(this%lkgwsink, this%nlakes, 'LKGWSINK', this%memoryPath)
498  call mem_allocate(this%ncncvr, this%nlakes, 'NCNCVR', this%memoryPath)
499  call mem_allocate(this%surfin, this%nlakes, 'SURFIN', this%memoryPath)
500  call mem_allocate(this%surfout, this%nlakes, 'SURFOUT', this%memoryPath)
501  call mem_allocate(this%surfout1, this%nlakes, 'SURFOUT1', this%memoryPath)
502  call mem_allocate(this%precip, this%nlakes, 'PRECIP', this%memoryPath)
503  call mem_allocate(this%precip1, this%nlakes, 'PRECIP1', this%memoryPath)
504  call mem_allocate(this%evap, this%nlakes, 'EVAP', this%memoryPath)
505  call mem_allocate(this%evap1, this%nlakes, 'EVAP1', this%memoryPath)
506  call mem_allocate(this%evapo, this%nlakes, 'EVAPO', this%memoryPath)
507  call mem_allocate(this%withr, this%nlakes, 'WITHR', this%memoryPath)
508  call mem_allocate(this%withr1, this%nlakes, 'WITHR1', this%memoryPath)
509  call mem_allocate(this%flwin, this%nlakes, 'FLWIN', this%memoryPath)
510  call mem_allocate(this%flwiter, this%nlakes, 'FLWITER', this%memoryPath)
511  call mem_allocate(this%flwiter1, this%nlakes, 'FLWITER1', this%memoryPath)
512  call mem_allocate(this%seep, this%nlakes, 'SEEP', this%memoryPath)
513  call mem_allocate(this%seep1, this%nlakes, 'SEEP1', this%memoryPath)
514  call mem_allocate(this%seep0, this%nlakes, 'SEEP0', this%memoryPath)
515  call mem_allocate(this%stageiter, this%nlakes, 'STAGEITER', this%memoryPath)
516  call mem_allocate(this%chterm, this%nlakes, 'CHTERM', this%memoryPath)
517  !
518  ! -- lake boundary and stages
519  call mem_allocate(this%iboundpak, this%nlakes, 'IBOUND', this%memoryPath)
520  call mem_allocate(this%xnewpak, this%nlakes, 'XNEWPAK', this%memoryPath)
521  call mem_allocate(this%xoldpak, this%nlakes, 'XOLDPAK', this%memoryPath)
522  !
523  ! -- lake iteration variables
524  call mem_allocate(this%iseepc, this%nlakes, 'ISEEPC', this%memoryPath)
525  call mem_allocate(this%idhc, this%nlakes, 'IDHC', this%memoryPath)
526  call mem_allocate(this%en1, this%nlakes, 'EN1', this%memoryPath)
527  call mem_allocate(this%en2, this%nlakes, 'EN2', this%memoryPath)
528  call mem_allocate(this%r1, this%nlakes, 'R1', this%memoryPath)
529  call mem_allocate(this%r2, this%nlakes, 'R2', this%memoryPath)
530  call mem_allocate(this%dh0, this%nlakes, 'DH0', this%memoryPath)
531  call mem_allocate(this%s0, this%nlakes, 'S0', this%memoryPath)
532  call mem_allocate(this%qgwf0, this%nlakes, 'QGWF0', this%memoryPath)
533  !
534  ! -- allocate character storage not managed by the memory manager
535  allocate (this%lakename(this%nlakes)) ! ditch after boundnames allocated??
536  allocate (this%status(this%nlakes))
537  !
538  do n = 1, this%nlakes
539  this%ntabrow(n) = 0
540  this%status(n) = 'ACTIVE'
541  this%laketop(n) = -dep20
542  this%lakebot(n) = dep20
543  this%sareamax(n) = dzero
544  this%iboundpak(n) = 1
545  this%xnewpak(n) = dep20
546  this%xoldpak(n) = dep20
547  !
548  ! -- initialize boundary values to zero
549  this%rainfall(n) = dzero
550  this%evaporation(n) = dzero
551  this%runoff(n) = dzero
552  this%inflow(n) = dzero
553  this%withdrawal(n) = dzero
554  end do
555  !
556  ! -- allocate local storage for aux variables
557  if (this%naux > 0) then
558  allocate (caux(this%naux))
559  end if
560  !
561  ! -- allocate and initialize temporary variables
562  allocate (nboundchk(this%nlakes))
563  do n = 1, this%nlakes
564  nboundchk(n) = 0
565  end do
566  !
567  ! -- read lake well data
568  ! -- get lakes block
569  call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, &
570  supportopenclose=.true.)
571  !
572  ! -- parse locations block if detected
573  if (isfound) then
574  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
575  ' PACKAGEDATA'
576  nlak = 0
577  nconn = 0
578  do
579  call this%parser%GetNextLine(endofblock)
580  if (endofblock) exit
581  n = this%parser%GetInteger()
582  !
583  if (n < 1 .or. n > this%nlakes) then
584  write (errmsg, '(a,1x,i0)') 'lakeno MUST BE > 0 and <= ', this%nlakes
585  call store_error(errmsg)
586  cycle
587  end if
588  !
589  ! -- increment nboundchk
590  nboundchk(n) = nboundchk(n) + 1
591  !
592  ! -- strt
593  this%strt(n) = this%parser%GetDouble()
594  !
595  ! nlakeconn
596  ival = this%parser%GetInteger()
597  !
598  if (ival < 0) then
599  write (errmsg, '(a,1x,i0)') 'nlakeconn MUST BE >= 0 for lake ', n
600  call store_error(errmsg)
601  end if
602  !
603  nconn = nconn + ival
604  this%nlakeconn(n) = ival
605  !
606  ! -- get aux data
607  do iaux = 1, this%naux
608  call this%parser%GetString(caux(iaux))
609  end do
610  !
611  ! -- set default bndName
612  write (cno, '(i9.9)') n
613  bndname = 'Lake'//cno
614  !
615  ! -- lakename
616  if (this%inamedbound /= 0) then
617  call this%parser%GetStringCaps(bndnametemp)
618  if (bndnametemp /= '') then
619  bndname = bndnametemp
620  end if
621  end if
622  this%lakename(n) = bndname
623  !
624  ! -- fill time series aware data
625  ! -- fill aux data
626  do jj = 1, this%naux
627  text = caux(jj)
628  ii = n
629  bndelem => this%lauxvar(jj, ii)
630  call read_value_or_time_series_adv(text, ii, jj, bndelem, &
631  this%packName, 'AUX', &
632  this%tsManager, this%iprpak, &
633  this%auxname(jj))
634  end do
635  !
636  nlak = nlak + 1
637  end do
638  !
639  ! -- check for duplicate or missing lakes
640  do n = 1, this%nlakes
641  if (nboundchk(n) == 0) then
642  write (errmsg, '(a,1x,i0)') 'NO DATA SPECIFIED FOR LAKE', n
643  call store_error(errmsg)
644  else if (nboundchk(n) > 1) then
645  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
646  'DATA FOR LAKE', n, 'SPECIFIED', nboundchk(n), 'TIMES'
647  call store_error(errmsg)
648  end if
649  end do
650  !
651  write (this%iout, '(1x,a)') 'END OF '//trim(adjustl(this%text))// &
652  ' PACKAGEDATA'
653  else
654  call store_error('REQUIRED PACKAGEDATA BLOCK NOT FOUND.')
655  end if
656  !
657  ! -- terminate if any errors were detected
658  if (count_errors() > 0) then
659  call this%parser%StoreErrorUnit()
660  end if
661  !
662  ! -- set MAXBOUND
663  this%MAXBOUND = nconn
664  write (this%iout, '(//4x,a,i7)') 'MAXBOUND = ', this%maxbound
665  !
666  ! -- set idxlakeconn
667  this%idxlakeconn(1) = 1
668  do n = 1, this%nlakes
669  this%idxlakeconn(n + 1) = this%idxlakeconn(n) + this%nlakeconn(n)
670  end do
671  !
672  ! -- deallocate local storage for aux variables
673  if (this%naux > 0) then
674  deallocate (caux)
675  end if
676  !
677  ! -- deallocate local storage for nboundchk
678  deallocate (nboundchk)
679  end subroutine lak_read_lakes
680 
681  !> @brief Read the lake connections for this package
682  !<
683  subroutine lak_read_lake_connections(this)
684  use constantsmodule, only: linelength, lenvarname
685  use simmodule, only: store_error, count_errors
686  ! -- dummy
687  class(laktype), intent(inout) :: this
688  ! -- local
689  character(len=LINELENGTH) :: keyword, cellid
690  integer(I4B) :: ierr, ival
691  logical(LGP) :: isfound, endOfBlock
692  logical(LGP) :: is_lake_bed
693  real(DP) :: rval
694  integer(I4B) :: j, n
695  integer(I4B) :: nn
696  integer(I4B) :: ipos, ipos0
697  integer(I4B) :: icellid, icellid0
698  real(DP) :: top
699  real(DP) :: bot
700  integer(I4B), dimension(:), pointer, contiguous :: nboundchk
701  character(len=LENVARNAME) :: ctypenm
702  !
703  ! -- allocate local storage
704  allocate (nboundchk(this%MAXBOUND))
705  do n = 1, this%MAXBOUND
706  nboundchk(n) = 0
707  end do
708  !
709  ! -- get connectiondata block
710  call this%parser%GetBlock('CONNECTIONDATA', isfound, ierr, &
711  supportopenclose=.true.)
712  !
713  ! -- parse connectiondata block if detected
714  if (isfound) then
715  ! -- allocate connection data using memory manager
716  call mem_allocate(this%imap, this%MAXBOUND, 'IMAP', this%memoryPath)
717  call mem_allocate(this%cellid, this%MAXBOUND, 'CELLID', this%memoryPath)
718  call mem_allocate(this%nodesontop, this%MAXBOUND, 'NODESONTOP', &
719  this%memoryPath)
720  call mem_allocate(this%ictype, this%MAXBOUND, 'ICTYPE', this%memoryPath)
721  call mem_allocate(this%bedleak, this%MAXBOUND, 'BEDLEAK', this%memoryPath) ! don't need to save this - use a temporary vector
722  call mem_allocate(this%belev, this%MAXBOUND, 'BELEV', this%memoryPath)
723  call mem_allocate(this%telev, this%MAXBOUND, 'TELEV', this%memoryPath)
724  call mem_allocate(this%connlength, this%MAXBOUND, 'CONNLENGTH', &
725  this%memoryPath)
726  call mem_allocate(this%connwidth, this%MAXBOUND, 'CONNWIDTH', &
727  this%memoryPath)
728  call mem_allocate(this%sarea, this%MAXBOUND, 'SAREA', this%memoryPath)
729  call mem_allocate(this%warea, this%MAXBOUND, 'WAREA', this%memoryPath)
730  call mem_allocate(this%satcond, this%MAXBOUND, 'SATCOND', this%memoryPath)
731  call mem_allocate(this%simcond, this%MAXBOUND, 'SIMCOND', this%memoryPath)
732  call mem_allocate(this%simlakgw, this%MAXBOUND, 'SIMLAKGW', this%memoryPath)
733  !
734  ! -- process the lake connection data
735  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
736  ' LAKE_CONNECTIONS'
737  do
738  call this%parser%GetNextLine(endofblock)
739  if (endofblock) exit
740  n = this%parser%GetInteger()
741  !
742  if (n < 1 .or. n > this%nlakes) then
743  write (errmsg, '(a,1x,i0)') 'lakeno MUST BE > 0 and <= ', this%nlakes
744  call store_error(errmsg)
745  cycle
746  end if
747  !
748  ! -- read connection number
749  ival = this%parser%GetInteger()
750  if (ival < 1 .or. ival > this%nlakeconn(n)) then
751  write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
752  'iconn FOR LAKE ', n, 'MUST BE > 1 and <= ', this%nlakeconn(n)
753  call store_error(errmsg)
754  cycle
755  end if
756  !
757  j = ival
758  ipos = this%idxlakeconn(n) + ival - 1
759  !
760  ! -- set imap
761  this%imap(ipos) = n
762  !
763  !
764  ! -- increment nboundchk
765  nboundchk(ipos) = nboundchk(ipos) + 1
766  !
767  ! -- read gwfnodes from the line
768  call this%parser%GetCellid(this%dis%ndim, cellid)
769  nn = this%dis%noder_from_cellid(cellid, &
770  this%parser%iuactive, this%iout)
771  !
772  ! -- determine if a valid cell location was provided
773  if (nn < 1) then
774  write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
775  'INVALID cellid FOR LAKE ', n, 'connection', j
776  call store_error(errmsg)
777  end if
778  !
779  ! -- set gwf cellid for connection
780  this%cellid(ipos) = nn
781  this%nodesontop(ipos) = nn
782  !
783  ! -- read ictype
784  call this%parser%GetStringCaps(keyword)
785  select case (keyword)
786  case ('VERTICAL')
787  this%ictype(ipos) = 0
788  case ('HORIZONTAL')
789  this%ictype(ipos) = 1
790  case ('EMBEDDEDH')
791  this%ictype(ipos) = 2
792  case ('EMBEDDEDV')
793  this%ictype(ipos) = 3
794  case default
795  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,a,a)') &
796  'UNKNOWN ctype FOR LAKE ', n, 'connection', j, &
797  '(', trim(keyword), ')'
798  call store_error(errmsg)
799  end select
800  write (ctypenm, '(a16)') keyword
801  !
802  ! -- bed leakance
803  !this%bedleak(ipos) = this%parser%GetDouble() !TODO: use this when NONE keyword deprecated
804  call this%parser%GetStringCaps(keyword)
805  select case (keyword)
806  case ('NONE')
807  is_lake_bed = .false.
808  this%bedleak(ipos) = dnodata
809  !
810  ! -- create warning message
811  write (warnmsg, '(2(a,1x,i0,1x),a,1pe8.1,a)') &
812  'BEDLEAK for connection', j, 'in lake', n, 'is specified to '// &
813  'be NONE. Lake connections where the lake-GWF connection '// &
814  'conductance is solely a function of aquifer properties '// &
815  'in the connected GWF cell should be specified with a '// &
816  'DNODATA (', dnodata, ') value.'
817  !
818  ! -- create deprecation warning
819  call deprecation_warning('CONNECTIONDATA', 'bedleak=NONE', '6.4.3', &
820  warnmsg, this%parser%GetUnit())
821  case default
822  read (keyword, *) rval
823  if (is_close(rval, dnodata)) then
824  is_lake_bed = .false.
825  else
826  is_lake_bed = .true.
827  end if
828  this%bedleak(ipos) = rval
829  end select
830  !
831  if (is_lake_bed .and. this%bedleak(ipos) < dzero) then
832  write (errmsg, '(a,1x,i0,1x,a)') 'bedleak FOR LAKE ', n, 'MUST BE >= 0'
833  call store_error(errmsg)
834  end if
835  !
836  ! -- belev
837  this%belev(ipos) = this%parser%GetDouble()
838  !
839  ! -- telev
840  this%telev(ipos) = this%parser%GetDouble()
841  !
842  ! -- connection length
843  rval = this%parser%GetDouble()
844  if (rval <= dzero) then
845  if (this%ictype(ipos) == 1 .or. this%ictype(ipos) == 2 .or. &
846  this%ictype(ipos) == 3) then
847  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,a,1x,a)') &
848  'connection length (connlen) FOR LAKE ', n, &
849  ', CONNECTION NO.', j, ', MUST BE > 0 FOR SPECIFIED ', &
850  'connection type (ctype)', ctypenm
851  call store_error(errmsg)
852  else
853  rval = dzero
854  end if
855  end if
856  this%connlength(ipos) = rval
857  !
858  ! -- connection width
859  rval = this%parser%GetDouble()
860  if (rval < dzero) then
861  if (this%ictype(ipos) == 1) then
862  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
863  'cell width (connwidth) FOR LAKE ', n, &
864  ' HORIZONTAL CONNECTION ', j, 'MUST BE >= 0'
865  call store_error(errmsg)
866  else
867  rval = dzero
868  end if
869  end if
870  this%connwidth(ipos) = rval
871  end do
872  write (this%iout, '(1x,a)') &
873  'END OF '//trim(adjustl(this%text))//' CONNECTIONDATA'
874  else
875  call store_error('REQUIRED CONNECTIONDATA BLOCK NOT FOUND.')
876  end if
877  !
878  ! -- terminate if any errors were detected
879  if (count_errors() > 0) then
880  call this%parser%StoreErrorUnit()
881  end if
882  !
883  ! -- check that embedded lakes have only one connection
884  do n = 1, this%nlakes
885  j = 0
886  do ipos = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
887  if (this%ictype(ipos) /= 2 .and. this%ictype(ipos) /= 3) cycle
888  j = j + 1
889  if (j > 1) then
890  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
891  'nlakeconn FOR LAKE', n, 'EMBEDDED CONNECTION', j, ' EXCEEDS 1.'
892  call store_error(errmsg)
893  end if
894  end do
895  end do
896  ! -- check that an embedded lake is not in the same cell as a lake
897  ! with a vertical connection
898  do n = 1, this%nlakes
899  ipos0 = this%idxlakeconn(n)
900  icellid0 = this%cellid(ipos0)
901  if (this%ictype(ipos0) /= 2 .and. this%ictype(ipos0) /= 3) cycle
902  do nn = 1, this%nlakes
903  if (nn == n) cycle
904  j = 0
905  do ipos = this%idxlakeconn(nn), this%idxlakeconn(nn + 1) - 1
906  j = j + 1
907  icellid = this%cellid(ipos)
908  if (icellid == icellid0) then
909  if (this%ictype(ipos) == 0) then
910  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
911  'EMBEDDED LAKE', n, &
912  'CANNOT COINCIDE WITH VERTICAL CONNECTION', j, &
913  'IN LAKE', nn, '.'
914  call store_error(errmsg)
915  end if
916  end if
917  end do
918  end do
919  end do
920  !
921  ! -- process the data
922  do n = 1, this%nlakes
923  j = 0
924  do ipos = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
925  j = j + 1
926  nn = this%cellid(ipos)
927  top = this%dis%top(nn)
928  bot = this%dis%bot(nn)
929  ! vertical connection
930  if (this%ictype(ipos) == 0) then
931  this%telev(ipos) = top + this%surfdep
932  this%belev(ipos) = top
933  this%lakebot(n) = min(this%belev(ipos), this%lakebot(n))
934  ! horizontal connection
935  else if (this%ictype(ipos) == 1) then
936  if (this%belev(ipos) == this%telev(ipos)) then
937  this%telev(ipos) = top
938  this%belev(ipos) = bot
939  else
940  if (this%belev(ipos) >= this%telev(ipos)) then
941  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
942  'telev FOR LAKE ', n, ' HORIZONTAL CONNECTION ', j, &
943  'MUST BE >= belev'
944  call store_error(errmsg)
945  else if (this%belev(ipos) < bot) then
946  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.7,1x,a)') &
947  'belev FOR LAKE ', n, ' HORIZONTAL CONNECTION ', j, &
948  'MUST BE >= cell bottom (', bot, ')'
949  call store_error(errmsg)
950  else if (this%telev(ipos) > top) then
951  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.7,1x,a)') &
952  'telev FOR LAKE ', n, ' HORIZONTAL CONNECTION ', j, &
953  'MUST BE <= cell top (', top, ')'
954  call store_error(errmsg)
955  end if
956  end if
957  this%laketop(n) = max(this%telev(ipos), this%laketop(n))
958  this%lakebot(n) = min(this%belev(ipos), this%lakebot(n))
959  ! embedded connections
960  else if (this%ictype(ipos) == 2 .or. this%ictype(ipos) == 3) then
961  this%telev(ipos) = top
962  this%belev(ipos) = bot
963  this%lakebot(n) = bot
964  end if
965  !
966  ! -- check for missing or duplicate lake connections
967  if (nboundchk(ipos) == 0) then
968  write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
969  'NO DATA SPECIFIED FOR LAKE', n, 'CONNECTION', j
970  call store_error(errmsg)
971  else if (nboundchk(ipos) > 1) then
972  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
973  'DATA FOR LAKE', n, 'CONNECTION', j, &
974  'SPECIFIED', nboundchk(ipos), 'TIMES'
975  call store_error(errmsg)
976  end if
977  !
978  ! -- set laketop if it has not been assigned
979  end do
980  if (this%laketop(n) == -dep20) then
981  this%laketop(n) = this%lakebot(n) + 100.
982  end if
983  end do
984  !
985  ! -- deallocate local variable
986  deallocate (nboundchk)
987  !
988  ! -- write summary of lake_connection error messages
989  if (count_errors() > 0) then
990  call this%parser%StoreErrorUnit()
991  end if
992  end subroutine lak_read_lake_connections
993 
994  !> @brief Read the lake tables for this package
995  !<
996  subroutine lak_read_tables(this)
997  use constantsmodule, only: linelength
998  use simmodule, only: store_error, count_errors
999  ! -- dummy
1000  class(laktype), intent(inout) :: this
1001  ! -- local
1002  type(laktabtype), dimension(:), allocatable :: laketables
1003  character(len=LINELENGTH) :: line
1004  character(len=LINELENGTH) :: keyword
1005  integer(I4B) :: ierr
1006  logical(LGP) :: isfound, endOfBlock
1007  integer(I4B) :: n
1008  integer(I4B) :: iconn
1009  integer(I4B) :: ntabs
1010  integer(I4B), dimension(:), pointer, contiguous :: nboundchk
1011  !
1012  ! -- skip of no outlets
1013  if (this%ntables < 1) return
1014  !
1015  ! -- allocate and initialize nboundchk
1016  allocate (nboundchk(this%nlakes))
1017  do n = 1, this%nlakes
1018  nboundchk(n) = 0
1019  end do
1020  !
1021  ! -- allocate derived type for table data
1022  allocate (laketables(this%nlakes))
1023  !
1024  ! -- get lake_tables block
1025  call this%parser%GetBlock('TABLES', isfound, ierr, &
1026  supportopenclose=.true.)
1027  !
1028  ! -- parse lake_tables block if detected
1029  if (isfound) then
1030  ntabs = 0
1031  ! -- process the lake table data
1032  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
1033  ' LAKE_TABLES'
1034  readtable: do
1035  call this%parser%GetNextLine(endofblock)
1036  if (endofblock) exit
1037  n = this%parser%GetInteger()
1038  !
1039  if (n < 1 .or. n > this%nlakes) then
1040  write (errmsg, '(a,1x,i0)') 'lakeno MUST BE > 0 and <= ', this%nlakes
1041  call store_error(errmsg)
1042  cycle readtable
1043  end if
1044  !
1045  ! -- increment ntab and nboundchk
1046  ntabs = ntabs + 1
1047  nboundchk(n) = nboundchk(n) + 1
1048  !
1049  ! -- read FILE keyword
1050  call this%parser%GetStringCaps(keyword)
1051  select case (keyword)
1052  case ('TAB6')
1053  call this%parser%GetStringCaps(keyword)
1054  if (trim(adjustl(keyword)) /= 'FILEIN') then
1055  errmsg = 'TAB6 keyword must be followed by "FILEIN" '// &
1056  'then by filename.'
1057  call store_error(errmsg)
1058  cycle readtable
1059  end if
1060  call this%parser%GetString(line)
1061  call this%lak_read_table(n, line, laketables(n))
1062  case default
1063  write (errmsg, '(a,1x,i0,1x,a)') &
1064  'LAKE TABLE ENTRY for LAKE ', n, 'MUST INCLUDE TAB6 KEYWORD'
1065  call store_error(errmsg)
1066  cycle readtable
1067  end select
1068  end do readtable
1069  !
1070  write (this%iout, '(1x,a)') &
1071  'END OF '//trim(adjustl(this%text))//' LAKE_TABLES'
1072  !
1073  ! -- check for missing or duplicate lake connections
1074  if (ntabs < this%ntables) then
1075  write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
1076  'TABLE DATA ARE SPECIFIED', ntabs, &
1077  'TIMES BUT NTABLES IS SET TO', this%ntables
1078  call store_error(errmsg)
1079  end if
1080  do n = 1, this%nlakes
1081  if (this%ntabrow(n) > 0 .and. nboundchk(n) > 1) then
1082  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
1083  'TABLE DATA FOR LAKE', n, 'SPECIFIED', nboundchk(n), 'TIMES'
1084  call store_error(errmsg)
1085  end if
1086  end do
1087  else
1088  call store_error('REQUIRED TABLES BLOCK NOT FOUND.')
1089  end if
1090  !
1091  ! -- deallocate local storage
1092  deallocate (nboundchk)
1093  !
1094  ! -- write summary of lake_table error messages
1095  if (count_errors() > 0) then
1096  call this%parser%StoreErrorUnit()
1097  end if
1098  !
1099  ! -- convert laketables to vectors
1100  call this%laktables_to_vectors(laketables)
1101  !
1102  ! -- destroy laketables
1103  do n = 1, this%nlakes
1104  if (this%ntabrow(n) > 0) then
1105  deallocate (laketables(n)%tabstage)
1106  deallocate (laketables(n)%tabvolume)
1107  deallocate (laketables(n)%tabsarea)
1108  iconn = this%idxlakeconn(n)
1109  if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
1110  deallocate (laketables(n)%tabwarea)
1111  end if
1112  end if
1113  end do
1114  deallocate (laketables)
1115  end subroutine lak_read_tables
1116 
1117  !> @brief Copy the laketables structure data into flattened vectors that are
1118  !! stored in the memory manager
1119  !<
1120  subroutine laktables_to_vectors(this, laketables)
1121  class(laktype), intent(inout) :: this
1122  type(laktabtype), intent(in), dimension(:), contiguous :: laketables
1123  integer(I4B) :: n
1124  integer(I4B) :: ntabrows
1125  integer(I4B) :: j
1126  integer(I4B) :: ipos
1127  integer(I4B) :: iconn
1128  !
1129  ! -- allocate index array for lak tables
1130  call mem_allocate(this%ialaktab, this%nlakes + 1, 'IALAKTAB', this%memoryPath)
1131  !
1132  ! -- Move the laktables structure information into flattened arrays
1133  this%ialaktab(1) = 1
1134  do n = 1, this%nlakes
1135  ! -- ialaktab contains a pointer into the flattened lak table data
1136  this%ialaktab(n + 1) = this%ialaktab(n) + this%ntabrow(n)
1137  end do
1138  !
1139  ! -- Allocate vectors for storing all lake table data
1140  ntabrows = this%ialaktab(this%nlakes + 1) - 1
1141  call mem_allocate(this%tabstage, ntabrows, 'TABSTAGE', this%memoryPath)
1142  call mem_allocate(this%tabvolume, ntabrows, 'TABVOLUME', this%memoryPath)
1143  call mem_allocate(this%tabsarea, ntabrows, 'TABSAREA', this%memoryPath)
1144  call mem_allocate(this%tabwarea, ntabrows, 'TABWAREA', this%memoryPath)
1145  !
1146  ! -- Copy data from laketables into vectors
1147  do n = 1, this%nlakes
1148  j = 1
1149  do ipos = this%ialaktab(n), this%ialaktab(n + 1) - 1
1150  this%tabstage(ipos) = laketables(n)%tabstage(j)
1151  this%tabvolume(ipos) = laketables(n)%tabvolume(j)
1152  this%tabsarea(ipos) = laketables(n)%tabsarea(j)
1153  iconn = this%idxlakeconn(n)
1154  if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
1155  !
1156  ! -- tabwarea only filled for ictype 2 and 3
1157  this%tabwarea(ipos) = laketables(n)%tabwarea(j)
1158  else
1159  this%tabwarea(ipos) = dzero
1160  end if
1161  j = j + 1
1162  end do
1163  end do
1164  end subroutine laktables_to_vectors
1165 
1166  !> @brief Read the lake table for this package
1167  !<
1168  subroutine lak_read_table(this, ilak, filename, laketable)
1169  use constantsmodule, only: linelength
1170  use inputoutputmodule, only: openfile
1171  use simmodule, only: store_error, count_errors
1172  ! -- dummy
1173  class(laktype), intent(inout) :: this
1174  integer(I4B), intent(in) :: ilak
1175  character(len=*), intent(in) :: filename
1176  type(laktabtype), intent(inout) :: laketable
1177  ! -- local
1178  character(len=LINELENGTH) :: keyword
1179  integer(I4B) :: ierr
1180  logical(LGP) :: isfound, endOfBlock
1181  integer(I4B) :: iu
1182  integer(I4B) :: n
1183  integer(I4B) :: ipos
1184  integer(I4B) :: j
1185  integer(I4B) :: jmin
1186  integer(I4B) :: iconn
1187  real(DP) :: vol
1188  real(DP) :: sa
1189  real(DP) :: wa
1190  real(DP) :: v
1191  real(DP) :: v0
1192  type(blockparsertype) :: parser
1193  ! -- formats
1194  character(len=*), parameter :: fmttaberr = &
1195  &'(a,1x,i0,1x,a,1x,g15.6,1x,a,1x,i0,1x,a,1x,i0,1x,a,1x,g15.6,1x,a)'
1196  !
1197  ! -- initialize locals
1198  n = 0
1199  j = 0
1200  !
1201  ! -- open the table file
1202  iu = 0
1203  call openfile(iu, this%iout, filename, 'LAKE TABLE')
1204  call parser%Initialize(iu, this%iout)
1205  !
1206  ! -- get dimensions block
1207  call parser%GetBlock('DIMENSIONS', isfound, ierr, supportopenclose=.true.)
1208  !
1209  ! -- parse lak table dimensions block if detected
1210  if (isfound) then
1211  ! -- process the lake table dimension data
1212  if (this%iprpak /= 0) then
1213  write (this%iout, '(/1x,a)') &
1214  'PROCESSING '//trim(adjustl(this%text))//' DIMENSIONS'
1215  end if
1216  readdims: do
1217  call parser%GetNextLine(endofblock)
1218  if (endofblock) exit
1219  call parser%GetStringCaps(keyword)
1220  select case (keyword)
1221  case ('NROW')
1222  n = parser%GetInteger()
1223 
1224  if (n < 1) then
1225  write (errmsg, '(a)') 'LAKE TABLE NROW MUST BE > 0'
1226  call store_error(errmsg)
1227  end if
1228  case ('NCOL')
1229  j = parser%GetInteger()
1230 
1231  if (this%ictype(ilak) == 2 .or. this%ictype(ilak) == 3) then
1232  jmin = 4
1233  else
1234  jmin = 3
1235  end if
1236  if (j < jmin) then
1237  write (errmsg, '(a,1x,i0)') 'LAKE TABLE NCOL MUST BE >= ', jmin
1238  call store_error(errmsg)
1239  end if
1240  !
1241  case default
1242  write (errmsg, '(a,a)') &
1243  'UNKNOWN '//trim(this%text)//' DIMENSIONS KEYWORD: ', trim(keyword)
1244  call store_error(errmsg)
1245  end select
1246  end do readdims
1247  if (this%iprpak /= 0) then
1248  write (this%iout, '(1x,a)') &
1249  'END OF '//trim(adjustl(this%text))//' DIMENSIONS'
1250  end if
1251  else
1252  call store_error('REQUIRED DIMENSIONS BLOCK NOT FOUND.')
1253  end if
1254  !
1255  ! -- check that ncol and nrow have been specified
1256  if (n < 1) then
1257  write (errmsg, '(a)') &
1258  'NROW NOT SPECIFIED IN THE LAKE TABLE DIMENSIONS BLOCK'
1259  call store_error(errmsg)
1260  end if
1261  if (j < 1) then
1262  write (errmsg, '(a)') &
1263  'NCOL NOT SPECIFIED IN THE LAKE TABLE DIMENSIONS BLOCK'
1264  call store_error(errmsg)
1265  end if
1266  !
1267  ! -- only read the lake table data if n and j are specified to be greater
1268  ! than zero
1269  if (n * j > 0) then
1270  !
1271  ! -- allocate space
1272  this%ntabrow(ilak) = n
1273  allocate (laketable%tabstage(n))
1274  allocate (laketable%tabvolume(n))
1275  allocate (laketable%tabsarea(n))
1276  ipos = this%idxlakeconn(ilak)
1277  if (this%ictype(ipos) == 2 .or. this%ictype(ipos) == 3) then
1278  allocate (laketable%tabwarea(n))
1279  end if
1280  !
1281  ! -- get table block
1282  call parser%GetBlock('TABLE', isfound, ierr, supportopenclose=.true.)
1283  !
1284  ! -- parse well_connections block if detected
1285  if (isfound) then
1286  !
1287  ! -- process the table data
1288  if (this%iprpak /= 0) then
1289  write (this%iout, '(/1x,a)') &
1290  'PROCESSING '//trim(adjustl(this%text))//' TABLE'
1291  end if
1292  iconn = this%idxlakeconn(ilak)
1293  ipos = 0
1294  readtabledata: do
1295  call parser%GetNextLine(endofblock)
1296  if (endofblock) exit
1297  ipos = ipos + 1
1298  if (ipos > this%ntabrow(ilak)) then
1299  cycle readtabledata
1300  end if
1301  laketable%tabstage(ipos) = parser%GetDouble()
1302  laketable%tabvolume(ipos) = parser%GetDouble()
1303  laketable%tabsarea(ipos) = parser%GetDouble()
1304  if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
1305  laketable%tabwarea(ipos) = parser%GetDouble()
1306  end if
1307  end do readtabledata
1308  !
1309  if (this%iprpak /= 0) then
1310  write (this%iout, '(1x,a)') &
1311  'END OF '//trim(adjustl(this%text))//' TABLE'
1312  end if
1313  else
1314  call store_error('REQUIRED TABLE BLOCK NOT FOUND.')
1315  end if
1316  !
1317  ! -- error condition if number of rows read are not equal to nrow
1318  if (ipos /= this%ntabrow(ilak)) then
1319  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
1320  'NROW SET TO', this%ntabrow(ilak), 'BUT', ipos, 'ROWS WERE READ'
1321  call store_error(errmsg)
1322  end if
1323  !
1324  ! -- set lake bottom based on table if it is an embedded lake
1325  iconn = this%idxlakeconn(ilak)
1326  if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
1327  do n = 1, this%ntabrow(ilak)
1328  vol = laketable%tabvolume(n)
1329  sa = laketable%tabsarea(n)
1330  wa = laketable%tabwarea(n)
1331  vol = vol * sa * wa
1332  ! -- check if all entries are zero
1333  if (vol > dzero) exit
1334  ! -- set lake bottom
1335  this%lakebot(ilak) = laketable%tabstage(n)
1336  this%belev(ilak) = laketable%tabstage(n)
1337  end do
1338  ! -- set maximum surface area for rainfall
1339  n = this%ntabrow(ilak)
1340  this%sareamax(ilak) = laketable%tabsarea(n)
1341  end if
1342  !
1343  ! -- verify the table data
1344  do n = 2, this%ntabrow(ilak)
1345  v = laketable%tabstage(n)
1346  v0 = laketable%tabstage(n - 1)
1347  if (v <= v0) then
1348  write (errmsg, fmttaberr) &
1349  'TABLE STAGE ENTRY', n, '(', laketable%tabstage(n), ') FOR LAKE ', &
1350  ilak, 'MUST BE GREATER THAN THE PREVIOUS STAGE ENTRY', &
1351  n - 1, '(', laketable%tabstage(n - 1), ')'
1352  call store_error(errmsg)
1353  end if
1354  v = laketable%tabvolume(n)
1355  v0 = laketable%tabvolume(n - 1)
1356  if (v <= v0) then
1357  write (errmsg, fmttaberr) &
1358  'TABLE VOLUME ENTRY', n, '(', laketable%tabvolume(n), &
1359  ') FOR LAKE ', &
1360  ilak, 'MUST BE GREATER THAN THE PREVIOUS VOLUME ENTRY', &
1361  n - 1, '(', laketable%tabvolume(n - 1), ')'
1362  call store_error(errmsg)
1363  end if
1364  v = laketable%tabsarea(n)
1365  v0 = laketable%tabsarea(n - 1)
1366  if (v < v0) then
1367  write (errmsg, fmttaberr) &
1368  'TABLE SURFACE AREA ENTRY', n, '(', &
1369  laketable%tabsarea(n), ') FOR LAKE ', ilak, &
1370  'MUST BE GREATER THAN OR EQUAL TO THE PREVIOUS SURFACE AREA ENTRY', &
1371  n - 1, '(', laketable%tabsarea(n - 1), ')'
1372  call store_error(errmsg)
1373  end if
1374  iconn = this%idxlakeconn(ilak)
1375  if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
1376  v = laketable%tabwarea(n)
1377  v0 = laketable%tabwarea(n - 1)
1378  if (v < v0) then
1379  write (errmsg, fmttaberr) &
1380  'TABLE EXCHANGE AREA ENTRY', n, '(', &
1381  laketable%tabwarea(n), ') FOR LAKE ', ilak, &
1382  'MUST BE GREATER THAN OR EQUAL TO THE PREVIOUS EXCHANGE AREA '// &
1383  'ENTRY', n - 1, '(', laketable%tabwarea(n - 1), ')'
1384  call store_error(errmsg)
1385  end if
1386  end if
1387  end do
1388  end if
1389  !
1390  ! -- write summary of lake table error messages
1391  if (count_errors() > 0) then
1392  call parser%StoreErrorUnit()
1393  end if
1394  !
1395  ! Close the table file and clear other parser members
1396  call parser%Clear()
1397  end subroutine lak_read_table
1398 
1399  !> @brief Read the lake outlets for this package
1400  !<
1401  subroutine lak_read_outlets(this)
1402  use constantsmodule, only: linelength
1403  use simmodule, only: store_error, count_errors
1404  use timeseriesmanagermodule, only: read_value_or_time_series_adv
1405  ! -- dummy
1406  class(laktype), intent(inout) :: this
1407  ! -- local
1408  character(len=LINELENGTH) :: text, keyword
1409  character(len=LENBOUNDNAME) :: bndName
1410  character(len=9) :: citem
1411  integer(I4B) :: ierr, ival
1412  logical(LGP) :: isfound, endOfBlock
1413  integer(I4B) :: n
1414  integer(I4B) :: jj
1415  integer(I4B), dimension(:), pointer, contiguous :: nboundchk
1416  real(DP), pointer :: bndElem => null()
1417  !
1418  ! -- get well_connections block
1419  call this%parser%GetBlock('OUTLETS', isfound, ierr, &
1420  supportopenclose=.true., blockrequired=.false.)
1421  !
1422  ! -- parse outlets block if detected
1423  if (isfound) then
1424  if (this%noutlets > 0) then
1425  !
1426  ! -- allocate and initialize local variables
1427  allocate (nboundchk(this%noutlets))
1428  do n = 1, this%noutlets
1429  nboundchk(n) = 0
1430  end do
1431  !
1432  ! -- allocate outlet data using memory manager
1433  call mem_allocate(this%lakein, this%NOUTLETS, 'LAKEIN', this%memoryPath)
1434  call mem_allocate(this%lakeout, this%NOUTLETS, 'LAKEOUT', this%memoryPath)
1435  call mem_allocate(this%iouttype, this%NOUTLETS, 'IOUTTYPE', &
1436  this%memoryPath)
1437  call mem_allocate(this%outrate, this%NOUTLETS, 'OUTRATE', this%memoryPath)
1438  call mem_allocate(this%outinvert, this%NOUTLETS, 'OUTINVERT', &
1439  this%memoryPath)
1440  call mem_allocate(this%outwidth, this%NOUTLETS, 'OUTWIDTH', &
1441  this%memoryPath)
1442  call mem_allocate(this%outrough, this%NOUTLETS, 'OUTROUGH', &
1443  this%memoryPath)
1444  call mem_allocate(this%outslope, this%NOUTLETS, 'OUTSLOPE', &
1445  this%memoryPath)
1446  call mem_allocate(this%simoutrate, this%NOUTLETS, 'SIMOUTRATE', &
1447  this%memoryPath)
1448  !
1449  ! -- initialize outlet rate
1450  do n = 1, this%noutlets
1451  this%outrate(n) = dzero
1452  end do
1453  !
1454  ! -- process the lake connection data
1455  write (this%iout, '(/1x,a)') &
1456  'PROCESSING '//trim(adjustl(this%text))//' OUTLETS'
1457  readoutlet: do
1458  call this%parser%GetNextLine(endofblock)
1459  if (endofblock) exit
1460  n = this%parser%GetInteger()
1461 
1462  if (n < 1 .or. n > this%noutlets) then
1463  write (errmsg, '(a,1x,i0)') &
1464  'outletno MUST BE > 0 and <= ', this%noutlets
1465  call store_error(errmsg)
1466  cycle readoutlet
1467  end if
1468  !
1469  ! -- increment nboundchk
1470  nboundchk(n) = nboundchk(n) + 1
1471  !
1472  ! -- read outlet lakein
1473  ival = this%parser%GetInteger()
1474  if (ival < 1 .or. ival > this%nlakes) then
1475  write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
1476  'lakein FOR OUTLET ', n, 'MUST BE > 0 and <= ', this%nlakes
1477  call store_error(errmsg)
1478  cycle readoutlet
1479  end if
1480  this%lakein(n) = ival
1481  !
1482  ! -- read outlet lakeout
1483  ival = this%parser%GetInteger()
1484  if (ival < 0 .or. ival > this%nlakes) then
1485  write (errmsg, '(a,1x,i0,1x,a,1x,i0)') &
1486  'lakeout FOR OUTLET ', n, 'MUST BE >= 0 and <= ', this%nlakes
1487  call store_error(errmsg)
1488  cycle readoutlet
1489  end if
1490  this%lakeout(n) = ival
1491  !
1492  ! -- read ictype
1493  call this%parser%GetStringCaps(keyword)
1494  select case (keyword)
1495  case ('SPECIFIED')
1496  this%iouttype(n) = 0
1497  case ('MANNING')
1498  this%iouttype(n) = 1
1499  case ('WEIR')
1500  this%iouttype(n) = 2
1501  case default
1502  write (errmsg, '(a,1x,i0,1x,a,a,a)') &
1503  'UNKNOWN couttype FOR OUTLET ', n, '(', trim(keyword), ')'
1504  call store_error(errmsg)
1505  cycle readoutlet
1506  end select
1507  !
1508  ! -- build bndname for outlet
1509  write (citem, '(i9.9)') n
1510  bndname = 'OUTLET'//citem
1511  !
1512  ! -- set a few variables for timeseries aware variables
1513  jj = 1
1514  !
1515  ! -- outlet invert
1516  call this%parser%GetString(text)
1517  bndelem => this%outinvert(n)
1518  call read_value_or_time_series_adv(text, n, jj, bndelem, &
1519  this%packName, 'BND', &
1520  this%tsManager, this%iprpak, &
1521  'INVERT')
1522  !
1523  ! -- outlet width
1524  call this%parser%GetString(text)
1525  bndelem => this%outwidth(n)
1526  call read_value_or_time_series_adv(text, n, jj, bndelem, &
1527  this%packName, 'BND', &
1528  this%tsManager, this%iprpak, 'WIDTH')
1529  !
1530  ! -- outlet roughness
1531  call this%parser%GetString(text)
1532  bndelem => this%outrough(n)
1533  call read_value_or_time_series_adv(text, n, jj, bndelem, &
1534  this%packName, 'BND', &
1535  this%tsManager, this%iprpak, 'ROUGH')
1536  !
1537  ! -- outlet slope
1538  call this%parser%GetString(text)
1539  bndelem => this%outslope(n)
1540  call read_value_or_time_series_adv(text, n, jj, bndelem, &
1541  this%packName, 'BND', &
1542  this%tsManager, this%iprpak, 'SLOPE')
1543  end do readoutlet
1544  write (this%iout, '(1x,a)') 'END OF '//trim(adjustl(this%text))// &
1545  ' OUTLETS'
1546  !
1547  ! -- check for duplicate or missing outlets
1548  do n = 1, this%noutlets
1549  if (nboundchk(n) == 0) then
1550  write (errmsg, '(a,1x,i0)') 'NO DATA SPECIFIED FOR OUTLET', n
1551  call store_error(errmsg)
1552  else if (nboundchk(n) > 1) then
1553  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
1554  'DATA FOR OUTLET', n, 'SPECIFIED', nboundchk(n), 'TIMES'
1555  call store_error(errmsg)
1556  end if
1557  end do
1558  !
1559  ! -- deallocate local storage
1560  deallocate (nboundchk)
1561  else
1562  write (errmsg, '(a,1x,a)') &
1563  'AN OUTLETS BLOCK SHOULD NOT BE SPECIFIED IF NOUTLETS IS NOT', &
1564  'SPECIFIED OR IS SPECIFIED TO BE 0.'
1565  call store_error(errmsg)
1566  end if
1567  !
1568  else
1569  if (this%noutlets > 0) then
1570  call store_error('REQUIRED OUTLETS BLOCK NOT FOUND.')
1571  end if
1572  end if
1573  !
1574  ! -- write summary of lake_connection error messages
1575  ierr = count_errors()
1576  if (ierr > 0) then
1577  call this%parser%StoreErrorUnit()
1578  end if
1579  end subroutine lak_read_outlets
1580 
1581  !> @brief Read the dimensions for this package
1582  !<
1583  subroutine lak_read_dimensions(this)
1584  use constantsmodule, only: linelength
1585  use simmodule, only: store_error, count_errors
1586  ! -- dummy
1587  class(laktype), intent(inout) :: this
1588  ! -- local
1589  character(len=LINELENGTH) :: keyword
1590  integer(I4B) :: ierr
1591  logical(LGP) :: isfound, endOfBlock
1592  !
1593  ! -- initialize dimensions to -1
1594  this%nlakes = -1
1595  this%maxbound = -1
1596  !
1597  ! -- get dimensions block
1598  call this%parser%GetBlock('DIMENSIONS', isfound, ierr, &
1599  supportopenclose=.true.)
1600  !
1601  ! -- parse dimensions block if detected
1602  if (isfound) then
1603  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
1604  ' DIMENSIONS'
1605  do
1606  call this%parser%GetNextLine(endofblock)
1607  if (endofblock) exit
1608  call this%parser%GetStringCaps(keyword)
1609  select case (keyword)
1610  case ('NLAKES')
1611  this%nlakes = this%parser%GetInteger()
1612  write (this%iout, '(4x,a,i7)') 'NLAKES = ', this%nlakes
1613  case ('NOUTLETS')
1614  this%noutlets = this%parser%GetInteger()
1615  write (this%iout, '(4x,a,i7)') 'NOUTLETS = ', this%noutlets
1616  case ('NTABLES')
1617  this%ntables = this%parser%GetInteger()
1618  write (this%iout, '(4x,a,i7)') 'NTABLES = ', this%ntables
1619  case default
1620  write (errmsg, '(a,a)') &
1621  'UNKNOWN '//trim(this%text)//' DIMENSION: ', trim(keyword)
1622  call store_error(errmsg)
1623  end select
1624  end do
1625  write (this%iout, '(1x,a)') &
1626  'END OF '//trim(adjustl(this%text))//' DIMENSIONS'
1627  else
1628  call store_error('REQUIRED DIMENSIONS BLOCK NOT FOUND.')
1629  end if
1630  !
1631  if (this%nlakes < 0) then
1632  write (errmsg, '(a)') &
1633  'NLAKES WAS NOT SPECIFIED OR WAS SPECIFIED INCORRECTLY.'
1634  call store_error(errmsg)
1635  end if
1636  !
1637  ! -- stop if errors were encountered in the DIMENSIONS block
1638  if (count_errors() > 0) then
1639  call this%parser%StoreErrorUnit()
1640  end if
1641  !
1642  ! -- read lakes block
1643  call this%lak_read_lakes()
1644  !
1645  ! -- read lake_connections block
1646  call this%lak_read_lake_connections()
1647  !
1648  ! -- read tables block
1649  call this%lak_read_tables()
1650  !
1651  ! -- read outlets block
1652  call this%lak_read_outlets()
1653  !
1654  ! -- Call define_listlabel to construct the list label that is written
1655  ! when PRINT_INPUT option is used.
1656  call this%define_listlabel()
1657  !
1658  ! -- setup the budget object
1659  call this%lak_setup_budobj()
1660  !
1661  ! -- setup the stage table object
1662  call this%lak_setup_tableobj()
1663  end subroutine lak_read_dimensions
1664 
1665  !> @brief Read the initial parameters for this package
1666  !<
1667  subroutine lak_read_initial_attr(this)
1668  use constantsmodule, only: linelength
1669  use memoryhelpermodule, only: create_mem_path
1670  use simmodule, only: store_error, count_errors
1671  use timeseriesmanagermodule, only: read_value_or_time_series_adv
1672  ! -- dummy
1673  class(laktype), intent(inout) :: this
1674  ! -- local
1675  character(len=LINELENGTH) :: text
1676  integer(I4B) :: j, jj, n
1677  integer(I4B) :: nn
1678  integer(I4B) :: idx
1679  real(DP) :: top
1680  real(DP) :: bot
1681  real(DP) :: k
1682  real(DP) :: area
1683  real(DP) :: length
1684  real(DP) :: s
1685  real(DP) :: dx
1686  real(DP) :: c
1687  real(DP) :: sa
1688  real(DP) :: wa
1689  real(DP) :: v
1690  real(DP) :: fact
1691  real(DP) :: c1
1692  real(DP) :: c2
1693  real(DP), allocatable, dimension(:) :: clb, caq
1694  character(len=14) :: cbedleak
1695  character(len=14) :: cbedcond
1696  character(len=10), dimension(0:3) :: ctype
1697  character(len=15) :: nodestr
1698  real(DP), pointer :: bndElem => null()
1699  ! -- data
1700  data ctype(0)/'VERTICAL '/
1701  data ctype(1)/'HORIZONTAL'/
1702  data ctype(2)/'EMBEDDEDH '/
1703  data ctype(3)/'EMBEDDEDV '/
1704  !
1705  ! -- initialize xnewpak and set stage
1706  do n = 1, this%nlakes
1707  this%xnewpak(n) = this%strt(n)
1708  write (text, '(g15.7)') this%strt(n)
1709  jj = 1 ! For STAGE
1710  bndelem => this%stage(n)
1711  call read_value_or_time_series_adv(text, n, jj, bndelem, this%packName, &
1712  'BND', this%tsManager, this%iprpak, &
1713  'STAGE')
1714  end do
1715  !
1716  ! -- initialize status (iboundpak) of lakes to active
1717  do n = 1, this%nlakes
1718  if (this%status(n) == 'CONSTANT') then
1719  this%iboundpak(n) = -1
1720  else if (this%status(n) == 'INACTIVE') then
1721  this%iboundpak(n) = 0
1722  else if (this%status(n) == 'ACTIVE ') then
1723  this%iboundpak(n) = 1
1724  end if
1725  end do
1726  !
1727  ! -- set boundname for each connection
1728  if (this%inamedbound /= 0) then
1729  do n = 1, this%nlakes
1730  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1731  this%boundname(j) = this%lakename(n)
1732  end do
1733  end do
1734  end if
1735  !
1736  ! -- copy boundname into boundname_cst
1737  call this%copy_boundname()
1738  !
1739  ! -- set pointer to gwf iss and gwf hk
1740  call mem_setptr(this%gwfiss, 'ISS', create_mem_path(this%name_model))
1741  call mem_setptr(this%gwfk11, 'K11', create_mem_path(this%name_model, 'NPF'))
1742  call mem_setptr(this%gwfk33, 'K33', create_mem_path(this%name_model, 'NPF'))
1743  call mem_setptr(this%gwfik33, 'IK33', create_mem_path(this%name_model, 'NPF'))
1744  call mem_setptr(this%gwfsat, 'SAT', create_mem_path(this%name_model, 'NPF'))
1745  !
1746  ! -- allocate temporary storage
1747  allocate (clb(this%MAXBOUND))
1748  allocate (caq(this%MAXBOUND))
1749  !
1750  ! -- calculate saturated conductance for each connection
1751  do n = 1, this%nlakes
1752  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1753  nn = this%cellid(j)
1754  top = this%dis%top(nn)
1755  bot = this%dis%bot(nn)
1756  ! vertical connection
1757  if (this%ictype(j) == 0) then
1758  area = this%dis%area(nn)
1759  this%sarea(j) = area
1760  this%warea(j) = area
1761  this%sareamax(n) = this%sareamax(n) + area
1762  if (this%gwfik33 == 0) then
1763  k = this%gwfk11(nn)
1764  else
1765  k = this%gwfk33(nn)
1766  end if
1767  length = dhalf * (top - bot)
1768  ! horizontal connection
1769  else if (this%ictype(j) == 1) then
1770  area = (this%telev(j) - this%belev(j)) * this%connwidth(j)
1771  ! -- recalculate area if connected cell is confined and lake
1772  ! connection top and bot are equal to the cell top and bot
1773  if (top == this%telev(j) .and. bot == this%belev(j)) then
1774  if (this%icelltype(nn) == 0) then
1775  area = this%gwfsat(nn) * (top - bot) * this%connwidth(j)
1776  end if
1777  end if
1778  this%sarea(j) = dzero
1779  this%warea(j) = area
1780  this%sareamax(n) = this%sareamax(n) + dzero
1781  k = this%gwfk11(nn)
1782  length = this%connlength(j)
1783  ! embedded horizontal connection
1784  else if (this%ictype(j) == 2) then
1785  area = done
1786  this%sarea(j) = dzero
1787  this%warea(j) = area
1788  this%sareamax(n) = this%sareamax(n) + dzero
1789  k = this%gwfk11(nn)
1790  length = this%connlength(j)
1791  ! embedded vertical connection
1792  else if (this%ictype(j) == 3) then
1793  area = done
1794  this%sarea(j) = dzero
1795  this%warea(j) = area
1796  this%sareamax(n) = this%sareamax(n) + dzero
1797  if (this%gwfik33 == 0) then
1798  k = this%gwfk11(nn)
1799  else
1800  k = this%gwfk33(nn)
1801  end if
1802  length = this%connlength(j)
1803  end if
1804  if (is_close(this%bedleak(j), dnodata)) then
1805  clb(j) = dnodata
1806  else if (this%bedleak(j) > dzero) then
1807  clb(j) = done / this%bedleak(j)
1808  else
1809  clb(j) = dzero
1810  end if
1811  if (k > dzero) then
1812  caq(j) = length / k
1813  else
1814  caq(j) = dzero
1815  end if
1816  if (is_close(this%bedleak(j), dnodata)) then
1817  this%satcond(j) = area / caq(j)
1818  else if (clb(j) * caq(j) > dzero) then
1819  this%satcond(j) = area / (clb(j) + caq(j))
1820  else
1821  this%satcond(j) = dzero
1822  end if
1823  end do
1824  end do
1825  !
1826  ! -- write a summary of the conductance
1827  if (this%iprpak > 0) then
1828  write (this%iout, '(//,29x,a,/)') &
1829  'INTERFACE CONDUCTANCE BETWEEN LAKE AND AQUIFER CELLS'
1830  write (this%iout, '(1x,a)') &
1831  & ' LAKE CONNECTION CONNECTION LAKEBED'// &
1832  & ' C O N D U C T A N C E S '
1833  write (this%iout, '(1x,a)') &
1834  & ' NUMBER NUMBER CELLID DIRECTION LEAKANCE'// &
1835  & ' LAKEBED AQUIFER COMBINED'
1836  write (this%iout, "(1x,108('-'))")
1837  do n = 1, this%nlakes
1838  idx = 0
1839  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
1840  idx = idx + 1
1841  fact = done
1842  if (this%ictype(j) == 1) then
1843  fact = this%telev(j) - this%belev(j)
1844  if (abs(fact) > dzero) then
1845  fact = done / fact
1846  end if
1847  end if
1848  nn = this%cellid(j)
1849  area = this%warea(j)
1850  c1 = dzero
1851  if (is_close(clb(j), dnodata)) then
1852  cbedleak = ' NONE '
1853  cbedcond = ' NONE '
1854  else if (clb(j) > dzero) then
1855  c1 = area * fact / clb(j)
1856  write (cbedleak, '(g14.5)') this%bedleak(j)
1857  write (cbedcond, '(g14.5)') c1
1858  else
1859  write (cbedleak, '(g14.5)') c1
1860  write (cbedcond, '(g14.5)') c1
1861  end if
1862  c2 = dzero
1863  if (caq(j) > dzero) then
1864  c2 = area * fact / caq(j)
1865  end if
1866  call this%dis%noder_to_string(nn, nodestr)
1867  write (this%iout, &
1868  '(1x,i10,1x,i10,1x,a15,1x,a10,2(1x,a14),2(1x,g14.5))') &
1869  n, idx, nodestr, ctype(this%ictype(j)), cbedleak, &
1870  cbedcond, c2, this%satcond(j) * fact
1871  end do
1872  end do
1873  write (this%iout, "(1x,108('-'))")
1874  write (this%iout, '(1x,a)') &
1875  'IF VERTICAL CONNECTION, CONDUCTANCE (L^2/T) IS &
1876  &BETWEEN AQUIFER CELL AND OVERLYING LAKE CELL.'
1877  write (this%iout, '(1x,a)') &
1878  'IF HORIZONTAL CONNECTION, CONDUCTANCES ARE PER &
1879  &UNIT SATURATED THICKNESS (L/T).'
1880  write (this%iout, '(1x,a)') &
1881  'IF EMBEDDED CONNECTION, CONDUCTANCES ARE PER &
1882  &UNIT EXCHANGE AREA (1/T).'
1883  !
1884  ! write(this%iout,*) n, idx, nodestr, this%sarea(j), this%warea(j)
1885  !
1886  ! -- calculate stage, surface area, wetted area, volume relation
1887  do n = 1, this%nlakes
1888  write (this%iout, '(//1x,a,1x,i10)') 'STAGE/VOLUME RELATION FOR LAKE ', n
1889  write (this%iout, '(/1x,5(a14))') ' STAGE', ' SURFACE AREA', &
1890  & ' WETTED AREA', ' CONDUCTANCE', &
1891  & ' VOLUME'
1892  write (this%iout, "(1x,70('-'))")
1893  dx = (this%laketop(n) - this%lakebot(n)) / 150.
1894  s = this%lakebot(n)
1895  do j = 1, 151
1896  call this%lak_calculate_conductance(n, s, c)
1897  call this%lak_calculate_sarea(n, s, sa)
1898  call this%lak_calculate_warea(n, s, wa, s)
1899  call this%lak_calculate_vol(n, s, v)
1900  write (this%iout, '(1x,5(E14.5))') s, sa, wa, c, v
1901  s = s + dx
1902  end do
1903  write (this%iout, "(1x,70('-'))")
1904  !
1905  write (this%iout, '(//1x,a,1x,i10)') 'STAGE/VOLUME RELATION FOR LAKE ', n
1906  write (this%iout, '(/1x,4(a14))') ' ', ' ', &
1907  & ' CALCULATED', ' STAGE'
1908  write (this%iout, '(1x,4(a14))') ' STAGE', ' VOLUME', &
1909  & ' STAGE', ' DIFFERENCE'
1910  write (this%iout, "(1x,56('-'))")
1911  s = this%lakebot(n) - dx
1912  do j = 1, 156
1913  call this%lak_calculate_vol(n, s, v)
1914  call this%lak_vol2stage(n, v, c)
1915  write (this%iout, '(1x,4(E14.5))') s, v, c, s - c
1916  s = s + dx
1917  end do
1918  write (this%iout, "(1x,56('-'))")
1919  end do
1920  end if
1921  !
1922  ! -- finished with pointer to gwf hydraulic conductivity
1923  this%gwfk11 => null()
1924  this%gwfk33 => null()
1925  this%gwfsat => null()
1926  this%gwfik33 => null()
1927  !
1928  ! -- deallocate temporary storage
1929  deallocate (clb)
1930  deallocate (caq)
1931  end subroutine lak_read_initial_attr
1932 
1933  !> @brief Perform linear interpolation of two vectors.
1934  !!
1935  !! Function assumes x data is sorted in ascending order
1936  !<
1937  subroutine lak_linear_interpolation(this, n, x, y, z, v)
1938  ! -- dummy
1939  class(laktype), intent(inout) :: this
1940  integer(I4B), intent(in) :: n
1941  real(DP), dimension(n), intent(in) :: x
1942  real(DP), dimension(n), intent(in) :: y
1943  real(DP), intent(in) :: z
1944  real(DP), intent(inout) :: v
1945  ! -- local
1946  integer(I4B) :: i
1947  real(DP) :: dx, dydx
1948  ! code
1949  v = dzero
1950  ! below bottom of range - set to lowest value
1951  if (z <= x(1)) then
1952  v = y(1)
1953  ! above highest value
1954  ! slope calculated from interval between n and n-1
1955  else if (z > x(n)) then
1956  dx = x(n) - x(n - 1)
1957  dydx = dzero
1958  if (abs(dx) > dzero) then
1959  dydx = (y(n) - y(n - 1)) / dx
1960  end if
1961  dx = (z - x(n))
1962  v = y(n) + dydx * dx
1963  ! between lowest and highest value in current interval
1964  else
1965  do i = 2, n
1966  dx = x(i) - x(i - 1)
1967  dydx = dzero
1968  if (z >= x(i - 1) .and. z <= x(i)) then
1969  if (abs(dx) > dzero) then
1970  dydx = (y(i) - y(i - 1)) / dx
1971  end if
1972  dx = (z - x(i - 1))
1973  v = y(i - 1) + dydx * dx
1974  exit
1975  end if
1976  end do
1977  end if
1978  end subroutine lak_linear_interpolation
1979 
1980  !> @brief Calculate the surface area of a lake at a given stage
1981  !<
1982  subroutine lak_calculate_sarea(this, ilak, stage, sarea)
1983  ! -- dummy
1984  class(laktype), intent(inout) :: this
1985  integer(I4B), intent(in) :: ilak
1986  real(DP), intent(in) :: stage
1987  real(DP), intent(inout) :: sarea
1988  ! -- local
1989  integer(I4B) :: i
1990  integer(I4B) :: ifirst
1991  integer(I4B) :: ilast
1992  real(DP) :: topl
1993  real(DP) :: botl
1994  real(DP) :: sat
1995  real(DP) :: sa
1996  !
1997  sarea = dzero
1998  i = this%ntabrow(ilak)
1999  if (i > 0) then
2000  ifirst = this%ialaktab(ilak)
2001  ilast = this%ialaktab(ilak + 1) - 1
2002  if (stage <= this%tabstage(ifirst)) then
2003  sarea = this%tabsarea(ifirst)
2004  else if (stage >= this%tabstage(ilast)) then
2005  sarea = this%tabsarea(ilast)
2006  else
2007  call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
2008  this%tabsarea(ifirst:ilast), &
2009  stage, sarea)
2010  end if
2011  else
2012  do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2013  topl = this%telev(i)
2014  botl = this%belev(i)
2015  sat = squadraticsaturation(topl, botl, stage)
2016  sa = sat * this%sarea(i)
2017  sarea = sarea + sa
2018  end do
2019  end if
2020  end subroutine lak_calculate_sarea
2021 
2022  !> @brief Calculate the wetted area of a lake at a given stage.
2023  !<
2024  subroutine lak_calculate_warea(this, ilak, stage, warea, hin)
2025  ! -- dummy
2026  class(laktype), intent(inout) :: this
2027  integer(I4B), intent(in) :: ilak
2028  real(DP), intent(in) :: stage
2029  real(DP), intent(inout) :: warea
2030  real(DP), optional, intent(inout) :: hin
2031  ! -- local
2032  integer(I4B) :: i
2033  integer(I4B) :: igwfnode
2034  real(DP) :: head
2035  real(DP) :: wa
2036  !
2037  warea = dzero
2038  do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2039  if (present(hin)) then
2040  head = hin
2041  else
2042  igwfnode = this%cellid(i)
2043  head = this%xnew(igwfnode)
2044  end if
2045  call this%lak_calculate_conn_warea(ilak, i, stage, head, wa)
2046  warea = warea + wa
2047  end do
2048  end subroutine lak_calculate_warea
2049 
2050  !> @brief Calculate the wetted area of a lake connection at a given stage
2051  !<
2052  subroutine lak_calculate_conn_warea(this, ilak, iconn, stage, head, wa)
2053  ! -- dummy
2054  class(laktype), intent(inout) :: this
2055  integer(I4B), intent(in) :: ilak
2056  integer(I4B), intent(in) :: iconn
2057  real(DP), intent(in) :: stage
2058  real(DP), intent(in) :: head
2059  real(DP), intent(inout) :: wa
2060  ! -- local
2061  integer(I4B) :: i
2062  integer(I4B) :: ifirst
2063  integer(I4B) :: ilast
2064  integer(I4B) :: node
2065  real(DP) :: topl
2066  real(DP) :: botl
2067  real(DP) :: vv
2068  real(DP) :: sat
2069  !
2070  wa = dzero
2071  topl = this%telev(iconn)
2072  botl = this%belev(iconn)
2073  call this%lak_calculate_cond_head(iconn, stage, head, vv)
2074  if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
2075  if (vv > topl) vv = topl
2076  i = this%ntabrow(ilak)
2077  ifirst = this%ialaktab(ilak)
2078  ilast = this%ialaktab(ilak + 1) - 1
2079  if (vv <= this%tabstage(ifirst)) then
2080  wa = this%tabwarea(ifirst)
2081  else if (vv >= this%tabstage(ilast)) then
2082  wa = this%tabwarea(ilast)
2083  else
2084  call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
2085  this%tabwarea(ifirst:ilast), &
2086  vv, wa)
2087  end if
2088  else
2089  node = this%cellid(iconn)
2090  ! -- confined cell
2091  if (this%icelltype(node) == 0) then
2092  sat = done
2093  ! -- convertible cell
2094  else
2095  sat = squadraticsaturation(topl, botl, vv)
2096  end if
2097  wa = sat * this%warea(iconn)
2098  end if
2099  end subroutine lak_calculate_conn_warea
2100 
2101  !> @brief Calculate the volume of a lake at a given stage
2102  !<
2103  subroutine lak_calculate_vol(this, ilak, stage, volume)
2104  ! -- dummy
2105  class(laktype), intent(inout) :: this
2106  integer(I4B), intent(in) :: ilak
2107  real(DP), intent(in) :: stage
2108  real(DP), intent(inout) :: volume
2109  ! -- local
2110  integer(I4B) :: i
2111  integer(I4B) :: ifirst
2112  integer(I4B) :: ilast
2113  real(DP) :: topl
2114  real(DP) :: botl
2115  real(DP) :: ds
2116  real(DP) :: sa
2117  real(DP) :: v
2118  real(DP) :: sat
2119  !
2120  volume = dzero
2121  i = this%ntabrow(ilak)
2122  if (i > 0) then
2123  ifirst = this%ialaktab(ilak)
2124  ilast = this%ialaktab(ilak + 1) - 1
2125  if (stage <= this%tabstage(ifirst)) then
2126  volume = this%tabvolume(ifirst)
2127  else if (stage >= this%tabstage(ilast)) then
2128  ds = stage - this%tabstage(ilast)
2129  sa = this%tabsarea(ilast)
2130  volume = this%tabvolume(ilast) + ds * sa
2131  else
2132  call this%lak_linear_interpolation(i, this%tabstage(ifirst:ilast), &
2133  this%tabvolume(ifirst:ilast), &
2134  stage, volume)
2135  end if
2136  else
2137  do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2138  topl = this%telev(i)
2139  botl = this%belev(i)
2140  sat = squadraticsaturation(topl, botl, stage)
2141  sa = sat * this%sarea(i)
2142  if (stage < botl) then
2143  v = dzero
2144  else if (stage > botl .and. stage < topl) then
2145  v = sa * (stage - botl)
2146  else
2147  v = sa * (topl - botl) + sa * (stage - topl)
2148  end if
2149  volume = volume + v
2150  end do
2151  end if
2152  end subroutine lak_calculate_vol
2153 
2154  !> @brief Calculate the total conductance for a lake at a provided stage
2155  !<
2156  subroutine lak_calculate_conductance(this, ilak, stage, conductance)
2157  ! -- dummy
2158  class(laktype), intent(inout) :: this
2159  integer(I4B), intent(in) :: ilak
2160  real(DP), intent(in) :: stage
2161  real(DP), intent(inout) :: conductance
2162  ! -- local
2163  integer(I4B) :: i
2164  real(DP) :: c
2165  !
2166  conductance = dzero
2167  do i = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2168  call this%lak_calculate_conn_conductance(ilak, i, stage, stage, c)
2169  conductance = conductance + c
2170  end do
2171  end subroutine lak_calculate_conductance
2172 
2173  !> @brief Calculate the controlling lake stage or groundwater head used to
2174  !! calculate the conductance for a lake connection from a provided stage and
2175  !! groundwater head
2176  !<
2177  subroutine lak_calculate_cond_head(this, iconn, stage, head, vv)
2178  ! -- dummy
2179  class(laktype), intent(inout) :: this
2180  integer(I4B), intent(in) :: iconn
2181  real(DP), intent(in) :: stage
2182  real(DP), intent(in) :: head
2183  real(DP), intent(inout) :: vv
2184  ! -- local
2185  real(DP) :: ss
2186  real(DP) :: hh
2187  real(DP) :: topl
2188  real(DP) :: botl
2189  !
2190  topl = this%telev(iconn)
2191  botl = this%belev(iconn)
2192  ss = min(stage, topl)
2193  hh = min(head, topl)
2194  if (this%igwhcopt > 0) then
2195  vv = hh
2196  else if (this%inewton > 0) then
2197  vv = max(ss, hh)
2198  else
2199  vv = dhalf * (ss + hh)
2200  end if
2201  end subroutine lak_calculate_cond_head
2202 
2203  !> @brief Calculate the conductance for a lake connection at a provided stage
2204  !! and groundwater head
2205  !<
2206  subroutine lak_calculate_conn_conductance(this, ilak, iconn, stage, head, cond)
2207  ! -- dummy
2208  class(laktype), intent(inout) :: this
2209  integer(I4B), intent(in) :: ilak
2210  integer(I4B), intent(in) :: iconn
2211  real(DP), intent(in) :: stage
2212  real(DP), intent(in) :: head
2213  real(DP), intent(inout) :: cond
2214  ! -- local
2215  integer(I4B) :: node
2216  !real(DP) :: ss
2217  !real(DP) :: hh
2218  real(DP) :: vv
2219  real(DP) :: topl
2220  real(DP) :: botl
2221  real(DP) :: sat
2222  real(DP) :: wa
2223  real(DP) :: vscratio
2224  !
2225  cond = dzero
2226  vscratio = done
2227  topl = this%telev(iconn)
2228  botl = this%belev(iconn)
2229  call this%lak_calculate_cond_head(iconn, stage, head, vv)
2230  sat = squadraticsaturation(topl, botl, vv)
2231  ! vertical connection
2232  ! use full saturated conductance if top and bottom of the lake connection
2233  ! are equal
2234  if (this%ictype(iconn) == 0) then
2235  if (abs(topl - botl) < dprec) then
2236  sat = done
2237  end if
2238  ! horizontal connection
2239  ! use full saturated conductance if the connected cell is not convertible
2240  else if (this%ictype(iconn) == 1) then
2241  node = this%cellid(iconn)
2242  if (this%icelltype(node) == 0) then
2243  sat = done
2244  end if
2245  ! embedded connection
2246  else if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
2247  node = this%cellid(iconn)
2248  if (this%icelltype(node) == 0) then
2249  vv = this%telev(iconn)
2250  call this%lak_calculate_conn_warea(ilak, iconn, vv, vv, wa)
2251  else
2252  call this%lak_calculate_conn_warea(ilak, iconn, stage, head, wa)
2253  end if
2254  sat = wa
2255  end if
2256  !
2257  ! -- account for viscosity effects (if vsc active)
2258  if (this%ivsc == 1) then
2259  ! flow from lake to aquifer
2260  if (stage > head) then
2261  vscratio = this%viscratios(1, iconn)
2262  ! flow from aquifer to lake
2263  else
2264  vscratio = this%viscratios(2, iconn)
2265  end if
2266  end if
2267  cond = sat * this%satcond(iconn) * vscratio
2268  end subroutine lak_calculate_conn_conductance
2269 
2270  !> @brief Calculate the total groundwater-lake flow at a provided stage
2271  !<
2272  subroutine lak_calculate_exchange(this, ilak, stage, totflow)
2273  ! -- dummy
2274  class(laktype), intent(inout) :: this
2275  integer(I4B), intent(in) :: ilak
2276  real(DP), intent(in) :: stage
2277  real(DP), intent(inout) :: totflow
2278  ! -- local
2279  integer(I4B) :: j
2280  integer(I4B) :: igwfnode
2281  real(DP) :: flow
2282  real(DP) :: hgwf
2283  !
2284  totflow = dzero
2285  do j = this%idxlakeconn(ilak), this%idxlakeconn(ilak + 1) - 1
2286  igwfnode = this%cellid(j)
2287  hgwf = this%xnew(igwfnode)
2288  call this%lak_calculate_conn_exchange(ilak, j, stage, hgwf, flow)
2289  totflow = totflow + flow
2290  end do
2291  end subroutine lak_calculate_exchange
2292 
2293  !> @brief Calculate the groundwater-lake flow at a provided stage and
2294  !! groundwater head
2295  !<
2296  subroutine lak_calculate_conn_exchange(this, ilak, iconn, stage, head, flow, &
2297  gwfhcof, gwfrhs)
2298  ! -- dummy
2299  class(laktype), intent(inout) :: this
2300  integer(I4B), intent(in) :: ilak
2301  integer(I4B), intent(in) :: iconn
2302  real(DP), intent(in) :: stage
2303  real(DP), intent(in) :: head
2304  real(DP), intent(inout) :: flow
2305  real(DP), intent(inout), optional :: gwfhcof
2306  real(DP), intent(inout), optional :: gwfrhs
2307  ! -- local
2308  real(DP) :: botl
2309  real(DP) :: cond
2310  real(DP) :: ss
2311  real(DP) :: hh
2312  real(DP) :: gwfhcof0
2313  real(DP) :: gwfrhs0
2314  !
2315  flow = dzero
2316  call this%lak_calculate_conn_conductance(ilak, iconn, stage, head, cond)
2317  botl = this%belev(iconn)
2318  !
2319  ! -- Set ss to stage or botl
2320  if (stage >= botl) then
2321  ss = stage
2322  else
2323  ss = botl
2324  end if
2325  !
2326  ! -- set hh to head or botl
2327  if (head >= botl) then
2328  hh = head
2329  else
2330  hh = botl
2331  end if
2332  !
2333  ! -- calculate flow, positive into lake
2334  flow = cond * (hh - ss)
2335  !
2336  ! -- Calculate gwfhcof and gwfrhs
2337  if (head >= botl) then
2338  gwfhcof0 = -cond
2339  gwfrhs0 = -cond * ss
2340  else
2341  gwfhcof0 = dzero
2342  gwfrhs0 = flow
2343  end if
2344  !
2345  ! Add density contributions, if active
2346  if (this%idense /= 0) then
2347  call this%lak_calculate_density_exchange(iconn, stage, head, cond, botl, &
2348  flow, gwfhcof0, gwfrhs0)
2349  end if
2350  !
2351  ! -- If present update gwfhcof and gwfrhs
2352  if (present(gwfhcof)) gwfhcof = gwfhcof0
2353  if (present(gwfrhs)) gwfrhs = gwfrhs0
2354  end subroutine lak_calculate_conn_exchange
2355 
2356  !> @brief Calculate the groundwater-lake flow at a provided stage and
2357  !! groundwater head
2358  !<
2359  subroutine lak_estimate_conn_exchange(this, iflag, ilak, iconn, idry, stage, &
2360  head, flow, source, gwfhcof, gwfrhs)
2361  ! -- dummy
2362  class(laktype), intent(inout) :: this
2363  integer(I4B), intent(in) :: iflag
2364  integer(I4B), intent(in) :: ilak
2365  integer(I4B), intent(in) :: iconn
2366  integer(I4B), intent(inout) :: idry
2367  real(DP), intent(in) :: stage
2368  real(DP), intent(in) :: head
2369  real(DP), intent(inout) :: flow
2370  real(DP), intent(inout) :: source
2371  real(DP), intent(inout), optional :: gwfhcof
2372  real(DP), intent(inout), optional :: gwfrhs
2373  ! -- local
2374  real(DP) :: gwfhcof0, gwfrhs0
2375  !
2376  flow = dzero
2377  idry = 0
2378  call this%lak_calculate_conn_exchange(ilak, iconn, stage, head, flow, &
2379  gwfhcof0, gwfrhs0)
2380  if (iflag == 1) then
2381  if (flow > dzero) then
2382  source = source + flow
2383  end if
2384  else if (iflag == 2) then
2385  if (-flow > source) then
2386  flow = -source
2387  source = dzero
2388  idry = 1
2389  else if (flow < dzero) then
2390  source = source + flow
2391  end if
2392  end if
2393  !
2394  ! -- Set gwfhcof and gwfrhs if present
2395  if (present(gwfhcof)) gwfhcof = gwfhcof0
2396  if (present(gwfrhs)) gwfrhs = gwfrhs0
2397  end subroutine lak_estimate_conn_exchange
2398 
2399  !> @brief Calculate the storage change in a lake based on provided stages
2400  !! and a passed delt
2401  !<
2402  subroutine lak_calculate_storagechange(this, ilak, stage, stage0, delt, dvr)
2403  ! -- dummy
2404  class(laktype), intent(inout) :: this
2405  integer(I4B), intent(in) :: ilak
2406  real(DP), intent(in) :: stage
2407  real(DP), intent(in) :: stage0
2408  real(DP), intent(in) :: delt
2409  real(DP), intent(inout) :: dvr
2410  ! -- local
2411  real(DP) :: v
2412  real(DP) :: v0
2413  !
2414  dvr = dzero
2415  if (this%gwfiss /= 1) then
2416  call this%lak_calculate_vol(ilak, stage, v)
2417  call this%lak_calculate_vol(ilak, stage0, v0)
2418  dvr = (v0 - v) / delt
2419  end if
2420  end subroutine lak_calculate_storagechange
2421 
2422  !> @brief Calculate the rainfall for a lake
2423  !<
2424  subroutine lak_calculate_rainfall(this, ilak, stage, ra)
2425  ! -- dummy
2426  class(laktype), intent(inout) :: this
2427  integer(I4B), intent(in) :: ilak
2428  real(DP), intent(in) :: stage
2429  real(DP), intent(inout) :: ra
2430  ! -- local
2431  integer(I4B) :: iconn
2432  real(DP) :: sa
2433  !
2434  ! -- rainfall
2435  iconn = this%idxlakeconn(ilak)
2436  if (this%ictype(iconn) == 2 .or. this%ictype(iconn) == 3) then
2437  sa = this%sareamax(ilak)
2438  else
2439  call this%lak_calculate_sarea(ilak, stage, sa)
2440  end if
2441  ra = this%rainfall(ilak) * sa
2442  end subroutine lak_calculate_rainfall
2443 
2444  !> @brief Calculate runoff to a lake
2445  !<
2446  subroutine lak_calculate_runoff(this, ilak, ro)
2447  ! -- dummy
2448  class(laktype), intent(inout) :: this
2449  integer(I4B), intent(in) :: ilak
2450  real(DP), intent(inout) :: ro
2451  !
2452  ! -- runoff
2453  ro = this%runoff(ilak)
2454  end subroutine lak_calculate_runoff
2455 
2456  !> @brief Calculate specified inflow to a lake
2457  !<
2458  subroutine lak_calculate_inflow(this, ilak, qin)
2459  ! -- dummy
2460  class(laktype), intent(inout) :: this
2461  integer(I4B), intent(in) :: ilak
2462  real(DP), intent(inout) :: qin
2463  !
2464  ! -- inflow to lake
2465  qin = this%inflow(ilak)
2466  end subroutine lak_calculate_inflow
2467 
2468  !> @brief Calculate the external flow terms to a lake
2469  !<
2470  subroutine lak_calculate_external(this, ilak, ex)
2471  ! -- dummy
2472  class(laktype), intent(inout) :: this
2473  integer(I4B), intent(in) :: ilak
2474  real(DP), intent(inout) :: ex
2475  !
2476  ! -- If mover is active, add receiver water to rhs and
2477  ! store available water (as positive value)
2478  ex = dzero
2479  if (this%imover == 1) then
2480  ex = this%pakmvrobj%get_qfrommvr(ilak)
2481  end if
2482  end subroutine lak_calculate_external
2483 
2484  !> @brief Calculate the withdrawal from a lake subject to an available volume
2485  !<
2486  subroutine lak_calculate_withdrawal(this, ilak, avail, wr)
2487  ! -- dummy
2488  class(laktype), intent(inout) :: this
2489  integer(I4B), intent(in) :: ilak
2490  real(DP), intent(inout) :: avail
2491  real(DP), intent(inout) :: wr
2492  !
2493  ! -- withdrawals - limit to sum of inflows and available volume
2494  wr = this%withdrawal(ilak)
2495  if (wr > avail) then
2496  wr = -avail
2497  else
2498  if (wr > dzero) then
2499  wr = -wr
2500  end if
2501  end if
2502  avail = avail + wr
2503  end subroutine lak_calculate_withdrawal
2504 
2505  !> @brief Calculate the evaporation from a lake at a provided stage subject
2506  !! to an available volume
2507  !<
2508  subroutine lak_calculate_evaporation(this, ilak, stage, avail, ev)
2509  ! -- dummy
2510  class(laktype), intent(inout) :: this
2511  integer(I4B), intent(in) :: ilak
2512  real(DP), intent(in) :: stage
2513  real(DP), intent(inout) :: avail
2514  real(DP), intent(inout) :: ev
2515  ! -- local
2516  real(DP) :: sa
2517  !
2518  ! -- evaporation - limit to sum of inflows and available volume
2519  call this%lak_calculate_sarea(ilak, stage, sa)
2520  ev = sa * this%evaporation(ilak)
2521  if (ev > avail) then
2522  if (is_close(avail, dprec)) then
2523  ev = dzero
2524  else
2525  ev = -avail
2526  end if
2527  else
2528  ev = -ev
2529  end if
2530  avail = avail + ev
2531  end subroutine lak_calculate_evaporation
2532 
2533  !> @brief Calculate the outlet inflow to a lake
2534  !<
2535  subroutine lak_calculate_outlet_inflow(this, ilak, outinf)
2536  ! -- dummy
2537  class(laktype), intent(inout) :: this
2538  integer(I4B), intent(in) :: ilak
2539  real(DP), intent(inout) :: outinf
2540  ! -- local
2541  integer(I4B) :: n
2542  !
2543  outinf = dzero
2544  do n = 1, this%noutlets
2545  if (this%lakeout(n) == ilak) then
2546  outinf = outinf - this%simoutrate(n)
2547  if (this%imover == 1) then
2548  outinf = outinf - this%pakmvrobj%get_qtomvr(n)
2549  end if
2550  end if
2551  end do
2552  end subroutine lak_calculate_outlet_inflow
2553 
2554  !> @brief Calculate the outlet outflow from a lake
2555  !<
2556  subroutine lak_calculate_outlet_outflow(this, ilak, stage, avail, outoutf)
2557  ! -- dummy
2558  class(laktype), intent(inout) :: this
2559  integer(I4B), intent(in) :: ilak
2560  real(DP), intent(in) :: stage
2561  real(DP), intent(inout) :: avail
2562  real(DP), intent(inout) :: outoutf
2563  ! -- local
2564  integer(I4B) :: n
2565  real(DP) :: g
2566  real(DP) :: d
2567  real(DP) :: c
2568  real(DP) :: gsm
2569  real(DP) :: rate
2570  !
2571  outoutf = dzero
2572  do n = 1, this%noutlets
2573  if (this%lakein(n) == ilak) then
2574  rate = dzero
2575  d = stage - this%outinvert(n)
2576  if (this%outdmax > dzero) then
2577  if (d > this%outdmax) d = this%outdmax
2578  end if
2579  g = dgravity * this%convlength * this%convtime * this%convtime
2580  select case (this%iouttype(n))
2581  ! specified rate
2582  case (0)
2583  rate = this%outrate(n)
2584  if (-rate > avail) then
2585  rate = -avail
2586  end if
2587  ! manning
2588  case (1)
2589  if (d > dzero) then
2590  c = (this%convlength**donethird) * this%convtime
2591  gsm = dzero
2592  if (this%outrough(n) > dzero) then
2593  gsm = done / this%outrough(n)
2594  end if
2595  rate = -c * gsm * this%outwidth(n) * (d**dfivethirds) * &
2596  sqrt(this%outslope(n))
2597  end if
2598  ! weir
2599  case (2)
2600  if (d > dzero) then
2601  rate = -dtwothirds * dcd * this%outwidth(n) * d * &
2602  sqrt(dtwo * g * d)
2603  end if
2604  end select
2605  this%simoutrate(n) = rate
2606  avail = avail + rate
2607  outoutf = outoutf + rate
2608  end if
2609  end do
2610  end subroutine lak_calculate_outlet_outflow
2611 
2612  !> @brief Get the outlet inflow to a lake from another lake
2613  !<
2614  subroutine lak_get_internal_inlet(this, ilak, outinf)
2615  ! -- dummy
2616  class(laktype), intent(inout) :: this
2617  integer(I4B), intent(in) :: ilak
2618  real(DP), intent(inout) :: outinf
2619  ! -- local
2620  integer(I4B) :: n
2621  !
2622  outinf = dzero
2623  do n = 1, this%noutlets
2624  if (this%lakeout(n) == ilak) then
2625  outinf = outinf - this%simoutrate(n)
2626  if (this%imover == 1) then
2627  outinf = outinf - this%pakmvrobj%get_qtomvr(n)
2628  end if
2629  end if
2630  end do
2631  end subroutine lak_get_internal_inlet
2632 
2633  !> @brief Get the outlet from a lake to another lake
2634  !<
2635  subroutine lak_get_internal_outlet(this, ilak, outoutf)
2636  ! -- dummy
2637  class(laktype), intent(inout) :: this
2638  integer(I4B), intent(in) :: ilak
2639  real(DP), intent(inout) :: outoutf
2640  ! -- local
2641  integer(I4B) :: n
2642  !
2643  outoutf = dzero
2644  do n = 1, this%noutlets
2645  if (this%lakein(n) == ilak) then
2646  if (this%lakeout(n) < 1) cycle
2647  outoutf = outoutf + this%simoutrate(n)
2648  end if
2649  end do
2650  end subroutine lak_get_internal_outlet
2651 
2652  !> @brief Get the outlet outflow from a lake to an external boundary
2653  !<
2654  subroutine lak_get_external_outlet(this, ilak, outoutf)
2655  ! -- dummy
2656  class(laktype), intent(inout) :: this
2657  integer(I4B), intent(in) :: ilak
2658  real(DP), intent(inout) :: outoutf
2659  ! -- local
2660  integer(I4B) :: n
2661  !
2662  outoutf = dzero
2663  do n = 1, this%noutlets
2664  if (this%lakein(n) == ilak) then
2665  if (this%lakeout(n) > 0) cycle
2666  outoutf = outoutf + this%simoutrate(n)
2667  end if
2668  end do
2669  end subroutine lak_get_external_outlet
2670 
2671  !> @brief Get the mover outflow from a lake to an external boundary
2672  !<
2673  subroutine lak_get_external_mover(this, ilak, outoutf)
2674  ! -- dummy
2675  class(laktype), intent(inout) :: this
2676  integer(I4B), intent(in) :: ilak
2677  real(DP), intent(inout) :: outoutf
2678  ! -- local
2679  integer(I4B) :: n
2680  !
2681  outoutf = dzero
2682  if (this%imover == 1) then
2683  do n = 1, this%noutlets
2684  if (this%lakein(n) == ilak) then
2685  if (this%lakeout(n) > 0) cycle
2686  outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
2687  end if
2688  end do
2689  end if
2690  end subroutine lak_get_external_mover
2691 
2692  !> @brief Get the mover outflow from a lake to another lake
2693  !<
2694  subroutine lak_get_internal_mover(this, ilak, outoutf)
2695  ! -- dummy
2696  class(laktype), intent(inout) :: this
2697  integer(I4B), intent(in) :: ilak
2698  real(DP), intent(inout) :: outoutf
2699  ! -- local
2700  integer(I4B) :: n
2701  !
2702  outoutf = dzero
2703  if (this%imover == 1) then
2704  do n = 1, this%noutlets
2705  if (this%lakein(n) == ilak) then
2706  if (this%lakeout(n) < 1) cycle
2707  outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
2708  end if
2709  end do
2710  end if
2711  end subroutine lak_get_internal_mover
2712 
2713  !> @brief Get the outlet to mover from a lake
2714  !<
2715  subroutine lak_get_outlet_tomover(this, ilak, outoutf)
2716  ! -- dummy
2717  class(laktype), intent(inout) :: this
2718  integer(I4B), intent(in) :: ilak
2719  real(DP), intent(inout) :: outoutf
2720  ! -- local
2721  integer(I4B) :: n
2722  !
2723  outoutf = dzero
2724  if (this%imover == 1) then
2725  do n = 1, this%noutlets
2726  if (this%lakein(n) == ilak) then
2727  outoutf = outoutf + this%pakmvrobj%get_qtomvr(n)
2728  end if
2729  end do
2730  end if
2731  end subroutine lak_get_outlet_tomover
2732 
2733  !> @brief Determine the stage from a provided volume
2734  !<
2735  subroutine lak_vol2stage(this, ilak, vol, stage)
2736  ! -- dummy
2737  class(laktype), intent(inout) :: this
2738  integer(I4B), intent(in) :: ilak
2739  real(DP), intent(in) :: vol
2740  real(DP), intent(inout) :: stage
2741  ! -- local
2742  integer(I4B) :: i
2743  integer(I4B) :: ibs
2744  real(DP) :: s0, s1, sm
2745  real(DP) :: v0, v1, vm
2746  real(DP) :: f0, f1, fm
2747  real(DP) :: sa
2748  real(DP) :: en0, en1
2749  real(DP) :: ds, ds0
2750  real(DP) :: denom
2751  !
2752  s0 = this%lakebot(ilak)
2753  call this%lak_calculate_vol(ilak, s0, v0)
2754  s1 = this%laketop(ilak)
2755  call this%lak_calculate_vol(ilak, s1, v1)
2756  ! -- zero volume
2757  if (vol <= v0) then
2758  stage = s0
2759  ! -- linear relation between stage and volume above top of lake
2760  else if (vol >= v1) then
2761  call this%lak_calculate_sarea(ilak, s1, sa)
2762  stage = s1 + (vol - v1) / sa
2763  ! -- use combination of secant and bisection
2764  else
2765  en0 = s0
2766  en1 = s1
2767  ! sm = s1 ! causes divide by zero in 1st line in secantbisection loop
2768  ! sm = s0 ! causes divide by zero in 1st line in secantbisection loop
2769  sm = dzero
2770  f0 = vol - v0
2771  f1 = vol - v1
2772  ibs = 0
2773  secantbisection: do i = 1, 150
2774  denom = f1 - f0
2775  if (denom /= dzero) then
2776  ds = f1 * (s1 - s0) / denom
2777  else
2778  ibs = 13
2779  end if
2780  if (i == 1) then
2781  ds0 = ds
2782  end if
2783  ! -- use bisection if end points are exceeded
2784  if (sm < en0 .or. sm > en1) ibs = 13
2785  ! -- use bisection if secant method stagnates or if
2786  ! ds exceeds previous ds - bisection would occur
2787  ! after conditions exceeded in 13 iterations
2788  if (ds * ds0 < dprec .or. abs(ds) > abs(ds0)) ibs = ibs + 1
2789  if (ibs > 12) then
2790  ds = dhalf * (s1 - s0)
2791  ibs = 0
2792  end if
2793  sm = s1 - ds
2794  if (abs(ds) < dem6) then
2795  exit secantbisection
2796  end if
2797  call this%lak_calculate_vol(ilak, sm, vm)
2798  fm = vol - vm
2799  s0 = s1
2800  f0 = f1
2801  s1 = sm
2802  f1 = fm
2803  ds0 = ds
2804  end do secantbisection
2805  stage = sm
2806  if (abs(ds) >= dem6) then
2807  write (this%iout, '(1x,a,1x,i0,4(1x,a,1x,g15.6))') &
2808  & 'LAK_VOL2STAGE failed for lake', ilak, 'volume error =', fm, &
2809  & 'finding stage (', stage, ') for volume =', vol, &
2810  & 'final change in stage =', ds
2811  end if
2812  end if
2813  end subroutine lak_vol2stage
2814 
2815  !> @brief Determine if a valid lake or outlet number has been specified
2816  function lak_check_valid(this, itemno) result(ierr)
2817  ! -- modules
2818  use simmodule, only: store_error
2819  ! -- return
2820  integer(I4B) :: ierr
2821  ! -- dummy
2822  class(laktype), intent(inout) :: this
2823  integer(I4B), intent(in) :: itemno
2824  ! -- local
2825  integer(I4B) :: ival
2826  !
2827  ierr = 0
2828  ival = abs(itemno)
2829  if (itemno > 0) then
2830  if (ival < 1 .or. ival > this%nlakes) then
2831  write (errmsg, '(a,1x,i0,1x,a,1x,i0,a)') &
2832  'LAKENO', itemno, 'must be greater than 0 and less than or equal to', &
2833  this%nlakes, '.'
2834  call store_error(errmsg)
2835  ierr = 1
2836  end if
2837  else
2838  if (ival < 1 .or. ival > this%noutlets) then
2839  write (errmsg, '(a,1x,i0,1x,a,1x,i0,a)') &
2840  'IOUTLET', itemno, 'must be greater than 0 and less than or equal to', &
2841  this%noutlets, '.'
2842  call store_error(errmsg)
2843  ierr = 1
2844  end if
2845  end if
2846  end function lak_check_valid
2847 
2848  !> @brief Set a stress period attribute for lakweslls(itemno) using keywords
2849  !<
2850  subroutine lak_set_stressperiod(this, itemno)
2851  ! -- modules
2852  use timeseriesmanagermodule, only: read_value_or_time_series_adv
2853  use simmodule, only: store_error
2854  ! -- dummy
2855  class(laktype), intent(inout) :: this
2856  integer(I4B), intent(in) :: itemno
2857  ! -- local
2858  character(len=LINELENGTH) :: text
2859  character(len=LINELENGTH) :: caux
2860  character(len=LINELENGTH) :: keyword
2861  integer(I4B) :: ierr
2862  integer(I4B) :: ii
2863  integer(I4B) :: jj
2864  real(DP), pointer :: bndElem => null()
2865  !
2866  ! -- read line
2867  call this%parser%GetStringCaps(keyword)
2868  select case (keyword)
2869  case ('STATUS')
2870  ierr = this%lak_check_valid(itemno)
2871  if (ierr /= 0) then
2872  goto 999
2873  end if
2874  call this%parser%GetStringCaps(text)
2875  this%status(itemno) = text(1:8)
2876  if (text == 'CONSTANT') then
2877  this%iboundpak(itemno) = -1
2878  else if (text == 'INACTIVE') then
2879  this%iboundpak(itemno) = 0
2880  else if (text == 'ACTIVE') then
2881  this%iboundpak(itemno) = 1
2882  else
2883  write (errmsg, '(a,a)') &
2884  'Unknown '//trim(this%text)//' lak status keyword: ', text//'.'
2885  call store_error(errmsg)
2886  end if
2887  case ('STAGE')
2888  ierr = this%lak_check_valid(itemno)
2889  if (ierr /= 0) then
2890  goto 999
2891  end if
2892  call this%parser%GetString(text)
2893  jj = 1 ! For STAGE
2894  bndelem => this%stage(itemno)
2895  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
2896  this%packName, 'BND', this%tsManager, &
2897  this%iprpak, 'STAGE')
2898  case ('RAINFALL')
2899  ierr = this%lak_check_valid(itemno)
2900  if (ierr /= 0) then
2901  goto 999
2902  end if
2903  call this%parser%GetString(text)
2904  jj = 1 ! For RAINFALL
2905  bndelem => this%rainfall(itemno)
2906  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
2907  this%packName, 'BND', this%tsManager, &
2908  this%iprpak, 'RAINFALL')
2909  if (this%rainfall(itemno) < dzero) then
2910  write (errmsg, '(a,i0,a,G0,a)') &
2911  'Lake ', itemno, ' was assigned a rainfall value of ', &
2912  this%rainfall(itemno), '. Rainfall must be positive.'
2913  call store_error(errmsg)
2914  end if
2915  case ('EVAPORATION')
2916  ierr = this%lak_check_valid(itemno)
2917  if (ierr /= 0) then
2918  goto 999
2919  end if
2920  call this%parser%GetString(text)
2921  jj = 1 ! For EVAPORATION
2922  bndelem => this%evaporation(itemno)
2923  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
2924  this%packName, 'BND', this%tsManager, &
2925  this%iprpak, 'EVAPORATION')
2926  if (this%evaporation(itemno) < dzero) then
2927  write (errmsg, '(a,i0,a,G0,a)') &
2928  'Lake ', itemno, ' was assigned an evaporation value of ', &
2929  this%evaporation(itemno), '. Evaporation must be positive.'
2930  call store_error(errmsg)
2931  end if
2932  case ('RUNOFF')
2933  ierr = this%lak_check_valid(itemno)
2934  if (ierr /= 0) then
2935  goto 999
2936  end if
2937  call this%parser%GetString(text)
2938  jj = 1 ! For RUNOFF
2939  bndelem => this%runoff(itemno)
2940  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
2941  this%packName, 'BND', this%tsManager, &
2942  this%iprpak, 'RUNOFF')
2943  if (this%runoff(itemno) < dzero) then
2944  write (errmsg, '(a,i0,a,G0,a)') &
2945  'Lake ', itemno, ' was assigned a runoff value of ', &
2946  this%runoff(itemno), '. Runoff must be positive.'
2947  call store_error(errmsg)
2948  end if
2949  case ('INFLOW')
2950  ierr = this%lak_check_valid(itemno)
2951  if (ierr /= 0) then
2952  goto 999
2953  end if
2954  call this%parser%GetString(text)
2955  jj = 1 ! For specified INFLOW
2956  bndelem => this%inflow(itemno)
2957  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
2958  this%packName, 'BND', this%tsManager, &
2959  this%iprpak, 'INFLOW')
2960  if (this%inflow(itemno) < dzero) then
2961  write (errmsg, '(a,i0,a,G0,a)') &
2962  'Lake ', itemno, ' was assigned an inflow value of ', &
2963  this%inflow(itemno), '. Inflow must be positive.'
2964  call store_error(errmsg)
2965  end if
2966  case ('WITHDRAWAL')
2967  ierr = this%lak_check_valid(itemno)
2968  if (ierr /= 0) then
2969  goto 999
2970  end if
2971  call this%parser%GetString(text)
2972  jj = 1 ! For specified WITHDRAWAL
2973  bndelem => this%withdrawal(itemno)
2974  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
2975  this%packName, 'BND', this%tsManager, &
2976  this%iprpak, 'WITHDRAWAL')
2977  if (this%withdrawal(itemno) < dzero) then
2978  write (errmsg, '(a,i0,a,G0,a)') &
2979  'Lake ', itemno, ' was assigned a withdrawal value of ', &
2980  this%withdrawal(itemno), '. Withdrawal must be positive.'
2981  call store_error(errmsg)
2982  end if
2983  case ('RATE')
2984  ierr = this%lak_check_valid(-itemno)
2985  if (ierr /= 0) then
2986  goto 999
2987  end if
2988  call this%parser%GetString(text)
2989  jj = 1 ! For specified OUTLET RATE
2990  bndelem => this%outrate(itemno)
2991  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
2992  this%packName, 'BND', this%tsManager, &
2993  this%iprpak, 'RATE')
2994  case ('INVERT')
2995  ierr = this%lak_check_valid(-itemno)
2996  if (ierr /= 0) then
2997  goto 999
2998  end if
2999  call this%parser%GetString(text)
3000  jj = 1 ! For OUTLET INVERT
3001  bndelem => this%outinvert(itemno)
3002  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
3003  this%packName, 'BND', this%tsManager, &
3004  this%iprpak, 'INVERT')
3005  case ('WIDTH')
3006  ierr = this%lak_check_valid(-itemno)
3007  if (ierr /= 0) then
3008  goto 999
3009  end if
3010  call this%parser%GetString(text)
3011  jj = 1 ! For OUTLET WIDTH
3012  bndelem => this%outwidth(itemno)
3013  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
3014  this%packName, 'BND', this%tsManager, &
3015  this%iprpak, 'WIDTH')
3016  case ('ROUGH')
3017  ierr = this%lak_check_valid(-itemno)
3018  if (ierr /= 0) then
3019  goto 999
3020  end if
3021  call this%parser%GetString(text)
3022  jj = 1 ! For OUTLET ROUGHNESS
3023  bndelem => this%outrough(itemno)
3024  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
3025  this%packName, 'BND', this%tsManager, &
3026  this%iprpak, 'ROUGH')
3027  case ('SLOPE')
3028  ierr = this%lak_check_valid(-itemno)
3029  if (ierr /= 0) then
3030  goto 999
3031  end if
3032  call this%parser%GetString(text)
3033  jj = 1 ! For OUTLET SLOPE
3034  bndelem => this%outslope(itemno)
3035  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
3036  this%packName, 'BND', this%tsManager, &
3037  this%iprpak, 'SLOPE')
3038  case ('AUXILIARY')
3039  ierr = this%lak_check_valid(itemno)
3040  if (ierr /= 0) then
3041  goto 999
3042  end if
3043  call this%parser%GetStringCaps(caux)
3044  do jj = 1, this%naux
3045  if (trim(adjustl(caux)) /= trim(adjustl(this%auxname(jj)))) cycle
3046  call this%parser%GetString(text)
3047  ii = itemno
3048  bndelem => this%lauxvar(jj, ii)
3049  call read_value_or_time_series_adv(text, itemno, jj, bndelem, &
3050  this%packName, 'AUX', &
3051  this%tsManager, this%iprpak, &
3052  this%auxname(jj))
3053  exit
3054  end do
3055  case default
3056  write (errmsg, '(2a)') &
3057  'Unknown '//trim(this%text)//' lak data keyword: ', &
3058  trim(keyword)//'.'
3059  end select
3060  !
3061  ! -- Return
3062 999 return
3063  end subroutine lak_set_stressperiod
3064 
3065  !> @brief Issue a parameter error for lakweslls(ilak)
3066  !!
3067  !! Read itmp and new boundaries if itmp > 0
3068  !<
3069  subroutine lak_set_attribute_error(this, ilak, keyword, msg)
3070  ! -- modules
3071  use simmodule, only: store_error
3072  ! -- dummy
3073  class(laktype), intent(inout) :: this
3074  integer(I4B), intent(in) :: ilak
3075  character(len=*), intent(in) :: keyword
3076  character(len=*), intent(in) :: msg
3077  !
3078  if (len(msg) == 0) then
3079  write (errmsg, '(a,1x,a,1x,i0,1x,a)') &
3080  keyword, ' for LAKE', ilak, 'has already been set.'
3081  else
3082  write (errmsg, '(a,1x,a,1x,i0,1x,a)') keyword, ' for LAKE', ilak, msg
3083  end if
3084  call store_error(errmsg)
3085  end subroutine lak_set_attribute_error
3086 
3087  !> @brief Set options specific to LakType
3088  !!
3089  !! lak_options overrides BndType%bnd_options
3090  !<
3091  subroutine lak_options(this, option, found)
3092  ! -- modules
3093  use constantsmodule, only: maxcharlen, dzero, mnormal
3094  use openspecmodule, only: access, form
3095  use simmodule, only: store_error
3096  use inputoutputmodule, only: urword, getunit, openfile
3097  ! -- dummy
3098  class(laktype), intent(inout) :: this
3099  character(len=*), intent(inout) :: option
3100  logical(LGP), intent(inout) :: found
3101  ! -- local
3102  character(len=MAXCHARLEN) :: fname, keyword
3103  real(DP) :: r
3104  ! -- formats
3105  character(len=*), parameter :: fmtlengthconv = &
3106  &"(4x, 'LENGTH CONVERSION VALUE (',g15.7,') SPECIFIED.')"
3107  character(len=*), parameter :: fmttimeconv = &
3108  &"(4x, 'TIME CONVERSION VALUE (',g15.7,') SPECIFIED.')"
3109  character(len=*), parameter :: fmtoutdmax = &
3110  &"(4x, 'MAXIMUM OUTLET WATER DEPTH (',g15.7,') SPECIFIED.')"
3111  character(len=*), parameter :: fmtlakeopt = &
3112  &"(4x, 'LAKE ', a, ' VALUE (',g15.7,') SPECIFIED.')"
3113  character(len=*), parameter :: fmtlakbin = &
3114  "(4x, 'LAK ', 1x, a, 1x, ' WILL BE SAVED TO FILE: ', &
3115  &a, /4x, 'OPENED ON UNIT: ', I0)"
3116  character(len=*), parameter :: fmtiter = &
3117  &"(4x, 'MAXIMUM LAK ITERATION VALUE (',i0,') SPECIFIED.')"
3118  character(len=*), parameter :: fmtdmaxchg = &
3119  &"(4x, 'MAXIMUM STAGE CHANGE VALUE (',g0,') SPECIFIED.')"
3120  !
3121  found = .true.
3122  select case (option)
3123  case ('PRINT_STAGE')
3124  this%iprhed = 1
3125  write (this%iout, '(4x,a)') trim(adjustl(this%text))// &
3126  ' STAGES WILL BE PRINTED TO LISTING FILE.'
3127  case ('STAGE')
3128  call this%parser%GetStringCaps(keyword)
3129  if (keyword == 'FILEOUT') then
3130  call this%parser%GetString(fname)
3131  this%istageout = getunit()
3132  call openfile(this%istageout, this%iout, fname, 'DATA(BINARY)', &
3133  form, access, 'REPLACE', mode_opt=mnormal)
3134  write (this%iout, fmtlakbin) 'STAGE', trim(adjustl(fname)), &
3135  this%istageout
3136  else
3137  call store_error('OPTIONAL STAGE KEYWORD MUST BE FOLLOWED BY FILEOUT')
3138  end if
3139  case ('BUDGET')
3140  call this%parser%GetStringCaps(keyword)
3141  if (keyword == 'FILEOUT') then
3142  call this%parser%GetString(fname)
3143  this%ibudgetout = getunit()
3144  call openfile(this%ibudgetout, this%iout, fname, 'DATA(BINARY)', &
3145  form, access, 'REPLACE', mode_opt=mnormal)
3146  write (this%iout, fmtlakbin) 'BUDGET', trim(adjustl(fname)), &
3147  this%ibudgetout
3148  else
3149  call store_error('OPTIONAL BUDGET KEYWORD MUST BE FOLLOWED BY FILEOUT')
3150  end if
3151  case ('BUDGETCSV')
3152  call this%parser%GetStringCaps(keyword)
3153  if (keyword == 'FILEOUT') then
3154  call this%parser%GetString(fname)
3155  this%ibudcsv = getunit()
3156  call openfile(this%ibudcsv, this%iout, fname, 'CSV', &
3157  filstat_opt='REPLACE')
3158  write (this%iout, fmtlakbin) 'BUDGET CSV', trim(adjustl(fname)), &
3159  this%ibudcsv
3160  else
3161  call store_error('OPTIONAL BUDGETCSV KEYWORD MUST BE FOLLOWED BY &
3162  &FILEOUT')
3163  end if
3164  case ('PACKAGE_CONVERGENCE')
3165  call this%parser%GetStringCaps(keyword)
3166  if (keyword == 'FILEOUT') then
3167  call this%parser%GetString(fname)
3168  this%ipakcsv = getunit()
3169  call openfile(this%ipakcsv, this%iout, fname, 'CSV', &
3170  filstat_opt='REPLACE', mode_opt=mnormal)
3171  write (this%iout, fmtlakbin) 'PACKAGE_CONVERGENCE', &
3172  trim(adjustl(fname)), this%ipakcsv
3173  else
3174  call store_error('OPTIONAL PACKAGE_CONVERGENCE KEYWORD MUST BE '// &
3175  'FOLLOWED BY FILEOUT')
3176  end if
3177  case ('MOVER')
3178  this%imover = 1
3179  write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
3180  case ('LENGTH_CONVERSION')
3181  this%convlength = this%parser%GetDouble()
3182  write (this%iout, fmtlengthconv) this%convlength
3183  case ('TIME_CONVERSION')
3184  this%convtime = this%parser%GetDouble()
3185  write (this%iout, fmttimeconv) this%convtime
3186  case ('SURFDEP')
3187  r = this%parser%GetDouble()
3188  if (r < dzero) then
3189  r = dzero
3190  end if
3191  this%surfdep = r
3192  write (this%iout, fmtlakeopt) 'SURFDEP', this%surfdep
3193  case ('MAXIMUM_ITERATIONS')
3194  this%maxlakit = this%parser%GetInteger()
3195  write (this%iout, fmtiter) this%maxlakit
3196  case ('MAXIMUM_STAGE_CHANGE')
3197  r = this%parser%GetDouble()
3198  this%dmaxchg = r
3199  this%delh = dp999 * r
3200  write (this%iout, fmtdmaxchg) this%dmaxchg
3201  !
3202  ! -- right now these are options that are only available in the
3203  ! development version and are not included in the documentation.
3204  ! These options are only available when IDEVELOPMODE in
3205  ! constants module is set to 1
3206  case ('DEV_GROUNDWATER_HEAD_CONDUCTANCE')
3207  call this%parser%DevOpt()
3208  this%igwhcopt = 1
3209  write (this%iout, '(4x,a)') &
3210  'CONDUCTANCE FOR HORIZONTAL CONNECTIONS WILL BE CALCULATED &
3211  &USING THE GROUNDWATER HEAD'
3212  case ('DEV_MAXIMUM_OUTLET_DEPTH')
3213  call this%parser%DevOpt()
3214  this%outdmax = this%parser%GetDouble()
3215  write (this%iout, fmtoutdmax) this%outdmax
3216  case ('DEV_NO_FINAL_CHECK')
3217  call this%parser%DevOpt()
3218  this%iconvchk = 0
3219  write (this%iout, '(4x,a)') &
3220  'A FINAL CONVERGENCE CHECK OF THE CHANGE IN LAKE STAGES &
3221  &WILL NOT BE MADE'
3222  case ('DEV_NO_FINAL_RESIDUAL_CHECK')
3223  call this%parser%DevOpt()
3224  this%iconvresidchk = 0
3225  write (this%iout, '(4x,a)') &
3226  'A FINAL CONVERGENCE CHECK OF THE CHANGE IN LAKE RESIDUALS &
3227  &WILL NOT BE MADE'
3228  case ('DEV_MAXIMUM_PERCENT_DIFFERENCE')
3229  call this%parser%DevOpt()
3230  r = this%parser%GetDouble()
3231  if (r < dzero) then
3232  r = dem1
3233  end if
3234  this%pdmax = r
3235  write (this%iout, fmtlakeopt) 'MAXIMUM_PERCENT_DIFFERENCE', this%pdmax
3236  case default
3237  !
3238  ! -- No options found
3239  found = .false.
3240  end select
3241  end subroutine lak_options
3242 
3243  !> @brief Allocate and Read
3244  !!
3245  !! Create new LAK package and point bndobj to the new package
3246  !<
3247  subroutine lak_ar(this)
3248  ! -- dummy
3249  class(laktype), intent(inout) :: this
3250  !
3251  call this%obs%obs_ar()
3252  !
3253  ! -- Allocate arrays in LAK and in package superclass
3254  call this%lak_allocate_arrays()
3255  !
3256  ! -- read optional initial package parameters
3257  call this%read_initial_attr()
3258  !
3259  ! -- setup pakmvrobj
3260  if (this%imover /= 0) then
3261  allocate (this%pakmvrobj)
3262  call this%pakmvrobj%ar(this%noutlets, this%nlakes, this%memoryPath)
3263  end if
3264  end subroutine lak_ar
3265 
3266  !> @brief Read and Prepare
3267  !!
3268  !! Read itmp and read new boundaries if itmp > 0
3269  !<
3270  subroutine lak_rp(this)
3271  ! -- modules
3272  use constantsmodule, only: linelength
3273  use tdismodule, only: kper, nper
3274  use simmodule, only: store_error, count_errors
3275  ! -- dummy
3276  class(laktype), intent(inout) :: this
3277  ! -- local
3278  character(len=LINELENGTH) :: title
3279  character(len=LINELENGTH) :: line
3280  character(len=LINELENGTH) :: text
3281  logical(LGP) :: isfound
3282  logical(LGP) :: endOfBlock
3283  integer(I4B) :: ierr
3284  integer(I4B) :: node
3285  integer(I4B) :: n
3286  integer(I4B) :: itemno
3287  integer(I4B) :: j
3288  ! -- formats
3289  character(len=*), parameter :: fmtblkerr = &
3290  &"('Looking for BEGIN PERIOD iper. Found ', a, ' instead.')"
3291  character(len=*), parameter :: fmtlsp = &
3292  &"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
3293  !
3294  ! -- set nbound to maxbound
3295  this%nbound = this%maxbound
3296  !
3297  ! -- Set ionper to the stress period number for which a new block of data
3298  ! will be read.
3299  if (this%inunit == 0) return
3300  !
3301  ! -- get stress period data
3302  if (this%ionper < kper) then
3303  !
3304  ! -- get period block
3305  call this%parser%GetBlock('PERIOD', isfound, ierr, &
3306  supportopenclose=.true., &
3307  blockrequired=.false.)
3308  if (isfound) then
3309  !
3310  ! -- read ionper and check for increasing period numbers
3311  call this%read_check_ionper()
3312  else
3313  !
3314  ! -- PERIOD block not found
3315  if (ierr < 0) then
3316  ! -- End of file found; data applies for remainder of simulation.
3317  this%ionper = nper + 1
3318  else
3319  ! -- Found invalid block
3320  call this%parser%GetCurrentLine(line)
3321  write (errmsg, fmtblkerr) adjustl(trim(line))
3322  call store_error(errmsg)
3323  call this%parser%StoreErrorUnit()
3324  end if
3325  end if
3326  end if
3327  !
3328  ! -- Read data if ionper == kper
3329  if (this%ionper == kper) then
3330  !
3331  ! -- setup table for period data
3332  if (this%iprpak /= 0) then
3333  !
3334  ! -- reset the input table object
3335  title = trim(adjustl(this%text))//' PACKAGE ('// &
3336  trim(adjustl(this%packName))//') DATA FOR PERIOD'
3337  write (title, '(a,1x,i6)') trim(adjustl(title)), kper
3338  call table_cr(this%inputtab, this%packName, title)
3339  call this%inputtab%table_df(1, 4, this%iout, finalize=.false.)
3340  text = 'NUMBER'
3341  call this%inputtab%initialize_column(text, 10, alignment=tabcenter)
3342  text = 'KEYWORD'
3343  call this%inputtab%initialize_column(text, 20, alignment=tableft)
3344  do n = 1, 2
3345  write (text, '(a,1x,i6)') 'VALUE', n
3346  call this%inputtab%initialize_column(text, 15, alignment=tabcenter)
3347  end do
3348  end if
3349  !
3350  ! -- read the data
3351  this%check_attr = 1
3352  stressperiod: do
3353  call this%parser%GetNextLine(endofblock)
3354  if (endofblock) exit
3355  !
3356  ! -- get lake or outlet number
3357  itemno = this%parser%GetInteger()
3358  !
3359  ! -- read data from the rest of the line
3360  call this%lak_set_stressperiod(itemno)
3361  !
3362  ! -- write line to table
3363  if (this%iprpak /= 0) then
3364  call this%parser%GetCurrentLine(line)
3365  call this%inputtab%line_to_columns(line)
3366  end if
3367  end do stressperiod
3368  !
3369  if (this%iprpak /= 0) then
3370  call this%inputtab%finalize_table()
3371  end if
3372  !
3373  ! -- using stress period data from the previous stress period
3374  else
3375  write (this%iout, fmtlsp) trim(this%filtyp)
3376  end if
3377  !
3378  ! -- write summary of lake stress period error messages
3379  if (count_errors() > 0) then
3380  call this%parser%StoreErrorUnit()
3381  end if
3382  !
3383  ! -- fill bound array with lake stage, conductance, and bottom elevation
3384  do n = 1, this%nlakes
3385  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3386  node = this%cellid(j)
3387  this%nodelist(j) = node
3388  this%bound(1, j) = this%xnewpak(n)
3389  this%bound(2, j) = this%satcond(j)
3390  this%bound(3, j) = this%belev(j)
3391  end do
3392  end do
3393  !
3394  ! -- copy lakein into iprmap so mvr budget contains lake instead of outlet
3395  if (this%imover == 1) then
3396  do n = 1, this%noutlets
3397  this%pakmvrobj%iprmap(n) = this%lakein(n)
3398  end do
3399  end if
3400  end subroutine lak_rp
3401 
3402  !> @brief Add package connection to matrix
3403  !<
3404  subroutine lak_ad(this)
3405  ! -- modules
3406  use simvariablesmodule, only: ifailedstepretry
3407  ! -- dummy
3408  class(laktype) :: this
3409  ! -- local
3410  integer(I4B) :: n
3411  integer(I4B) :: j
3412  integer(I4B) :: iaux
3413  !
3414  ! -- Advance the time series
3415  call this%TsManager%ad()
3416  !
3417  ! -- update auxiliary variables by copying from the derived-type time
3418  ! series variable into the bndpackage auxvar variable so that this
3419  ! information is properly written to the GWF budget file
3420  if (this%naux > 0) then
3421  do n = 1, this%nlakes
3422  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3423  do iaux = 1, this%naux
3424  if (this%noupdateauxvar(iaux) /= 0) cycle
3425  this%auxvar(iaux, j) = this%lauxvar(iaux, n)
3426  end do
3427  end do
3428  end do
3429  end if
3430  !
3431  ! -- Update or restore state
3432  if (ifailedstepretry == 0) then
3433  !
3434  ! -- copy xnew into xold and set xnewpak to stage for
3435  ! constant stage lakes
3436  do n = 1, this%nlakes
3437  this%xoldpak(n) = this%xnewpak(n)
3438  this%stageiter(n) = this%xnewpak(n)
3439  if (this%iboundpak(n) < 0) then
3440  this%xnewpak(n) = this%stage(n)
3441  end if
3442  this%seep0(n) = dzero
3443  end do
3444  else
3445  !
3446  ! -- copy xold back into xnew as this is a
3447  ! retry of this time step
3448  do n = 1, this%nlakes
3449  this%xnewpak(n) = this%xoldpak(n)
3450  this%stageiter(n) = this%xnewpak(n)
3451  if (this%iboundpak(n) < 0) then
3452  this%xnewpak(n) = this%stage(n)
3453  end if
3454  this%seep0(n) = dzero
3455  end do
3456  end if
3457  !
3458  ! -- pakmvrobj ad
3459  if (this%imover == 1) then
3460  call this%pakmvrobj%ad()
3461  end if
3462  !
3463  ! -- For each observation, push simulated value and corresponding
3464  ! simulation time from "current" to "preceding" and reset
3465  ! "current" value.
3466  call this%obs%obs_ad()
3467  end subroutine lak_ad
3468 
3469  !> @brief Formulate the HCOF and RHS terms
3470  !!
3471  !! Skip if no lakes, otherwise calculate hcof and rhs
3472  !<
3473  subroutine lak_cf(this)
3474  ! -- dummy
3475  class(laktype) :: this
3476  ! -- local
3477  integer(I4B) :: j, n
3478  integer(I4B) :: igwfnode
3479  real(DP) :: hlak, bottom_lake
3480  !
3481  ! -- save groundwater seepage for lake solution
3482  do n = 1, this%nlakes
3483  this%seep0(n) = this%seep(n)
3484  end do
3485  !
3486  ! -- save variables for convergence check
3487  do n = 1, this%nlakes
3488  this%s0(n) = this%xnewpak(n)
3489  call this%lak_calculate_exchange(n, this%s0(n), this%qgwf0(n))
3490  end do
3491  !
3492  ! -- find highest active cell
3493  do n = 1, this%nlakes
3494  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3495  ! -- skip horizontal connections
3496  if (this%ictype(j) /= 0) then
3497  cycle
3498  end if
3499  igwfnode = this%nodesontop(j)
3500  if (this%ibound(igwfnode) == 0) then
3501  call this%dis%highest_active(igwfnode, this%ibound)
3502  end if
3503  this%nodelist(j) = igwfnode
3504  this%cellid(j) = igwfnode
3505  end do
3506  end do
3507  !
3508  ! -- reset ibound for cells where lake stage is above the bottom
3509  ! of the lake in the cell or the lake is inactive - only applied to
3510  ! vertical connections
3511  do n = 1, this%nlakes
3512  !
3513  hlak = this%xnewpak(n)
3514  !
3515  ! -- Go through lake connections
3516  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3517  !
3518  ! -- assign gwf node number
3519  igwfnode = this%cellid(j)
3520  !
3521  ! -- skip inactive or constant head GWF cells
3522  if (this%ibound(igwfnode) < 1) then
3523  cycle
3524  end if
3525  !
3526  ! -- skip horizontal connections
3527  if (this%ictype(j) /= 0) then
3528  cycle
3529  end if
3530  !
3531  ! -- skip embedded lakes
3532  if (this%ictype(j) == 2 .or. this%ictype(j) == 3) then
3533  cycle
3534  end if
3535  !
3536  ! -- Mark ibound for wet lakes or inactive lakes; reset to 1 otherwise
3537  bottom_lake = this%belev(j)
3538  if (hlak > bottom_lake .or. this%iboundpak(n) == 0) then
3539  this%ibound(igwfnode) = iwetlake
3540  else
3541  this%ibound(igwfnode) = 1
3542  end if
3543  end do
3544  !
3545  end do
3546  !
3547  ! -- Store the lake stage and cond in bound array for other
3548  ! packages, such as the BUY package
3549  call this%lak_bound_update()
3550  end subroutine lak_cf
3551 
3552  !> @brief Copy rhs and hcof into solution rhs and amat
3553  !<
3554  subroutine lak_fc(this, rhs, ia, idxglo, matrix_sln)
3555  ! -- dummy
3556  class(laktype) :: this
3557  real(DP), dimension(:), intent(inout) :: rhs
3558  integer(I4B), dimension(:), intent(in) :: ia
3559  integer(I4B), dimension(:), intent(in) :: idxglo
3560  class(matrixbasetype), pointer :: matrix_sln
3561  ! -- local
3562  integer(I4B) :: j, n
3563  integer(I4B) :: igwfnode
3564  integer(I4B) :: ipossymd
3565  !
3566  ! -- pakmvrobj fc
3567  if (this%imover == 1) then
3568  call this%pakmvrobj%fc()
3569  end if
3570  !
3571  ! -- make a stab at a solution
3572  call this%lak_solve()
3573  !
3574  ! -- add terms to the gwf matrix
3575  do n = 1, this%nlakes
3576  if (this%iboundpak(n) == 0) cycle
3577  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3578  igwfnode = this%cellid(j)
3579  if (this%ibound(igwfnode) < 1) cycle
3580  ipossymd = idxglo(ia(igwfnode))
3581  call matrix_sln%add_value_pos(ipossymd, this%hcof(j))
3582  rhs(igwfnode) = rhs(igwfnode) + this%rhs(j)
3583  end do
3584  end do
3585  end subroutine lak_fc
3586 
3587  !> @brief Fill newton terms
3588  !<
3589  subroutine lak_fn(this, rhs, ia, idxglo, matrix_sln)
3590  ! -- dummy
3591  class(laktype) :: this
3592  real(DP), dimension(:), intent(inout) :: rhs
3593  integer(I4B), dimension(:), intent(in) :: ia
3594  integer(I4B), dimension(:), intent(in) :: idxglo
3595  class(matrixbasetype), pointer :: matrix_sln
3596  ! -- local
3597  integer(I4B) :: j, n
3598  integer(I4B) :: ipos
3599  integer(I4B) :: igwfnode
3600  integer(I4B) :: idry
3601  real(DP) :: hlak
3602  real(DP) :: avail
3603  real(DP) :: ra
3604  real(DP) :: ro
3605  real(DP) :: qinf
3606  real(DP) :: ex
3607  real(DP) :: head
3608  real(DP) :: q
3609  real(DP) :: q1
3610  real(DP) :: rterm
3611  real(DP) :: drterm
3612  !
3613  do n = 1, this%nlakes
3614  if (this%iboundpak(n) == 0) cycle
3615  hlak = this%xnewpak(n)
3616  call this%lak_calculate_available(n, hlak, avail, &
3617  ra, ro, qinf, ex, this%delh)
3618  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
3619  igwfnode = this%cellid(j)
3620  ipos = ia(igwfnode)
3621  head = this%xnew(igwfnode)
3622  if (-this%hcof(j) > dzero) then
3623  if (this%ibound(igwfnode) > 0) then
3624  ! -- estimate lake-aquifer exchange with perturbed groundwater head
3625  ! exchange is relative to the lake
3626  !avail = DEP20
3627  call this%lak_estimate_conn_exchange(2, n, j, idry, hlak, &
3628  head + this%delh, q1, avail)
3629  q1 = -q1
3630  ! -- calculate unperturbed lake-aquifer exchange
3631  q = this%hcof(j) * head - this%rhs(j)
3632  ! -- calculate rterm
3633  rterm = this%hcof(j) * head
3634  ! -- calculate derivative
3635  drterm = (q1 - q) / this%delh
3636  ! -- add terms to convert conductance formulation into
3637  ! newton-raphson formulation
3638  call matrix_sln%add_value_pos(idxglo(ipos), drterm - this%hcof(j))
3639  rhs(igwfnode) = rhs(igwfnode) - rterm + drterm * head
3640  end if
3641  end if
3642  end do
3643  end do
3644  end subroutine lak_fn
3645 
3646  !> @brief Final convergence check for package
3647  !<
3648  subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
3649  ! -- modules
3650  use tdismodule, only: totim, kstp, kper, delt
3651  ! -- dummy
3652  class(laktype), intent(inout) :: this
3653  integer(I4B), intent(in) :: innertot
3654  integer(I4B), intent(in) :: kiter
3655  integer(I4B), intent(in) :: iend
3656  integer(I4B), intent(in) :: icnvgmod
3657  character(len=LENPAKLOC), intent(inout) :: cpak
3658  integer(I4B), intent(inout) :: ipak
3659  real(DP), intent(inout) :: dpak
3660  ! -- local
3661  character(len=LENPAKLOC) :: cloc
3662  character(len=LINELENGTH) :: tag
3663  integer(I4B) :: icheck
3664  integer(I4B) :: ipakfail
3665  integer(I4B) :: locdhmax
3666  integer(I4B) :: locresidmax
3667  integer(I4B) :: locdgwfmax
3668  integer(I4B) :: locdqoutmax
3669  integer(I4B) :: locdqfrommvrmax
3670  integer(I4B) :: ntabrows
3671  integer(I4B) :: ntabcols
3672  integer(I4B) :: n
3673  real(DP) :: q
3674  real(DP) :: q0
3675  real(DP) :: qtolfact
3676  real(DP) :: area
3677  real(DP) :: gwf0
3678  real(DP) :: gwf
3679  real(DP) :: dh
3680  real(DP) :: resid
3681  real(DP) :: dgwf
3682  real(DP) :: hlak0
3683  real(DP) :: hlak
3684  real(DP) :: qout0
3685  real(DP) :: qout
3686  real(DP) :: dqout
3687  real(DP) :: inf
3688  real(DP) :: ra
3689  real(DP) :: ro
3690  real(DP) :: qinf
3691  real(DP) :: ex
3692  real(DP) :: dhmax
3693  real(DP) :: residmax
3694  real(DP) :: dgwfmax
3695  real(DP) :: dqoutmax
3696  real(DP) :: dqfrommvr
3697  real(DP) :: dqfrommvrmax
3698  !
3699  ! -- initialize local variables
3700  icheck = this%iconvchk
3701  ipakfail = 0
3702  locdhmax = 0
3703  locresidmax = 0
3704  locdgwfmax = 0
3705  locdqoutmax = 0
3706  locdqfrommvrmax = 0
3707  dhmax = dzero
3708  residmax = dzero
3709  dgwfmax = dzero
3710  dqoutmax = dzero
3711  dqfrommvrmax = dzero
3712  !
3713  ! -- if not saving package convergence data on check convergence if
3714  ! the model is considered converged
3715  if (this%ipakcsv == 0) then
3716  if (icnvgmod == 0) then
3717  icheck = 0
3718  end if
3719  !
3720  ! -- saving package convergence data
3721  else
3722  !
3723  ! -- header for package csv
3724  if (.not. associated(this%pakcsvtab)) then
3725  !
3726  ! -- determine the number of columns and rows
3727  ntabrows = 1
3728  ntabcols = 11
3729  if (this%noutlets > 0) then
3730  ntabcols = ntabcols + 2
3731  end if
3732  if (this%imover == 1) then
3733  ntabcols = ntabcols + 2
3734  end if
3735  !
3736  ! -- setup table
3737  call table_cr(this%pakcsvtab, this%packName, '')
3738  call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, &
3739  lineseparator=.false., separator=',', &
3740  finalize=.false.)
3741  !
3742  ! -- add columns to package csv
3743  tag = 'total_inner_iterations'
3744  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
3745  tag = 'totim'
3746  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
3747  tag = 'kper'
3748  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
3749  tag = 'kstp'
3750  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
3751  tag = 'nouter'
3752  call this%pakcsvtab%initialize_column(tag, 10, alignment=tableft)
3753  tag = 'dvmax'
3754  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
3755  tag = 'dvmax_loc'
3756  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
3757  tag = 'residmax'
3758  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
3759  tag = 'residmax_loc'
3760  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
3761  tag = 'dgwfmax'
3762  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
3763  tag = 'dgwfmax_loc'
3764  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
3765  if (this%noutlets > 0) then
3766  tag = 'dqoutmax'
3767  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
3768  tag = 'dqoutmax_loc'
3769  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
3770  end if
3771  if (this%imover == 1) then
3772  tag = 'dqfrommvrmax'
3773  call this%pakcsvtab%initialize_column(tag, 15, alignment=tableft)
3774  tag = 'dqfrommvrmax_loc'
3775  call this%pakcsvtab%initialize_column(tag, 16, alignment=tableft)
3776  end if
3777  end if
3778  end if
3779  !
3780  ! -- perform package convergence check
3781  if (icheck /= 0) then
3782  final_check: do n = 1, this%nlakes
3783  if (this%iboundpak(n) < 1) cycle
3784  !
3785  ! -- set previous and current lake stage
3786  hlak0 = this%s0(n)
3787  hlak = this%xnewpak(n)
3788  !
3789  ! -- stage difference
3790  dh = hlak0 - hlak
3791  !
3792  ! -- calculate surface area
3793  call this%lak_calculate_sarea(n, hlak, area)
3794  !
3795  ! -- set the Q to length factor
3796  if (area > dzero) then
3797  qtolfact = delt / area
3798  else
3799  qtolfact = dzero
3800  end if
3801  !
3802  ! -- difference in the residual
3803  call this%lak_calculate_residual(n, hlak, resid)
3804  resid = resid * qtolfact
3805  !
3806  ! -- change in gwf exchange
3807  dgwf = dzero
3808  if (area > dzero) then
3809  gwf0 = this%qgwf0(n)
3810  call this%lak_calculate_exchange(n, hlak, gwf)
3811  dgwf = (gwf0 - gwf) * qtolfact
3812  end if
3813  !
3814  ! -- change in outflows
3815  dqout = dzero
3816  if (this%noutlets > 0) then
3817  if (area > dzero) then
3818  call this%lak_calculate_available(n, hlak0, inf, ra, ro, qinf, ex)
3819  call this%lak_calculate_outlet_outflow(n, hlak0, inf, qout0)
3820  call this%lak_calculate_available(n, hlak, inf, ra, ro, qinf, ex)
3821  call this%lak_calculate_outlet_outflow(n, hlak, inf, qout)
3822  dqout = (qout0 - qout) * qtolfact
3823  end if
3824  end if
3825  !
3826  ! -- q from mvr
3827  dqfrommvr = dzero
3828  if (this%imover == 1) then
3829  q = this%pakmvrobj%get_qfrommvr(n)
3830  q0 = this%pakmvrobj%get_qfrommvr0(n)
3831  dqfrommvr = qtolfact * (q0 - q)
3832  end if
3833  !
3834  ! -- evaluate magnitude of differences
3835  if (n == 1) then
3836  locdhmax = n
3837  dhmax = dh
3838  locdgwfmax = n
3839  residmax = resid
3840  locresidmax = n
3841  dgwfmax = dgwf
3842  locdqoutmax = n
3843  dqoutmax = dqout
3844  dqfrommvrmax = dqfrommvr
3845  locdqfrommvrmax = n
3846  else
3847  if (abs(dh) > abs(dhmax)) then
3848  locdhmax = n
3849  dhmax = dh
3850  end if
3851  if (abs(resid) > abs(residmax)) then
3852  locresidmax = n
3853  residmax = resid
3854  end if
3855  if (abs(dgwf) > abs(dgwfmax)) then
3856  locdgwfmax = n
3857  dgwfmax = dgwf
3858  end if
3859  if (abs(dqout) > abs(dqoutmax)) then
3860  locdqoutmax = n
3861  dqoutmax = dqout
3862  end if
3863  if (abs(dqfrommvr) > abs(dqfrommvrmax)) then
3864  dqfrommvrmax = dqfrommvr
3865  locdqfrommvrmax = n
3866  end if
3867  end if
3868  end do final_check
3869  !
3870  ! -- set dpak and cpak
3871  if (abs(dhmax) > abs(dpak)) then
3872  ipak = locdhmax
3873  dpak = dhmax
3874  write (cloc, "(a,'-',a)") &
3875  trim(this%packName), 'stage'
3876  cpak = trim(cloc)
3877  end if
3878  if (abs(residmax) > abs(dpak)) then
3879  ipak = locresidmax
3880  dpak = residmax
3881  write (cloc, "(a,'-',a)") &
3882  trim(this%packName), 'residual'
3883  cpak = trim(cloc)
3884  end if
3885  if (abs(dgwfmax) > abs(dpak)) then
3886  ipak = locdgwfmax
3887  dpak = dgwfmax
3888  write (cloc, "(a,'-',a)") &
3889  trim(this%packName), 'gwf'
3890  cpak = trim(cloc)
3891  end if
3892  if (this%noutlets > 0) then
3893  if (abs(dqoutmax) > abs(dpak)) then
3894  ipak = locdqoutmax
3895  dpak = dqoutmax
3896  write (cloc, "(a,'-',a)") &
3897  trim(this%packName), 'outlet'
3898  cpak = trim(cloc)
3899  end if
3900  end if
3901  if (this%imover == 1) then
3902  if (abs(dqfrommvrmax) > abs(dpak)) then
3903  ipak = locdqfrommvrmax
3904  dpak = dqfrommvrmax
3905  write (cloc, "(a,'-',a)") trim(this%packName), 'qfrommvr'
3906  cpak = trim(cloc)
3907  end if
3908  end if
3909  !
3910  ! -- write convergence data to package csv
3911  if (this%ipakcsv /= 0) then
3912  !
3913  ! -- write the data
3914  call this%pakcsvtab%add_term(innertot)
3915  call this%pakcsvtab%add_term(totim)
3916  call this%pakcsvtab%add_term(kper)
3917  call this%pakcsvtab%add_term(kstp)
3918  call this%pakcsvtab%add_term(kiter)
3919  call this%pakcsvtab%add_term(dhmax)
3920  call this%pakcsvtab%add_term(locdhmax)
3921  call this%pakcsvtab%add_term(residmax)
3922  call this%pakcsvtab%add_term(locresidmax)
3923  call this%pakcsvtab%add_term(dgwfmax)
3924  call this%pakcsvtab%add_term(locdgwfmax)
3925  if (this%noutlets > 0) then
3926  call this%pakcsvtab%add_term(dqoutmax)
3927  call this%pakcsvtab%add_term(locdqoutmax)
3928  end if
3929  if (this%imover == 1) then
3930  call this%pakcsvtab%add_term(dqfrommvrmax)
3931  call this%pakcsvtab%add_term(locdqfrommvrmax)
3932  end if
3933  !
3934  ! -- finalize the package csv
3935  if (iend == 1) then
3936  call this%pakcsvtab%finalize_table()
3937  end if
3938  end if
3939  end if
3940  end subroutine lak_cc
3941 
3942  !> @brief Calculate flows
3943  !<
3944  subroutine lak_cq(this, x, flowja, iadv)
3945  ! -- modules
3946  use tdismodule, only: delt
3947  ! -- dummy
3948  class(laktype), intent(inout) :: this
3949  real(DP), dimension(:), intent(in) :: x
3950  real(DP), dimension(:), contiguous, intent(inout) :: flowja
3951  integer(I4B), optional, intent(in) :: iadv
3952  ! -- local
3953  real(DP) :: rrate
3954  real(DP) :: chratin, chratout
3955  ! -- for budget
3956  integer(I4B) :: j, n
3957  real(DP) :: hlak
3958  real(DP) :: v0, v1
3959  !
3960  call this%lak_solve(update=.false.)
3961  !
3962  ! -- call base functionality in bnd_cq. This will calculate lake-gwf flows
3963  ! and put them into this%simvals
3964  call this%BndType%bnd_cq(x, flowja, iadv=1)
3965  !
3966  ! -- calculate several budget terms
3967  chratin = dzero
3968  chratout = dzero
3969  do n = 1, this%nlakes
3970  this%chterm(n) = dzero
3971  if (this%iboundpak(n) == 0) cycle
3972  hlak = this%xnewpak(n)
3973  call this%lak_calculate_vol(n, hlak, v1)
3974  !
3975  ! -- add budget terms for active lakes
3976  if (this%iboundpak(n) /= 0) then
3977  !
3978  ! -- rainfall
3979  rrate = this%precip(n)
3980  call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
3981  !
3982  ! -- evaporation
3983  rrate = this%evap(n)
3984  call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
3985  !
3986  ! -- runoff
3987  rrate = this%runoff(n)
3988  call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
3989  !
3990  ! -- inflow
3991  rrate = this%inflow(n)
3992  call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
3993  !
3994  ! -- withdrawals
3995  rrate = this%withr(n)
3996  call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
3997  !
3998  ! -- add lake storage changes
3999  rrate = dzero
4000  if (this%iboundpak(n) > 0) then
4001  if (this%gwfiss /= 1) then
4002  call this%lak_calculate_vol(n, this%xoldpak(n), v0)
4003  rrate = -(v1 - v0) / delt
4004  call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4005  end if
4006  end if
4007  this%qsto(n) = rrate
4008  !
4009  ! -- add external outlets
4010  call this%lak_get_external_outlet(n, rrate)
4011  call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4012  !
4013  ! -- add mover terms
4014  if (this%imover == 1) then
4015  if (this%iboundpak(n) /= 0) then
4016  rrate = this%pakmvrobj%get_qfrommvr(n)
4017  else
4018  rrate = dzero
4019  end if
4020  call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4021  end if
4022  end if
4023  end do
4024  !
4025  ! -- gwf flow and constant flow to lake
4026  do n = 1, this%nlakes
4027  rrate = dzero
4028  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4029  ! simvals is from aquifer perspective, and so it is positive
4030  ! for flow into the aquifer. Need to switch sign for lake
4031  ! perspective.
4032  rrate = -this%simvals(j)
4033  this%qleak(j) = rrate
4034  if (this%iboundpak(n) /= 0) then
4035  call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
4036  end if
4037  end do
4038  end do
4039  !
4040  ! -- fill the budget object
4041  call this%lak_fill_budobj()
4042  end subroutine lak_cq
4043 
4044  !> @brief Output LAK package flow terms
4045  !<
4046  subroutine lak_ot_package_flows(this, icbcfl, ibudfl)
4047  use tdismodule, only: kstp, kper, delt, pertim, totim
4048  class(laktype) :: this
4049  integer(I4B), intent(in) :: icbcfl
4050  integer(I4B), intent(in) :: ibudfl
4051  integer(I4B) :: ibinun
4052  !
4053  ! -- write the flows from the budobj
4054  ibinun = 0
4055  if (this%ibudgetout /= 0) then
4056  ibinun = this%ibudgetout
4057  end if
4058  if (icbcfl == 0) ibinun = 0
4059  if (ibinun > 0) then
4060  call this%budobj%save_flows(this%dis, ibinun, kstp, kper, delt, &
4061  pertim, totim, this%iout)
4062  end if
4063  !
4064  ! -- Print lake flows table
4065  if (ibudfl /= 0 .and. this%iprflow /= 0) then
4066  call this%budobj%write_flowtable(this%dis, kstp, kper)
4067  end if
4068  end subroutine lak_ot_package_flows
4069 
4070  !> @brief Write flows to binary file and/or print flows to budget
4071  !<
4072  subroutine lak_ot_model_flows(this, icbcfl, ibudfl, icbcun, imap)
4073  class(laktype) :: this
4074  integer(I4B), intent(in) :: icbcfl
4075  integer(I4B), intent(in) :: ibudfl
4076  integer(I4B), intent(in) :: icbcun
4077  integer(I4B), dimension(:), optional, intent(in) :: imap
4078  !
4079  ! -- write the flows from the budobj
4080  call this%BndType%bnd_ot_model_flows(icbcfl, ibudfl, icbcun, this%imap)
4081  end subroutine lak_ot_model_flows
4082 
4083  !> @brief Save LAK-calculated values to binary file
4084  !<
4085  subroutine lak_ot_dv(this, idvsave, idvprint)
4086  use tdismodule, only: kstp, kper, pertim, totim
4087  use constantsmodule, only: dhnoflo, dhdry
4088  use inputoutputmodule, only: ulasav
4089  class(laktype) :: this
4090  integer(I4B), intent(in) :: idvsave
4091  integer(I4B), intent(in) :: idvprint
4092  integer(I4B) :: ibinun
4093  integer(I4B) :: n
4094  real(DP) :: v
4095  real(DP) :: d
4096  real(DP) :: stage
4097  real(DP) :: sa
4098  real(DP) :: wa
4099  !
4100  ! -- set unit number for binary dependent variable output
4101  ibinun = 0
4102  if (this%istageout /= 0) then
4103  ibinun = this%istageout
4104  end if
4105  if (idvsave == 0) ibinun = 0
4106  !
4107  ! -- write lake binary output
4108  if (ibinun > 0) then
4109  do n = 1, this%nlakes
4110  v = this%xnewpak(n)
4111  d = v - this%lakebot(n)
4112  if (this%iboundpak(n) == 0) then
4113  v = dhnoflo
4114  else if (d <= dzero) then
4115  v = dhdry
4116  end if
4117  this%dbuff(n) = v
4118  end do
4119  call ulasav(this%dbuff, ' STAGE', kstp, kper, pertim, totim, &
4120  this%nlakes, 1, 1, ibinun)
4121  end if
4122  !
4123  ! -- Print lake stage table
4124  if (idvprint /= 0 .and. this%iprhed /= 0) then
4125  !
4126  ! -- set table kstp and kper
4127  call this%stagetab%set_kstpkper(kstp, kper)
4128  !
4129  ! -- write data
4130  do n = 1, this%nlakes
4131  if (this%iboundpak(n) == 0) then
4132  stage = dhnoflo
4133  sa = dhnoflo
4134  wa = dhnoflo
4135  v = dhnoflo
4136  else
4137  stage = this%xnewpak(n)
4138  call this%lak_calculate_sarea(n, stage, sa)
4139  call this%lak_calculate_warea(n, stage, wa)
4140  call this%lak_calculate_vol(n, stage, v)
4141  end if
4142  if (this%inamedbound == 1) then
4143  call this%stagetab%add_term(this%lakename(n))
4144  end if
4145  call this%stagetab%add_term(n)
4146  call this%stagetab%add_term(stage)
4147  call this%stagetab%add_term(sa)
4148  call this%stagetab%add_term(wa)
4149  call this%stagetab%add_term(v)
4150  end do
4151  end if
4152  end subroutine lak_ot_dv
4153 
4154  !> @brief Write LAK budget to listing file
4155  !<
4156  subroutine lak_ot_bdsummary(this, kstp, kper, iout, ibudfl)
4157  ! -- module
4158  use tdismodule, only: totim, delt
4159  ! -- dummy
4160  class(laktype) :: this !< LakType object
4161  integer(I4B), intent(in) :: kstp !< time step number
4162  integer(I4B), intent(in) :: kper !< period number
4163  integer(I4B), intent(in) :: iout !< flag and unit number for the model listing file
4164  integer(I4B), intent(in) :: ibudfl !< flag indicating budget should be written
4165  !
4166  call this%budobj%write_budtable(kstp, kper, iout, ibudfl, totim, delt)
4167  end subroutine lak_ot_bdsummary
4168 
4169  !> @brief Deallocate objects
4170  !<
4171  subroutine lak_da(this)
4172  ! -- modules
4173  use memorymanagermodule, only: mem_deallocate
4174  ! -- dummy
4175  class(laktype) :: this
4176  !
4177  ! -- arrays
4178  deallocate (this%lakename)
4179  deallocate (this%status)
4180  deallocate (this%clakbudget)
4181  call mem_deallocate(this%dbuff)
4182  deallocate (this%cauxcbc)
4183  call mem_deallocate(this%qauxcbc)
4184  call mem_deallocate(this%qleak)
4185  call mem_deallocate(this%qsto)
4186  call mem_deallocate(this%denseterms)
4187  call mem_deallocate(this%viscratios)
4188  !
4189  ! -- tables
4190  if (this%ntables > 0) then
4191  call mem_deallocate(this%ialaktab)
4192  call mem_deallocate(this%tabstage)
4193  call mem_deallocate(this%tabvolume)
4194  call mem_deallocate(this%tabsarea)
4195  call mem_deallocate(this%tabwarea)
4196  end if
4197  !
4198  ! -- budobj
4199  call this%budobj%budgetobject_da()
4200  deallocate (this%budobj)
4201  nullify (this%budobj)
4202  !
4203  ! -- outlets
4204  if (this%noutlets > 0) then
4205  call mem_deallocate(this%lakein)
4206  call mem_deallocate(this%lakeout)
4207  call mem_deallocate(this%iouttype)
4208  call mem_deallocate(this%outrate)
4209  call mem_deallocate(this%outinvert)
4210  call mem_deallocate(this%outwidth)
4211  call mem_deallocate(this%outrough)
4212  call mem_deallocate(this%outslope)
4213  call mem_deallocate(this%simoutrate)
4214  end if
4215  !
4216  ! -- stage table
4217  if (this%iprhed > 0) then
4218  call this%stagetab%table_da()
4219  deallocate (this%stagetab)
4220  nullify (this%stagetab)
4221  end if
4222  !
4223  ! -- package csv table
4224  if (this%ipakcsv > 0) then
4225  if (associated(this%pakcsvtab)) then
4226  call this%pakcsvtab%table_da()
4227  deallocate (this%pakcsvtab)
4228  nullify (this%pakcsvtab)
4229  end if
4230  end if
4231  !
4232  ! -- scalars
4233  call mem_deallocate(this%iprhed)
4234  call mem_deallocate(this%istageout)
4235  call mem_deallocate(this%ibudgetout)
4236  call mem_deallocate(this%ibudcsv)
4237  call mem_deallocate(this%ipakcsv)
4238  call mem_deallocate(this%nlakes)
4239  call mem_deallocate(this%noutlets)
4240  call mem_deallocate(this%ntables)
4241  call mem_deallocate(this%convlength)
4242  call mem_deallocate(this%convtime)
4243  call mem_deallocate(this%outdmax)
4244  call mem_deallocate(this%igwhcopt)
4245  call mem_deallocate(this%iconvchk)
4246  call mem_deallocate(this%iconvresidchk)
4247  call mem_deallocate(this%maxlakit)
4248  call mem_deallocate(this%surfdep)
4249  call mem_deallocate(this%dmaxchg)
4250  call mem_deallocate(this%delh)
4251  call mem_deallocate(this%pdmax)
4252  call mem_deallocate(this%check_attr)
4253  call mem_deallocate(this%bditems)
4254  call mem_deallocate(this%cbcauxitems)
4255  call mem_deallocate(this%idense)
4256  !
4257  call mem_deallocate(this%nlakeconn)
4258  call mem_deallocate(this%idxlakeconn)
4259  call mem_deallocate(this%ntabrow)
4260  call mem_deallocate(this%strt)
4261  call mem_deallocate(this%laketop)
4262  call mem_deallocate(this%lakebot)
4263  call mem_deallocate(this%sareamax)
4264  call mem_deallocate(this%stage)
4265  call mem_deallocate(this%rainfall)
4266  call mem_deallocate(this%evaporation)
4267  call mem_deallocate(this%runoff)
4268  call mem_deallocate(this%inflow)
4269  call mem_deallocate(this%withdrawal)
4270  call mem_deallocate(this%lauxvar)
4271  call mem_deallocate(this%avail)
4272  call mem_deallocate(this%lkgwsink)
4273  call mem_deallocate(this%ncncvr)
4274  call mem_deallocate(this%surfin)
4275  call mem_deallocate(this%surfout)
4276  call mem_deallocate(this%surfout1)
4277  call mem_deallocate(this%precip)
4278  call mem_deallocate(this%precip1)
4279  call mem_deallocate(this%evap)
4280  call mem_deallocate(this%evap1)
4281  call mem_deallocate(this%evapo)
4282  call mem_deallocate(this%withr)
4283  call mem_deallocate(this%withr1)
4284  call mem_deallocate(this%flwin)
4285  call mem_deallocate(this%flwiter)
4286  call mem_deallocate(this%flwiter1)
4287  call mem_deallocate(this%seep)
4288  call mem_deallocate(this%seep1)
4289  call mem_deallocate(this%seep0)
4290  call mem_deallocate(this%stageiter)
4291  call mem_deallocate(this%chterm)
4292  !
4293  ! -- lake boundary and stages
4294  call mem_deallocate(this%iboundpak)
4295  call mem_deallocate(this%xnewpak)
4296  call mem_deallocate(this%xoldpak)
4297  !
4298  ! -- lake iteration variables
4299  call mem_deallocate(this%iseepc)
4300  call mem_deallocate(this%idhc)
4301  call mem_deallocate(this%en1)
4302  call mem_deallocate(this%en2)
4303  call mem_deallocate(this%r1)
4304  call mem_deallocate(this%r2)
4305  call mem_deallocate(this%dh0)
4306  call mem_deallocate(this%s0)
4307  call mem_deallocate(this%qgwf0)
4308  !
4309  ! -- lake connection variables
4310  call mem_deallocate(this%imap)
4311  call mem_deallocate(this%cellid)
4312  call mem_deallocate(this%nodesontop)
4313  call mem_deallocate(this%ictype)
4314  call mem_deallocate(this%bedleak)
4315  call mem_deallocate(this%belev)
4316  call mem_deallocate(this%telev)
4317  call mem_deallocate(this%connlength)
4318  call mem_deallocate(this%connwidth)
4319  call mem_deallocate(this%sarea)
4320  call mem_deallocate(this%warea)
4321  call mem_deallocate(this%satcond)
4322  call mem_deallocate(this%simcond)
4323  call mem_deallocate(this%simlakgw)
4324  !
4325  ! -- pointers to gwf variables
4326  nullify (this%gwfiss)
4327  !
4328  ! -- Parent object
4329  call this%BndType%bnd_da()
4330  end subroutine lak_da
4331 
4332  !> @brief Define the list heading that is written to iout when PRINT_INPUT
4333  !! option is used
4334  !<
4335  subroutine define_listlabel(this)
4336  ! -- modules
4337  class(laktype), intent(inout) :: this
4338  !
4339  ! -- create the header list label
4340  this%listlabel = trim(this%filtyp)//' NO.'
4341  if (this%dis%ndim == 3) then
4342  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
4343  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW'
4344  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL'
4345  elseif (this%dis%ndim == 2) then
4346  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
4347  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D'
4348  else
4349  write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE'
4350  end if
4351  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'STRESS RATE'
4352  if (this%inamedbound == 1) then
4353  write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
4354  end if
4355  end subroutine define_listlabel
4356 
4357  !> @brief Set pointers to model arrays and variables so that a package has
4358  !! access to these things
4359  !<
4360  subroutine lak_set_pointers(this, neq, ibound, xnew, xold, flowja)
4361  ! -- dummy
4362  class(laktype) :: this
4363  integer(I4B), pointer :: neq
4364  integer(I4B), dimension(:), pointer, contiguous :: ibound
4365  real(DP), dimension(:), pointer, contiguous :: xnew
4366  real(DP), dimension(:), pointer, contiguous :: xold
4367  real(DP), dimension(:), pointer, contiguous :: flowja
4368  !
4369  ! -- call base BndType set_pointers
4370  call this%BndType%set_pointers(neq, ibound, xnew, xold, flowja)
4371  !
4372  ! -- Set the LAK pointers
4373  !
4374  ! -- set package pointers
4375  !istart = this%dis%nodes + this%ioffset + 1
4376  !iend = istart + this%nlakes - 1
4377  !this%iboundpak => this%ibound(istart:iend)
4378  !this%xnewpak => this%xnew(istart:iend)
4379  !
4380  ! -- initialize xnewpak
4381  !do n = 1, this%nlakes
4382  ! this%xnewpak(n) = DEP20
4383  !end do
4384  end subroutine lak_set_pointers
4385 
4386  !> @brief Procedures related to observations (type-bound)
4387  !!
4388  !! Return true because LAK package supports observations. Overrides
4389  !! BndType%bnd_obs_supported()
4390  !<
4391  logical function lak_obs_supported(this)
4392  ! -- dummy
4393  class(laktype) :: this
4394  !
4395  lak_obs_supported = .true.
4396  end function lak_obs_supported
4397 
4398  !> @brief Store observation type supported by LAK package. Overrides
4399  !! BndType%bnd_df_obs
4400  !<
4401  subroutine lak_df_obs(this)
4402  ! -- dummy
4403  class(laktype) :: this
4404  ! -- local
4405  integer(I4B) :: indx
4406  !
4407  ! -- Store obs type and assign procedure pointer
4408  ! for stage observation type.
4409  call this%obs%StoreObsType('stage', .false., indx)
4410  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4411  !
4412  ! -- Store obs type and assign procedure pointer
4413  ! for ext-inflow observation type.
4414  call this%obs%StoreObsType('ext-inflow', .true., indx)
4415  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4416  !
4417  ! -- Store obs type and assign procedure pointer
4418  ! for outlet-inflow observation type.
4419  call this%obs%StoreObsType('outlet-inflow', .true., indx)
4420  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4421  !
4422  ! -- Store obs type and assign procedure pointer
4423  ! for inflow observation type.
4424  call this%obs%StoreObsType('inflow', .true., indx)
4425  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4426  !
4427  ! -- Store obs type and assign procedure pointer
4428  ! for from-mvr observation type.
4429  call this%obs%StoreObsType('from-mvr', .true., indx)
4430  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4431  !
4432  ! -- Store obs type and assign procedure pointer
4433  ! for rainfall observation type.
4434  call this%obs%StoreObsType('rainfall', .true., indx)
4435  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4436  !
4437  ! -- Store obs type and assign procedure pointer
4438  ! for runoff observation type.
4439  call this%obs%StoreObsType('runoff', .true., indx)
4440  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4441  !
4442  ! -- Store obs type and assign procedure pointer
4443  ! for lak observation type.
4444  call this%obs%StoreObsType('lak', .true., indx)
4445  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4446  !
4447  ! -- Store obs type and assign procedure pointer
4448  ! for evaporation observation type.
4449  call this%obs%StoreObsType('evaporation', .true., indx)
4450  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4451  !
4452  ! -- Store obs type and assign procedure pointer
4453  ! for withdrawal observation type.
4454  call this%obs%StoreObsType('withdrawal', .true., indx)
4455  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4456  !
4457  ! -- Store obs type and assign procedure pointer
4458  ! for ext-outflow observation type.
4459  call this%obs%StoreObsType('ext-outflow', .true., indx)
4460  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4461  !
4462  ! -- Store obs type and assign procedure pointer
4463  ! for to-mvr observation type.
4464  call this%obs%StoreObsType('to-mvr', .true., indx)
4465  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4466  !
4467  ! -- Store obs type and assign procedure pointer
4468  ! for storage observation type.
4469  call this%obs%StoreObsType('storage', .true., indx)
4470  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4471  !
4472  ! -- Store obs type and assign procedure pointer
4473  ! for constant observation type.
4474  call this%obs%StoreObsType('constant', .true., indx)
4475  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4476  !
4477  ! -- Store obs type and assign procedure pointer
4478  ! for outlet observation type.
4479  call this%obs%StoreObsType('outlet', .true., indx)
4480  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4481  !
4482  ! -- Store obs type and assign procedure pointer
4483  ! for volume observation type.
4484  call this%obs%StoreObsType('volume', .true., indx)
4485  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4486  !
4487  ! -- Store obs type and assign procedure pointer
4488  ! for surface-area observation type.
4489  call this%obs%StoreObsType('surface-area', .true., indx)
4490  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4491  !
4492  ! -- Store obs type and assign procedure pointer
4493  ! for wetted-area observation type.
4494  call this%obs%StoreObsType('wetted-area', .true., indx)
4495  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4496  !
4497  ! -- Store obs type and assign procedure pointer
4498  ! for conductance observation type.
4499  call this%obs%StoreObsType('conductance', .true., indx)
4500  this%obs%obsData(indx)%ProcessIdPtr => lak_process_obsid
4501  end subroutine lak_df_obs
4502 
4503  !> @brief Calculate observations this time step and call ObsType%SaveOneSimval
4504  !! for each LakType observation.
4505  !<
4506  subroutine lak_bd_obs(this)
4507  ! -- dummy
4508  class(laktype) :: this
4509  ! -- local
4510  integer(I4B) :: i
4511  integer(I4B) :: igwfnode
4512  integer(I4B) :: j
4513  integer(I4B) :: jj
4514  integer(I4B) :: n
4515  real(DP) :: hgwf
4516  real(DP) :: hlak
4517  real(DP) :: v
4518  real(DP) :: v2
4519  type(observetype), pointer :: obsrv => null()
4520  !
4521  ! Write simulated values for all LAK observations
4522  if (this%obs%npakobs > 0) then
4523  call this%obs%obs_bd_clear()
4524  do i = 1, this%obs%npakobs
4525  obsrv => this%obs%pakobs(i)%obsrv
4526  do j = 1, obsrv%indxbnds_count
4527  v = dnodata
4528  jj = obsrv%indxbnds(j)
4529  select case (obsrv%ObsTypeId)
4530  case ('STAGE')
4531  if (this%iboundpak(jj) /= 0) then
4532  v = this%xnewpak(jj)
4533  end if
4534  case ('EXT-INFLOW')
4535  if (this%iboundpak(jj) /= 0) then
4536  call this%lak_calculate_inflow(jj, v)
4537  end if
4538  case ('OUTLET-INFLOW')
4539  if (this%iboundpak(jj) /= 0) then
4540  call this%lak_calculate_outlet_inflow(jj, v)
4541  end if
4542  case ('INFLOW')
4543  if (this%iboundpak(jj) /= 0) then
4544  call this%lak_calculate_inflow(jj, v)
4545  call this%lak_calculate_outlet_inflow(jj, v2)
4546  v = v + v2
4547  end if
4548  case ('FROM-MVR')
4549  if (this%iboundpak(jj) /= 0) then
4550  if (this%imover == 1) then
4551  v = this%pakmvrobj%get_qfrommvr(jj)
4552  end if
4553  end if
4554  case ('RAINFALL')
4555  if (this%iboundpak(jj) /= 0) then
4556  v = this%precip(jj)
4557  end if
4558  case ('RUNOFF')
4559  if (this%iboundpak(jj) /= 0) then
4560  v = this%runoff(jj)
4561  end if
4562  case ('LAK')
4563  n = this%imap(jj)
4564  if (this%iboundpak(n) /= 0) then
4565  igwfnode = this%cellid(jj)
4566  hgwf = this%xnew(igwfnode)
4567  if (this%hcof(jj) /= dzero) then
4568  v = -(this%hcof(jj) * (this%xnewpak(n) - hgwf))
4569  else
4570  v = -this%rhs(jj)
4571  end if
4572  end if
4573  case ('EVAPORATION')
4574  if (this%iboundpak(jj) /= 0) then
4575  v = this%evap(jj)
4576  end if
4577  case ('WITHDRAWAL')
4578  if (this%iboundpak(jj) /= 0) then
4579  v = this%withr(jj)
4580  end if
4581  case ('EXT-OUTFLOW')
4582  n = this%lakein(jj)
4583  if (this%iboundpak(n) /= 0) then
4584  if (this%lakeout(jj) == 0) then
4585  v = this%simoutrate(jj)
4586  if (v < dzero) then
4587  if (this%imover == 1) then
4588  v = v + this%pakmvrobj%get_qtomvr(jj)
4589  end if
4590  end if
4591  end if
4592  end if
4593  case ('TO-MVR')
4594  n = this%lakein(jj)
4595  if (this%iboundpak(n) /= 0) then
4596  if (this%imover == 1) then
4597  v = this%pakmvrobj%get_qtomvr(jj)
4598  if (v > dzero) then
4599  v = -v
4600  end if
4601  end if
4602  end if
4603  case ('STORAGE')
4604  if (this%iboundpak(jj) /= 0) then
4605  v = this%qsto(jj)
4606  end if
4607  case ('CONSTANT')
4608  if (this%iboundpak(jj) /= 0) then
4609  v = this%chterm(jj)
4610  end if
4611  case ('OUTLET')
4612  n = this%lakein(jj)
4613  if (this%iboundpak(n) /= 0) then
4614  v = this%simoutrate(jj)
4615  end if
4616  case ('VOLUME')
4617  if (this%iboundpak(jj) /= 0) then
4618  call this%lak_calculate_vol(jj, this%xnewpak(jj), v)
4619  end if
4620  case ('SURFACE-AREA')
4621  if (this%iboundpak(jj) /= 0) then
4622  hlak = this%xnewpak(jj)
4623  call this%lak_calculate_sarea(jj, hlak, v)
4624  end if
4625  case ('WETTED-AREA')
4626  n = this%imap(jj)
4627  if (this%iboundpak(n) /= 0) then
4628  hlak = this%xnewpak(n)
4629  igwfnode = this%cellid(jj)
4630  hgwf = this%xnew(igwfnode)
4631  call this%lak_calculate_conn_warea(n, jj, hlak, hgwf, v)
4632  end if
4633  case ('CONDUCTANCE')
4634  n = this%imap(jj)
4635  if (this%iboundpak(n) /= 0) then
4636  hlak = this%xnewpak(n)
4637  igwfnode = this%cellid(jj)
4638  hgwf = this%xnew(igwfnode)
4639  call this%lak_calculate_conn_conductance(n, jj, hlak, hgwf, v)
4640  end if
4641  case default
4642  errmsg = 'Unrecognized observation type: '//trim(obsrv%ObsTypeId)
4643  call store_error(errmsg)
4644  end select
4645  call this%obs%SaveOneSimval(obsrv, v)
4646  end do
4647  end do
4648  !
4649  ! -- write summary of error messages
4650  if (count_errors() > 0) then
4651  call store_error_unit(this%inunit)
4652  end if
4653  end if
4654  end subroutine lak_bd_obs
4655 
4656  !> @brief Process each observation
4657  !!
4658  !! Only done the first stress period since boundaries are fixed for the
4659  !! simulation
4660  !<
4661  subroutine lak_rp_obs(this)
4662  use tdismodule, only: kper
4663  ! -- dummy
4664  class(laktype), intent(inout) :: this
4665  ! -- local
4666  integer(I4B) :: i
4667  integer(I4B) :: j
4668  integer(I4B) :: nn1
4669  integer(I4B) :: nn2
4670  integer(I4B) :: jj
4671  character(len=LENBOUNDNAME) :: bname
4672  logical(LGP) :: jfound
4673  class(observetype), pointer :: obsrv => null()
4674  ! -- formats
4675 10 format('Boundary "', a, '" for observation "', a, &
4676  '" is invalid in package "', a, '"')
4677  !
4678  ! -- process each package observation
4679  ! only done the first stress period since boundaries are fixed
4680  ! for the simulation
4681  if (kper == 1) then
4682  do i = 1, this%obs%npakobs
4683  obsrv => this%obs%pakobs(i)%obsrv
4684  !
4685  ! -- get node number 1
4686  nn1 = obsrv%NodeNumber
4687  if (nn1 == namedboundflag) then
4688  bname = obsrv%FeatureName
4689  if (bname /= '') then
4690  ! -- Observation lake is based on a boundary name.
4691  ! Iterate through all lakes to identify and store
4692  ! corresponding index in bound array.
4693  jfound = .false.
4694  if (obsrv%ObsTypeId == 'LAK' .or. &
4695  obsrv%ObsTypeId == 'CONDUCTANCE' .or. &
4696  obsrv%ObsTypeId == 'WETTED-AREA') then
4697  do j = 1, this%nlakes
4698  do jj = this%idxlakeconn(j), this%idxlakeconn(j + 1) - 1
4699  if (this%boundname(jj) == bname) then
4700  jfound = .true.
4701  call obsrv%AddObsIndex(jj)
4702  end if
4703  end do
4704  end do
4705  else if (obsrv%ObsTypeId == 'EXT-OUTFLOW' .or. &
4706  obsrv%ObsTypeId == 'TO-MVR' .or. &
4707  obsrv%ObsTypeId == 'OUTLET') then
4708  do j = 1, this%noutlets
4709  jj = this%lakein(j)
4710  if (this%lakename(jj) == bname) then
4711  jfound = .true.
4712  call obsrv%AddObsIndex(j)
4713  end if
4714  end do
4715  else
4716  do j = 1, this%nlakes
4717  if (this%lakename(j) == bname) then
4718  jfound = .true.
4719  call obsrv%AddObsIndex(j)
4720  end if
4721  end do
4722  end if
4723  if (.not. jfound) then
4724  write (errmsg, 10) &
4725  trim(bname), trim(obsrv%Name), trim(this%packName)
4726  call store_error(errmsg)
4727  end if
4728  end if
4729  else
4730  if (obsrv%indxbnds_count == 0) then
4731  if (obsrv%ObsTypeId == 'LAK' .or. &
4732  obsrv%ObsTypeId == 'CONDUCTANCE' .or. &
4733  obsrv%ObsTypeId == 'WETTED-AREA') then
4734  nn2 = obsrv%NodeNumber2
4735  j = this%idxlakeconn(nn1) + nn2 - 1
4736  call obsrv%AddObsIndex(j)
4737  else
4738  call obsrv%AddObsIndex(nn1)
4739  end if
4740  else
4741  errmsg = 'Programming error in lak_rp_obs'
4742  call store_error(errmsg)
4743  end if
4744  end if
4745  !
4746  ! -- catch non-cumulative observation assigned to observation defined
4747  ! by a boundname that is assigned to more than one element
4748  if (obsrv%ObsTypeId == 'STAGE') then
4749  if (obsrv%indxbnds_count > 1) then
4750  write (errmsg, '(a,3(1x,a))') &
4751  trim(adjustl(obsrv%ObsTypeId)), &
4752  'for observation', trim(adjustl(obsrv%Name)), &
4753  ' must be assigned to a lake with a unique boundname.'
4754  call store_error(errmsg)
4755  end if
4756  end if
4757  !
4758  ! -- check that index values are valid
4759  if (obsrv%ObsTypeId == 'TO-MVR' .or. &
4760  obsrv%ObsTypeId == 'EXT-OUTFLOW' .or. &
4761  obsrv%ObsTypeId == 'OUTLET') then
4762  do j = 1, obsrv%indxbnds_count
4763  nn1 = obsrv%indxbnds(j)
4764  if (nn1 < 1 .or. nn1 > this%noutlets) then
4765  write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
4766  trim(adjustl(obsrv%ObsTypeId)), &
4767  ' outlet must be > 0 and <=', this%noutlets, &
4768  '(specified value is ', nn1, ')'
4769  call store_error(errmsg)
4770  end if
4771  end do
4772  else if (obsrv%ObsTypeId == 'LAK' .or. &
4773  obsrv%ObsTypeId == 'CONDUCTANCE' .or. &
4774  obsrv%ObsTypeId == 'WETTED-AREA') then
4775  do j = 1, obsrv%indxbnds_count
4776  nn1 = obsrv%indxbnds(j)
4777  if (nn1 < 1 .or. nn1 > this%maxbound) then
4778  write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
4779  trim(adjustl(obsrv%ObsTypeId)), &
4780  'lake connection number must be > 0 and <=', this%maxbound, &
4781  '(specified value is ', nn1, ')'
4782  call store_error(errmsg)
4783  end if
4784  end do
4785  else
4786  do j = 1, obsrv%indxbnds_count
4787  nn1 = obsrv%indxbnds(j)
4788  if (nn1 < 1 .or. nn1 > this%nlakes) then
4789  write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
4790  trim(adjustl(obsrv%ObsTypeId)), &
4791  ' lake must be > 0 and <=', this%nlakes, &
4792  '(specified value is ', nn1, ')'
4793  call store_error(errmsg)
4794  end if
4795  end do
4796  end if
4797  end do
4798  !
4799  ! -- evaluate if there are any observation errors
4800  if (count_errors() > 0) then
4801  call store_error_unit(this%inunit)
4802  end if
4803  end if
4804  end subroutine lak_rp_obs
4805 
4806  !
4807  ! -- Procedures related to observations (NOT type-bound)
4808 
4809  !> @brief This procedure is pointed to by ObsDataType%ProcesssIdPtr. It
4810  !! processes the ID string of an observation definition for LAK package
4811  !! observations.
4812  !<
4813  subroutine lak_process_obsid(obsrv, dis, inunitobs, iout)
4814  ! -- dummy
4815  type(observetype), intent(inout) :: obsrv
4816  class(disbasetype), intent(in) :: dis
4817  integer(I4B), intent(in) :: inunitobs
4818  integer(I4B), intent(in) :: iout
4819  ! -- local
4820  integer(I4B) :: nn1, nn2
4821  integer(I4B) :: icol, istart, istop
4822  character(len=LINELENGTH) :: string
4823  character(len=LENBOUNDNAME) :: bndname
4824  !
4825  string = obsrv%IDstring
4826  ! -- Extract lake number from string and store it.
4827  ! If 1st item is not an integer(I4B), it should be a
4828  ! lake name--deal with it.
4829  icol = 1
4830  ! -- get lake number or boundary name
4831  call extract_idnum_or_bndname(string, icol, istart, istop, nn1, bndname)
4832  if (nn1 == namedboundflag) then
4833  obsrv%FeatureName = bndname
4834  else
4835  if (obsrv%ObsTypeId == 'LAK' .or. obsrv%ObsTypeId == 'CONDUCTANCE' .or. &
4836  obsrv%ObsTypeId == 'WETTED-AREA') then
4837  call extract_idnum_or_bndname(string, icol, istart, istop, nn2, bndname)
4838  if (len_trim(bndname) < 1 .and. nn2 < 0) then
4839  write (errmsg, '(a,1x,a,a,1x,a,1x,a)') &
4840  'For observation type', trim(adjustl(obsrv%ObsTypeId)), &
4841  ', ID given as an integer and not as boundname,', &
4842  'but ID2 (iconn) is missing. Either change ID to valid', &
4843  'boundname or supply valid entry for ID2.'
4844  call store_error(errmsg)
4845  end if
4846  if (nn2 == namedboundflag) then
4847  obsrv%FeatureName = bndname
4848  ! -- reset nn1
4849  nn1 = nn2
4850  else
4851  obsrv%NodeNumber2 = nn2
4852  end if
4853  end if
4854  end if
4855  ! -- store lake number (NodeNumber)
4856  obsrv%NodeNumber = nn1
4857  end subroutine lak_process_obsid
4858 
4859  !
4860  ! -- private LAK methods
4861  !
4862 
4863  !> @brief Accumulate constant head terms for budget
4864  !<
4865  subroutine lak_accumulate_chterm(this, ilak, rrate, chratin, chratout)
4866  ! -- dummy
4867  class(laktype) :: this
4868  integer(I4B), intent(in) :: ilak
4869  real(DP), intent(in) :: rrate
4870  real(DP), intent(inout) :: chratin
4871  real(DP), intent(inout) :: chratout
4872  ! -- locals
4873  real(DP) :: q
4874  !
4875  ! code
4876  if (this%iboundpak(ilak) < 0) then
4877  q = -rrate
4878  this%chterm(ilak) = this%chterm(ilak) + q
4879  !
4880  ! -- See if flow is into lake or out of lake.
4881  if (q < dzero) then
4882  !
4883  ! -- Flow is out of lake subtract rate from ratout.
4884  chratout = chratout - q
4885  else
4886  !
4887  ! -- Flow is into lake; add rate to ratin.
4888  chratin = chratin + q
4889  end if
4890  end if
4891  end subroutine lak_accumulate_chterm
4892 
4893  !> @brief Store the lake head and connection conductance in the bound array
4894  !<
4895  subroutine lak_bound_update(this)
4896  ! -- dummy
4897  class(laktype), intent(inout) :: this
4898  ! -- local
4899  integer(I4B) :: j, n, node
4900  real(DP) :: hlak, head, clak
4901  !
4902  ! -- Return if no lak lakes
4903  if (this%nbound == 0) return
4904  !
4905  ! -- Calculate hcof and rhs for each lak entry
4906  do n = 1, this%nlakes
4907  hlak = this%xnewpak(n)
4908  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
4909  node = this%cellid(j)
4910  head = this%xnew(node)
4911  call this%lak_calculate_conn_conductance(n, j, hlak, head, clak)
4912  this%bound(1, j) = hlak
4913  this%bound(2, j) = clak
4914  end do
4915  end do
4916  end subroutine lak_bound_update
4917 
4918  !> @brief Solve for lake stage
4919  !<
4920  subroutine lak_solve(this, update)
4921  ! -- modules
4922  use tdismodule, only: delt
4923  ! -- dummy
4924  class(laktype), intent(inout) :: this
4925  logical(LGP), intent(in), optional :: update
4926  ! -- local
4927  logical(LGP) :: lupdate
4928  integer(I4B) :: i
4929  integer(I4B) :: j
4930  integer(I4B) :: n
4931  integer(I4B) :: iicnvg
4932  integer(I4B) :: iter
4933  integer(I4B) :: maxiter
4934  integer(I4B) :: ncnv
4935  integer(I4B) :: idry
4936  integer(I4B) :: idry1
4937  integer(I4B) :: igwfnode
4938  integer(I4B) :: ibflg
4939  integer(I4B) :: idhp
4940  real(DP) :: hlak
4941  real(DP) :: hlak0
4942  real(DP) :: v0
4943  real(DP) :: v1
4944  real(DP) :: head
4945  real(DP) :: ra
4946  real(DP) :: ro
4947  real(DP) :: qinf
4948  real(DP) :: ex
4949  real(DP) :: ev
4950  real(DP) :: outinf
4951  real(DP) :: qlakgw
4952  real(DP) :: qlakgw1
4953  real(DP) :: gwfhcof
4954  real(DP) :: gwfrhs
4955  real(DP) :: avail
4956  real(DP) :: resid
4957  real(DP) :: resid1
4958  real(DP) :: residb
4959  real(DP) :: wr
4960  real(DP) :: derv
4961  real(DP) :: dh
4962  real(DP) :: adh
4963  real(DP) :: adh0
4964  real(DP) :: delh
4965  real(DP) :: ts
4966  real(DP) :: area
4967  real(DP) :: qtolfact
4968  !
4969  ! -- set lupdate
4970  if (present(update)) then
4971  lupdate = update
4972  else
4973  lupdate = .true.
4974  end if
4975  !
4976  ! -- initialize
4977  avail = dzero
4978  delh = this%delh
4979  !
4980  ! -- initialize
4981  do n = 1, this%nlakes
4982  this%ncncvr(n) = 0
4983  this%surfin(n) = dzero
4984  this%surfout(n) = dzero
4985  this%surfout1(n) = dzero
4986  if (this%xnewpak(n) < this%lakebot(n)) then
4987  this%xnewpak(n) = this%lakebot(n)
4988  end if
4989  if (this%gwfiss /= 0) then
4990  this%xoldpak(n) = this%xnewpak(n)
4991  end if
4992  ! -- lake iteration items
4993  this%iseepc(n) = 0
4994  this%idhc(n) = 0
4995  this%en1(n) = this%lakebot(n)
4996  call this%lak_calculate_residual(n, this%en1(n), this%r1(n))
4997  this%en2(n) = this%laketop(n)
4998  call this%lak_calculate_residual(n, this%en2(n), this%r2(n))
4999  end do
5000  do n = 1, this%noutlets
5001  this%simoutrate(n) = dzero
5002  end do
5003  !
5004  ! -- sum up inflows from mover inflows
5005  do n = 1, this%nlakes
5006  call this%lak_calculate_outlet_inflow(n, this%surfin(n))
5007  end do
5008  !
5009  ! -- sum up overland runoff, inflows, and external flows into lake
5010  ! (includes maximum lake volume)
5011  do n = 1, this%nlakes
5012  hlak0 = this%xoldpak(n)
5013  hlak = this%xnewpak(n)
5014  call this%lak_calculate_runoff(n, ro)
5015  call this%lak_calculate_inflow(n, qinf)
5016  call this%lak_calculate_external(n, ex)
5017  call this%lak_calculate_vol(n, hlak0, v0)
5018  call this%lak_calculate_vol(n, hlak, v1)
5019  this%flwin(n) = this%surfin(n) + ro + qinf + ex + &
5020  max(v0, v1) / delt
5021  end do
5022  !
5023  ! -- sum up inflows from upstream outlets
5024  do n = 1, this%nlakes
5025  call this%lak_calculate_outlet_inflow(n, outinf)
5026  this%flwin(n) = this%flwin(n) + outinf
5027  end do
5028  !
5029  iicnvg = 0
5030  maxiter = this%maxlakit
5031  !
5032  ! -- outer loop
5033  converge: do iter = 1, maxiter
5034  ncnv = 0
5035  do n = 1, this%nlakes
5036  if (this%ncncvr(n) == 0) ncnv = 1
5037  end do
5038  if (iter == maxiter) ncnv = 0
5039  if (ncnv == 0) iicnvg = 1
5040  !
5041  ! -- initialize variables
5042  do n = 1, this%nlakes
5043  this%evap(n) = dzero
5044  this%precip(n) = dzero
5045  this%precip1(n) = dzero
5046  this%seep(n) = dzero
5047  this%seep1(n) = dzero
5048  this%evap(n) = dzero
5049  this%evap1(n) = dzero
5050  this%evapo(n) = dzero
5051  this%withr(n) = dzero
5052  this%withr1(n) = dzero
5053  this%flwiter(n) = this%flwin(n)
5054  this%flwiter1(n) = this%flwin(n)
5055  if (this%gwfiss /= 0) then
5056  this%flwiter(n) = dep20
5057  this%flwiter1(n) = dep20
5058  end if
5059  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5060  this%hcof(j) = dzero
5061  this%rhs(j) = dzero
5062  end do
5063  end do
5064  !
5065  estseep: do i = 1, 2
5066  lakseep: do n = 1, this%nlakes
5067  ! -- skip inactive lakes
5068  if (this%iboundpak(n) == 0) then
5069  cycle lakseep
5070  end if
5071  ! - set xoldpak to xnewpak if steady-state
5072  if (this%gwfiss /= 0) then
5073  this%xoldpak(n) = this%xnewpak(n)
5074  end if
5075  hlak = this%xnewpak(n)
5076  calcconnseep: do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5077  igwfnode = this%cellid(j)
5078  head = this%xnew(igwfnode)
5079  if (this%ncncvr(n) /= 2) then
5080  if (this%ibound(igwfnode) > 0) then
5081  call this%lak_estimate_conn_exchange(i, n, j, idry, hlak, &
5082  head, qlakgw, &
5083  this%flwiter(n), &
5084  gwfhcof, gwfrhs)
5085  call this%lak_estimate_conn_exchange(i, n, j, idry1, &
5086  hlak + delh, head, qlakgw1, &
5087  this%flwiter1(n))
5088  !
5089  ! -- add to gwf matrix
5090  if (ncnv == 0 .and. i == 2) then
5091  if (j == this%maxbound) then
5092  this%ncncvr(n) = 2
5093  end if
5094  if (idry /= 1) then
5095  this%hcof(j) = gwfhcof
5096  this%rhs(j) = gwfrhs
5097  else
5098  this%hcof(j) = dzero
5099  this%rhs(j) = qlakgw
5100  end if
5101  end if
5102  if (i == 2) then
5103  this%seep(n) = this%seep(n) + qlakgw
5104  this%seep1(n) = this%seep1(n) + qlakgw1
5105  end if
5106  end if
5107  end if
5108  !
5109  end do calcconnseep
5110  end do lakseep
5111  end do estseep
5112  !
5113  laklevel: do n = 1, this%nlakes
5114  ! -- skip inactive lakes
5115  if (this%iboundpak(n) == 0) then
5116  this%ncncvr(n) = 1
5117  cycle laklevel
5118  end if
5119  ibflg = 0
5120  hlak = this%xnewpak(n)
5121  if (iter < maxiter) then
5122  this%stageiter(n) = this%xnewpak(n)
5123  end if
5124  call this%lak_calculate_rainfall(n, hlak, ra)
5125  this%precip(n) = ra
5126  this%flwiter(n) = this%flwiter(n) + ra
5127  call this%lak_calculate_rainfall(n, hlak + delh, ra)
5128  this%precip1(n) = ra
5129  this%flwiter1(n) = this%flwiter1(n) + ra
5130  !
5131  ! -- limit withdrawals to lake inflows and lake storage
5132  call this%lak_calculate_withdrawal(n, this%flwiter(n), wr)
5133  this%withr(n) = wr
5134  call this%lak_calculate_withdrawal(n, this%flwiter1(n), wr)
5135  this%withr1(n) = wr
5136  !
5137  ! -- limit evaporation to lake inflows and lake storage
5138  call this%lak_calculate_evaporation(n, hlak, this%flwiter(n), ev)
5139  this%evap(n) = ev
5140  call this%lak_calculate_evaporation(n, hlak + delh, this%flwiter1(n), ev)
5141  this%evap1(n) = ev
5142  !
5143  ! -- no outlet flow if evaporation consumes all water
5144  call this%lak_calculate_outlet_outflow(n, hlak + delh, &
5145  this%flwiter1(n), &
5146  this%surfout1(n))
5147  call this%lak_calculate_outlet_outflow(n, hlak, this%flwiter(n), &
5148  this%surfout(n))
5149  !
5150  ! -- update the surface inflow values
5151  call this%lak_calculate_outlet_inflow(n, this%surfin(n))
5152  !
5153  !
5154  if (ncnv == 1) then
5155  if (this%iboundpak(n) > 0 .and. lupdate .eqv. .true.) then
5156  !
5157  ! -- recalculate flwin
5158  hlak0 = this%xoldpak(n)
5159  hlak = this%xnewpak(n)
5160  call this%lak_calculate_vol(n, hlak0, v0)
5161  call this%lak_calculate_vol(n, hlak, v1)
5162  call this%lak_calculate_runoff(n, ro)
5163  call this%lak_calculate_inflow(n, qinf)
5164  call this%lak_calculate_external(n, ex)
5165  this%flwin(n) = this%surfin(n) + ro + qinf + ex + &
5166  max(v0, v1) / delt
5167  !
5168  ! -- compute new lake stage using Newton's method
5169  resid = this%precip(n) + this%evap(n) + this%withr(n) + ro + &
5170  qinf + ex + this%surfin(n) + &
5171  this%surfout(n) + this%seep(n)
5172  resid1 = this%precip1(n) + this%evap1(n) + this%withr1(n) + ro + &
5173  qinf + ex + this%surfin(n) + &
5174  this%surfout1(n) + this%seep1(n)
5175  !
5176  ! -- add storage changes for transient stress periods
5177  hlak = this%xnewpak(n)
5178  if (this%gwfiss /= 1) then
5179  call this%lak_calculate_vol(n, hlak, v1)
5180  resid = resid + (v0 - v1) / delt
5181  call this%lak_calculate_vol(n, hlak + delh, v1)
5182  resid1 = resid1 + (v0 - v1) / delt
5183  end if
5184  !
5185  ! -- determine the derivative and the stage change
5186  if (abs(resid1 - resid) > dzero) then
5187  derv = (resid1 - resid) / delh
5188  dh = dzero
5189  if (abs(derv) > dprec) then
5190  dh = resid / derv
5191  end if
5192  else
5193  if (resid < dzero) then
5194  resid = dzero
5195  end if
5196  call this%lak_vol2stage(n, resid, dh)
5197  dh = hlak - dh
5198  this%ncncvr(n) = 1
5199  end if
5200  !
5201  ! -- determine if the updated stage is outside the endpoints
5202  ts = hlak - dh
5203  if (iter == 1) this%dh0(n) = dh
5204  adh = abs(dh)
5205  adh0 = abs(this%dh0(n))
5206  if ((ts >= this%en2(n)) .or. (ts < this%en1(n))) then
5207  ! -- use bisection if dh is increasing or updated stage is below the
5208  ! bottom of the lake
5209  if ((adh > adh0) .or. (ts - this%lakebot(n)) < dprec) then
5210  residb = resid
5211  call this%lak_bisection(n, ibflg, hlak, ts, dh, residb)
5212  end if
5213  end if
5214  !
5215  ! -- set seep0 on the first lake iteration
5216  if (iter == 1) then
5217  this%seep0(n) = this%seep(n)
5218  end if
5219  !
5220  ! -- check for slow convergence
5221  if (this%seep(n) * this%seep0(n) < dprec) then
5222  this%iseepc(n) = this%iseepc(n) + 1
5223  else
5224  this%iseepc(n) = 0
5225  end if
5226  ! -- determine of convergence is slow and oscillating
5227  idhp = 0
5228  if (dh * this%dh0(n) < dprec) idhp = 1
5229  ! -- determine if stage change is increasing
5230  adh = abs(dh)
5231  if (adh > adh0) idhp = 1
5232  ! -- increment idhc convergence flag
5233  if (idhp == 1) then
5234  this%idhc(n) = this%idhc(n) + 1
5235  end if
5236  !
5237  ! -- switch to bisection when the Newton-Raphson method oscillates
5238  ! or when convergence is slow
5239  if (ibflg == 1) then
5240  if (this%iseepc(n) > 7 .or. this%idhc(n) > 12) then
5241  call this%lak_bisection(n, ibflg, hlak, ts, dh, residb)
5242  end if
5243  end if
5244  else
5245  dh = dzero
5246  end if
5247  !
5248  ! -- update lake stage
5249  hlak = hlak - dh
5250  if (hlak < this%lakebot(n)) then
5251  hlak = this%lakebot(n)
5252  end if
5253  !
5254  ! -- calculate surface area
5255  call this%lak_calculate_sarea(n, hlak, area)
5256  !
5257  ! -- set the Q to length factor
5258  if (area > dzero) then
5259  qtolfact = delt / area
5260  else
5261  qtolfact = dzero
5262  end if
5263  !
5264  ! -- recalculate the residual
5265  call this%lak_calculate_residual(n, hlak, resid)
5266  !
5267  ! -- evaluate convergence
5268  !if (ABS(dh) < delh) then
5269  if (abs(dh) < delh .and. abs(resid) * qtolfact < this%dmaxchg) then
5270  this%ncncvr(n) = 1
5271  end if
5272  this%xnewpak(n) = hlak
5273  !
5274  ! -- save iterates for lake
5275  this%seep0(n) = this%seep(n)
5276  this%dh0(n) = dh
5277  end if
5278  end do laklevel
5279  !
5280  if (iicnvg == 1) exit converge
5281  !
5282  end do converge
5283  !
5284  ! -- Mover terms: store outflow after diversion loss
5285  ! as qformvr and reduce outflow (qd)
5286  ! by how much was actually sent to the mover
5287  if (this%imover == 1) then
5288  do n = 1, this%noutlets
5289  call this%pakmvrobj%accumulate_qformvr(n, -this%simoutrate(n))
5290  end do
5291  end if
5292  end subroutine lak_solve
5293 
5294  !> @ brief Lake package bisection method
5295  !!
5296  !! Use bisection method to find lake stage that reduces the residual
5297  !<
5298  subroutine lak_bisection(this, n, ibflg, hlak, temporary_stage, dh, residual)
5299  ! -- dummy
5300  class(laktype), intent(inout) :: this
5301  integer(I4B), intent(in) :: n !< lake number
5302  integer(I4B), intent(inout) :: ibflg !< bisection flag
5303  real(DP), intent(in) :: hlak !< lake stage
5304  real(DP), intent(inout) :: temporary_stage !< temporary lake stage
5305  real(DP), intent(inout) :: dh !< lake stage change
5306  real(DP), intent(inout) :: residual !< lake residual
5307  ! -- local
5308  integer(I4B) :: i
5309  real(DP) :: temporary_stage0
5310  real(DP) :: residuala
5311  real(DP) :: endpoint1
5312  real(DP) :: endpoint2
5313  ! -- code
5314  ibflg = 1
5315  temporary_stage0 = hlak
5316  endpoint1 = this%en1(n)
5317  endpoint2 = this%en2(n)
5318  call this%lak_calculate_residual(n, temporary_stage, residuala)
5319  if (hlak > endpoint1 .and. hlak < endpoint2) then
5320  endpoint2 = hlak
5321  end if
5322  do i = 1, this%maxlakit
5323  temporary_stage = dhalf * (endpoint1 + endpoint2)
5324  call this%lak_calculate_residual(n, temporary_stage, residual)
5325  if (abs(residual) == dzero .or. &
5326  abs(temporary_stage0 - temporary_stage) < this%dmaxchg) then
5327  exit
5328  end if
5329  call this%lak_calculate_residual(n, endpoint1, residuala)
5330  ! -- change end points
5331  ! -- root is between temporary_stage and endpoint2
5332  if (sign(done, residuala) == sign(done, residual)) then
5333  endpoint1 = temporary_stage
5334  ! -- root is between endpoint1 and temporary_stage
5335  else
5336  endpoint2 = temporary_stage
5337  end if
5338  temporary_stage0 = temporary_stage
5339  end do
5340  dh = hlak - temporary_stage
5341  end subroutine lak_bisection
5342 
5343  !> @brief Calculate the available volumetric rate for a lake given a passed
5344  !! stage
5345  !<
5346  subroutine lak_calculate_available(this, n, hlak, avail, &
5347  ra, ro, qinf, ex, headp)
5348  ! -- modules
5349  use tdismodule, only: delt
5350  ! -- dummy
5351  class(laktype), intent(inout) :: this
5352  integer(I4B), intent(in) :: n
5353  real(DP), intent(in) :: hlak
5354  real(DP), intent(inout) :: avail
5355  real(DP), intent(inout) :: ra
5356  real(DP), intent(inout) :: ro
5357  real(DP), intent(inout) :: qinf
5358  real(DP), intent(inout) :: ex
5359  real(DP), intent(in), optional :: headp
5360  ! -- local
5361  integer(I4B) :: j
5362  integer(I4B) :: idry
5363  integer(I4B) :: igwfnode
5364  real(DP) :: hp
5365  real(DP) :: head
5366  real(DP) :: qlakgw
5367  real(DP) :: v0
5368  !
5369  ! -- set hp
5370  if (present(headp)) then
5371  hp = headp
5372  else
5373  hp = dzero
5374  end if
5375  !
5376  ! -- initialize
5377  avail = dzero
5378  !
5379  ! -- calculate the aquifer sources to the lake
5380  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5381  igwfnode = this%cellid(j)
5382  if (this%ibound(igwfnode) == 0) cycle
5383  head = this%xnew(igwfnode) + hp
5384  call this%lak_estimate_conn_exchange(1, n, j, idry, hlak, head, qlakgw, &
5385  avail)
5386  end do
5387  !
5388  ! -- add rainfall
5389  call this%lak_calculate_rainfall(n, hlak, ra)
5390  avail = avail + ra
5391  !
5392  ! -- calculate runoff
5393  call this%lak_calculate_runoff(n, ro)
5394  avail = avail + ro
5395  !
5396  ! -- calculate inflow
5397  call this%lak_calculate_inflow(n, qinf)
5398  avail = avail + qinf
5399  !
5400  ! -- calculate external flow terms
5401  call this%lak_calculate_external(n, ex)
5402  avail = avail + ex
5403  !
5404  ! -- calculate volume available in storage
5405  call this%lak_calculate_vol(n, this%xoldpak(n), v0)
5406  avail = avail + v0 / delt
5407  end subroutine lak_calculate_available
5408 
5409  !> @brief Calculate the residual for a lake given a passed stage
5410  !<
5411  subroutine lak_calculate_residual(this, n, hlak, resid, headp)
5412  ! -- modules
5413  use tdismodule, only: delt
5414  ! -- dummy
5415  class(laktype), intent(inout) :: this
5416  integer(I4B), intent(in) :: n
5417  real(DP), intent(in) :: hlak
5418  real(DP), intent(inout) :: resid
5419  real(DP), intent(in), optional :: headp
5420  ! -- local
5421  integer(I4B) :: j
5422  integer(I4B) :: idry
5423  integer(I4B) :: igwfnode
5424  real(DP) :: hp
5425  real(DP) :: avail
5426  real(DP) :: head
5427  real(DP) :: ra
5428  real(DP) :: ro
5429  real(DP) :: qinf
5430  real(DP) :: ex
5431  real(DP) :: ev
5432  real(DP) :: wr
5433  real(DP) :: sout
5434  real(DP) :: sin
5435  real(DP) :: qlakgw
5436  real(DP) :: seep
5437  real(DP) :: hlak0
5438  real(DP) :: v0
5439  real(DP) :: v1
5440  !
5441  ! -- set hp
5442  if (present(headp)) then
5443  hp = headp
5444  else
5445  hp = dzero
5446  end if
5447  !
5448  ! -- initialize
5449  resid = dzero
5450  avail = dzero
5451  seep = dzero
5452  !
5453  ! -- calculate the available water
5454  call this%lak_calculate_available(n, hlak, avail, &
5455  ra, ro, qinf, ex, hp)
5456  !
5457  ! -- calculate groundwater seepage
5458  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5459  igwfnode = this%cellid(j)
5460  if (this%ibound(igwfnode) == 0) cycle
5461  head = this%xnew(igwfnode) + hp
5462  call this%lak_estimate_conn_exchange(2, n, j, idry, hlak, head, qlakgw, &
5463  avail)
5464  seep = seep + qlakgw
5465  end do
5466  !
5467  ! -- limit withdrawals to lake inflows and lake storage
5468  call this%lak_calculate_withdrawal(n, avail, wr)
5469  !
5470  ! -- limit evaporation to lake inflows and lake storage
5471  call this%lak_calculate_evaporation(n, hlak, avail, ev)
5472  !
5473  ! -- no outlet flow if evaporation consumes all water
5474  call this%lak_calculate_outlet_outflow(n, hlak, avail, sout)
5475  !
5476  ! -- update the surface inflow values
5477  call this%lak_calculate_outlet_inflow(n, sin)
5478  !
5479  ! -- calculate residual
5480  resid = ra + ev + wr + ro + qinf + ex + sin + sout + seep
5481  !
5482  ! -- include storage
5483  if (this%gwfiss /= 1) then
5484  hlak0 = this%xoldpak(n)
5485  call this%lak_calculate_vol(n, hlak0, v0)
5486  call this%lak_calculate_vol(n, hlak, v1)
5487  resid = resid + (v0 - v1) / delt
5488  end if
5489  end subroutine lak_calculate_residual
5490 
5491  !> @brief Set up the budget object that stores all the lake flows
5492  !<
5493  subroutine lak_setup_budobj(this)
5494  ! -- modules
5495  use constantsmodule, only: lenbudtxt
5496  ! -- dummy
5497  class(laktype) :: this
5498  ! -- local
5499  integer(I4B) :: nbudterm
5500  integer(I4B) :: nlen
5501  integer(I4B) :: j, n, n1, n2
5502  integer(I4B) :: maxlist, naux
5503  integer(I4B) :: idx
5504  real(DP) :: q
5505  character(len=LENBUDTXT) :: text
5506  character(len=LENBUDTXT), dimension(1) :: auxtxt
5507  !
5508  ! -- Determine the number of lake budget terms. These are fixed for
5509  ! the simulation and cannot change
5510  nbudterm = 9
5511  nlen = 0
5512  do n = 1, this%noutlets
5513  if (this%lakein(n) > 0 .and. this%lakeout(n) > 0) then
5514  nlen = nlen + 1
5515  end if
5516  end do
5517  if (nlen > 0) nbudterm = nbudterm + 1
5518  if (this%imover == 1) nbudterm = nbudterm + 2
5519  if (this%naux > 0) nbudterm = nbudterm + 1
5520  !
5521  ! -- set up budobj
5522  call budgetobject_cr(this%budobj, this%packName)
5523  call this%budobj%budgetobject_df(this%nlakes, nbudterm, 0, 0, &
5524  ibudcsv=this%ibudcsv)
5525  idx = 0
5526  !
5527  ! -- Go through and set up each budget term. nlen is the number
5528  ! of outlets that discharge into another lake
5529  if (nlen > 0) then
5530  text = ' FLOW-JA-FACE'
5531  idx = idx + 1
5532  maxlist = 2 * nlen
5533  naux = 0
5534  call this%budobj%budterm(idx)%initialize(text, &
5535  this%name_model, &
5536  this%packName, &
5537  this%name_model, &
5538  this%packName, &
5539  maxlist, .false., .false., &
5540  naux, ordered_id1=.false.)
5541  !
5542  ! -- store connectivity
5543  call this%budobj%budterm(idx)%reset(2 * nlen)
5544  q = dzero
5545  do n = 1, this%noutlets
5546  n1 = this%lakein(n)
5547  n2 = this%lakeout(n)
5548  if (n1 > 0 .and. n2 > 0) then
5549  call this%budobj%budterm(idx)%update_term(n1, n2, q)
5550  call this%budobj%budterm(idx)%update_term(n2, n1, -q)
5551  end if
5552  end do
5553  end if
5554  !
5555  ! --
5556  text = ' GWF'
5557  idx = idx + 1
5558  maxlist = this%maxbound
5559  naux = 1
5560  auxtxt(1) = ' FLOW-AREA'
5561  call this%budobj%budterm(idx)%initialize(text, &
5562  this%name_model, &
5563  this%packName, &
5564  this%name_model, &
5565  this%name_model, &
5566  maxlist, .false., .true., &
5567  naux, auxtxt)
5568  call this%budobj%budterm(idx)%reset(this%maxbound)
5569  q = dzero
5570  do n = 1, this%nlakes
5571  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5572  n2 = this%cellid(j)
5573  call this%budobj%budterm(idx)%update_term(n, n2, q)
5574  end do
5575  end do
5576  !
5577  ! --
5578  text = ' RAINFALL'
5579  idx = idx + 1
5580  maxlist = this%nlakes
5581  naux = 0
5582  call this%budobj%budterm(idx)%initialize(text, &
5583  this%name_model, &
5584  this%packName, &
5585  this%name_model, &
5586  this%packName, &
5587  maxlist, .false., .false., &
5588  naux)
5589  !
5590  ! --
5591  text = ' EVAPORATION'
5592  idx = idx + 1
5593  maxlist = this%nlakes
5594  naux = 0
5595  call this%budobj%budterm(idx)%initialize(text, &
5596  this%name_model, &
5597  this%packName, &
5598  this%name_model, &
5599  this%packName, &
5600  maxlist, .false., .false., &
5601  naux)
5602  !
5603  ! --
5604  text = ' RUNOFF'
5605  idx = idx + 1
5606  maxlist = this%nlakes
5607  naux = 0
5608  call this%budobj%budterm(idx)%initialize(text, &
5609  this%name_model, &
5610  this%packName, &
5611  this%name_model, &
5612  this%packName, &
5613  maxlist, .false., .false., &
5614  naux)
5615  !
5616  ! --
5617  text = ' EXT-INFLOW'
5618  idx = idx + 1
5619  maxlist = this%nlakes
5620  naux = 0
5621  call this%budobj%budterm(idx)%initialize(text, &
5622  this%name_model, &
5623  this%packName, &
5624  this%name_model, &
5625  this%packName, &
5626  maxlist, .false., .false., &
5627  naux)
5628  !
5629  ! --
5630  text = ' WITHDRAWAL'
5631  idx = idx + 1
5632  maxlist = this%nlakes
5633  naux = 0
5634  call this%budobj%budterm(idx)%initialize(text, &
5635  this%name_model, &
5636  this%packName, &
5637  this%name_model, &
5638  this%packName, &
5639  maxlist, .false., .false., &
5640  naux)
5641  !
5642  ! --
5643  text = ' EXT-OUTFLOW'
5644  idx = idx + 1
5645  maxlist = this%nlakes
5646  naux = 0
5647  call this%budobj%budterm(idx)%initialize(text, &
5648  this%name_model, &
5649  this%packName, &
5650  this%name_model, &
5651  this%packName, &
5652  maxlist, .false., .false., &
5653  naux)
5654  !
5655  ! --
5656  text = ' STORAGE'
5657  idx = idx + 1
5658  maxlist = this%nlakes
5659  naux = 1
5660  auxtxt(1) = ' VOLUME'
5661  call this%budobj%budterm(idx)%initialize(text, &
5662  this%name_model, &
5663  this%packName, &
5664  this%name_model, &
5665  this%packName, &
5666  maxlist, .false., .false., &
5667  naux, auxtxt)
5668  !
5669  ! --
5670  text = ' CONSTANT'
5671  idx = idx + 1
5672  maxlist = this%nlakes
5673  naux = 0
5674  call this%budobj%budterm(idx)%initialize(text, &
5675  this%name_model, &
5676  this%packName, &
5677  this%name_model, &
5678  this%packName, &
5679  maxlist, .false., .false., &
5680  naux)
5681  !
5682  ! --
5683  if (this%imover == 1) then
5684  !
5685  ! --
5686  text = ' FROM-MVR'
5687  idx = idx + 1
5688  maxlist = this%nlakes
5689  naux = 0
5690  call this%budobj%budterm(idx)%initialize(text, &
5691  this%name_model, &
5692  this%packName, &
5693  this%name_model, &
5694  this%packName, &
5695  maxlist, .false., .false., &
5696  naux)
5697  !
5698  ! --
5699  text = ' TO-MVR'
5700  idx = idx + 1
5701  maxlist = this%noutlets
5702  naux = 0
5703  call this%budobj%budterm(idx)%initialize(text, &
5704  this%name_model, &
5705  this%packName, &
5706  this%name_model, &
5707  this%packName, &
5708  maxlist, .false., .false., &
5709  naux, ordered_id1=.false.)
5710  !
5711  ! -- store to-mvr connection information
5712  call this%budobj%budterm(idx)%reset(this%noutlets)
5713  q = dzero
5714  do n = 1, this%noutlets
5715  n1 = this%lakein(n)
5716  call this%budobj%budterm(idx)%update_term(n1, n1, q)
5717  end do
5718  end if
5719  !
5720  ! --
5721  naux = this%naux
5722  if (naux > 0) then
5723  !
5724  ! --
5725  text = ' AUXILIARY'
5726  idx = idx + 1
5727  maxlist = this%nlakes
5728  call this%budobj%budterm(idx)%initialize(text, &
5729  this%name_model, &
5730  this%packName, &
5731  this%name_model, &
5732  this%packName, &
5733  maxlist, .false., .false., &
5734  naux, this%auxname)
5735  end if
5736  !
5737  ! -- if lake flow for each reach are written to the listing file
5738  if (this%iprflow /= 0) then
5739  call this%budobj%flowtable_df(this%iout)
5740  end if
5741  end subroutine lak_setup_budobj
5742 
5743  !> @brief Copy flow terms into this%budobj
5744  !<
5745  subroutine lak_fill_budobj(this)
5746  ! -- dummy
5747  class(laktype) :: this
5748  ! -- local
5749  integer(I4B) :: naux
5750  real(DP), dimension(:), allocatable :: auxvartmp
5751  !integer(I4B) :: i
5752  integer(I4B) :: j
5753  integer(I4B) :: n
5754  integer(I4B) :: n1
5755  integer(I4B) :: n2
5756  integer(I4B) :: ii
5757  integer(I4B) :: jj
5758  integer(I4B) :: idx
5759  integer(I4B) :: nlen
5760  real(DP) :: v, v1
5761  real(DP) :: q
5762  real(DP) :: lkstg, gwhead, wa
5763  !
5764  ! -- initialize counter
5765  idx = 0
5766 
5767  ! -- FLOW JA FACE
5768  nlen = 0
5769  do n = 1, this%noutlets
5770  if (this%lakein(n) > 0 .and. this%lakeout(n) > 0) then
5771  nlen = nlen + 1
5772  end if
5773  end do
5774  if (nlen > 0) then
5775  idx = idx + 1
5776  call this%budobj%budterm(idx)%reset(2 * nlen)
5777  do n = 1, this%noutlets
5778  n1 = this%lakein(n)
5779  n2 = this%lakeout(n)
5780  if (n1 > 0 .and. n2 > 0) then
5781  q = this%simoutrate(n)
5782  if (this%imover == 1) then
5783  q = q + this%pakmvrobj%get_qtomvr(n)
5784  end if
5785  call this%budobj%budterm(idx)%update_term(n1, n2, q)
5786  call this%budobj%budterm(idx)%update_term(n2, n1, -q)
5787  end if
5788  end do
5789  end if
5790  !
5791  ! -- GWF (LEAKAGE)
5792  idx = idx + 1
5793  call this%budobj%budterm(idx)%reset(this%maxbound)
5794  do n = 1, this%nlakes
5795  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
5796  n2 = this%cellid(j)
5797  q = this%qleak(j)
5798  lkstg = this%xnewpak(n)
5799  ! -- For the case when the lak stage is exactly equal
5800  ! to the lake bottom, the wetted area is not returned
5801  ! equal to 0.0
5802  gwhead = this%xnew(n2)
5803  call this%lak_calculate_conn_warea(n, j, lkstg, gwhead, wa)
5804  ! -- For thermal conduction between a lake and a gw cell,
5805  ! the shared wetted area should be reset to zero when the lake
5806  ! stage is below the cell bottom
5807  if (this%belev(j) > lkstg) wa = dzero
5808  this%qauxcbc(1) = wa
5809  call this%budobj%budterm(idx)%update_term(n, n2, q, this%qauxcbc)
5810  end do
5811  end do
5812  !
5813  ! -- RAIN
5814  idx = idx + 1
5815  call this%budobj%budterm(idx)%reset(this%nlakes)
5816  do n = 1, this%nlakes
5817  q = this%precip(n)
5818  call this%budobj%budterm(idx)%update_term(n, n, q)
5819  end do
5820  !
5821  ! -- EVAPORATION
5822  idx = idx + 1
5823  call this%budobj%budterm(idx)%reset(this%nlakes)
5824  do n = 1, this%nlakes
5825  q = this%evap(n)
5826  call this%budobj%budterm(idx)%update_term(n, n, q)
5827  end do
5828  !
5829  ! -- RUNOFF
5830  idx = idx + 1
5831  call this%budobj%budterm(idx)%reset(this%nlakes)
5832  do n = 1, this%nlakes
5833  q = this%runoff(n)
5834  call this%budobj%budterm(idx)%update_term(n, n, q)
5835  end do
5836  !
5837  ! -- INFLOW
5838  idx = idx + 1
5839  call this%budobj%budterm(idx)%reset(this%nlakes)
5840  do n = 1, this%nlakes
5841  q = this%inflow(n)
5842  call this%budobj%budterm(idx)%update_term(n, n, q)
5843  end do
5844  !
5845  ! -- WITHDRAWAL
5846  idx = idx + 1
5847  call this%budobj%budterm(idx)%reset(this%nlakes)
5848  do n = 1, this%nlakes
5849  q = this%withr(n)
5850  call this%budobj%budterm(idx)%update_term(n, n, q)
5851  end do
5852  !
5853  ! -- EXTERNAL OUTFLOW
5854  idx = idx + 1
5855  call this%budobj%budterm(idx)%reset(this%nlakes)
5856  do n = 1, this%nlakes
5857  call this%lak_get_external_outlet(n, q)
5858  ! subtract tomover from external outflow
5859  call this%lak_get_external_mover(n, v)
5860  q = q + v
5861  call this%budobj%budterm(idx)%update_term(n, n, q)
5862  end do
5863  !
5864  ! -- STORAGE
5865  idx = idx + 1
5866  call this%budobj%budterm(idx)%reset(this%nlakes)
5867  do n = 1, this%nlakes
5868  call this%lak_calculate_vol(n, this%xnewpak(n), v1)
5869  q = this%qsto(n)
5870  this%qauxcbc(1) = v1
5871  call this%budobj%budterm(idx)%update_term(n, n, q, this%qauxcbc)
5872  end do
5873  !
5874  ! -- CONSTANT FLOW
5875  idx = idx + 1
5876  call this%budobj%budterm(idx)%reset(this%nlakes)
5877  do n = 1, this%nlakes
5878  q = this%chterm(n)
5879  call this%budobj%budterm(idx)%update_term(n, n, q)
5880  end do
5881  !
5882  ! -- MOVER
5883  if (this%imover == 1) then
5884  !
5885  ! -- FROM MOVER
5886  idx = idx + 1
5887  call this%budobj%budterm(idx)%reset(this%nlakes)
5888  do n = 1, this%nlakes
5889  q = this%pakmvrobj%get_qfrommvr(n)
5890  call this%budobj%budterm(idx)%update_term(n, n, q)
5891  end do
5892  !
5893  ! -- TO MOVER
5894  idx = idx + 1
5895  call this%budobj%budterm(idx)%reset(this%noutlets)
5896  do n = 1, this%noutlets
5897  n1 = this%lakein(n)
5898  q = this%pakmvrobj%get_qtomvr(n)
5899  if (q > dzero) then
5900  q = -q
5901  end if
5902  call this%budobj%budterm(idx)%update_term(n1, n1, q)
5903  end do
5904  !
5905  end if
5906  !
5907  ! -- AUXILIARY VARIABLES
5908  naux = this%naux
5909  if (naux > 0) then
5910  idx = idx + 1
5911  allocate (auxvartmp(naux))
5912  call this%budobj%budterm(idx)%reset(this%nlakes)
5913  do n = 1, this%nlakes
5914  q = dzero
5915  do jj = 1, naux
5916  ii = n
5917  auxvartmp(jj) = this%lauxvar(jj, ii)
5918  end do
5919  call this%budobj%budterm(idx)%update_term(n, n, q, auxvartmp)
5920  end do
5921  deallocate (auxvartmp)
5922  end if
5923  !
5924  ! --Terms are filled, now accumulate them for this time step
5925  call this%budobj%accumulate_terms()
5926  end subroutine lak_fill_budobj
5927 
5928  !> @brief Set up the table object that is used to write the lak stage data
5929  !!
5930  !! The terms listed here must correspond in number and order to the ones
5931  !! written to the stage table in the lak_ot method
5932  !<
5933  subroutine lak_setup_tableobj(this)
5934  ! -- modules
5935  use constantsmodule, only: linelength, lenbudtxt
5936  ! -- dummy
5937  class(laktype) :: this
5938  ! -- local
5939  integer(I4B) :: nterms
5940  character(len=LINELENGTH) :: title
5941  character(len=LINELENGTH) :: text
5942  !
5943  ! -- setup stage table
5944  if (this%iprhed > 0) then
5945  !
5946  ! -- Determine the number of lake stage terms. These are fixed for
5947  ! the simulation and cannot change. This includes FLOW-JA-FACE
5948  ! so they can be written to the binary budget files, but these internal
5949  ! flows are not included as part of the budget table.
5950  nterms = 5
5951  if (this%inamedbound == 1) then
5952  nterms = nterms + 1
5953  end if
5954  !
5955  ! -- set up table title
5956  title = trim(adjustl(this%text))//' PACKAGE ('// &
5957  trim(adjustl(this%packName))//') STAGES FOR EACH CONTROL VOLUME'
5958  !
5959  ! -- set up stage tableobj
5960  call table_cr(this%stagetab, this%packName, title)
5961  call this%stagetab%table_df(this%nlakes, nterms, this%iout, &
5962  transient=.true.)
5963  !
5964  ! -- Go through and set up table budget term
5965  if (this%inamedbound == 1) then
5966  text = 'NAME'
5967  call this%stagetab%initialize_column(text, 20, alignment=tableft)
5968  end if
5969  !
5970  ! -- lake number
5971  text = 'NUMBER'
5972  call this%stagetab%initialize_column(text, 10, alignment=tabcenter)
5973  !
5974  ! -- lake stage
5975  text = 'STAGE'
5976  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5977  !
5978  ! -- lake surface area
5979  text = 'SURFACE AREA'
5980  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5981  !
5982  ! -- lake wetted area
5983  text = 'WETTED AREA'
5984  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5985  !
5986  ! -- lake volume
5987  text = 'VOLUME'
5988  call this%stagetab%initialize_column(text, 12, alignment=tabcenter)
5989  end if
5990  end subroutine lak_setup_tableobj
5991 
5992  !> @brief Activate addition of density terms
5993  !<
5994  subroutine lak_activate_density(this)
5995  ! -- dummy
5996  class(laktype), intent(inout) :: this
5997  ! -- local
5998  integer(I4B) :: i, j
5999  !
6000  ! -- Set idense and reallocate denseterms to be of size MAXBOUND
6001  this%idense = 1
6002  call mem_reallocate(this%denseterms, 3, this%MAXBOUND, 'DENSETERMS', &
6003  this%memoryPath)
6004  do i = 1, this%maxbound
6005  do j = 1, 3
6006  this%denseterms(j, i) = dzero
6007  end do
6008  end do
6009  write (this%iout, '(/1x,a)') 'DENSITY TERMS HAVE BEEN ACTIVATED FOR LAKE &
6010  &PACKAGE: '//trim(adjustl(this%packName))
6011  end subroutine lak_activate_density
6012 
6013  !> @brief Activate viscosity terms
6014  !!
6015  !! Method to activate addition of viscosity terms for a LAK package reach.
6016  !<
6017  subroutine lak_activate_viscosity(this)
6018  ! -- modules
6019  use memorymanagermodule, only: mem_reallocate
6020  ! -- dummy variables
6021  class(laktype), intent(inout) :: this !< LakType object
6022  ! -- local variables
6023  integer(I4B) :: i
6024  integer(I4B) :: j
6025  !
6026  ! -- Set ivsc and reallocate viscratios to be of size MAXBOUND
6027  this%ivsc = 1
6028  call mem_reallocate(this%viscratios, 2, this%MAXBOUND, 'VISCRATIOS', &
6029  this%memoryPath)
6030  do i = 1, this%maxbound
6031  do j = 1, 2
6032  this%viscratios(j, i) = done
6033  end do
6034  end do
6035  write (this%iout, '(/1x,a)') 'VISCOSITY HAS BEEN ACTIVATED FOR LAK &
6036  &PACKAGE: '//trim(adjustl(this%packName))
6037  end subroutine lak_activate_viscosity
6038 
6039  !> @brief Calculate the groundwater-lake density exchange terms
6040  !!
6041  !! Arguments are as follows:
6042  !! iconn : lak-gwf connection number
6043  !! stage : lake stage
6044  !! head : gwf head
6045  !! cond : conductance
6046  !! botl : bottom elevation of this connection
6047  !! flow : calculated flow, updated here with density terms
6048  !! gwfhcof : gwf head coefficient, updated here with density terms
6049  !! gwfrhs : gwf right-hand-side value, updated here with density terms
6050  !!
6051  !! Member variable used here
6052  !! denseterms : shape (3, MAXBOUND), filled by buoyancy package
6053  !! col 1 is relative density of lake (denselak / denseref)
6054  !! col 2 is relative density of gwf cell (densegwf / denseref)
6055  !! col 3 is elevation of gwf cell
6056  !<
6057  subroutine lak_calculate_density_exchange(this, iconn, stage, head, cond, &
6058  botl, flow, gwfhcof, gwfrhs)
6059  ! -- dummy
6060  class(laktype), intent(inout) :: this
6061  integer(I4B), intent(in) :: iconn
6062  real(DP), intent(in) :: stage
6063  real(DP), intent(in) :: head
6064  real(DP), intent(in) :: cond
6065  real(DP), intent(in) :: botl
6066  real(DP), intent(inout) :: flow
6067  real(DP), intent(inout) :: gwfhcof
6068  real(DP), intent(inout) :: gwfrhs
6069  ! -- local
6070  real(DP) :: ss
6071  real(DP) :: hh
6072  real(DP) :: havg
6073  real(DP) :: rdenselak
6074  real(DP) :: rdensegwf
6075  real(DP) :: rdenseavg
6076  real(DP) :: elevlak
6077  real(DP) :: elevgwf
6078  real(DP) :: elevavg
6079  real(DP) :: d1
6080  real(DP) :: d2
6081  logical(LGP) :: stage_below_bot
6082  logical(LGP) :: head_below_bot
6083  !
6084  ! -- Set lak density to lak density or gwf density
6085  if (stage >= botl) then
6086  ss = stage
6087  stage_below_bot = .false.
6088  rdenselak = this%denseterms(1, iconn) ! lak rel density
6089  else
6090  ss = botl
6091  stage_below_bot = .true.
6092  rdenselak = this%denseterms(2, iconn) ! gwf rel density
6093  end if
6094  !
6095  ! -- set hh to head or botl
6096  if (head >= botl) then
6097  hh = head
6098  head_below_bot = .false.
6099  rdensegwf = this%denseterms(2, iconn) ! gwf rel density
6100  else
6101  hh = botl
6102  head_below_bot = .true.
6103  rdensegwf = this%denseterms(1, iconn) ! lak rel density
6104  end if
6105  !
6106  ! -- todo: hack because denseterms not updated in a cf calculation
6107  if (rdensegwf == dzero) return
6108  !
6109  ! -- Update flow
6110  if (stage_below_bot .and. head_below_bot) then
6111  !
6112  ! -- flow is zero, so no terms are updated
6113  !
6114  else
6115  !
6116  ! -- calculate average relative density
6117  rdenseavg = dhalf * (rdenselak + rdensegwf)
6118  !
6119  ! -- Add contribution of first density term:
6120  ! cond * (denseavg/denseref - 1) * (hgwf - hlak)
6121  d1 = cond * (rdenseavg - done)
6122  gwfhcof = gwfhcof - d1
6123  gwfrhs = gwfrhs - d1 * ss
6124  d1 = d1 * (hh - ss)
6125  flow = flow + d1
6126  !
6127  ! -- Add second density term if stage and head not below bottom
6128  if (.not. stage_below_bot .and. .not. head_below_bot) then
6129  !
6130  ! -- Add contribution of second density term:
6131  ! cond * (havg - elevavg) * (densegwf - denselak) / denseref
6132  elevgwf = this%denseterms(3, iconn)
6133  if (this%ictype(iconn) == 0 .or. this%ictype(iconn) == 3) then
6134  ! -- vertical or embedded vertical connection
6135  elevlak = botl
6136  else
6137  ! -- horizontal or embedded horizontal connection
6138  elevlak = elevgwf
6139  end if
6140  elevavg = dhalf * (elevlak + elevgwf)
6141  havg = dhalf * (hh + ss)
6142  d2 = cond * (havg - elevavg) * (rdensegwf - rdenselak)
6143  gwfrhs = gwfrhs + d2
6144  flow = flow + d2
6145  end if
6146  end if
6147  end subroutine lak_calculate_density_exchange
6148 
6149 end module lakmodule
memorymanagermodule::mem_allocate
Definition: MemoryManager.f90:54
memorymanagermodule::mem_deallocate
Definition: MemoryManager.f90:122
memorymanagermodule::mem_reallocate
Definition: MemoryManager.f90:79
memorymanagermodule::mem_setptr
Definition: MemoryManager.f90:89
blockparsermodule
This module contains block parser methods.
Definition: BlockParser.f90:7
bndmodule
This module contains the base boundary package.
Definition: BoundaryPackage.f90:8
budgetobjectmodule::budgetobject_cr
subroutine, public budgetobject_cr(this, name)
Create a new budget object.
Definition: BudgetObject.f90:83
constantsmodule
This module contains simulation constants.
Definition: Constants.f90:9
constantsmodule::linelength
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
constantsmodule::dhdry
real(dp), parameter dhdry
real dry cell constant
Definition: Constants.f90:94
constantsmodule::tabcenter
@ tabcenter
centered table column
Definition: Constants.f90:172
constantsmodule::tabright
@ tabright
right justified table column
Definition: Constants.f90:173
constantsmodule::tableft
@ tableft
left justified table column
Definition: Constants.f90:171
constantsmodule::mnormal
@ mnormal
normal output mode
Definition: Constants.f90:206
constantsmodule::dtwothirds
real(dp), parameter dtwothirds
real constant 2/3
Definition: Constants.f90:70
constantsmodule::tabucstring
@ tabucstring
upper case string table data
Definition: Constants.f90:180
constantsmodule::tabstring
@ tabstring
string table data
Definition: Constants.f90:179
constantsmodule::tabreal
@ tabreal
real table data
Definition: Constants.f90:182
constantsmodule::tabinteger
@ tabinteger
integer table data
Definition: Constants.f90:181
constantsmodule::lenpackagename
integer(i4b), parameter lenpackagename
maximum length of the package name
Definition: Constants.f90:23
constantsmodule::iwetlake
integer(i4b), parameter iwetlake
integer constant for a dry lake
Definition: Constants.f90:52
constantsmodule::deight
real(dp), parameter deight
real constant 8
Definition: Constants.f90:83
constantsmodule::dfivethirds
real(dp), parameter dfivethirds
real constant 5/3
Definition: Constants.f90:78
constantsmodule::dp999
real(dp), parameter dp999
real constant 999/1000
Definition: Constants.f90:74
constantsmodule::namedboundflag
integer(i4b), parameter namedboundflag
named bound flag
Definition: Constants.f90:49
constantsmodule::donethird
real(dp), parameter donethird
real constant 1/3
Definition: Constants.f90:67
constantsmodule::dnodata
real(dp), parameter dnodata
real no data constant
Definition: Constants.f90:95
constantsmodule::dhnoflo
real(dp), parameter dhnoflo
real no flow constant
Definition: Constants.f90:93
constantsmodule::dhundred
real(dp), parameter dhundred
real constant 100
Definition: Constants.f90:86
constantsmodule::lenpakloc
integer(i4b), parameter lenpakloc
maximum length of a package location
Definition: Constants.f90:50
constantsmodule::lentimeseriesname
integer(i4b), parameter lentimeseriesname
maximum length of a time series name
Definition: Constants.f90:42
constantsmodule::dep20
real(dp), parameter dep20
real constant 1e20
Definition: Constants.f90:91
constantsmodule::dem1
real(dp), parameter dem1
real constant 1e-1
Definition: Constants.f90:103
constantsmodule::maxadpit
integer(i4b), parameter maxadpit
maximum advanced package Newton-Raphson iterations
Definition: Constants.f90:53
constantsmodule::lenvarname
integer(i4b), parameter lenvarname
maximum length of a variable name
Definition: Constants.f90:17
constantsmodule::dhalf
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
constantsmodule::lenftype
integer(i4b), parameter lenftype
maximum length of a package type (DIS, WEL, OC, etc.)
Definition: Constants.f90:39
constantsmodule::dgravity
real(dp), parameter dgravity
real constant gravitational acceleration (m/(s s))
Definition: Constants.f90:132
constantsmodule::dpi
real(dp), parameter dpi
real constant
Definition: Constants.f90:128
constantsmodule::lenboundname
integer(i4b), parameter lenboundname
maximum length of a bound name
Definition: Constants.f90:36
constantsmodule::dem4
real(dp), parameter dem4
real constant 1e-4
Definition: Constants.f90:107
constantsmodule::dem30
real(dp), parameter dem30
real constant 1e-30
Definition: Constants.f90:118
constantsmodule::dem6
real(dp), parameter dem6
real constant 1e-6
Definition: Constants.f90:109
constantsmodule::dzero
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
constantsmodule::dcd
real(dp), parameter dcd
real constant weir coefficient in SI units
Definition: Constants.f90:133
constantsmodule::dem5
real(dp), parameter dem5
real constant 1e-5
Definition: Constants.f90:108
constantsmodule::dten
real(dp), parameter dten
real constant 10
Definition: Constants.f90:84
constantsmodule::dprec
real(dp), parameter dprec
real constant machine precision
Definition: Constants.f90:120
constantsmodule::maxcharlen
integer(i4b), parameter maxcharlen
maximum length of char string
Definition: Constants.f90:47
constantsmodule::dp7
real(dp), parameter dp7
real constant 7/10
Definition: Constants.f90:71
constantsmodule::dem9
real(dp), parameter dem9
real constant 1e-9
Definition: Constants.f90:112
constantsmodule::dem2
real(dp), parameter dem2
real constant 1e-2
Definition: Constants.f90:105
constantsmodule::dtwo
real(dp), parameter dtwo
real constant 2
Definition: Constants.f90:79
constantsmodule::lenbudtxt
integer(i4b), parameter lenbudtxt
maximum length of a budget component names
Definition: Constants.f90:37
constantsmodule::dthree
real(dp), parameter dthree
real constant 3
Definition: Constants.f90:80
constantsmodule::done
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
geomutilmodule::get_node
integer(i4b) function, public get_node(ilay, irow, icol, nlay, nrow, ncol)
Get node number, given layer, row, and column indices for a structured grid. If any argument is inval...
Definition: GeomUtil.f90:83
inputoutputmodule::extract_idnum_or_bndname
subroutine, public extract_idnum_or_bndname(line, icol, istart, istop, idnum, bndname)
Starting at position icol, define string as line(istart:istop).
Definition: InputOutput.f90:1292
inputoutputmodule::getunit
integer(i4b) function, public getunit()
Get a free unit number.
Definition: InputOutput.f90:170
inputoutputmodule::ulasav
subroutine, public ulasav(buf, text, kstp, kper, pertim, totim, ncol, nrow, ilay, ichn)
Save 1 layer array on disk.
Definition: InputOutput.f90:903
inputoutputmodule::openfile
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
inputoutputmodule::urword
subroutine, public urword(line, icol, istart, istop, ncode, n, r, iout, in)
Extract a word from a string.
Definition: InputOutput.f90:406
kindmodule
This module defines variable data types.
Definition: kind.f90:8
lakmodule::lak_activate_density
subroutine lak_activate_density(this)
Activate addition of density terms.
Definition: gwf-lak.f90:5995
lakmodule::lak_cq
subroutine lak_cq(this, x, flowja, iadv)
Calculate flows.
Definition: gwf-lak.f90:3945
lakmodule::lak_read_outlets
subroutine lak_read_outlets(this)
Read the lake outlets for this package.
Definition: gwf-lak.f90:1402
lakmodule::lak_vol2stage
subroutine lak_vol2stage(this, ilak, vol, stage)
Determine the stage from a provided volume.
Definition: gwf-lak.f90:2736
lakmodule::lak_calculate_outlet_outflow
subroutine lak_calculate_outlet_outflow(this, ilak, stage, avail, outoutf)
Calculate the outlet outflow from a lake.
Definition: gwf-lak.f90:2557
lakmodule::lak_read_lake_connections
subroutine lak_read_lake_connections(this)
Read the lake connections for this package.
Definition: gwf-lak.f90:684
lakmodule::lak_calculate_sarea
subroutine lak_calculate_sarea(this, ilak, stage, sarea)
Calculate the surface area of a lake at a given stage.
Definition: gwf-lak.f90:1983
lakmodule::lak_fn
subroutine lak_fn(this, rhs, ia, idxglo, matrix_sln)
Fill newton terms.
Definition: gwf-lak.f90:3590
lakmodule::lak_read_tables
subroutine lak_read_tables(this)
Read the lake tables for this package.
Definition: gwf-lak.f90:997
lakmodule::lak_set_stressperiod
subroutine lak_set_stressperiod(this, itemno)
Set a stress period attribute for lakweslls(itemno) using keywords.
Definition: gwf-lak.f90:2851
lakmodule::lak_ot_package_flows
subroutine lak_ot_package_flows(this, icbcfl, ibudfl)
Output LAK package flow terms.
Definition: gwf-lak.f90:4047
lakmodule::lak_accumulate_chterm
subroutine lak_accumulate_chterm(this, ilak, rrate, chratin, chratout)
Accumulate constant head terms for budget.
Definition: gwf-lak.f90:4866
lakmodule::lak_process_obsid
subroutine lak_process_obsid(obsrv, dis, inunitobs, iout)
This procedure is pointed to by ObsDataTypeProcesssIdPtr. It processes the ID string of an observatio...
Definition: gwf-lak.f90:4814
lakmodule::lak_get_external_mover
subroutine lak_get_external_mover(this, ilak, outoutf)
Get the mover outflow from a lake to an external boundary.
Definition: gwf-lak.f90:2674
lakmodule::lak_cf
subroutine lak_cf(this)
Formulate the HCOF and RHS terms.
Definition: gwf-lak.f90:3474
lakmodule::lak_calculate_density_exchange
subroutine lak_calculate_density_exchange(this, iconn, stage, head, cond, botl, flow, gwfhcof, gwfrhs)
Calculate the groundwater-lake density exchange terms.
Definition: gwf-lak.f90:6059
lakmodule::lak_calculate_evaporation
subroutine lak_calculate_evaporation(this, ilak, stage, avail, ev)
Calculate the evaporation from a lake at a provided stage subject to an available volume.
Definition: gwf-lak.f90:2509
lakmodule::lak_read_lakes
subroutine lak_read_lakes(this)
Read the dimensions for this package.
Definition: gwf-lak.f90:453
lakmodule::lak_setup_budobj
subroutine lak_setup_budobj(this)
Set up the budget object that stores all the lake flows.
Definition: gwf-lak.f90:5494
lakmodule::lak_calculate_external
subroutine lak_calculate_external(this, ilak, ex)
Calculate the external flow terms to a lake.
Definition: gwf-lak.f90:2471
lakmodule::lak_calculate_warea
subroutine lak_calculate_warea(this, ilak, stage, warea, hin)
Calculate the wetted area of a lake at a given stage.
Definition: gwf-lak.f90:2025
lakmodule::lak_calculate_vol
subroutine lak_calculate_vol(this, ilak, stage, volume)
Calculate the volume of a lake at a given stage.
Definition: gwf-lak.f90:2104
lakmodule::lak_da
subroutine lak_da(this)
Deallocate objects.
Definition: gwf-lak.f90:4172
lakmodule::lak_calculate_storagechange
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.
Definition: gwf-lak.f90:2403
lakmodule::lak_get_internal_outlet
subroutine lak_get_internal_outlet(this, ilak, outoutf)
Get the outlet from a lake to another lake.
Definition: gwf-lak.f90:2636
lakmodule::lak_calculate_conductance
subroutine lak_calculate_conductance(this, ilak, stage, conductance)
Calculate the total conductance for a lake at a provided stage.
Definition: gwf-lak.f90:2157
lakmodule::laktables_to_vectors
subroutine laktables_to_vectors(this, laketables)
Copy the laketables structure data into flattened vectors that are stored in the memory manager.
Definition: gwf-lak.f90:1121
lakmodule::lak_calculate_residual
subroutine lak_calculate_residual(this, n, hlak, resid, headp)
Calculate the residual for a lake given a passed stage.
Definition: gwf-lak.f90:5412
lakmodule::lak_get_outlet_tomover
subroutine lak_get_outlet_tomover(this, ilak, outoutf)
Get the outlet to mover from a lake.
Definition: gwf-lak.f90:2716
lakmodule::lak_allocate_arrays
subroutine lak_allocate_arrays(this)
Allocate scalar members.
Definition: gwf-lak.f90:389
lakmodule::lak_calculate_available
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.
Definition: gwf-lak.f90:5348
lakmodule::lak_obs_supported
logical function lak_obs_supported(this)
Procedures related to observations (type-bound)
Definition: gwf-lak.f90:4392
lakmodule::lak_set_pointers
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.
Definition: gwf-lak.f90:4361
lakmodule::lak_calculate_inflow
subroutine lak_calculate_inflow(this, ilak, qin)
Calculate specified inflow to a lake.
Definition: gwf-lak.f90:2459
lakmodule::text
character(len=lenpackagename) text
Definition: gwf-lak.f90:44
lakmodule::lak_ot_model_flows
subroutine lak_ot_model_flows(this, icbcfl, ibudfl, icbcun, imap)
Write flows to binary file and/or print flows to budget.
Definition: gwf-lak.f90:4073
lakmodule::lak_rp_obs
subroutine lak_rp_obs(this)
Process each observation.
Definition: gwf-lak.f90:4662
lakmodule::lak_calculate_conn_conductance
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.
Definition: gwf-lak.f90:2207
lakmodule::lak_calculate_exchange
subroutine lak_calculate_exchange(this, ilak, stage, totflow)
Calculate the total groundwater-lake flow at a provided stage.
Definition: gwf-lak.f90:2273
lakmodule::lak_linear_interpolation
subroutine lak_linear_interpolation(this, n, x, y, z, v)
Perform linear interpolation of two vectors.
Definition: gwf-lak.f90:1938
lakmodule::lak_setup_tableobj
subroutine lak_setup_tableobj(this)
Set up the table object that is used to write the lak stage data.
Definition: gwf-lak.f90:5934
lakmodule::lak_ot_dv
subroutine lak_ot_dv(this, idvsave, idvprint)
Save LAK-calculated values to binary file.
Definition: gwf-lak.f90:4086
lakmodule::lak_options
subroutine lak_options(this, option, found)
Set options specific to LakType.
Definition: gwf-lak.f90:3092
lakmodule::lak_get_external_outlet
subroutine lak_get_external_outlet(this, ilak, outoutf)
Get the outlet outflow from a lake to an external boundary.
Definition: gwf-lak.f90:2655
lakmodule::lak_calculate_rainfall
subroutine lak_calculate_rainfall(this, ilak, stage, ra)
Calculate the rainfall for a lake.
Definition: gwf-lak.f90:2425
lakmodule::lak_get_internal_inlet
subroutine lak_get_internal_inlet(this, ilak, outinf)
Get the outlet inflow to a lake from another lake.
Definition: gwf-lak.f90:2615
lakmodule::lak_cc
subroutine lak_cc(this, innertot, kiter, iend, icnvgmod, cpak, ipak, dpak)
Final convergence check for package.
Definition: gwf-lak.f90:3649
lakmodule::lak_rp
subroutine lak_rp(this)
Read and Prepare.
Definition: gwf-lak.f90:3271
lakmodule::lak_read_dimensions
subroutine lak_read_dimensions(this)
Read the dimensions for this package.
Definition: gwf-lak.f90:1584
lakmodule::lak_activate_viscosity
subroutine lak_activate_viscosity(this)
Activate viscosity terms.
Definition: gwf-lak.f90:6018
lakmodule::lak_read_initial_attr
subroutine lak_read_initial_attr(this)
Read the initial parameters for this package.
Definition: gwf-lak.f90:1668
lakmodule::lak_check_valid
integer(i4b) function lak_check_valid(this, itemno)
Determine if a valid lake or outlet number has been specified.
Definition: gwf-lak.f90:2817
lakmodule::lak_bisection
subroutine lak_bisection(this, n, ibflg, hlak, temporary_stage, dh, residual)
@ brief Lake package bisection method
Definition: gwf-lak.f90:5299
lakmodule::lak_ad
subroutine lak_ad(this)
Add package connection to matrix.
Definition: gwf-lak.f90:3405
lakmodule::lak_solve
subroutine lak_solve(this, update)
Solve for lake stage.
Definition: gwf-lak.f90:4921
lakmodule::ftype
character(len=lenftype) ftype
Definition: gwf-lak.f90:43
lakmodule::lak_fc
subroutine lak_fc(this, rhs, ia, idxglo, matrix_sln)
Copy rhs and hcof into solution rhs and amat.
Definition: gwf-lak.f90:3555
lakmodule::lak_set_attribute_error
subroutine lak_set_attribute_error(this, ilak, keyword, msg)
Issue a parameter error for lakweslls(ilak)
Definition: gwf-lak.f90:3070
lakmodule::lak_calculate_runoff
subroutine lak_calculate_runoff(this, ilak, ro)
Calculate runoff to a lake.
Definition: gwf-lak.f90:2447
lakmodule::define_listlabel
subroutine define_listlabel(this)
Define the list heading that is written to iout when PRINT_INPUT option is used.
Definition: gwf-lak.f90:4336
lakmodule::lak_read_table
subroutine lak_read_table(this, ilak, filename, laketable)
Read the lake table for this package.
Definition: gwf-lak.f90:1169
lakmodule::lak_calculate_withdrawal
subroutine lak_calculate_withdrawal(this, ilak, avail, wr)
Calculate the withdrawal from a lake subject to an available volume.
Definition: gwf-lak.f90:2487
lakmodule::lak_calculate_conn_warea
subroutine lak_calculate_conn_warea(this, ilak, iconn, stage, head, wa)
Calculate the wetted area of a lake connection at a given stage.
Definition: gwf-lak.f90:2053
lakmodule::lak_df_obs
subroutine lak_df_obs(this)
Store observation type supported by LAK package. Overrides BndTypebnd_df_obs.
Definition: gwf-lak.f90:4402
lakmodule::lak_allocate_scalars
subroutine lak_allocate_scalars(this)
Allocate scalar members.
Definition: gwf-lak.f90:329
lakmodule::lak_bound_update
subroutine lak_bound_update(this)
Store the lake head and connection conductance in the bound array.
Definition: gwf-lak.f90:4896
lakmodule::lak_calculate_cond_head
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...
Definition: gwf-lak.f90:2178
lakmodule::lak_ot_bdsummary
subroutine lak_ot_bdsummary(this, kstp, kper, iout, ibudfl)
Write LAK budget to listing file.
Definition: gwf-lak.f90:4157
lakmodule::lak_create
subroutine, public lak_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
Create a new LAK Package and point bndobj to the new package.
Definition: gwf-lak.f90:291
lakmodule::lak_bd_obs
subroutine lak_bd_obs(this)
Calculate observations this time step and call ObsTypeSaveOneSimval for each LakType observation.
Definition: gwf-lak.f90:4507
lakmodule::lak_estimate_conn_exchange
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.
Definition: gwf-lak.f90:2361
lakmodule::lak_fill_budobj
subroutine lak_fill_budobj(this)
Copy flow terms into thisbudobj.
Definition: gwf-lak.f90:5746
lakmodule::lak_get_internal_mover
subroutine lak_get_internal_mover(this, ilak, outoutf)
Get the mover outflow from a lake to another lake.
Definition: gwf-lak.f90:2695
lakmodule::lak_calculate_outlet_inflow
subroutine lak_calculate_outlet_inflow(this, ilak, outinf)
Calculate the outlet inflow to a lake.
Definition: gwf-lak.f90:2536
lakmodule::lak_calculate_conn_exchange
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.
Definition: gwf-lak.f90:2298
lakmodule::lak_ar
subroutine lak_ar(this)
Allocate and Read.
Definition: gwf-lak.f90:3248
mathutilmodule::is_close
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
Definition: MathUtil.f90:46
memoryhelpermodule::create_mem_path
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
Definition: MemoryHelper.f90:22
observemodule
This module contains the derived types ObserveType and ObsDataType.
Definition: Observe.f90:15
obsmodule
This module contains the derived type ObsType.
Definition: Obs.f90:127
openspecmodule::access
character(len=20) access
Definition: OpenSpec.f90:7
openspecmodule::form
character(len=20) form
Definition: OpenSpec.f90:7
simmodule
This module contains simulation methods.
Definition: Sim.f90:10
simmodule::store_error
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
simmodule::count_errors
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
simmodule::deprecation_warning
subroutine, public deprecation_warning(cblock, cvar, cver, endmsg, iunit)
Store deprecation warning message.
Definition: Sim.f90:256
simmodule::store_error_unit
subroutine, public store_error_unit(iunit, terminate)
Store the file unit number.
Definition: Sim.f90:168
simvariablesmodule
This module contains simulation variables.
Definition: SimVariables.f90:9
simvariablesmodule::errmsg
character(len=maxcharlen) errmsg
error message string
Definition: SimVariables.f90:28
simvariablesmodule::ifailedstepretry
integer(i4b) ifailedstepretry
current retry for this time step
Definition: SimVariables.f90:42
simvariablesmodule::warnmsg
character(len=maxcharlen) warnmsg
warning message string
Definition: SimVariables.f90:29
smoothingmodule::squadraticsaturation
real(dp) function squadraticsaturation(top, bot, x, eps)
@ brief sQuadraticSaturation
Definition: SmoothingFunctions.f90:276
smoothingmodule::squadraticsaturationderivative
real(dp) function squadraticsaturationderivative(top, bot, x, eps)
@ brief Derivative of the quadratic saturation function
Definition: SmoothingFunctions.f90:365
smoothingmodule::sqsaturationderivative
real(dp) function sqsaturationderivative(top, bot, x, c1, c2)
@ brief sQSaturationDerivative
Definition: SmoothingFunctions.f90:467
smoothingmodule::sqsaturation
real(dp) function sqsaturation(top, bot, x, c1, c2)
@ brief sQSaturation
Definition: SmoothingFunctions.f90:413
tablemodule::table_cr
subroutine, public table_cr(this, name, title)
Definition: Table.f90:87
tdismodule::pertim
real(dp), pointer, public pertim
time relative to start of stress period
Definition: tdis.f90:30
tdismodule::totim
real(dp), pointer, public totim
time relative to start of simulation
Definition: tdis.f90:32
tdismodule::kstp
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
tdismodule::kper
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
tdismodule::delt
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:29
tdismodule::nper
integer(i4b), pointer, public nper
number of stress period
Definition: tdis.f90:21
timeseriesmanagermodule::read_value_or_time_series_adv
subroutine, public read_value_or_time_series_adv(textInput, ii, jj, bndElem, pkgName, auxOrBnd, tsManager, iprpak, varName)
Call this subroutine from advanced packages to define timeseries link for a variable (varName).
Definition: TimeSeriesManager.f90:586
bndmodule::bndtype
@ brief BndType
Definition: BoundaryPackage.f90:48