MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
UzfCellGroup.f90
Go to the documentation of this file.
1 module uzfcellgroupmodule
2 
3  use kindmodule, only: dp, i4b
4  use constantsmodule, only: dzero, dem30, dem20, dem15, dem14, dem12, dem10, &
5  dem9, dem7, dem6, dem5, dem4, dem3, dhalf, done, &
6  dtwo, dthree, dep20
7  use smoothingmodule
8  use tdismodule, only: itmuni, delt, kper
9  use uzfetutilmodule, only: etfunc_lin, etfunc_nlin
10 
11  implicit none
12  private
13  public :: uzfcellgrouptype
14 
15  type :: uzfcellgrouptype
16 
17  integer(I4B) :: imem_manager
18  real(dp), pointer, dimension(:), contiguous :: thtr => null()
19  real(dp), pointer, dimension(:), contiguous :: thts => null()
20  real(dp), pointer, dimension(:), contiguous :: thti => null()
21  real(dp), pointer, dimension(:), contiguous :: eps => null()
22  real(dp), pointer, dimension(:), contiguous :: extwc => null()
23  real(dp), pointer, dimension(:), contiguous :: ha => null()
24  real(dp), pointer, dimension(:), contiguous :: hroot => null()
25  real(dp), pointer, dimension(:), contiguous :: rootact => null()
26  real(dp), pointer, dimension(:), contiguous :: etact => null()
27  real(dp), dimension(:, :), pointer, contiguous :: uzspst => null()
28  real(dp), dimension(:, :), pointer, contiguous :: uzthst => null()
29  real(dp), dimension(:, :), pointer, contiguous :: uzflst => null()
30  real(dp), dimension(:, :), pointer, contiguous :: uzdpst => null()
31  integer(I4B), pointer, dimension(:), contiguous :: nwavst => null()
32  real(dp), pointer, dimension(:), contiguous :: totflux => null()
33  integer(I4B), pointer, dimension(:), contiguous :: nwav => null()
34  integer(I4B), pointer, dimension(:), contiguous :: ntrail => null()
35  real(dp), pointer, dimension(:), contiguous :: sinf => null()
36  real(dp), pointer, dimension(:), contiguous :: finf => null()
37  real(dp), pointer, dimension(:), contiguous :: pet => null()
38  real(dp), pointer, dimension(:), contiguous :: petmax => null()
39  real(dp), pointer, dimension(:), contiguous :: extdp => null()
40  real(dp), pointer, dimension(:), contiguous :: extdpuz => null()
41  real(dp), pointer, dimension(:), contiguous :: finf_rej => null()
42  real(dp), pointer, dimension(:), contiguous :: gwet => null()
43  real(dp), pointer, dimension(:), contiguous :: uzfarea => null()
44  real(dp), pointer, dimension(:), contiguous :: cellarea => null()
45  real(dp), pointer, dimension(:), contiguous :: celtop => null()
46  real(dp), pointer, dimension(:), contiguous :: celbot => null()
47  real(dp), pointer, dimension(:), contiguous :: landtop => null()
48  real(dp), pointer, dimension(:), contiguous :: watab => null()
49  real(dp), pointer, dimension(:), contiguous :: watabold => null()
50  real(dp), pointer, dimension(:), contiguous :: vks => null()
51  real(dp), pointer, dimension(:), contiguous :: surfdep => null()
52  real(dp), pointer, dimension(:), contiguous :: surflux => null()
53  real(dp), pointer, dimension(:), contiguous :: surfluxbelow => null()
54  real(dp), pointer, dimension(:), contiguous :: surfseep => null()
55  real(dp), pointer, dimension(:), contiguous :: gwpet => null()
56  integer(I4B), pointer, dimension(:), contiguous :: landflag => null()
57  integer(I4B), pointer, dimension(:), contiguous :: ivertcon => null()
58 
59  contains
60 
61  procedure :: init
62  procedure :: setdata
63  procedure :: sethead
64  procedure :: setdatauzfarea
65  procedure :: setdatafinf
66  procedure :: setdataet
67  procedure :: setdataetwc
68  procedure :: setdataetha
69  procedure :: setwaves
70  procedure :: wave_shift
71  procedure :: routewaves
72  procedure :: uzflow
73  procedure :: addrech
74  procedure :: trailwav
75  procedure :: leadwav
76  procedure :: advance
77  procedure :: solve
78  procedure :: unsat_stor
79  procedure :: update_wav
80  procedure :: simgwet
81  procedure :: caph
82  procedure :: rate_et_z
83  procedure :: uzet
84  procedure :: uz_rise
85  procedure :: rejfinf
86  procedure :: gwseep
87  procedure :: setbelowpet
88  procedure :: setgwpet
89  procedure :: dealloc
90  procedure :: get_water_content_at_depth
91  procedure :: get_wcnew
92  end type uzfcellgrouptype
93 
94 contains
95 
96  !> @brief Allocate and set uzf object variables
97  !<
98  subroutine init(this, ncells, nwav, memory_path)
99  ! -- modules
100  use memorymanagermodule, only: mem_allocate
101  ! -- dummy
102  class(uzfcellgrouptype) :: this
103  integer(I4B), intent(in) :: nwav
104  integer(I4B), intent(in) :: ncells
105  character(len=*), intent(in), optional :: memory_path
106  ! -- local
107  integer(I4B) :: icell
108  !
109  ! -- Use mem_allocate if memory path is passed in, otherwise it's a temp object
110  if (present(memory_path)) then
111  this%imem_manager = 1
112  call mem_allocate(this%uzdpst, nwav, ncells, 'UZDPST', memory_path)
113  call mem_allocate(this%uzthst, nwav, ncells, 'UZTHST', memory_path)
114  call mem_allocate(this%uzflst, nwav, ncells, 'UZFLST', memory_path)
115  call mem_allocate(this%uzspst, nwav, ncells, 'UZSPST', memory_path)
116  call mem_allocate(this%nwavst, ncells, 'NWAVST', memory_path)
117  call mem_allocate(this%thtr, ncells, 'THTR', memory_path)
118  call mem_allocate(this%thts, ncells, 'THTS', memory_path)
119  call mem_allocate(this%thti, ncells, 'THTI', memory_path)
120  call mem_allocate(this%eps, ncells, 'EPS', memory_path)
121  call mem_allocate(this%ha, ncells, 'HA', memory_path)
122  call mem_allocate(this%hroot, ncells, 'HROOT', memory_path)
123  call mem_allocate(this%rootact, ncells, 'ROOTACT', memory_path)
124  call mem_allocate(this%extwc, ncells, 'EXTWC', memory_path)
125  call mem_allocate(this%etact, ncells, 'ETACT', memory_path)
126  call mem_allocate(this%nwav, ncells, 'NWAV', memory_path)
127  call mem_allocate(this%ntrail, ncells, 'NTRAIL', memory_path)
128  call mem_allocate(this%totflux, ncells, 'TOTFLUX', memory_path)
129  call mem_allocate(this%sinf, ncells, 'SINF', memory_path)
130  call mem_allocate(this%finf, ncells, 'FINF', memory_path)
131  call mem_allocate(this%finf_rej, ncells, 'FINF_REJ', memory_path)
132  call mem_allocate(this%gwet, ncells, 'GWET', memory_path)
133  call mem_allocate(this%uzfarea, ncells, 'UZFAREA', memory_path)
134  call mem_allocate(this%cellarea, ncells, 'CELLAREA', memory_path)
135  call mem_allocate(this%celtop, ncells, 'CELTOP', memory_path)
136  call mem_allocate(this%celbot, ncells, 'CELBOT', memory_path)
137  call mem_allocate(this%landtop, ncells, 'LANDTOP', memory_path)
138  call mem_allocate(this%watab, ncells, 'WATAB', memory_path)
139  call mem_allocate(this%watabold, ncells, 'WATABOLD', memory_path)
140  call mem_allocate(this%surfdep, ncells, 'SURFDEP', memory_path)
141  call mem_allocate(this%vks, ncells, 'VKS', memory_path)
142  call mem_allocate(this%surflux, ncells, 'SURFLUX', memory_path)
143  call mem_allocate(this%surfluxbelow, ncells, 'SURFLUXBELOW', memory_path)
144  call mem_allocate(this%surfseep, ncells, 'SURFSEEP', memory_path)
145  call mem_allocate(this%gwpet, ncells, 'GWPET', memory_path)
146  call mem_allocate(this%pet, ncells, 'PET', memory_path)
147  call mem_allocate(this%petmax, ncells, 'PETMAX', memory_path)
148  call mem_allocate(this%extdp, ncells, 'EXTDP', memory_path)
149  call mem_allocate(this%extdpuz, ncells, 'EXTDPUZ', memory_path)
150  call mem_allocate(this%landflag, ncells, 'LANDFLAG', memory_path)
151  call mem_allocate(this%ivertcon, ncells, 'IVERTCON', memory_path)
152  else
153  this%imem_manager = 0
154  allocate (this%uzdpst(nwav, ncells))
155  allocate (this%uzthst(nwav, ncells))
156  allocate (this%uzflst(nwav, ncells))
157  allocate (this%uzspst(nwav, ncells))
158  allocate (this%nwavst(ncells))
159  allocate (this%thtr(ncells))
160  allocate (this%thts(ncells))
161  allocate (this%thti(ncells))
162  allocate (this%eps(ncells))
163  allocate (this%ha(ncells))
164  allocate (this%hroot(ncells))
165  allocate (this%rootact(ncells))
166  allocate (this%extwc(ncells))
167  allocate (this%etact(ncells))
168  allocate (this%nwav(ncells))
169  allocate (this%ntrail(ncells))
170  allocate (this%totflux(ncells))
171  allocate (this%sinf(ncells))
172  allocate (this%finf(ncells))
173  allocate (this%finf_rej(ncells))
174  allocate (this%gwet(ncells))
175  allocate (this%uzfarea(ncells))
176  allocate (this%cellarea(ncells))
177  allocate (this%celtop(ncells))
178  allocate (this%celbot(ncells))
179  allocate (this%landtop(ncells))
180  allocate (this%watab(ncells))
181  allocate (this%watabold(ncells))
182  allocate (this%surfdep(ncells))
183  allocate (this%vks(ncells))
184  allocate (this%surflux(ncells))
185  allocate (this%surfluxbelow(ncells))
186  allocate (this%surfseep(ncells))
187  allocate (this%gwpet(ncells))
188  allocate (this%pet(ncells))
189  allocate (this%petmax(ncells))
190  allocate (this%extdp(ncells))
191  allocate (this%extdpuz(ncells))
192  allocate (this%landflag(ncells))
193  allocate (this%ivertcon(ncells))
194  end if
195  do icell = 1, ncells
196  this%uzdpst(:, icell) = dzero
197  this%uzthst(:, icell) = dzero
198  this%uzflst(:, icell) = dzero
199  this%uzspst(:, icell) = dzero
200  this%nwavst(icell) = 1
201  this%thtr(icell) = dzero
202  this%thts(icell) = dzero
203  this%thti(icell) = dzero
204  this%eps(icell) = dzero
205  this%ha(icell) = dzero
206  this%hroot(icell) = dzero
207  this%rootact(icell) = dzero
208  this%extwc(icell) = dzero
209  this%etact(icell) = dzero
210  this%nwav(icell) = nwav
211  this%ntrail(icell) = 0
212  this%totflux(icell) = dzero
213  this%sinf(icell) = dzero
214  this%finf(icell) = dzero
215  this%finf_rej(icell) = dzero
216  this%gwet(icell) = dzero
217  this%uzfarea(icell) = dzero
218  this%cellarea(icell) = dzero
219  this%celtop(icell) = dzero
220  this%celbot(icell) = dzero
221  this%landtop(icell) = dzero
222  this%watab(icell) = dzero
223  this%watabold(icell) = dzero
224  this%surfdep(icell) = dzero
225  this%vks(icell) = dzero
226  this%surflux(icell) = dzero
227  this%surfluxbelow(icell) = dzero
228  this%surfseep(icell) = dzero
229  this%gwpet(icell) = dzero
230  this%pet(icell) = dzero
231  this%petmax(icell) = dzero
232  this%extdp(icell) = dzero
233  this%extdpuz(icell) = dzero
234  this%landflag(icell) = 0
235  this%ivertcon(icell) = 0
236  end do
237  end subroutine init
238 
239  !> @brief Deallocate uzf object variables
240  !<
241  subroutine dealloc(this)
242  ! -- modules
243  use memorymanagermodule, only: mem_deallocate
244  ! -- dummy
245  class(uzfcellgrouptype) :: this
246  !
247  ! -- deallocate based on whether or not memory manager was used
248  if (this%imem_manager == 0) then
249  deallocate (this%uzdpst)
250  deallocate (this%uzthst)
251  deallocate (this%uzflst)
252  deallocate (this%uzspst)
253  deallocate (this%nwavst)
254  deallocate (this%thtr)
255  deallocate (this%thts)
256  deallocate (this%thti)
257  deallocate (this%eps)
258  deallocate (this%ha)
259  deallocate (this%hroot)
260  deallocate (this%rootact)
261  deallocate (this%extwc)
262  deallocate (this%etact)
263  deallocate (this%nwav)
264  deallocate (this%ntrail)
265  deallocate (this%totflux)
266  deallocate (this%sinf)
267  deallocate (this%finf)
268  deallocate (this%finf_rej)
269  deallocate (this%gwet)
270  deallocate (this%uzfarea)
271  deallocate (this%cellarea)
272  deallocate (this%celtop)
273  deallocate (this%celbot)
274  deallocate (this%landtop)
275  deallocate (this%watab)
276  deallocate (this%watabold)
277  deallocate (this%surfdep)
278  deallocate (this%vks)
279  deallocate (this%surflux)
280  deallocate (this%surfluxbelow)
281  deallocate (this%surfseep)
282  deallocate (this%gwpet)
283  deallocate (this%pet)
284  deallocate (this%petmax)
285  deallocate (this%extdp)
286  deallocate (this%extdpuz)
287  deallocate (this%landflag)
288  deallocate (this%ivertcon)
289  else
290  call mem_deallocate(this%uzdpst)
291  call mem_deallocate(this%uzthst)
292  call mem_deallocate(this%uzflst)
293  call mem_deallocate(this%uzspst)
294  call mem_deallocate(this%nwavst)
295  call mem_deallocate(this%thtr)
296  call mem_deallocate(this%thts)
297  call mem_deallocate(this%thti)
298  call mem_deallocate(this%eps)
299  call mem_deallocate(this%ha)
300  call mem_deallocate(this%hroot)
301  call mem_deallocate(this%rootact)
302  call mem_deallocate(this%extwc)
303  call mem_deallocate(this%etact)
304  call mem_deallocate(this%nwav)
305  call mem_deallocate(this%ntrail)
306  call mem_deallocate(this%totflux)
307  call mem_deallocate(this%sinf)
308  call mem_deallocate(this%finf)
309  call mem_deallocate(this%finf_rej)
310  call mem_deallocate(this%gwet)
311  call mem_deallocate(this%uzfarea)
312  call mem_deallocate(this%cellarea)
313  call mem_deallocate(this%celtop)
314  call mem_deallocate(this%celbot)
315  call mem_deallocate(this%landtop)
316  call mem_deallocate(this%watab)
317  call mem_deallocate(this%watabold)
318  call mem_deallocate(this%surfdep)
319  call mem_deallocate(this%vks)
320  call mem_deallocate(this%surflux)
321  call mem_deallocate(this%surfluxbelow)
322  call mem_deallocate(this%surfseep)
323  call mem_deallocate(this%gwpet)
324  call mem_deallocate(this%pet)
325  call mem_deallocate(this%petmax)
326  call mem_deallocate(this%extdp)
327  call mem_deallocate(this%extdpuz)
328  call mem_deallocate(this%landflag)
329  call mem_deallocate(this%ivertcon)
330  end if
331  end subroutine dealloc
332 
333  !> @brief Set uzf object material properties
334  !<
335  subroutine setdata(this, icell, area, top, bot, surfdep, vks, thtr, thts, &
336  thti, eps, ntrail, landflag, ivertcon)
337  ! -- dummy
338  class(uzfcellgrouptype) :: this
339  integer(I4B), intent(in) :: icell
340  real(DP), intent(in) :: area
341  real(DP), intent(in) :: top
342  real(DP), intent(in) :: bot
343  real(DP), intent(in) :: surfdep
344  real(DP), intent(in) :: vks
345  real(DP), intent(in) :: thtr
346  real(DP), intent(in) :: thts
347  real(DP), intent(in) :: thti
348  real(DP), intent(in) :: eps
349  integer(I4B), intent(in) :: ntrail
350  integer(I4B), intent(in) :: landflag
351  integer(I4B), intent(in) :: ivertcon
352  !
353  ! -- set the values for uzf cell icell
354  this%landflag(icell) = landflag
355  this%ivertcon(icell) = ivertcon
356  this%surfdep(icell) = surfdep
357  this%uzfarea(icell) = area
358  this%cellarea(icell) = area
359  if (this%landflag(icell) == 1) then
360  this%celtop(icell) = top - dhalf * this%surfdep(icell)
361  else
362  this%celtop(icell) = top
363  end if
364  this%celbot(icell) = bot
365  this%vks(icell) = vks
366  this%thtr(icell) = thtr
367  this%thts(icell) = thts
368  this%thti(icell) = thti
369  this%eps(icell) = eps
370  this%ntrail(icell) = ntrail
371  this%pet(icell) = dzero
372  this%extdp(icell) = dzero
373  this%extwc(icell) = dzero
374  this%ha(icell) = dzero
375  this%hroot(icell) = dzero
376  end subroutine setdata
377 
378  !> @brief Set initial head for uzf object
379  !<
380  subroutine sethead(this, icell, hgwf)
381  ! -- dummy
382  class(uzfcellgrouptype) :: this
383  integer(I4B), intent(in) :: icell
384  real(DP), intent(in) :: hgwf
385  !
386  ! -- set initial head
387  this%watab(icell) = this%celbot(icell)
388  if (hgwf > this%celbot(icell)) this%watab(icell) = hgwf
389  if (this%watab(icell) > this%celtop(icell)) &
390  this%watab(icell) = this%celtop(icell)
391  this%watabold(icell) = this%watab(icell)
392  end subroutine sethead
393 
394  !> @brief Set infiltration
395  !<
396  subroutine setdatafinf(this, icell, finf)
397  ! -- dummy
398  class(uzfcellgrouptype) :: this
399  integer(I4B), intent(in) :: icell
400  real(DP), intent(in) :: finf
401  !
402  if (this%landflag(icell) == 1) then
403  this%sinf(icell) = finf
404  this%finf(icell) = finf
405  else
406  this%sinf(icell) = dzero
407  this%finf(icell) = dzero
408  end if
409  this%finf_rej(icell) = dzero
410  this%surflux(icell) = dzero
411  this%surfluxbelow(icell) = dzero
412  end subroutine setdatafinf
413 
414  !> @brief Set uzfarea using cellarea and areamult
415  !<
416  subroutine setdatauzfarea(this, icell, areamult)
417  ! -- dummy
418  class(uzfcellgrouptype) :: this
419  integer(I4B), intent(in) :: icell
420  real(DP), intent(in) :: areamult
421  !
422  ! -- set uzf area
423  this%uzfarea(icell) = this%cellarea(icell) * areamult
424  end subroutine setdatauzfarea
425 
426  !> @brief Set unsaturated ET-related variables
427  !<
428  subroutine setdataet(this, icell, jbelow, pet, extdp)
429  ! -- dummy
430  class(uzfcellgrouptype) :: this
431  integer(I4B), intent(in) :: icell
432  integer(I4B), intent(in) :: jbelow
433  real(DP), intent(in) :: pet
434  real(DP), intent(in) :: extdp
435  ! -- local
436  real(DP) :: thick
437  !
438  if (this%landflag(icell) == 1) then
439  this%pet(icell) = pet
440  this%gwpet(icell) = pet
441  else
442  this%pet(icell) = dzero
443  this%gwpet(icell) = dzero
444  end if
445  thick = this%celtop(icell) - this%celbot(icell)
446  this%extdp(icell) = extdp
447  if (this%landflag(icell) > 0) then
448  this%landtop(icell) = this%celtop(icell)
449  this%petmax(icell) = this%pet(icell)
450  end if
451  !
452  ! -- set uz extinction depth
453  if (this%landtop(icell) - this%extdp(icell) < this%celbot(icell)) then
454  this%extdpuz(icell) = thick
455  else
456  this%extdpuz(icell) = this%celtop(icell) - &
457  (this%landtop(icell) - this%extdp(icell))
458  end if
459  if (this%extdpuz(icell) < dzero) this%extdpuz(icell) = dzero
460  if (this%extdpuz(icell) > dem7 .and. this%extdp(icell) < dem7) &
461  this%extdp(icell) = this%extdpuz(icell)
462  !
463  ! -- set pet for underlying cell
464  if (jbelow > 0) then
465  this%landtop(jbelow) = this%landtop(icell)
466  this%petmax(jbelow) = this%petmax(icell)
467  end if
468  end subroutine setdataet
469 
470  !> @brief Subtract aet from pet to calculate residual et for gw
471  !<
472  subroutine setgwpet(this, icell)
473  ! -- modules
474  use tdismodule, only: delt
475  ! -- dummy
476  class(uzfcellgrouptype) :: this
477  integer(I4B), intent(in) :: icell
478  ! -- dummy
479  real(DP) :: pet
480  !
481  pet = dzero
482  !
483  ! -- reduce pet for gw by uzet
484  pet = this%pet(icell) - this%etact(icell) / delt
485  if (pet < dzero) pet = dzero
486  this%gwpet(icell) = pet
487  end subroutine setgwpet
488 
489  !> @brief Subtract aet from pet to calculate residual et for deeper cells
490  !<
491  subroutine setbelowpet(this, icell, jbelow)
492  ! -- modules
493  use tdismodule, only: delt
494  ! -- dummy
495  class(uzfcellgrouptype) :: this
496  integer(I4B), intent(in) :: icell
497  integer(I4B), intent(in) :: jbelow
498  ! -- dummy
499  real(DP) :: pet
500  !
501  pet = dzero
502  !
503  ! -- transfer unmet pet to lower cell
504  !
505  if (this%extdpuz(jbelow) > dem3) then
506  pet = this%pet(icell) - this%etact(icell) / delt - &
507  this%gwet(icell) / this%uzfarea(icell)
508  if (pet < dzero) pet = dzero
509  end if
510  this%pet(jbelow) = pet
511  end subroutine setbelowpet
512 
513  !> @brief Set extinction water content
514  !<
515  subroutine setdataetwc(this, icell, jbelow, extwc)
516  ! -- dummy
517  class(uzfcellgrouptype) :: this
518  integer(I4B), intent(in) :: icell
519  integer(I4B), intent(in) :: jbelow
520  real(DP), intent(in) :: extwc
521  !
522  ! -- set extinction water content
523  this%extwc(icell) = extwc
524  if (jbelow > 0) this%extwc(jbelow) = extwc
525  end subroutine setdataetwc
526 
527  !> @brief Set variables for head-based unsaturated flow
528  !<
529  subroutine setdataetha(this, icell, jbelow, ha, hroot, rootact)
530  ! -- dummy
531  class(uzfcellgrouptype) :: this
532  integer(I4B), intent(in) :: icell
533  integer(I4B), intent(in) :: jbelow
534  real(DP), intent(in) :: ha
535  real(DP), intent(in) :: hroot
536  real(DP), intent(in) :: rootact
537  !
538  ! -- set variables
539  this%ha(icell) = ha
540  this%hroot(icell) = hroot
541  this%rootact(icell) = rootact
542  if (jbelow > 0) then
543  this%ha(jbelow) = ha
544  this%hroot(jbelow) = hroot
545  this%rootact(jbelow) = rootact
546  end if
547  end subroutine setdataetha
548 
549  !> @brief Set variables to advance to new time step. nothing yet.
550  !<
551  subroutine advance(this, icell)
552  ! -- dummy
553  class(uzfcellgrouptype) :: this
554  integer(I4B), intent(in) :: icell
555  !
556  ! -- set variables
557  this%surfseep(icell) = dzero
558  end subroutine advance
559 
560  !> @brief Formulate the unsaturated flow object, calculate terms for gwf
561  !! equation
562  !<
563  subroutine solve(this, thiswork, jbelow, icell, totfluxtot, ietflag, &
564  issflag, iseepflag, hgwf, qfrommvr, ierr, &
565  reset_state, trhs, thcof, deriv, watercontent)
566  ! -- modules
567  use tdismodule, only: delt
568  ! -- dummy
569  class(uzfcellgrouptype) :: this
570  type(uzfcellgrouptype) :: thiswork !< work object for resetting wave state
571  integer(I4B), intent(in) :: jbelow !< number of underlying uzf object or 0 if none
572  integer(I4B), intent(in) :: icell !< number of this uzf object
573  real(DP), intent(inout) :: totfluxtot !<
574  integer(I4B), intent(in) :: ietflag !< et is off (0) or based one water content (1) or pressure (2)
575  integer(I4B), intent(in) :: issflag !< steady state flag
576  integer(I4B), intent(in) :: iseepflag !< discharge to land is active (1) or not (0)
577  real(DP), intent(in) :: hgwf !< head for cell icell
578  real(DP), intent(in) :: qfrommvr !< water inflow from mover
579  integer(I4B), intent(inout) :: ierr !< flag indicating not enough waves
580  logical, intent(in) :: reset_state !< flag indicating that waves should be reset after solution
581  real(DP), intent(inout), optional :: trhs !< total uzf rhs contribution to GWF model
582  real(DP), intent(inout), optional :: thcof !< total uzf hcof contribution to GWF model
583  real(DP), intent(inout), optional :: deriv !< derivate term for contribution to GWF model
584  real(DP), intent(inout), optional :: watercontent !< calculated water content
585  ! -- local
586  real(DP) :: test
587  real(DP) :: scale
588  real(DP) :: seep
589  real(DP) :: finfact
590  real(DP) :: derivfinf
591  real(DP) :: trhsfinf
592  real(DP) :: thcoffinf
593  real(DP) :: trhsseep
594  real(DP) :: thcofseep
595  real(DP) :: deriv1
596  real(DP) :: deriv2
597  !
598  ! -- initialize variables
599  totfluxtot = dzero
600  trhsfinf = dzero
601  thcoffinf = dzero
602  trhsseep = dzero
603  thcofseep = dzero
604  this%finf_rej(icell) = dzero
605  this%surflux(icell) = this%finf(icell) + qfrommvr / this%uzfarea(icell)
606  this%watab(icell) = hgwf
607  this%surfseep(icell) = dzero
608  seep = dzero
609  finfact = dzero
610  this%etact(icell) = dzero
611  this%surfluxbelow(icell) = dzero
612  if (this%ivertcon(icell) > 0) then
613  this%finf(jbelow) = dzero
614  end if
615  if (this%watab(icell) < this%celbot(icell)) &
616  this%watab(icell) = this%celbot(icell)
617  !
618  ! -- initialize derivative variables
619  deriv1 = dzero
620  deriv2 = dzero
621  derivfinf = dzero
622  !
623  ! -- save wave states for resetting after iteration.
624  if (reset_state) then
625  call thiswork%wave_shift(this, 1, icell, 0, 1, this%nwavst(icell), 1)
626  end if
627  !
628  if (this%watab(icell) > this%celtop(icell)) &
629  this%watab(icell) = this%celtop(icell)
630  !
631  ! -- add water from mover to applied infiltration.
632  if (this%surflux(icell) > this%vks(icell)) then
633  this%surflux(icell) = this%vks(icell)
634  end if
635  !
636  ! -- saturation excess rejected infiltration
637  if (this%landflag(icell) == 1) then
638  call this%rejfinf(icell, deriv1, hgwf, trhsfinf, thcoffinf, finfact)
639  this%surflux(icell) = finfact
640  end if
641  !
642  ! -- calculate rejected infiltration
643  this%finf_rej(icell) = this%finf(icell) + &
644  (qfrommvr / this%uzfarea(icell)) - this%surflux(icell)
645  !
646  ! -- calculate groundwater discharge
647  if (iseepflag > 0 .and. this%landflag(icell) == 1) then
648  call this%gwseep(icell, deriv2, scale, hgwf, trhsseep, thcofseep, seep)
649  this%surfseep(icell) = seep
650  end if
651  !
652  ! -- route water through unsat zone, calc. storage change and recharge
653  test = this%watab(icell)
654  if (this%watabold(icell) - test < -dem15) test = this%watabold(icell)
655  if (this%celtop(icell) - test > dem15) then
656  if (issflag == 0) then
657  call this%routewaves(totfluxtot, delt, ietflag, icell, ierr)
658  if (ierr > 0) return
659  call this%uz_rise(icell, totfluxtot)
660  this%totflux(icell) = totfluxtot
661  if (this%ivertcon(icell) > 0) then
662  call this%addrech(icell, jbelow, hgwf, trhsfinf, thcoffinf, &
663  derivfinf, delt)
664  end if
665  else
666  this%totflux(icell) = this%surflux(icell) * delt
667  totfluxtot = this%surflux(icell) * delt
668  end if
669  thcoffinf = dzero
670  trhsfinf = this%totflux(icell) * this%uzfarea(icell) / delt
671  if (.not. reset_state) then
672  call this%update_wav(icell, delt, issflag, 0)
673  end if
674  else
675  this%totflux(icell) = this%surflux(icell) * delt
676  totfluxtot = this%surflux(icell) * delt
677  if (.not. reset_state) then
678  call this%update_wav(icell, delt, issflag, 1)
679  end if
680  end if
681  !
682  ! -- If formulating, then these variables will be present
683  if (present(deriv)) deriv = deriv1 + deriv2 + derivfinf
684  if (present(trhs)) trhs = trhsfinf + trhsseep
685  if (present(thcof)) thcof = thcoffinf + thcofseep
686  !
687  ! -- Assign water content prior to resetting waves
688  if (present(watercontent)) then
689  watercontent = this%get_wcnew(icell)
690  end if
691  !
692  ! -- reset waves to previous state for next iteration
693  if (reset_state) then
694  call this%wave_shift(thiswork, icell, 1, 0, 1, thiswork%nwavst(1), 1)
695  end if
696  end subroutine solve
697 
698  !> @brief Add recharge or infiltration to cells
699  !<
700  subroutine addrech(this, icell, jbelow, hgwf, trhs, thcof, deriv, delt)
701  ! -- dummy
702  class(uzfcellgrouptype) :: this
703  integer(I4B), intent(in) :: icell
704  integer(I4B), intent(in) :: jbelow
705  real(DP), intent(inout) :: trhs
706  real(DP), intent(inout) :: thcof
707  real(DP), intent(inout) :: deriv
708  real(DP), intent(in) :: delt
709  real(DP), intent(in) :: hgwf
710  ! -- local
711  real(DP) :: fcheck
712  real(DP) :: x, scale, range
713  !
714  ! -- initialize
715  range = dem5
716  deriv = dzero
717  thcof = dzero
718  trhs = this%uzfarea(icell) * this%totflux(icell) / delt
719  if (this%totflux(icell) < dem14) return
720  scale = done
721  !
722  ! -- smoothly reduce flow between cells when head close to cell top
723  x = hgwf - (this%celbot(icell) - range)
724  call sscurve(x, range, deriv, scale)
725  deriv = this%uzfarea(icell) * deriv * this%totflux(icell) / delt
726  this%finf(jbelow) = (done - scale) * this%totflux(icell) / delt
727  fcheck = this%finf(jbelow) - this%vks(jbelow)
728  !
729  ! -- reduce flow between cells when vks is too small
730  if (fcheck < dem14) fcheck = dzero
731  this%finf(jbelow) = this%finf(jbelow) - fcheck
732  this%surfluxbelow(icell) = this%finf(jbelow)
733  this%totflux(icell) = scale * this%totflux(icell) + fcheck * delt
734  trhs = this%uzfarea(icell) * this%totflux(icell) / delt
735  end subroutine addrech
736 
737  !> @brief Reject applied infiltration due to low vks
738  !<
739  subroutine rejfinf(this, icell, deriv, hgwf, trhs, thcof, finfact)
740  ! -- dummy
741  class(uzfcellgrouptype) :: this
742  integer(I4B), intent(in) :: icell
743  real(DP), intent(inout) :: deriv
744  real(DP), intent(inout) :: finfact
745  real(DP), intent(inout) :: thcof
746  real(DP), intent(inout) :: trhs
747  real(DP), intent(in) :: hgwf
748  ! -- local
749  real(DP) :: x, range, scale, q
750  !
751  range = this%surfdep(icell)
752  q = this%surflux(icell)
753  finfact = q
754  trhs = finfact * this%uzfarea(icell)
755  x = this%celtop(icell) - hgwf
756  call slinear(x, range, deriv, scale)
757  deriv = -q * deriv * this%uzfarea(icell) * scale
758  if (scale < done) then
759  finfact = q * scale
760  trhs = finfact * this%uzfarea(icell) * this%celtop(icell) / range
761  thcof = finfact * this%uzfarea(icell) / range
762  end if
763  end subroutine rejfinf
764 
765  !> @brief Calculate groundwater discharge to land surface
766  !<
767  subroutine gwseep(this, icell, deriv, scale, hgwf, trhs, thcof, seep)
768  ! -- dummy
769  class(uzfcellgrouptype) :: this
770  integer(I4B), intent(in) :: icell
771  real(DP), intent(inout) :: deriv
772  real(DP), intent(inout) :: trhs
773  real(DP), intent(inout) :: thcof
774  real(DP), intent(inout) :: seep
775  real(DP), intent(out) :: scale
776  real(DP), intent(in) :: hgwf
777  ! -- local
778  real(DP) :: x, range, y, deriv1, d1, d2, Q
779  !
780  seep = dzero
781  deriv = dzero
782  deriv1 = dzero
783  d1 = dzero
784  d2 = dzero
785  scale = dzero
786  q = this%uzfarea(icell) * this%vks(icell)
787  range = this%surfdep(icell)
788  x = hgwf - this%celtop(icell)
789  call scubiclinear(x, range, deriv1, y)
790  scale = y
791  seep = scale * q * (hgwf - this%celtop(icell)) / range
792  trhs = scale * q * this%celtop(icell) / range
793  thcof = -scale * q / range
794  d1 = -deriv1 * q * x / range
795  d2 = -scale * q / range
796  deriv = d1 + d2
797  if (seep < dzero) then
798  seep = dzero
799  deriv = dzero
800  trhs = dzero
801  thcof = dzero
802  end if
803  end subroutine gwseep
804 
805  !> @brief Calculate gwf et using residual uzf pet
806  !<
807  subroutine simgwet(this, igwetflag, icell, hgwf, trhs, thcof, det)
808  ! -- dummy
809  class(uzfcellgrouptype) :: this
810  integer(I4B), intent(in) :: igwetflag
811  integer(I4B), intent(in) :: icell
812  real(DP), intent(in) :: hgwf
813  real(DP), intent(inout) :: trhs
814  real(DP), intent(inout) :: thcof
815  real(DP), intent(inout) :: det
816  ! -- local
817  real(DP) :: s, x, c, b, et
818  !
819  this%gwet(icell) = dzero
820  trhs = dzero
821  thcof = dzero
822  det = dzero
823  s = this%landtop(icell)
824  x = this%extdp(icell)
825  c = this%gwpet(icell)
826  b = this%celbot(icell)
827  if (b > hgwf) return
828  if (x < dem6) return
829  if (igwetflag == 1) then
830  et = etfunc_lin(s, x, c, det, trhs, thcof, hgwf, &
831  this%celtop(icell), this%celbot(icell))
832  else if (igwetflag == 2) then
833  et = etfunc_nlin(s, x, c, det, trhs, thcof, hgwf)
834  end if
835  ! this%gwet(icell) = et * this%uzfarea(icell)
836  trhs = trhs * this%uzfarea(icell)
837  thcof = thcof * this%uzfarea(icell)
838  this%gwet(icell) = trhs - (thcof * hgwf)
839  ! write(99,*)'in group', icell, this%gwet(icell)
840  end subroutine simgwet
841 
842  !> @brief Calculate recharge due to a rise in the gwf head
843  !<
844  subroutine uz_rise(this, icell, totfluxtot)
845  ! -- dummy
846  class(uzfcellgrouptype) :: this
847  integer(I4B), intent(in) :: icell
848  real(DP), intent(inout) :: totfluxtot
849  ! -- local
850  real(DP) :: fm1, fm2, d1
851  !
852  ! -- additional recharge from a rising water table
853  if (this%watab(icell) - this%watabold(icell) > dem30) then
854  d1 = this%celtop(icell) - this%watabold(icell)
855  fm1 = this%unsat_stor(icell, d1)
856  d1 = this%celtop(icell) - this%watab(icell)
857  fm2 = this%unsat_stor(icell, d1)
858  totfluxtot = totfluxtot + (fm1 - fm2)
859  end if
860  end subroutine uz_rise
861 
862  !> @brief Reset waves to default values at start of simulation
863  !<
864  subroutine setwaves(this, icell)
865  ! -- dummy
866  class(uzfcellgrouptype) :: this
867  ! -- local
868  integer(I4B), intent(in) :: icell
869  real(DP) :: bottom, top
870  integer(I4B) :: jk
871  real(DP) :: thick
872  !
873  ! -- initialize
874  this%totflux(icell) = dzero
875  this%nwavst(icell) = 1
876  this%uzdpst(:, icell) = dzero
877  thick = this%celtop(icell) - this%watab(icell)
878  do jk = 1, this%nwav(icell)
879  this%uzthst(jk, icell) = this%thtr(icell)
880  end do
881  !
882  ! -- initialize waves for first stress period
883  if (thick > dzero) then
884  this%uzdpst(1, icell) = thick
885  this%uzthst(1, icell) = this%thti(icell)
886  top = this%uzthst(1, icell) - this%thtr(icell)
887  if (top < dzero) top = dzero
888  bottom = this%thts(icell) - this%thtr(icell)
889  if (bottom < dzero) bottom = dzero
890  this%uzflst(1, icell) = this%vks(icell) * (top / bottom)**this%eps(icell)
891  if (this%uzthst(1, icell) < this%thtr(icell)) &
892  this%uzthst(1, icell) = this%thtr(icell)
893  !
894  ! -- calculate water stored in the unsaturated zone
895  if (top > dzero) then
896  this%uzspst(1, icell) = dzero
897  else
898  this%uzflst(1, icell) = dzero
899  this%uzspst(1, icell) = dzero
900  end if
901  !
902  ! no unsaturated zone
903  else
904  this%uzflst(1, icell) = dzero
905  this%uzdpst(1, icell) = dzero
906  this%uzspst(1, icell) = dzero
907  this%uzthst(1, icell) = this%thtr(icell)
908  end if
909  end subroutine
910 
911  !> @brief Prepare and route waves over time step
912  !<
913  subroutine routewaves(this, totfluxtot, delt, ietflag, icell, ierr)
914  ! -- dummy
915  class(uzfcellgrouptype) :: this
916  real(DP), intent(inout) :: totfluxtot
917  real(DP), intent(in) :: delt
918  integer(I4B), intent(in) :: ietflag
919  integer(I4B), intent(in) :: icell
920  integer(I4B), intent(inout) :: ierr
921  ! -- local
922  real(DP) :: thick, thickold
923  integer(I4B) :: idelt, iwav, ik
924  !
925  ! -- initialize
926  this%totflux(icell) = dzero
927  this%etact(icell) = dzero
928  thick = this%celtop(icell) - this%watab(icell)
929  thickold = this%celtop(icell) - this%watabold(icell)
930  !
931  ! -- no uz, clear waves
932  if (thickold < dzero) then
933  do iwav = 1, 6
934  this%uzthst(iwav, icell) = this%thtr(icell)
935  this%uzdpst(iwav, icell) = dzero
936  this%uzspst(iwav, icell) = dzero
937  this%uzflst(iwav, icell) = dzero
938  this%nwavst(icell) = 1
939  end do
940  end if
941  idelt = 1
942  do ik = 1, idelt
943  call this%uzflow(thick, thickold, delt, ietflag, icell, ierr)
944  if (ierr > 0) return
945  totfluxtot = totfluxtot + this%totflux(icell)
946  end do
947  end subroutine routewaves
948 
949  !> @brief Copy waves or shift waves in arrays
950  !<
951  subroutine wave_shift(this, this2, icell, icell2, shft, strt, stp, cntr)
952  ! -- dummy
953  class(uzfcellgrouptype) :: this
954  type(uzfcellgrouptype) :: this2
955  integer(I4B), intent(in) :: icell
956  integer(I4B), intent(in) :: icell2
957  integer(I4B), intent(in) :: shft
958  integer(I4B), intent(in) :: strt
959  integer(I4B), intent(in) :: stp
960  integer(I4B), intent(in) :: cntr
961  ! -- local
962  integer(I4B) :: j
963  !
964  ! -- copy waves from one uzf cell group to another
965  do j = strt, stp, cntr
966  this%uzthst(j, icell) = this2%uzthst(j + shft, icell2)
967  this%uzdpst(j, icell) = this2%uzdpst(j + shft, icell2)
968  this%uzflst(j, icell) = this2%uzflst(j + shft, icell2)
969  this%uzspst(j, icell) = this2%uzspst(j + shft, icell2)
970  end do
971  this%nwavst(icell) = this2%nwavst(icell2)
972  end subroutine
973 
974  !> @brief Method of Characteristics solution for kinematic wave equation
975  !<
976  subroutine uzflow(this, thick, thickold, delt, ietflag, icell, ierr)
977  ! -- dummy
978  class(uzfcellgrouptype) :: this
979  real(DP), intent(inout) :: thickold
980  real(DP), intent(inout) :: thick
981  real(DP), intent(in) :: delt
982  integer(I4B), intent(in) :: ietflag
983  integer(I4B), intent(in) :: icell
984  integer(I4B), intent(inout) :: ierr
985  ! -- local
986  real(DP) :: ffcheck, time, feps1, feps2
987  real(DP) :: thetadif, thetab, fluxb, oldsflx
988  integer(I4B) :: itrailflg, itester
989  !
990  time = dzero
991  this%totflux(icell) = dzero
992  itrailflg = 0
993  oldsflx = this%uzflst(this%nwavst(icell), icell)
994  call factors(feps1, feps2)
995  !
996  ! -- check for falling or rising water table
997  if ((thick - thickold) > feps1) then
998  thetadif = abs(this%uzthst(1, icell) - this%thtr(icell))
999  if (thetadif > dem6) then
1000  call this%wave_shift(this, icell, icell, -1, &
1001  this%nwavst(icell) + 1, 2, -1)
1002  if (this%uzdpst(2, icell) < dem30) &
1003  this%uzdpst(2, icell) = (this%ntrail(icell) + dtwo) * dem6
1004  if (this%uzthst(2, icell) > this%thtr(icell)) then
1005  this%uzspst(2, icell) = this%uzflst(2, icell) / &
1006  (this%uzthst(2, icell) - this%thtr(icell))
1007  else
1008  this%uzspst(2, icell) = dzero
1009  end if
1010  this%uzthst(1, icell) = this%thtr(icell)
1011  this%uzflst(1, icell) = dzero
1012  this%uzspst(1, icell) = dzero
1013  this%uzdpst(1, icell) = thick
1014  this%nwavst(icell) = this%nwavst(icell) + 1
1015  if (this%nwavst(icell) >= this%nwav(icell)) then
1016  ! -- too many waves error
1017  ierr = 1
1018  return
1019  end if
1020  else
1021  this%uzdpst(1, icell) = thick
1022  end if
1023  end if
1024  thetab = this%uzthst(1, icell)
1025  fluxb = this%uzflst(1, icell)
1026  this%totflux(icell) = dzero
1027  itester = 0
1028  ffcheck = (this%surflux(icell) - this%uzflst(this%nwavst(icell), icell))
1029  !
1030  ! -- increase new waves in infiltration changes
1031  if (ffcheck > feps2 .OR. ffcheck < -feps2) then
1032  this%nwavst(icell) = this%nwavst(icell) + 1
1033  if (this%nwavst(icell) >= this%nwav(icell)) then
1034  !
1035  ! -- too many waves error
1036  ierr = 1
1037  return
1038  end if
1039  else if (this%nwavst(icell) == 1) then
1040  itester = 1
1041  end if
1042  if (this%nwavst(icell) > 1) then
1043  if (ffcheck < -feps2) then
1044  call this%trailwav(icell, ierr)
1045  if (ierr > 0) return
1046  itrailflg = 1
1047  end if
1048  call this%leadwav(time, itester, itrailflg, thetab, fluxb, ffcheck, &
1049  feps2, delt, icell)
1050  end if
1051  if (itester == 1) then
1052  this%totflux(icell) = this%totflux(icell) + &
1053  (delt - time) * this%uzflst(1, icell)
1054  time = dzero
1055  itester = 0
1056  end if
1057  !
1058  ! -- simulate et
1059  if (ietflag > 0) call this%uzet(icell, delt, ietflag, ierr)
1060  if (ierr > 0) return
1061  end subroutine uzflow
1062 
1063  !> @brief Calculate unit specific tolerances
1064  !<
1065  subroutine factors(feps1, feps2)
1066  ! -- dummy
1067  real(DP), intent(out) :: feps1
1068  real(DP), intent(out) :: feps2
1069  real(DP) :: factor1
1070  real(DP) :: factor2
1071  !
1072  ! calculate constants for uzflow
1073  factor1 = done
1074  factor2 = done
1075  feps1 = dem9
1076  feps2 = dem9
1077  if (itmuni == 1) then
1078  factor1 = done / 86400.d0
1079  else if (itmuni == 2) then
1080  factor1 = done / 1440.d0
1081  else if (itmuni == 3) then
1082  factor1 = done / 24.0d0
1083  else if (itmuni == 5) then
1084  factor1 = 365.0d0
1085  end if
1086  factor2 = done / 0.3048
1087  feps1 = feps1 * factor1 * factor2
1088  feps2 = feps2 * factor1 * factor2
1089  end subroutine factors
1090 
1091  !> @brief Create and set trail waves
1092  !<
1093  subroutine trailwav(this, icell, ierr)
1094  ! -- dummy
1095  class(uzfcellgrouptype) :: this
1096  integer(I4B), intent(in) :: icell
1097  integer(I4B), intent(inout) :: ierr
1098  ! -- local
1099  real(DP) :: smoist, smoistinc, ftrail, eps_m1
1100  real(DP) :: thtsrinv
1101  real(DP) :: flux1, flux2, theta1, theta2
1102  real(DP) :: fnuminc
1103  integer(I4B) :: j, jj, jk, nwavstm1
1104  !
1105  ! -- initialize
1106  eps_m1 = dble(this%eps(icell)) - done
1107  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1108  nwavstm1 = this%nwavst(icell) - 1
1109  !
1110  ! -- initialize trailwaves
1111  smoist = (((this%surflux(icell) / this%vks(icell))** &
1112  (done / this%eps(icell))) * &
1113  (this%thts(icell) - this%thtr(icell))) + this%thtr(icell)
1114  if (this%uzthst(nwavstm1, icell) - smoist > dem9) then
1115  fnuminc = dzero
1116  do jk = 1, this%ntrail(icell)
1117  fnuminc = fnuminc + float(jk)
1118  end do
1119  smoistinc = (this%uzthst(nwavstm1, icell) - smoist) / (fnuminc - done)
1120  jj = this%ntrail(icell)
1121  ftrail = dble(this%ntrail(icell)) + done
1122  do j = this%nwavst(icell), this%nwavst(icell) + this%ntrail(icell) - 1
1123  if (j > this%nwav(icell)) then
1124  ! -- too many waves error
1125  ierr = 1
1126  return
1127  end if
1128  if (j > this%nwavst(icell)) then
1129  this%uzthst(j, icell) = this%uzthst(j - 1, icell) &
1130  - ((ftrail - float(jj)) * smoistinc)
1131  else
1132  this%uzthst(j, icell) = this%uzthst(j - 1, icell) - dem9
1133  end if
1134  jj = jj - 1
1135  if (this%uzthst(j, icell) <= this%thtr(icell) + dem9) &
1136  this%uzthst(j, icell) = this%thtr(icell) + dem9
1137  this%uzflst(j, icell) = &
1138  this%vks(icell) * (((this%uzthst(j, icell) - this%thtr(icell)) * &
1139  thtsrinv)**this%eps(icell))
1140  theta2 = this%uzthst(j - 1, icell)
1141  flux2 = this%uzflst(j - 1, icell)
1142  flux1 = this%uzflst(j, icell)
1143  theta1 = this%uzthst(j, icell)
1144  this%uzspst(j, icell) = leadspeed(theta1, theta2, flux1, flux2, &
1145  this%thts(icell), this%thtr(icell), &
1146  this%eps(icell), this%vks(icell))
1147  this%uzdpst(j, icell) = dzero
1148  if (j == this%nwavst(icell)) then
1149  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1150  (this%ntrail(icell) + 1) * dem9
1151  else
1152  this%uzdpst(j, icell) = this%uzdpst(j - 1, icell) - dem9
1153  end if
1154  end do
1155  this%nwavst(icell) = this%nwavst(icell) + this%ntrail(icell) - 1
1156  if (this%nwavst(icell) >= this%nwav(icell)) then
1157  ! -- too many waves error
1158  ierr = 1
1159  return
1160  end if
1161  else
1162  this%uzdpst(this%nwavst, icell) = dzero
1163  this%uzflst(this%nwavst, icell) = &
1164  this%vks(icell) * (((this%uzthst(this%nwavst, icell) - &
1165  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1166  this%uzthst(this%nwavst, icell) = smoist
1167  theta2 = this%uzthst(this%nwavst(icell) - 1, icell)
1168  flux2 = this%uzflst(this%nwavst(icell) - 1, icell)
1169  flux1 = this%uzflst(this%nwavst(icell), icell)
1170  theta1 = this%uzthst(this%nwavst(icell), icell)
1171  this%uzspst(this%nwavst(icell), icell) = &
1172  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1173  this%thtr(icell), this%eps(icell), this%vks(icell))
1174  end if
1175  end subroutine trailwav
1176 
1177  !> @brief Create a lead wave and route over time step
1178  !<
1179  subroutine leadwav(this, time, itester, itrailflg, thetab, fluxb, &
1180  ffcheck, feps2, delt, icell)
1181  ! -- dummy
1182  class(uzfcellgrouptype) :: this
1183  real(DP), intent(inout) :: thetab
1184  real(DP), intent(inout) :: fluxb
1185  real(DP), intent(in) :: feps2
1186  real(DP), intent(inout) :: time
1187  integer(I4B), intent(inout) :: itester
1188  integer(I4B), intent(inout) :: itrailflg
1189  real(DP), intent(inout) :: ffcheck
1190  real(DP), intent(in) :: delt
1191  integer(I4B), intent(in) :: icell
1192  ! -- local
1193  real(DP) :: bottomtime, shortest, fcheck
1194  real(DP) :: eps_m1, timenew, bottom, timedt
1195  real(DP) :: thtsrinv, diff, fluxhld2
1196  real(DP) :: flux1, flux2, theta1, theta2, ftest
1197  real(DP), allocatable, dimension(:) :: checktime
1198  integer(I4B) :: iflx, iremove, j, l
1199  integer(I4B) :: nwavp1, jshort
1200  integer(I4B), allocatable, dimension(:) :: more
1201  !
1202  allocate (checktime(this%nwavst(icell)))
1203  allocate (more(this%nwavst(icell)))
1204  ftest = dzero
1205  eps_m1 = dble(this%eps(icell)) - done
1206  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1207  !
1208  ! -- initialize new wave
1209  if (itrailflg == 0) then
1210  if (ffcheck > feps2) then
1211  this%uzflst(this%nwavst(icell), icell) = this%surflux(icell)
1212  if (this%uzflst(this%nwavst(icell), icell) < dem30) &
1213  this%uzflst(this%nwavst(icell), icell) = dzero
1214  this%uzthst(this%nwavst(icell), icell) = &
1215  (((this%uzflst(this%nwavst(icell), icell) / this%vks(icell))** &
1216  (done / this%eps(icell))) * (this%thts(icell) - this%thtr(icell))) &
1217  + this%thtr(icell)
1218  theta2 = this%uzthst(this%nwavst(icell), icell)
1219  flux2 = this%uzflst(this%nwavst(icell), icell)
1220  flux1 = this%uzflst(this%nwavst(icell) - 1, icell)
1221  theta1 = this%uzthst(this%nwavst(icell) - 1, icell)
1222  this%uzspst(this%nwavst(icell), icell) = &
1223  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1224  this%thtr(icell), this%eps(icell), this%vks(icell))
1225  this%uzdpst(this%nwavst(icell), icell) = dzero
1226  end if
1227  end if
1228  !
1229  ! -- route all waves and interception of waves over times step
1230  diff = done
1231  timedt = dzero
1232  iflx = 0
1233  fluxhld2 = this%uzflst(1, icell)
1234  if (this%nwavst(icell) == 0) itester = 1
1235  if (itester /= 1) then
1236  do while (diff > dem6)
1237  nwavp1 = this%nwavst(icell) + 1
1238  timedt = delt - time
1239  do j = 1, this%nwavst(icell)
1240  checktime(j) = dep20
1241  more(j) = 0
1242  end do
1243  shortest = timedt
1244  if (this%nwavst(icell) > 2) then
1245  j = 2
1246  !
1247  ! -- calculate time until wave overtakes wave ahead
1248  nwavp1 = this%nwavst(icell) + 1
1249  do while (j < nwavp1)
1250  ftest = this%uzspst(j - 1, icell) - this%uzspst(j, icell)
1251  if (abs(ftest) > dem30) then
1252  checktime(j) = (this%uzdpst(j, icell) - &
1253  this%uzdpst(j - 1, icell)) / ftest
1254  if (checktime(j) < dem30) checktime(j) = dep20
1255  end if
1256  j = j + 1
1257  end do
1258  end if
1259  !
1260  ! - calc time until wave reaches bottom of cell
1261  bottomtime = dep20
1262  if (this%nwavst(icell) > 1) then
1263  if (this%uzspst(2, icell) > dzero) then
1264  bottom = this%uzspst(2, icell)
1265  if (bottom < dem15) bottom = dem15
1266  bottomtime = (this%uzdpst(1, icell) - this%uzdpst(2, icell)) / bottom
1267  if (bottomtime < dzero) bottomtime = dem12
1268  end if
1269  end if
1270  !
1271  ! -- calc time for wave interception
1272  jshort = 0
1273  do j = this%nwavst(icell), 3, -1
1274  if (shortest - checktime(j) > -dem9) then
1275  more(j) = 1
1276  jshort = j
1277  shortest = checktime(j)
1278  end if
1279  end do
1280  do j = 3, this%nwavst(icell)
1281  if (shortest - checktime(j) < dem9) then
1282  if (j /= jshort) more(j) = 0
1283  end if
1284  end do
1285  !
1286  ! -- what happens first, waves hits bottom or interception
1287  iremove = 0
1288  timenew = time
1289  fcheck = (time + shortest) - delt
1290  if (shortest < dem7) fcheck = -done
1291  if (bottomtime < shortest .AND. time + bottomtime < delt) then
1292  j = 2
1293  do while (j < nwavp1)
1294  !
1295  ! -- route waves
1296  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1297  this%uzspst(j, icell) * bottomtime
1298  j = j + 1
1299  end do
1300  fluxb = this%uzflst(2, icell)
1301  thetab = this%uzthst(2, icell)
1302  iflx = 1
1303  call this%wave_shift(this, icell, icell, 1, 1, &
1304  this%nwavst(icell) - 1, 1)
1305  iremove = 1
1306  timenew = time + bottomtime
1307  this%uzspst(1, icell) = dzero
1308  !
1309  ! -- do waves intercept before end of time step
1310  else if (fcheck < dzero .AND. this%nwavst(icell) > 2) then
1311  j = 2
1312  do while (j < nwavp1)
1313  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1314  this%uzspst(j, icell) * shortest
1315  j = j + 1
1316  end do
1317  !
1318  ! -- combine waves that intercept, remove a wave
1319  j = 3
1320  l = j
1321  do while (j < this%nwavst(icell) + 1)
1322  if (more(j) == 1) then
1323  l = j
1324  theta2 = this%uzthst(j, icell)
1325  flux2 = this%uzflst(j, icell)
1326  if (j == 3) then
1327  flux1 = fluxb
1328  theta1 = thetab
1329  else
1330  flux1 = this%uzflst(j - 2, icell)
1331  theta1 = this%uzthst(j - 2, icell)
1332  end if
1333  this%uzspst(j, icell) = leadspeed(theta1, theta2, flux1, flux2, &
1334  this%thts(icell), &
1335  this%thtr(icell), &
1336  this%eps(icell), this%vks(icell))
1337  !
1338  ! -- update waves
1339  call this%wave_shift(this, icell, icell, 1, l - 1, &
1340  this%nwavst(icell) - 1, 1)
1341  l = this%nwavst(icell) + 1
1342  iremove = iremove + 1
1343  end if
1344  j = j + 1
1345  end do
1346  timenew = timenew + shortest
1347  !
1348  ! -- calc. total flux to bottom during remaining time in step
1349  else
1350  j = 2
1351  do while (j < nwavp1)
1352  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1353  this%uzspst(j, icell) * timedt
1354  j = j + 1
1355  end do
1356  timenew = delt
1357  end if
1358  this%totflux(icell) = this%totflux(icell) + fluxhld2 * (timenew - time)
1359  if (iflx == 1) then
1360  fluxhld2 = this%uzflst(1, icell)
1361  iflx = 0
1362  end if
1363  !
1364  ! -- remove dead waves
1365  this%nwavst(icell) = this%nwavst(icell) - iremove
1366  time = timenew
1367  diff = delt - time
1368  if (this%nwavst(icell) == 1) then
1369  itester = 1
1370  exit
1371  end if
1372  end do
1373  end if
1374  deallocate (checktime)
1375  deallocate (more)
1376  end subroutine leadwav
1377 
1378  !> @brief Calculates waves speed from dflux/dtheta
1379  !<
1380  function leadspeed(theta1, theta2, flux1, flux2, thts, thtr, eps, vks)
1381  ! -- Return
1382  real(dp) :: leadspeed
1383  ! -- dummy
1384  real(dp), intent(in) :: theta1
1385  real(dp), intent(in) :: theta2
1386  real(dp), intent(in) :: flux1
1387  real(dp), intent(inout) :: flux2
1388  real(dp), intent(in) :: thts
1389  real(dp), intent(in) :: thtr
1390  real(dp), intent(in) :: eps
1391  real(dp), intent(in) :: vks
1392  ! -- local
1393  real(dp) :: comp1, comp2, thsrinv, epsfksths
1394  real(dp) :: eps_m1, fhold, comp3
1395  !
1396  eps_m1 = eps - done
1397  thsrinv = done / (thts - thtr)
1398  epsfksths = eps * vks * thsrinv
1399  comp1 = theta2 - theta1
1400  comp2 = abs(flux2 - flux1)
1401  comp3 = theta1 - thtr
1402  if (comp2 < dem15) flux2 = flux1 + dem15
1403  if (abs(comp1) < dem30) then
1404  if (comp3 > dem30) fhold = (comp3 * thsrinv)**eps
1405  if (fhold < dem30) fhold = dem30
1406  leadspeed = epsfksths * (fhold**eps_m1)
1407  else
1408  leadspeed = (flux2 - flux1) / (theta2 - theta1)
1409  end if
1410  if (leadspeed < dem30) leadspeed = dem30
1411  end function leadspeed
1412 
1413  !> @brief Sums up mobile water over depth interval
1414  !<
1415  function unsat_stor(this, icell, d1)
1416  ! -- Return
1417  real(dp) :: unsat_stor
1418  ! -- dummy
1419  class(uzfcellgrouptype) :: this
1420  integer(I4B), intent(in) :: icell
1421  real(dp), intent(inout) :: d1
1422  ! -- local
1423  real(dp) :: fm
1424  integer(I4B) :: j, k, nwavm1, jj
1425  !
1426  fm = dzero
1427  j = this%nwavst(icell) + 1
1428  k = this%nwavst(icell)
1429  nwavm1 = k - 1
1430  if (d1 > this%uzdpst(1, icell)) d1 = this%uzdpst(1, icell)
1431  !
1432  ! -- find deepest wave above depth d1, counter held as j
1433  do while (k > 0)
1434  if (this%uzdpst(k, icell) - d1 < -dem30) j = k
1435  k = k - 1
1436  end do
1437  if (j > this%nwavst(icell)) then
1438  fm = fm + (this%uzthst(this%nwavst(icell), icell) - this%thtr(icell)) * d1
1439  elseif (this%nwavst(icell) > 1) then
1440  if (j > 1) then
1441  fm = fm + (this%uzthst(j - 1, icell) - this%thtr(icell)) &
1442  * (d1 - this%uzdpst(j, icell))
1443  end if
1444  do jj = j, nwavm1
1445  fm = fm + (this%uzthst(jj, icell) - this%thtr(icell)) &
1446  * (this%uzdpst(jj, icell) &
1447  - this%uzdpst(jj + 1, icell))
1448  end do
1449  fm = fm + (this%uzthst(this%nwavst(icell), icell) - this%thtr(icell)) &
1450  * (this%uzdpst(this%nwavst(icell), icell))
1451  else
1452  fm = fm + (this%uzthst(1, icell) - this%thtr(icell)) * d1
1453  end if
1454  unsat_stor = fm
1455  end function unsat_stor
1456 
1457  !> @brief Update to new state of uz at end of time step
1458  !<
1459  subroutine update_wav(this, icell, delt, iss, itest)
1460  ! -- dummy
1461  class(uzfcellgrouptype) :: this
1462  integer(I4B), intent(in) :: icell
1463  integer(I4B), intent(in) :: itest
1464  integer(I4B), intent(in) :: iss
1465  real(DP), intent(in) :: delt
1466  ! -- local
1467  real(DP) :: bot, depthsave, top
1468  real(DP) :: thick, thtsrinv
1469  integer(I4B) :: nwavhld, k, j
1470  !
1471  bot = this%watab(icell)
1472  top = this%celtop(icell)
1473  thick = top - bot
1474  nwavhld = this%nwavst(icell)
1475  if (itest == 1) then
1476  this%uzflst(1, icell) = dzero
1477  this%uzthst(1, icell) = this%thtr(icell)
1478  return
1479  end if
1480  if (iss == 1) then
1481  if (this%thts(icell) - this%thtr(icell) < dem7) then
1482  thtsrinv = done / dem7
1483  else
1484  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1485  end if
1486  this%totflux(icell) = this%surflux(icell) * delt
1487  this%watabold(icell) = this%watab(icell)
1488  this%uzthst(1, icell) = this%thti(icell)
1489  this%uzflst(1, icell) = &
1490  this%vks(icell) * (((this%uzthst(1, icell) - this%thtr(icell)) &
1491  * thtsrinv)**this%eps(icell))
1492  this%uzdpst(1, icell) = thick
1493  this%uzspst(1, icell) = thick
1494  this%nwavst(icell) = 1
1495  else
1496  !
1497  ! -- water table rises through waves
1498  if (this%watab(icell) - this%watabold(icell) > dem30) then
1499  depthsave = this%uzdpst(1, icell)
1500  j = 0
1501  k = this%nwavst(icell)
1502  do while (k > 0)
1503  if (this%uzdpst(k, icell) - thick < -dem30) j = k
1504  k = k - 1
1505  end do
1506  this%uzdpst(1, icell) = thick
1507  if (j > 1) then
1508  this%uzspst(1, icell) = dzero
1509  this%nwavst(icell) = this%nwavst(icell) - j + 2
1510  this%uzthst(1, icell) = this%uzthst(j - 1, icell)
1511  this%uzflst(1, icell) = this%uzflst(j - 1, icell)
1512  if (j > 2) call this%wave_shift(this, icell, icell, j - 2, 2, &
1513  nwavhld - (j - 2), 1)
1514  elseif (j == 0) then
1515  this%uzspst(1, icell) = dzero
1516  this%uzthst(1, icell) = this%uzthst(this%nwavst(icell), icell)
1517  this%uzflst(1, icell) = this%uzflst(this%nwavst(icell), icell)
1518  this%nwavst(icell) = 1
1519  end if
1520  end if
1521  !
1522  ! -- calculate new unsat. storage
1523  if (thick <= dzero) then
1524  this%uzspst(1, icell) = dzero
1525  this%nwavst(icell) = 1
1526  this%uzthst(1, icell) = this%thtr(icell)
1527  this%uzflst(1, icell) = dzero
1528  end if
1529  this%watabold(icell) = this%watab(icell)
1530  end if
1531  end subroutine update_wav
1532 
1533  !> @brief Remove water from uz due to et
1534  !<
1535  subroutine uzet(this, icell, delt, ietflag, ierr)
1536  ! -- dummy
1537  class(uzfcellgrouptype) :: this
1538  integer(I4B), intent(in) :: icell
1539  real(DP), intent(in) :: delt
1540  integer(I4B), intent(in) :: ietflag
1541  integer(I4B), intent(inout) :: ierr
1542  ! -- local
1543  type(uzfcellgrouptype) :: uzfktemp
1544  real(DP) :: diff
1545  real(DP) :: thetaout
1546  real(DP) :: fm
1547  real(DP) :: st
1548  real(DP) :: thtsrinv
1549  real(DP) :: epsfksthts
1550  real(DP) :: fmp
1551  real(DP) :: fktho
1552  real(DP) :: theta1
1553  real(DP) :: theta2
1554  real(DP) :: flux1
1555  real(DP) :: flux2
1556  real(DP) :: hcap
1557  real(DP) :: factor
1558  real(DP) :: tho
1559  real(DP) :: depth
1560  real(DP) :: extwc1
1561  real(DP) :: petsub
1562  integer(I4B) :: i
1563  integer(I4B) :: j
1564  integer(I4B) :: jhold
1565  integer(I4B) :: jk
1566  integer(I4B) :: kj
1567  integer(I4B) :: kk
1568  integer(I4B) :: numadd
1569  integer(I4B) :: k
1570  integer(I4B) :: nwv
1571  integer(I4B) :: itest
1572  !
1573  ! -- initialize
1574  this%etact(icell) = dzero
1575  if (this%extdpuz(icell) < dem7) return
1576  petsub = this%rootact(icell) * this%pet(icell) * &
1577  this%extdpuz(icell) / this%extdp(icell)
1578  thetaout = delt * petsub / this%extdp(icell)
1579  if (ietflag == 1) thetaout = delt * this%pet(icell) / this%extdp(icell)
1580  if (thetaout < dem10) return
1581  depth = this%uzdpst(1, icell)
1582  st = this%unsat_stor(icell, depth)
1583  if (st < dem4) return
1584  !
1585  ! -- allocate temporary wave storage.
1586  nwv = this%nwavst(icell)
1587  itest = 0
1588  call uzfktemp%init(1, nwv)
1589  !
1590  ! store original wave characteristics
1591  call uzfktemp%wave_shift(this, 1, icell, 0, 1, nwv, 1)
1592  factor = done
1593  this%etact(icell) = dzero
1594  if (this%thts(icell) - this%thtr(icell) < dem7) then
1595  thtsrinv = 1.0 / dem7
1596  else
1597  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1598  end if
1599  epsfksthts = this%eps(icell) * this%vks(icell) * thtsrinv
1600  this%etact(icell) = dzero
1601  fmp = dzero
1602  extwc1 = this%extwc(icell) - this%thtr(icell)
1603  if (extwc1 < dem6) extwc1 = dem7
1604  numadd = 0
1605  fm = st
1606  k = 0
1607  !
1608  ! -- loop for reducing aet to pet when et is head dependent
1609  do while (itest == 0)
1610  k = k + 1
1611  if (k > 1 .AND. abs(fmp - petsub) > dem5 * petsub) then
1612  factor = factor / (fm / petsub)
1613  end if
1614  !
1615  ! -- one wave shallower than extdp
1616  if (this%nwavst(icell) == 1 .AND. &
1617  this%uzdpst(1, icell) <= this%extdpuz(icell)) then
1618  if (ietflag == 2) then
1619  tho = this%uzthst(1, icell)
1620  fktho = this%uzflst(1, icell)
1621  hcap = this%caph(icell, tho)
1622  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1623  end if
1624  if ((this%uzthst(1, icell) - thetaout) > this%thtr(icell) + extwc1) then
1625  this%uzthst(1, icell) = this%uzthst(1, icell) - thetaout
1626  this%uzflst(1, icell) = &
1627  this%vks(icell) * (((this%uzthst(1, icell) - &
1628  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1629  else if (this%uzthst(1, icell) > this%thtr(icell) + extwc1) then
1630  this%uzthst(1, icell) = this%thtr(icell) + extwc1
1631  this%uzflst(1, icell) = &
1632  this%vks(icell) * (((this%uzthst(1, icell) - &
1633  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1634  end if
1635  !
1636  ! -- all waves shallower than extinction depth
1637  else if (this%nwavst(icell) > 1 .AND. &
1638  this%uzdpst(this%nwavst(icell), icell) > this%extdpuz(icell)) then
1639  if (ietflag == 2) then
1640  tho = this%uzthst(this%nwavst(icell), icell)
1641  fktho = this%uzflst(this%nwavst(icell), icell)
1642  hcap = this%caph(icell, tho)
1643  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1644  end if
1645  if (this%uzthst(this%nwavst(icell), icell) - thetaout > &
1646  this%thtr(icell) + extwc1) then
1647  this%uzthst(this%nwavst(icell) + 1, icell) = &
1648  this%uzthst(this%nwavst(icell), icell) - thetaout
1649  numadd = 1
1650  else if (this%uzthst(this%nwavst(icell), icell) > &
1651  this%thtr(icell) + extwc1) then
1652  this%uzthst(this%nwavst(icell) + 1, icell) = this%thtr(icell) + extwc1
1653  numadd = 1
1654  end if
1655  if (numadd == 1) then
1656  this%uzflst(this%nwavst(icell) + 1, icell) = &
1657  this%vks(icell) * &
1658  (((this%uzthst(this%nwavst(icell) + 1, icell) - &
1659  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1660  theta2 = this%uzthst(this%nwavst(icell) + 1, icell)
1661  flux2 = this%uzflst(this%nwavst(icell) + 1, icell)
1662  flux1 = this%uzflst(this%nwavst(icell), icell)
1663  theta1 = this%uzthst(this%nwavst(icell), icell)
1664  this%uzspst(this%nwavst(icell) + 1, icell) = &
1665  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1666  this%thtr(icell), this%eps(icell), this%vks(icell))
1667  this%uzdpst(this%nwavst(icell) + 1, icell) = this%extdpuz(icell)
1668  this%nwavst(icell) = this%nwavst(icell) + 1
1669  if (this%nwavst(icell) > this%nwav(icell)) then
1670  !
1671  ! -- too many waves error, deallocate temp arrays and return
1672  ierr = 1
1673  goto 500
1674  end if
1675  else
1676  numadd = 0
1677  end if
1678  !
1679  ! -- one wave below extinction depth
1680  else if (this%nwavst(icell) == 1) then
1681  if (ietflag == 2) then
1682  tho = this%uzthst(1, icell)
1683  fktho = this%uzflst(1, icell)
1684  hcap = this%caph(icell, tho)
1685  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1686  end if
1687  if ((this%uzthst(1, icell) - thetaout) > this%thtr(icell) + extwc1) then
1688  if (thetaout > dem30) then
1689  this%uzthst(2, icell) = this%uzthst(1, icell) - thetaout
1690  this%uzflst(2, icell) = &
1691  this%vks(icell) * (((this%uzthst(2, icell) - this%thtr(icell)) * &
1692  thtsrinv)**this%eps(icell))
1693  this%uzdpst(2, icell) = this%extdpuz(icell)
1694  theta2 = this%uzthst(2, icell)
1695  flux2 = this%uzflst(2, icell)
1696  flux1 = this%uzflst(1, icell)
1697  theta1 = this%uzthst(1, icell)
1698  this%uzspst(2, icell) = &
1699  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1700  this%thtr(icell), this%eps(icell), this%vks(icell))
1701  this%nwavst(icell) = this%nwavst(icell) + 1
1702  if (this%nwavst(icell) > this%nwav(icell)) then
1703  !
1704  ! -- too many waves error
1705  ierr = 1
1706  goto 500
1707  end if
1708  end if
1709  else if (this%uzthst(1, icell) > this%thtr(icell) + extwc1) then
1710  if (thetaout > dem30) then
1711  this%uzthst(2, icell) = this%thtr(icell) + extwc1
1712  this%uzflst(2, icell) = &
1713  this%vks(icell) * (((this%uzthst(2, icell) - &
1714  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1715  this%uzdpst(2, icell) = this%extdpuz(icell)
1716  theta2 = this%uzthst(2, icell)
1717  flux2 = this%uzflst(2, icell)
1718  flux1 = this%uzflst(1, icell)
1719  theta1 = this%uzthst(1, icell)
1720  this%uzspst(2, icell) = &
1721  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1722  this%thtr(icell), this%eps(icell), this%vks(icell))
1723  this%nwavst(icell) = this%nwavst(icell) + 1
1724  if (this%nwavst(icell) > this%nwav(icell)) then
1725  !
1726  ! -- too many waves error
1727  ierr = 1
1728  goto 500
1729  end if
1730  end if
1731  end if
1732  else
1733  !
1734  ! -- extinction depth splits waves
1735  if (this%uzdpst(1, icell) - this%extdpuz(icell) > dem7) then
1736  j = 2
1737  jk = 0
1738  !
1739  ! -- locate extinction depth between waves
1740  do while (jk == 0)
1741  diff = this%uzdpst(j, icell) - this%extdpuz(icell)
1742  if (diff > dzero) then
1743  j = j + 1
1744  else
1745  jk = 1
1746  end if
1747  end do
1748  kk = j
1749  if (this%uzthst(j, icell) > this%thtr(icell) + extwc1) then
1750  !
1751  ! -- create a wave at extinction depth
1752  if (abs(diff) > dem5) then
1753  call this%wave_shift(this, icell, icell, -1, &
1754  this%nwavst(icell) + 1, j, -1)
1755  this%uzdpst(j, icell) = this%extdpuz(icell)
1756  this%nwavst(icell) = this%nwavst(icell) + 1
1757  if (this%nwavst(icell) > this%nwav(icell)) then
1758  !
1759  ! -- too many waves error
1760  ierr = 1
1761  goto 500
1762  end if
1763  end if
1764  kk = j
1765  else
1766  jhold = this%nwavst(icell)
1767  i = j + 1
1768  do while (i < this%nwavst(icell))
1769  if (this%uzthst(i, icell) > this%thtr(icell) + extwc1) then
1770  jhold = i
1771  i = this%nwavst(icell) + 1
1772  end if
1773  i = i + 1
1774  end do
1775  j = jhold
1776  kk = jhold
1777  end if
1778  else
1779  kk = 1
1780  end if
1781  !
1782  ! -- all waves above extinction depth
1783  do while (kk <= this%nwavst(icell))
1784  if (ietflag == 2) then
1785  tho = this%uzthst(kk, icell)
1786  fktho = this%uzflst(kk, icell)
1787  hcap = this%caph(icell, tho)
1788  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1789  end if
1790  if (this%uzthst(kk, icell) > this%thtr(icell) + extwc1) then
1791  if (this%uzthst(kk, icell) - thetaout > &
1792  this%thtr(icell) + extwc1) then
1793  this%uzthst(kk, icell) = this%uzthst(kk, icell) - thetaout
1794  else if (this%uzthst(kk, icell) > this%thtr(icell) + extwc1) then
1795  this%uzthst(kk, icell) = this%thtr(icell) + extwc1
1796  end if
1797  if (kk == 1) then
1798  this%uzflst(kk, icell) = &
1799  this%vks(icell) * &
1800  (((this%uzthst(kk, icell) - &
1801  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1802  end if
1803  if (kk > 1) then
1804  flux1 = &
1805  this%vks(icell) * ((this%uzthst(kk - 1, icell) - &
1806  this%thtr(icell)) * thtsrinv)**this%eps(icell)
1807  flux2 = &
1808  this%vks(icell) * ((this%uzthst(kk, icell) - &
1809  this%thtr(icell)) * thtsrinv)**this%eps(icell)
1810  this%uzflst(kk, icell) = flux2
1811  theta2 = this%uzthst(kk, icell)
1812  theta1 = this%uzthst(kk - 1, icell)
1813  this%uzspst(kk, icell) = leadspeed(theta1, theta2, flux1, flux2, &
1814  this%thts(icell), &
1815  this%thtr(icell), &
1816  this%eps(icell), this%vks(icell))
1817  end if
1818  end if
1819  kk = kk + 1
1820  end do
1821  end if
1822  !
1823  ! -- calculate aet
1824  kj = 1
1825  do while (kj <= this%nwavst(icell) - 1)
1826  if (abs(this%uzthst(kj, icell) - this%uzthst(kj + 1, icell)) < dem6) then
1827  call this%wave_shift(this, icell, icell, 1, kj + 1, &
1828  this%nwavst(icell) - 1, 1)
1829  kj = kj - 1
1830  this%nwavst(icell) = this%nwavst(icell) - 1
1831  end if
1832  kj = kj + 1
1833  end do
1834  depth = this%uzdpst(1, icell)
1835  fm = this%unsat_stor(icell, depth)
1836  this%etact(icell) = st - fm
1837  fm = this%etact(icell) / delt
1838  if (this%etact(icell) < dzero) then
1839  call this%wave_shift(uzfktemp, icell, 1, 0, 1, nwv, 1)
1840  this%nwavst(icell) = nwv
1841  this%etact(icell) = dzero
1842  elseif (petsub - fm < -dem15 .AND. ietflag == 2) then
1843  !
1844  ! -- aet greater than pet, reset and try again
1845  call this%wave_shift(uzfktemp, icell, 1, 0, 1, nwv, 1)
1846  this%nwavst(icell) = nwv
1847  this%etact(icell) = dzero
1848  else
1849  itest = 1
1850  end if
1851  !
1852  ! -- end aet-pet loop for head dependent et
1853  fmp = fm
1854  if (k > 100) then
1855  itest = 1
1856  elseif (ietflag < 2) then
1857  fmp = petsub
1858  itest = 1
1859  end if
1860  end do
1861 500 continue
1862  !
1863  ! -- deallocate temporary worker
1864  call uzfktemp%dealloc()
1865  end subroutine uzet
1866 
1867  !> @brief Calculate capillary pressure head from B-C equation
1868  !<
1869  function caph(this, icell, tho)
1870  ! -- dummy
1871  class(uzfcellgrouptype) :: this
1872  integer(I4B), intent(in) :: icell
1873  real(dp), intent(in) :: tho
1874  ! -- local
1875  real(dp) :: caph, lambda, star
1876  !
1877  caph = -dem6
1878  star = (tho - this%thtr(icell)) / (this%thts(icell) - this%thtr(icell))
1879  if (star < dem15) star = dem15
1880  lambda = dtwo / (this%eps(icell) - dthree)
1881  if (star > dem15) then
1882  if (tho - this%thts(icell) < dem15) then
1883  caph = this%ha(icell) * star**(-done / lambda)
1884  else
1885  caph = dzero
1886  end if
1887  end if
1888  end function caph
1889 
1890  !> @brief Calculate capillary pressure-based uz et
1891  function rate_et_z(this, icell, factor, fktho, h)
1892  ! -- Return
1893  real(dp) :: rate_et_z
1894  ! -- dummy
1895  class(uzfcellgrouptype) :: this
1896  integer(I4B), intent(in) :: icell
1897  real(dp), intent(in) :: factor, fktho, h
1898  !
1899  rate_et_z = factor * fktho * (h - this%hroot(icell))
1900  if (rate_et_z < dzero) rate_et_z = dzero
1901  end function rate_et_z
1902 
1903  !> @brief Determine the water content at a specific depth
1904  !!
1905  !! Because UZF-calculated waves are internal to UZF objects, different water
1906  !! contents exists at different depths.
1907  !<
1908  function get_water_content_at_depth(this, icell, depth) result(theta_at_depth)
1909  ! -- dummy
1910  class(uzfcellgrouptype) :: this
1911  integer(I4B), intent(in) :: icell !< uzf cell containing depth
1912  real(dp), intent(in) :: depth !< depth within the cell
1913  ! -- return
1914  real(dp) :: theta_at_depth
1915  ! -- local
1916  real(dp) :: d1
1917  real(dp) :: d2
1918  real(dp) :: f1
1919  real(dp) :: f2
1920  !
1921  if (this%watab(icell) < this%celtop(icell)) then
1922  if (this%celtop(icell) - depth > this%watab(icell)) then
1923  d1 = depth - dem3
1924  d2 = depth + dem3
1925  f1 = this%unsat_stor(icell, d1)
1926  f2 = this%unsat_stor(icell, d2)
1927  theta_at_depth = this%thtr(icell) + (f2 - f1) / (d2 - d1)
1928  else
1929  theta_at_depth = this%thts(icell)
1930  end if
1931  else
1932  theta_at_depth = this%thts(icell)
1933  end if
1934  end function get_water_content_at_depth
1935 
1936  !> @brief Calculate and return the cell-based water content value
1937  !<
1938  function get_wcnew(this, icell) result(watercontent)
1939  ! -- dummy
1940  class(uzfcellgrouptype) :: this
1941  integer(I4B), intent(in) :: icell !< uzf cell containing depth
1942  ! -- return
1943  real(dp) :: watercontent
1944  ! -- local
1945  real(dp) :: top
1946  real(dp) :: bot
1947  real(dp) :: theta_r
1948  real(dp) :: thk
1949  real(dp) :: hgwf
1950  real(dp) :: fm
1951  real(dp) :: d
1952  !
1953  hgwf = this%watab(icell)
1954  top = this%celtop(icell)
1955  bot = this%celbot(icell)
1956  thk = top - max(bot, hgwf)
1957  if (thk > dzero) then
1958  theta_r = this%thtr(icell)
1959  d = thk
1960  fm = this%unsat_stor(icell, d)
1961  watercontent = fm / thk
1962  watercontent = watercontent + theta_r
1963  else
1964  watercontent = dzero
1965  end if
1966  end function get_wcnew
1967 
1968 end module uzfcellgroupmodule
init
subroutine init()
Definition: GridSorting.f90:24
memorymanagermodule::mem_allocate
Definition: MemoryManager.f90:54
memorymanagermodule::mem_deallocate
Definition: MemoryManager.f90:122
constantsmodule
This module contains simulation constants.
Definition: Constants.f90:9
constantsmodule::dem12
real(dp), parameter dem12
real constant 1e-12
Definition: Constants.f90:114
constantsmodule::dem20
real(dp), parameter dem20
real constant 1e-20
Definition: Constants.f90:117
constantsmodule::dem10
real(dp), parameter dem10
real constant 1e-10
Definition: Constants.f90:113
constantsmodule::dem7
real(dp), parameter dem7
real constant 1e-7
Definition: Constants.f90:110
constantsmodule::dep20
real(dp), parameter dep20
real constant 1e20
Definition: Constants.f90:91
constantsmodule::dem14
real(dp), parameter dem14
real constant 1e-14
Definition: Constants.f90:115
constantsmodule::dhalf
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
constantsmodule::dem3
real(dp), parameter dem3
real constant 1e-3
Definition: Constants.f90:106
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::dem5
real(dp), parameter dem5
real constant 1e-5
Definition: Constants.f90:108
constantsmodule::dem9
real(dp), parameter dem9
real constant 1e-9
Definition: Constants.f90:112
constantsmodule::dem15
real(dp), parameter dem15
real constant 1e-15
Definition: Constants.f90:116
constantsmodule::dtwo
real(dp), parameter dtwo
real constant 2
Definition: Constants.f90:79
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
kindmodule
This module defines variable data types.
Definition: kind.f90:8
smoothingmodule::slinear
subroutine slinear(x, range, dydx, y)
@ brief sLinear
Definition: SmoothingFunctions.f90:111
smoothingmodule::scubiclinear
subroutine scubiclinear(x, range, dydx, y)
@ brief sCubicLinear
Definition: SmoothingFunctions.f90:46
smoothingmodule::sscurve
subroutine sscurve(x, range, dydx, y)
@ brief SCurve
Definition: SmoothingFunctions.f90:15
tdismodule::itmuni
integer(i4b), pointer, public itmuni
flag indicating time units
Definition: tdis.f90:22
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
uzfcellgroupmodule::update_wav
subroutine update_wav(this, icell, delt, iss, itest)
Update to new state of uz at end of time step.
Definition: UzfCellGroup.f90:1460
uzfcellgroupmodule::factors
subroutine factors(feps1, feps2)
Calculate unit specific tolerances.
Definition: UzfCellGroup.f90:1066
uzfcellgroupmodule::uzflow
subroutine uzflow(this, thick, thickold, delt, ietflag, icell, ierr)
Method of Characteristics solution for kinematic wave equation.
Definition: UzfCellGroup.f90:977
uzfcellgroupmodule::advance
subroutine advance(this, icell)
Set variables to advance to new time step. nothing yet.
Definition: UzfCellGroup.f90:552
uzfcellgroupmodule::setbelowpet
subroutine setbelowpet(this, icell, jbelow)
Subtract aet from pet to calculate residual et for deeper cells.
Definition: UzfCellGroup.f90:492
uzfcellgroupmodule::gwseep
subroutine gwseep(this, icell, deriv, scale, hgwf, trhs, thcof, seep)
Calculate groundwater discharge to land surface.
Definition: UzfCellGroup.f90:768
uzfcellgroupmodule::sethead
subroutine sethead(this, icell, hgwf)
Set initial head for uzf object.
Definition: UzfCellGroup.f90:381
uzfcellgroupmodule::leadwav
subroutine leadwav(this, time, itester, itrailflg, thetab, fluxb, ffcheck, feps2, delt, icell)
Create a lead wave and route over time step.
Definition: UzfCellGroup.f90:1181
uzfcellgroupmodule::caph
real(dp) function caph(this, icell, tho)
Calculate capillary pressure head from B-C equation.
Definition: UzfCellGroup.f90:1870
uzfcellgroupmodule::setdatauzfarea
subroutine setdatauzfarea(this, icell, areamult)
Set uzfarea using cellarea and areamult.
Definition: UzfCellGroup.f90:417
uzfcellgroupmodule::setdataetha
subroutine setdataetha(this, icell, jbelow, ha, hroot, rootact)
Set variables for head-based unsaturated flow.
Definition: UzfCellGroup.f90:530
uzfcellgroupmodule::addrech
subroutine addrech(this, icell, jbelow, hgwf, trhs, thcof, deriv, delt)
Add recharge or infiltration to cells.
Definition: UzfCellGroup.f90:701
uzfcellgroupmodule::setwaves
subroutine setwaves(this, icell)
Reset waves to default values at start of simulation.
Definition: UzfCellGroup.f90:865
uzfcellgroupmodule::solve
subroutine solve(this, thiswork, jbelow, icell, totfluxtot, ietflag, issflag, iseepflag, hgwf, qfrommvr, ierr, reset_state, trhs, thcof, deriv, watercontent)
Formulate the unsaturated flow object, calculate terms for gwf equation.
Definition: UzfCellGroup.f90:566
uzfcellgroupmodule::setdata
subroutine setdata(this, icell, area, top, bot, surfdep, vks, thtr, thts, thti, eps, ntrail, landflag, ivertcon)
Set uzf object material properties.
Definition: UzfCellGroup.f90:337
uzfcellgroupmodule::leadspeed
real(dp) function leadspeed(theta1, theta2, flux1, flux2, thts, thtr, eps, vks)
Calculates waves speed from dflux/dtheta.
Definition: UzfCellGroup.f90:1381
uzfcellgroupmodule::get_water_content_at_depth
real(dp) function get_water_content_at_depth(this, icell, depth)
Determine the water content at a specific depth.
Definition: UzfCellGroup.f90:1909
uzfcellgroupmodule::wave_shift
subroutine wave_shift(this, this2, icell, icell2, shft, strt, stp, cntr)
Copy waves or shift waves in arrays.
Definition: UzfCellGroup.f90:952
uzfcellgroupmodule::rejfinf
subroutine rejfinf(this, icell, deriv, hgwf, trhs, thcof, finfact)
Reject applied infiltration due to low vks.
Definition: UzfCellGroup.f90:740
uzfcellgroupmodule::rate_et_z
real(dp) function rate_et_z(this, icell, factor, fktho, h)
Calculate capillary pressure-based uz et.
Definition: UzfCellGroup.f90:1892
uzfcellgroupmodule::simgwet
subroutine simgwet(this, igwetflag, icell, hgwf, trhs, thcof, det)
Calculate gwf et using residual uzf pet.
Definition: UzfCellGroup.f90:808
uzfcellgroupmodule::get_wcnew
real(dp) function get_wcnew(this, icell)
Calculate and return the cell-based water content value.
Definition: UzfCellGroup.f90:1939
uzfcellgroupmodule::uz_rise
subroutine uz_rise(this, icell, totfluxtot)
Calculate recharge due to a rise in the gwf head.
Definition: UzfCellGroup.f90:845
uzfcellgroupmodule::uzet
subroutine uzet(this, icell, delt, ietflag, ierr)
Remove water from uz due to et.
Definition: UzfCellGroup.f90:1536
uzfcellgroupmodule::setgwpet
subroutine setgwpet(this, icell)
Subtract aet from pet to calculate residual et for gw.
Definition: UzfCellGroup.f90:473
uzfcellgroupmodule::routewaves
subroutine routewaves(this, totfluxtot, delt, ietflag, icell, ierr)
Prepare and route waves over time step.
Definition: UzfCellGroup.f90:914
uzfcellgroupmodule::setdataet
subroutine setdataet(this, icell, jbelow, pet, extdp)
Set unsaturated ET-related variables.
Definition: UzfCellGroup.f90:429
uzfcellgroupmodule::dealloc
subroutine dealloc(this)
Deallocate uzf object variables.
Definition: UzfCellGroup.f90:242
uzfcellgroupmodule::unsat_stor
real(dp) function unsat_stor(this, icell, d1)
Sums up mobile water over depth interval.
Definition: UzfCellGroup.f90:1416
uzfcellgroupmodule::trailwav
subroutine trailwav(this, icell, ierr)
Create and set trail waves.
Definition: UzfCellGroup.f90:1094
uzfcellgroupmodule::setdataetwc
subroutine setdataetwc(this, icell, jbelow, extwc)
Set extinction water content.
Definition: UzfCellGroup.f90:516
uzfcellgroupmodule::setdatafinf
subroutine setdatafinf(this, icell, finf)
Set infiltration.
Definition: UzfCellGroup.f90:397
uzfetutilmodule::etfunc_lin
real(dp) function, public etfunc_lin(efflndsrf, extdp, resid_pet, deriv_et, trhs, thcof, hgwf, celtop, celbot)
Calculate gwf ET using linear decay ET function from mf-2005.
Definition: UzfEtUtil.f90:18
uzfetutilmodule::etfunc_nlin
real(dp) function, public etfunc_nlin(efflndsrf, extdp, resid_pet, deriv_et, trhs, thcof, hgwf)
Calculate gwf ET using a square decay ET function with smoothing at the specified extinction depth.
Definition: UzfEtUtil.f90:79