MODFLOW 6  version 6.8.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) then
659  if (reset_state) &
660  call this%wave_shift(thiswork, icell, 1, 0, 1, thiswork%nwavst(1), 1)
661  return
662  end if
663  call this%uz_rise(icell, totfluxtot)
664  this%totflux(icell) = totfluxtot
665  if (this%ivertcon(icell) > 0) then
666  call this%addrech(icell, jbelow, hgwf, trhsfinf, thcoffinf, &
667  derivfinf, delt)
668  end if
669  else
670  this%totflux(icell) = this%surflux(icell) * delt
671  totfluxtot = this%surflux(icell) * delt
672  end if
673  thcoffinf = dzero
674  trhsfinf = this%totflux(icell) * this%uzfarea(icell) / delt
675  if (.not. reset_state) then
676  call this%update_wav(icell, delt, issflag, 0)
677  end if
678  else
679  this%totflux(icell) = this%surflux(icell) * delt
680  totfluxtot = this%surflux(icell) * delt
681  if (.not. reset_state) then
682  call this%update_wav(icell, delt, issflag, 1)
683  end if
684  end if
685  !
686  ! -- If formulating, then these variables will be present
687  if (present(deriv)) deriv = deriv1 + deriv2 + derivfinf
688  if (present(trhs)) trhs = trhsfinf + trhsseep
689  if (present(thcof)) thcof = thcoffinf + thcofseep
690  !
691  ! -- Assign water content prior to resetting waves
692  if (present(watercontent)) then
693  watercontent = this%get_wcnew(icell)
694  end if
695  !
696  ! -- reset waves to previous state for next iteration
697  if (reset_state) then
698  call this%wave_shift(thiswork, icell, 1, 0, 1, thiswork%nwavst(1), 1)
699  end if
700  end subroutine solve
701 
702  !> @brief Add recharge or infiltration to cells
703  !<
704  subroutine addrech(this, icell, jbelow, hgwf, trhs, thcof, deriv, delt)
705  ! -- dummy
706  class(uzfcellgrouptype) :: this
707  integer(I4B), intent(in) :: icell
708  integer(I4B), intent(in) :: jbelow
709  real(DP), intent(inout) :: trhs
710  real(DP), intent(inout) :: thcof
711  real(DP), intent(inout) :: deriv
712  real(DP), intent(in) :: delt
713  real(DP), intent(in) :: hgwf
714  ! -- local
715  real(DP) :: fcheck
716  real(DP) :: x, scale, range
717  !
718  ! -- initialize
719  range = dem5
720  deriv = dzero
721  thcof = dzero
722  trhs = this%uzfarea(icell) * this%totflux(icell) / delt
723  if (this%totflux(icell) < dem14) return
724  scale = done
725  !
726  ! -- smoothly reduce flow between cells when head close to cell top
727  x = hgwf - (this%celbot(icell) - range)
728  call sscurve(x, range, deriv, scale)
729  deriv = this%uzfarea(icell) * deriv * this%totflux(icell) / delt
730  this%finf(jbelow) = (done - scale) * this%totflux(icell) / delt
731  fcheck = this%finf(jbelow) - this%vks(jbelow)
732  !
733  ! -- reduce flow between cells when vks is too small
734  if (fcheck < dem14) fcheck = dzero
735  this%finf(jbelow) = this%finf(jbelow) - fcheck
736  this%surfluxbelow(icell) = this%finf(jbelow)
737  this%totflux(icell) = scale * this%totflux(icell) + fcheck * delt
738  trhs = this%uzfarea(icell) * this%totflux(icell) / delt
739  end subroutine addrech
740 
741  !> @brief Reject applied infiltration due to low vks
742  !<
743  subroutine rejfinf(this, icell, deriv, hgwf, trhs, thcof, finfact)
744  ! -- dummy
745  class(uzfcellgrouptype) :: this
746  integer(I4B), intent(in) :: icell
747  real(DP), intent(inout) :: deriv
748  real(DP), intent(inout) :: finfact
749  real(DP), intent(inout) :: thcof
750  real(DP), intent(inout) :: trhs
751  real(DP), intent(in) :: hgwf
752  ! -- local
753  real(DP) :: x, range, scale, q
754  !
755  range = this%surfdep(icell)
756  q = this%surflux(icell)
757  finfact = q
758  trhs = finfact * this%uzfarea(icell)
759  x = this%celtop(icell) - hgwf
760  call slinear(x, range, deriv, scale)
761  deriv = -q * deriv * this%uzfarea(icell) * scale
762  if (scale < done) then
763  finfact = q * scale
764  trhs = finfact * this%uzfarea(icell) * this%celtop(icell) / range
765  thcof = finfact * this%uzfarea(icell) / range
766  end if
767  end subroutine rejfinf
768 
769  !> @brief Calculate groundwater discharge to land surface
770  !<
771  subroutine gwseep(this, icell, deriv, scale, hgwf, trhs, thcof, seep)
772  ! -- dummy
773  class(uzfcellgrouptype) :: this
774  integer(I4B), intent(in) :: icell
775  real(DP), intent(inout) :: deriv
776  real(DP), intent(inout) :: trhs
777  real(DP), intent(inout) :: thcof
778  real(DP), intent(inout) :: seep
779  real(DP), intent(out) :: scale
780  real(DP), intent(in) :: hgwf
781  ! -- local
782  real(DP) :: x, range, y, deriv1, d1, d2, Q
783  !
784  seep = dzero
785  deriv = dzero
786  deriv1 = dzero
787  d1 = dzero
788  d2 = dzero
789  scale = dzero
790  q = this%uzfarea(icell) * this%vks(icell)
791  range = this%surfdep(icell)
792  x = hgwf - this%celtop(icell)
793  call scubiclinear(x, range, deriv1, y)
794  scale = y
795  seep = scale * q * (hgwf - this%celtop(icell)) / range
796  trhs = scale * q * this%celtop(icell) / range
797  thcof = -scale * q / range
798  d1 = -deriv1 * q * x / range
799  d2 = -scale * q / range
800  deriv = d1 + d2
801  if (seep < dzero) then
802  seep = dzero
803  deriv = dzero
804  trhs = dzero
805  thcof = dzero
806  end if
807  end subroutine gwseep
808 
809  !> @brief Calculate gwf et using residual uzf pet
810  !<
811  subroutine simgwet(this, igwetflag, icell, hgwf, trhs, thcof, det)
812  ! -- dummy
813  class(uzfcellgrouptype) :: this
814  integer(I4B), intent(in) :: igwetflag
815  integer(I4B), intent(in) :: icell
816  real(DP), intent(in) :: hgwf
817  real(DP), intent(inout) :: trhs
818  real(DP), intent(inout) :: thcof
819  real(DP), intent(inout) :: det
820  ! -- local
821  real(DP) :: s, x, c, b, et
822  !
823  this%gwet(icell) = dzero
824  trhs = dzero
825  thcof = dzero
826  det = dzero
827  s = this%landtop(icell)
828  x = this%extdp(icell)
829  c = this%gwpet(icell)
830  b = this%celbot(icell)
831  if (b > hgwf) return
832  if (x < dem6) return
833  if (igwetflag == 1) then
834  et = etfunc_lin(s, x, c, det, trhs, thcof, hgwf, &
835  this%celtop(icell), this%celbot(icell))
836  else if (igwetflag == 2) then
837  et = etfunc_nlin(s, x, c, det, trhs, thcof, hgwf)
838  end if
839  ! this%gwet(icell) = et * this%uzfarea(icell)
840  trhs = trhs * this%uzfarea(icell)
841  thcof = thcof * this%uzfarea(icell)
842  this%gwet(icell) = trhs - (thcof * hgwf)
843  ! write(99,*)'in group', icell, this%gwet(icell)
844  end subroutine simgwet
845 
846  !> @brief Calculate recharge due to a rise in the gwf head
847  !<
848  subroutine uz_rise(this, icell, totfluxtot)
849  ! -- dummy
850  class(uzfcellgrouptype) :: this
851  integer(I4B), intent(in) :: icell
852  real(DP), intent(inout) :: totfluxtot
853  ! -- local
854  real(DP) :: fm1, fm2, d1
855  !
856  ! -- additional recharge from a rising water table
857  if (this%watab(icell) - this%watabold(icell) > dem30) then
858  d1 = this%celtop(icell) - this%watabold(icell)
859  fm1 = this%unsat_stor(icell, d1)
860  d1 = this%celtop(icell) - this%watab(icell)
861  fm2 = this%unsat_stor(icell, d1)
862  totfluxtot = totfluxtot + (fm1 - fm2)
863  end if
864  end subroutine uz_rise
865 
866  !> @brief Reset waves to default values at start of simulation
867  !<
868  subroutine setwaves(this, icell)
869  ! -- dummy
870  class(uzfcellgrouptype) :: this
871  ! -- local
872  integer(I4B), intent(in) :: icell
873  real(DP) :: bottom, top
874  integer(I4B) :: jk
875  real(DP) :: thick
876  !
877  ! -- initialize
878  this%totflux(icell) = dzero
879  this%nwavst(icell) = 1
880  this%uzdpst(:, icell) = dzero
881  thick = this%celtop(icell) - this%watab(icell)
882  do jk = 1, this%nwav(icell)
883  this%uzthst(jk, icell) = this%thtr(icell)
884  end do
885  !
886  ! -- initialize waves for first stress period
887  if (thick > dzero) then
888  this%uzdpst(1, icell) = thick
889  this%uzthst(1, icell) = this%thti(icell)
890  top = this%uzthst(1, icell) - this%thtr(icell)
891  if (top < dzero) top = dzero
892  bottom = this%thts(icell) - this%thtr(icell)
893  if (bottom < dzero) bottom = dzero
894  this%uzflst(1, icell) = this%vks(icell) * (top / bottom)**this%eps(icell)
895  if (this%uzthst(1, icell) < this%thtr(icell)) &
896  this%uzthst(1, icell) = this%thtr(icell)
897  !
898  ! -- calculate water stored in the unsaturated zone
899  if (top > dzero) then
900  this%uzspst(1, icell) = dzero
901  else
902  this%uzflst(1, icell) = dzero
903  this%uzspst(1, icell) = dzero
904  end if
905  !
906  ! no unsaturated zone
907  else
908  this%uzflst(1, icell) = dzero
909  this%uzdpst(1, icell) = dzero
910  this%uzspst(1, icell) = dzero
911  this%uzthst(1, icell) = this%thtr(icell)
912  end if
913  end subroutine
914 
915  !> @brief Prepare and route waves over time step
916  !<
917  subroutine routewaves(this, totfluxtot, delt, ietflag, icell, ierr)
918  ! -- dummy
919  class(uzfcellgrouptype) :: this
920  real(DP), intent(inout) :: totfluxtot
921  real(DP), intent(in) :: delt
922  integer(I4B), intent(in) :: ietflag
923  integer(I4B), intent(in) :: icell
924  integer(I4B), intent(inout) :: ierr
925  ! -- local
926  real(DP) :: thick, thickold
927  integer(I4B) :: idelt, iwav, ik
928  !
929  ! -- initialize
930  this%totflux(icell) = dzero
931  this%etact(icell) = dzero
932  thick = this%celtop(icell) - this%watab(icell)
933  thickold = this%celtop(icell) - this%watabold(icell)
934  !
935  ! -- no uz, clear waves
936  if (thickold < dzero) then
937  do iwav = 1, this%nwavst(icell)
938  this%uzthst(iwav, icell) = this%thtr(icell)
939  this%uzdpst(iwav, icell) = dzero
940  this%uzspst(iwav, icell) = dzero
941  this%uzflst(iwav, icell) = dzero
942  end do
943  this%nwavst(icell) = 1
944  end if
945  idelt = 1
946  do ik = 1, idelt
947  call this%uzflow(thick, thickold, delt, ietflag, icell, ierr)
948  if (ierr > 0) return
949  totfluxtot = totfluxtot + this%totflux(icell)
950  end do
951  end subroutine routewaves
952 
953  !> @brief Copy waves or shift waves in arrays
954  !<
955  subroutine wave_shift(this, this2, icell, icell2, shft, strt, stp, cntr)
956  ! -- dummy
957  class(uzfcellgrouptype) :: this
958  type(uzfcellgrouptype) :: this2
959  integer(I4B), intent(in) :: icell
960  integer(I4B), intent(in) :: icell2
961  integer(I4B), intent(in) :: shft
962  integer(I4B), intent(in) :: strt
963  integer(I4B), intent(in) :: stp
964  integer(I4B), intent(in) :: cntr
965  ! -- local
966  integer(I4B) :: j
967  !
968  ! -- copy waves from one uzf cell group to another
969  do j = strt, stp, cntr
970  this%uzthst(j, icell) = this2%uzthst(j + shft, icell2)
971  this%uzdpst(j, icell) = this2%uzdpst(j + shft, icell2)
972  this%uzflst(j, icell) = this2%uzflst(j + shft, icell2)
973  this%uzspst(j, icell) = this2%uzspst(j + shft, icell2)
974  end do
975  this%nwavst(icell) = this2%nwavst(icell2)
976  end subroutine
977 
978  !> @brief Method of Characteristics solution for kinematic wave equation
979  !<
980  subroutine uzflow(this, thick, thickold, delt, ietflag, icell, ierr)
981  ! -- dummy
982  class(uzfcellgrouptype) :: this
983  real(DP), intent(inout) :: thickold
984  real(DP), intent(inout) :: thick
985  real(DP), intent(in) :: delt
986  integer(I4B), intent(in) :: ietflag
987  integer(I4B), intent(in) :: icell
988  integer(I4B), intent(inout) :: ierr
989  ! -- local
990  real(DP) :: ffcheck, time, feps1, feps2
991  real(DP) :: thetadif, thetab, fluxb, oldsflx
992  integer(I4B) :: itrailflg, itester
993  !
994  time = dzero
995  this%totflux(icell) = dzero
996  itrailflg = 0
997  oldsflx = this%uzflst(this%nwavst(icell), icell)
998  call factors(feps1, feps2)
999  !
1000  ! -- check for falling or rising water table
1001  if ((thick - thickold) > feps1) then
1002  thetadif = abs(this%uzthst(1, icell) - this%thtr(icell))
1003  if (thetadif > dem6) then
1004  call this%wave_shift(this, icell, icell, -1, &
1005  this%nwavst(icell) + 1, 2, -1)
1006  if (this%uzdpst(2, icell) < dem30) &
1007  this%uzdpst(2, icell) = (this%ntrail(icell) + dtwo) * dem6
1008  if (this%uzthst(2, icell) > this%thtr(icell)) then
1009  this%uzspst(2, icell) = this%uzflst(2, icell) / &
1010  (this%uzthst(2, icell) - this%thtr(icell))
1011  else
1012  this%uzspst(2, icell) = dzero
1013  end if
1014  this%uzthst(1, icell) = this%thtr(icell)
1015  this%uzflst(1, icell) = dzero
1016  this%uzspst(1, icell) = dzero
1017  this%uzdpst(1, icell) = thick
1018  this%nwavst(icell) = this%nwavst(icell) + 1
1019  if (this%nwavst(icell) >= this%nwav(icell)) then
1020  ! -- too many waves error
1021  ierr = 1
1022  return
1023  end if
1024  else
1025  this%uzdpst(1, icell) = thick
1026  end if
1027  end if
1028  thetab = this%uzthst(1, icell)
1029  fluxb = this%uzflst(1, icell)
1030  this%totflux(icell) = dzero
1031  itester = 0
1032  ffcheck = (this%surflux(icell) - this%uzflst(this%nwavst(icell), icell))
1033  !
1034  ! -- increase new waves in infiltration changes
1035  if (ffcheck > feps2 .OR. ffcheck < -feps2) then
1036  this%nwavst(icell) = this%nwavst(icell) + 1
1037  if (this%nwavst(icell) >= this%nwav(icell)) then
1038  !
1039  ! -- too many waves error
1040  ierr = 1
1041  return
1042  end if
1043  else if (this%nwavst(icell) == 1) then
1044  itester = 1
1045  end if
1046  if (this%nwavst(icell) > 1) then
1047  if (ffcheck < -feps2) then
1048  call this%trailwav(icell, ierr)
1049  if (ierr > 0) return
1050  itrailflg = 1
1051  end if
1052  call this%leadwav(time, itester, itrailflg, thetab, fluxb, ffcheck, &
1053  feps2, delt, icell)
1054  end if
1055  if (itester == 1) then
1056  this%totflux(icell) = this%totflux(icell) + &
1057  (delt - time) * this%uzflst(1, icell)
1058  time = dzero
1059  itester = 0
1060  end if
1061  !
1062  ! -- simulate et
1063  if (ietflag > 0) call this%uzet(icell, delt, ietflag, ierr)
1064  if (ierr > 0) return
1065  end subroutine uzflow
1066 
1067  !> @brief Calculate unit specific tolerances
1068  !<
1069  subroutine factors(feps1, feps2)
1070  ! -- dummy
1071  real(DP), intent(out) :: feps1
1072  real(DP), intent(out) :: feps2
1073  real(DP) :: factor1
1074  real(DP) :: factor2
1075  !
1076  ! calculate constants for uzflow
1077  factor1 = done
1078  factor2 = done
1079  feps1 = dem9
1080  feps2 = dem9
1081  if (itmuni == 1) then
1082  factor1 = done / 86400.d0
1083  else if (itmuni == 2) then
1084  factor1 = done / 1440.d0
1085  else if (itmuni == 3) then
1086  factor1 = done / 24.0d0
1087  else if (itmuni == 5) then
1088  factor1 = 365.0d0
1089  end if
1090  factor2 = done / 0.3048
1091  feps1 = feps1 * factor1 * factor2
1092  feps2 = feps2 * factor1 * factor2
1093  end subroutine factors
1094 
1095  !> @brief Create and set trail waves
1096  !<
1097  subroutine trailwav(this, icell, ierr)
1098  ! -- dummy
1099  class(uzfcellgrouptype) :: this
1100  integer(I4B), intent(in) :: icell
1101  integer(I4B), intent(inout) :: ierr
1102  ! -- local
1103  real(DP) :: smoist, smoistinc, ftrail, eps_m1
1104  real(DP) :: thtsrinv
1105  real(DP) :: flux1, flux2, theta1, theta2
1106  real(DP) :: fnuminc
1107  integer(I4B) :: j, jj, jk, nwavstm1
1108  !
1109  ! -- initialize
1110  eps_m1 = dble(this%eps(icell)) - done
1111  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1112  nwavstm1 = this%nwavst(icell) - 1
1113  !
1114  ! -- initialize trailwaves
1115  smoist = (((this%surflux(icell) / this%vks(icell))** &
1116  (done / this%eps(icell))) * &
1117  (this%thts(icell) - this%thtr(icell))) + this%thtr(icell)
1118  if (this%uzthst(nwavstm1, icell) - smoist > dem9) then
1119  fnuminc = dzero
1120  do jk = 1, this%ntrail(icell)
1121  fnuminc = fnuminc + float(jk)
1122  end do
1123  smoistinc = (this%uzthst(nwavstm1, icell) - smoist) / (fnuminc - done)
1124  jj = this%ntrail(icell)
1125  ftrail = dble(this%ntrail(icell)) + done
1126  do j = this%nwavst(icell), this%nwavst(icell) + this%ntrail(icell) - 1
1127  if (j > this%nwav(icell)) then
1128  ! -- too many waves error
1129  ierr = 1
1130  return
1131  end if
1132  if (j > this%nwavst(icell)) then
1133  this%uzthst(j, icell) = this%uzthst(j - 1, icell) &
1134  - ((ftrail - float(jj)) * smoistinc)
1135  else
1136  this%uzthst(j, icell) = this%uzthst(j - 1, icell) - dem9
1137  end if
1138  jj = jj - 1
1139  if (this%uzthst(j, icell) <= this%thtr(icell) + dem9) &
1140  this%uzthst(j, icell) = this%thtr(icell) + dem9
1141  this%uzflst(j, icell) = &
1142  this%vks(icell) * (((this%uzthst(j, icell) - this%thtr(icell)) * &
1143  thtsrinv)**this%eps(icell))
1144  theta2 = this%uzthst(j - 1, icell)
1145  flux2 = this%uzflst(j - 1, icell)
1146  flux1 = this%uzflst(j, icell)
1147  theta1 = this%uzthst(j, icell)
1148  this%uzspst(j, icell) = leadspeed(theta1, theta2, flux1, flux2, &
1149  this%thts(icell), this%thtr(icell), &
1150  this%eps(icell), this%vks(icell))
1151  this%uzdpst(j, icell) = dzero
1152  if (j == this%nwavst(icell)) then
1153  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1154  (this%ntrail(icell) + 1) * dem9
1155  else
1156  this%uzdpst(j, icell) = this%uzdpst(j - 1, icell) - dem9
1157  end if
1158  end do
1159  this%nwavst(icell) = this%nwavst(icell) + this%ntrail(icell) - 1
1160  if (this%nwavst(icell) >= this%nwav(icell)) then
1161  ! -- too many waves error
1162  ierr = 1
1163  return
1164  end if
1165  else
1166  this%uzdpst(this%nwavst(icell), icell) = dzero
1167  this%uzflst(this%nwavst(icell), icell) = &
1168  this%vks(icell) * (((this%uzthst(this%nwavst(icell), icell) - &
1169  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1170  this%uzthst(this%nwavst(icell), icell) = smoist
1171  theta2 = this%uzthst(this%nwavst(icell) - 1, icell)
1172  flux2 = this%uzflst(this%nwavst(icell) - 1, icell)
1173  flux1 = this%uzflst(this%nwavst(icell), icell)
1174  theta1 = this%uzthst(this%nwavst(icell), icell)
1175  this%uzspst(this%nwavst(icell), icell) = &
1176  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1177  this%thtr(icell), this%eps(icell), this%vks(icell))
1178  end if
1179  end subroutine trailwav
1180 
1181  !> @brief Create a lead wave and route over time step
1182  !<
1183  subroutine leadwav(this, time, itester, itrailflg, thetab, fluxb, &
1184  ffcheck, feps2, delt, icell)
1185  ! -- dummy
1186  class(uzfcellgrouptype) :: this
1187  real(DP), intent(inout) :: thetab
1188  real(DP), intent(inout) :: fluxb
1189  real(DP), intent(in) :: feps2
1190  real(DP), intent(inout) :: time
1191  integer(I4B), intent(inout) :: itester
1192  integer(I4B), intent(inout) :: itrailflg
1193  real(DP), intent(inout) :: ffcheck
1194  real(DP), intent(in) :: delt
1195  integer(I4B), intent(in) :: icell
1196  ! -- local
1197  real(DP) :: bottomtime, shortest, fcheck
1198  real(DP) :: eps_m1, timenew, bottom, timedt
1199  real(DP) :: thtsrinv, diff, fluxhld2
1200  real(DP) :: flux1, flux2, theta1, theta2, ftest
1201  real(DP), allocatable, dimension(:) :: checktime
1202  integer(I4B) :: iflx, iremove, j, l
1203  integer(I4B) :: nwavp1, jshort
1204  integer(I4B), allocatable, dimension(:) :: more
1205  !
1206  allocate (checktime(this%nwavst(icell)))
1207  allocate (more(this%nwavst(icell)))
1208  ftest = dzero
1209  eps_m1 = dble(this%eps(icell)) - done
1210  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1211  !
1212  ! -- initialize new wave
1213  if (itrailflg == 0) then
1214  if (ffcheck > feps2) then
1215  this%uzflst(this%nwavst(icell), icell) = this%surflux(icell)
1216  if (this%uzflst(this%nwavst(icell), icell) < dem30) &
1217  this%uzflst(this%nwavst(icell), icell) = dzero
1218  this%uzthst(this%nwavst(icell), icell) = &
1219  (((this%uzflst(this%nwavst(icell), icell) / this%vks(icell))** &
1220  (done / this%eps(icell))) * (this%thts(icell) - this%thtr(icell))) &
1221  + this%thtr(icell)
1222  theta2 = this%uzthst(this%nwavst(icell), icell)
1223  flux2 = this%uzflst(this%nwavst(icell), icell)
1224  flux1 = this%uzflst(this%nwavst(icell) - 1, icell)
1225  theta1 = this%uzthst(this%nwavst(icell) - 1, icell)
1226  this%uzspst(this%nwavst(icell), icell) = &
1227  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1228  this%thtr(icell), this%eps(icell), this%vks(icell))
1229  this%uzdpst(this%nwavst(icell), icell) = dzero
1230  end if
1231  end if
1232  !
1233  ! -- route all waves and interception of waves over times step
1234  diff = done
1235  timedt = dzero
1236  iflx = 0
1237  fluxhld2 = this%uzflst(1, icell)
1238  if (this%nwavst(icell) == 0) itester = 1
1239  if (itester /= 1) then
1240  do while (diff > dem6)
1241  nwavp1 = this%nwavst(icell) + 1
1242  timedt = delt - time
1243  do j = 1, this%nwavst(icell)
1244  checktime(j) = dep20
1245  more(j) = 0
1246  end do
1247  shortest = timedt
1248  if (this%nwavst(icell) > 2) then
1249  j = 2
1250  !
1251  ! -- calculate time until wave overtakes wave ahead
1252  nwavp1 = this%nwavst(icell) + 1
1253  do while (j < nwavp1)
1254  ftest = this%uzspst(j - 1, icell) - this%uzspst(j, icell)
1255  if (abs(ftest) > dem30) then
1256  checktime(j) = (this%uzdpst(j, icell) - &
1257  this%uzdpst(j - 1, icell)) / ftest
1258  if (checktime(j) < dem30) checktime(j) = dep20
1259  end if
1260  j = j + 1
1261  end do
1262  end if
1263  !
1264  ! - calc time until wave reaches bottom of cell
1265  bottomtime = dep20
1266  if (this%nwavst(icell) > 1) then
1267  if (this%uzspst(2, icell) > dzero) then
1268  bottom = this%uzspst(2, icell)
1269  if (bottom < dem15) bottom = dem15
1270  bottomtime = (this%uzdpst(1, icell) - this%uzdpst(2, icell)) / bottom
1271  if (bottomtime < dzero) bottomtime = dem12
1272  end if
1273  end if
1274  !
1275  ! -- calc time for wave interception
1276  jshort = 0
1277  do j = this%nwavst(icell), 3, -1
1278  if (shortest - checktime(j) > -dem9) then
1279  more(j) = 1
1280  jshort = j
1281  shortest = checktime(j)
1282  end if
1283  end do
1284  do j = 3, this%nwavst(icell)
1285  if (shortest - checktime(j) < dem9) then
1286  if (j /= jshort) more(j) = 0
1287  end if
1288  end do
1289  !
1290  ! -- what happens first, waves hits bottom or interception
1291  iremove = 0
1292  timenew = time
1293  fcheck = (time + shortest) - delt
1294  if (shortest < dem7) fcheck = -done
1295  if (bottomtime < shortest .AND. time + bottomtime < delt) then
1296  j = 2
1297  do while (j < nwavp1)
1298  !
1299  ! -- route waves
1300  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1301  this%uzspst(j, icell) * bottomtime
1302  j = j + 1
1303  end do
1304  fluxb = this%uzflst(2, icell)
1305  thetab = this%uzthst(2, icell)
1306  iflx = 1
1307  call this%wave_shift(this, icell, icell, 1, 1, &
1308  this%nwavst(icell) - 1, 1)
1309  iremove = 1
1310  timenew = time + bottomtime
1311  this%uzspst(1, icell) = dzero
1312  !
1313  ! -- do waves intercept before end of time step
1314  else if (fcheck < dzero .AND. this%nwavst(icell) > 2) then
1315  j = 2
1316  do while (j < nwavp1)
1317  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1318  this%uzspst(j, icell) * shortest
1319  j = j + 1
1320  end do
1321  !
1322  ! -- combine waves that intercept, remove a wave
1323  j = 3
1324  l = j
1325  do while (j < this%nwavst(icell) + 1)
1326  if (more(j) == 1) then
1327  l = j
1328  theta2 = this%uzthst(j, icell)
1329  flux2 = this%uzflst(j, icell)
1330  if (j == 3) then
1331  flux1 = fluxb
1332  theta1 = thetab
1333  else
1334  flux1 = this%uzflst(j - 2, icell)
1335  theta1 = this%uzthst(j - 2, icell)
1336  end if
1337  this%uzspst(j, icell) = leadspeed(theta1, theta2, flux1, flux2, &
1338  this%thts(icell), &
1339  this%thtr(icell), &
1340  this%eps(icell), this%vks(icell))
1341  !
1342  ! -- update waves
1343  call this%wave_shift(this, icell, icell, 1, l - 1, &
1344  this%nwavst(icell) - 1, 1)
1345  l = this%nwavst(icell) + 1
1346  iremove = iremove + 1
1347  end if
1348  j = j + 1
1349  end do
1350  timenew = timenew + shortest
1351  !
1352  ! -- calc. total flux to bottom during remaining time in step
1353  else
1354  j = 2
1355  do while (j < nwavp1)
1356  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1357  this%uzspst(j, icell) * timedt
1358  j = j + 1
1359  end do
1360  timenew = delt
1361  end if
1362  this%totflux(icell) = this%totflux(icell) + fluxhld2 * (timenew - time)
1363  if (iflx == 1) then
1364  fluxhld2 = this%uzflst(1, icell)
1365  iflx = 0
1366  end if
1367  !
1368  ! -- remove dead waves
1369  this%nwavst(icell) = this%nwavst(icell) - iremove
1370  time = timenew
1371  diff = delt - time
1372  if (this%nwavst(icell) == 1) then
1373  itester = 1
1374  exit
1375  end if
1376  end do
1377  end if
1378  deallocate (checktime)
1379  deallocate (more)
1380  end subroutine leadwav
1381 
1382  !> @brief Calculates waves speed from dflux/dtheta
1383  !<
1384  function leadspeed(theta1, theta2, flux1, flux2, thts, thtr, eps, vks)
1385  ! -- Return
1386  real(dp) :: leadspeed
1387  ! -- dummy
1388  real(dp), intent(in) :: theta1
1389  real(dp), intent(in) :: theta2
1390  real(dp), intent(in) :: flux1
1391  real(dp), intent(inout) :: flux2
1392  real(dp), intent(in) :: thts
1393  real(dp), intent(in) :: thtr
1394  real(dp), intent(in) :: eps
1395  real(dp), intent(in) :: vks
1396  ! -- local
1397  real(dp) :: comp1, comp2, thsrinv, epsfksths
1398  real(dp) :: eps_m1, fhold, comp3
1399  !
1400  eps_m1 = eps - done
1401  thsrinv = done / (thts - thtr)
1402  epsfksths = eps * vks * thsrinv
1403  comp1 = theta2 - theta1
1404  comp2 = abs(flux2 - flux1)
1405  comp3 = theta1 - thtr
1406  if (comp2 < dem15) flux2 = flux1 + dem15
1407  if (abs(comp1) < dem30) then
1408  fhold = dem30
1409  if (comp3 > dem30) fhold = (comp3 * thsrinv)**eps
1410  if (fhold < dem30) fhold = dem30
1411  leadspeed = epsfksths * (fhold**eps_m1)
1412  else
1413  leadspeed = (flux2 - flux1) / (theta2 - theta1)
1414  end if
1415  if (leadspeed < dem30) leadspeed = dem30
1416  end function leadspeed
1417 
1418  !> @brief Sums up mobile water over depth interval
1419  !<
1420  function unsat_stor(this, icell, d1)
1421  ! -- Return
1422  real(dp) :: unsat_stor
1423  ! -- dummy
1424  class(uzfcellgrouptype) :: this
1425  integer(I4B), intent(in) :: icell
1426  real(dp), intent(inout) :: d1
1427  ! -- local
1428  real(dp) :: fm
1429  integer(I4B) :: j, k, nwavm1, jj
1430  !
1431  fm = dzero
1432  j = this%nwavst(icell) + 1
1433  k = this%nwavst(icell)
1434  nwavm1 = k - 1
1435  if (d1 > this%uzdpst(1, icell)) d1 = this%uzdpst(1, icell)
1436  !
1437  ! -- find deepest wave above depth d1, counter held as j
1438  do while (k > 0)
1439  if (this%uzdpst(k, icell) - d1 < -dem30) j = k
1440  k = k - 1
1441  end do
1442  if (j > this%nwavst(icell)) then
1443  fm = fm + (this%uzthst(this%nwavst(icell), icell) - this%thtr(icell)) * d1
1444  elseif (this%nwavst(icell) > 1) then
1445  if (j > 1) then
1446  fm = fm + (this%uzthst(j - 1, icell) - this%thtr(icell)) &
1447  * (d1 - this%uzdpst(j, icell))
1448  end if
1449  do jj = j, nwavm1
1450  fm = fm + (this%uzthst(jj, icell) - this%thtr(icell)) &
1451  * (this%uzdpst(jj, icell) &
1452  - this%uzdpst(jj + 1, icell))
1453  end do
1454  fm = fm + (this%uzthst(this%nwavst(icell), icell) - this%thtr(icell)) &
1455  * (this%uzdpst(this%nwavst(icell), icell))
1456  else
1457  fm = fm + (this%uzthst(1, icell) - this%thtr(icell)) * d1
1458  end if
1459  unsat_stor = fm
1460  end function unsat_stor
1461 
1462  !> @brief Update to new state of uz at end of time step
1463  !<
1464  subroutine update_wav(this, icell, delt, iss, itest)
1465  ! -- dummy
1466  class(uzfcellgrouptype) :: this
1467  integer(I4B), intent(in) :: icell
1468  integer(I4B), intent(in) :: itest
1469  integer(I4B), intent(in) :: iss
1470  real(DP), intent(in) :: delt
1471  ! -- local
1472  real(DP) :: bot, depthsave, top
1473  real(DP) :: thick, thtsrinv
1474  integer(I4B) :: nwavhld, k, j
1475  !
1476  bot = this%watab(icell)
1477  top = this%celtop(icell)
1478  thick = top - bot
1479  nwavhld = this%nwavst(icell)
1480  if (itest == 1) then
1481  this%uzflst(1, icell) = dzero
1482  this%uzthst(1, icell) = this%thtr(icell)
1483  return
1484  end if
1485  if (iss == 1) then
1486  if (this%thts(icell) - this%thtr(icell) < dem7) then
1487  thtsrinv = done / dem7
1488  else
1489  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1490  end if
1491  this%totflux(icell) = this%surflux(icell) * delt
1492  this%watabold(icell) = this%watab(icell)
1493  this%uzthst(1, icell) = this%thti(icell)
1494  this%uzflst(1, icell) = &
1495  this%vks(icell) * (((this%uzthst(1, icell) - this%thtr(icell)) &
1496  * thtsrinv)**this%eps(icell))
1497  this%uzdpst(1, icell) = thick
1498  this%uzspst(1, icell) = thick
1499  this%nwavst(icell) = 1
1500  else
1501  !
1502  ! -- water table rises through waves
1503  if (this%watab(icell) - this%watabold(icell) > dem30) then
1504  depthsave = this%uzdpst(1, icell)
1505  j = 0
1506  k = this%nwavst(icell)
1507  do while (k > 0)
1508  if (this%uzdpst(k, icell) - thick < -dem30) j = k
1509  k = k - 1
1510  end do
1511  this%uzdpst(1, icell) = thick
1512  if (j > 1) then
1513  this%uzspst(1, icell) = dzero
1514  this%nwavst(icell) = this%nwavst(icell) - j + 2
1515  this%uzthst(1, icell) = this%uzthst(j - 1, icell)
1516  this%uzflst(1, icell) = this%uzflst(j - 1, icell)
1517  if (j > 2) call this%wave_shift(this, icell, icell, j - 2, 2, &
1518  nwavhld - (j - 2), 1)
1519  elseif (j == 0) then
1520  this%uzspst(1, icell) = dzero
1521  this%uzthst(1, icell) = this%uzthst(this%nwavst(icell), icell)
1522  this%uzflst(1, icell) = this%uzflst(this%nwavst(icell), icell)
1523  this%nwavst(icell) = 1
1524  end if
1525  end if
1526  !
1527  ! -- calculate new unsat. storage
1528  if (thick <= dzero) then
1529  this%uzspst(1, icell) = dzero
1530  this%nwavst(icell) = 1
1531  this%uzthst(1, icell) = this%thtr(icell)
1532  this%uzflst(1, icell) = dzero
1533  end if
1534  this%watabold(icell) = this%watab(icell)
1535  end if
1536  end subroutine update_wav
1537 
1538  !> @brief Remove water from uz due to et
1539  !<
1540  subroutine uzet(this, icell, delt, ietflag, ierr)
1541  ! -- dummy
1542  class(uzfcellgrouptype) :: this
1543  integer(I4B), intent(in) :: icell
1544  real(DP), intent(in) :: delt
1545  integer(I4B), intent(in) :: ietflag
1546  integer(I4B), intent(inout) :: ierr
1547  ! -- local
1548  type(uzfcellgrouptype) :: uzfktemp
1549  real(DP) :: diff
1550  real(DP) :: thetaout
1551  real(DP) :: fm
1552  real(DP) :: st
1553  real(DP) :: thtsrinv
1554  real(DP) :: epsfksthts
1555  real(DP) :: fmp
1556  real(DP) :: fktho
1557  real(DP) :: theta1
1558  real(DP) :: theta2
1559  real(DP) :: flux1
1560  real(DP) :: flux2
1561  real(DP) :: hcap
1562  real(DP) :: factor
1563  real(DP) :: tho
1564  real(DP) :: depth
1565  real(DP) :: extwc1
1566  real(DP) :: petsub
1567  integer(I4B) :: i
1568  integer(I4B) :: j
1569  integer(I4B) :: jhold
1570  integer(I4B) :: jk
1571  integer(I4B) :: kj
1572  integer(I4B) :: kk
1573  integer(I4B) :: numadd
1574  integer(I4B) :: k
1575  integer(I4B) :: nwv
1576  integer(I4B) :: itest
1577  !
1578  ! -- initialize
1579  this%etact(icell) = dzero
1580  if (this%extdpuz(icell) < dem7) return
1581  petsub = this%rootact(icell) * this%pet(icell) * &
1582  this%extdpuz(icell) / this%extdp(icell)
1583  thetaout = delt * petsub / this%extdp(icell)
1584  if (ietflag == 1) thetaout = delt * this%pet(icell) / this%extdp(icell)
1585  if (thetaout < dem10) return
1586  depth = this%uzdpst(1, icell)
1587  st = this%unsat_stor(icell, depth)
1588  if (st < dem4) return
1589  !
1590  ! -- allocate temporary wave storage.
1591  nwv = this%nwavst(icell)
1592  itest = 0
1593  call uzfktemp%init(1, nwv)
1594  !
1595  ! store original wave characteristics
1596  call uzfktemp%wave_shift(this, 1, icell, 0, 1, nwv, 1)
1597  factor = done
1598  this%etact(icell) = dzero
1599  if (this%thts(icell) - this%thtr(icell) < dem7) then
1600  thtsrinv = 1.0 / dem7
1601  else
1602  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1603  end if
1604  epsfksthts = this%eps(icell) * this%vks(icell) * thtsrinv
1605  this%etact(icell) = dzero
1606  fmp = dzero
1607  extwc1 = this%extwc(icell) - this%thtr(icell)
1608  if (extwc1 < dem6) extwc1 = dem7
1609  numadd = 0
1610  fm = st
1611  k = 0
1612  !
1613  ! -- loop for reducing aet to pet when et is head dependent
1614  do while (itest == 0)
1615  k = k + 1
1616  if (k > 1 .AND. abs(fmp - petsub) > dem5 * petsub) then
1617  factor = factor / (fm / petsub)
1618  end if
1619  !
1620  ! -- one wave shallower than extdp
1621  if (this%nwavst(icell) == 1 .AND. &
1622  this%uzdpst(1, icell) <= this%extdpuz(icell)) then
1623  if (ietflag == 2) then
1624  tho = this%uzthst(1, icell)
1625  fktho = this%uzflst(1, icell)
1626  hcap = this%caph(icell, tho)
1627  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1628  end if
1629  if ((this%uzthst(1, icell) - thetaout) > this%thtr(icell) + extwc1) then
1630  this%uzthst(1, icell) = this%uzthst(1, icell) - thetaout
1631  this%uzflst(1, icell) = &
1632  this%vks(icell) * (((this%uzthst(1, icell) - &
1633  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1634  else if (this%uzthst(1, icell) > this%thtr(icell) + extwc1) then
1635  this%uzthst(1, icell) = this%thtr(icell) + extwc1
1636  this%uzflst(1, icell) = &
1637  this%vks(icell) * (((this%uzthst(1, icell) - &
1638  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1639  end if
1640  !
1641  ! -- all waves shallower than extinction depth
1642  else if (this%nwavst(icell) > 1 .AND. &
1643  this%uzdpst(this%nwavst(icell), icell) > this%extdpuz(icell)) then
1644  if (ietflag == 2) then
1645  tho = this%uzthst(this%nwavst(icell), icell)
1646  fktho = this%uzflst(this%nwavst(icell), icell)
1647  hcap = this%caph(icell, tho)
1648  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1649  end if
1650  if (this%nwavst(icell) + 1 > this%nwav(icell)) then
1651  !
1652  ! -- too many waves error
1653  ierr = 1
1654  goto 500
1655  end if
1656  if (this%uzthst(this%nwavst(icell), icell) - thetaout > &
1657  this%thtr(icell) + extwc1) then
1658  this%uzthst(this%nwavst(icell) + 1, icell) = &
1659  this%uzthst(this%nwavst(icell), icell) - thetaout
1660  numadd = 1
1661  else if (this%uzthst(this%nwavst(icell), icell) > &
1662  this%thtr(icell) + extwc1) then
1663  this%uzthst(this%nwavst(icell) + 1, icell) = this%thtr(icell) + extwc1
1664  numadd = 1
1665  end if
1666  if (numadd == 1) then
1667  this%uzflst(this%nwavst(icell) + 1, icell) = &
1668  this%vks(icell) * &
1669  (((this%uzthst(this%nwavst(icell) + 1, icell) - &
1670  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1671  theta2 = this%uzthst(this%nwavst(icell) + 1, icell)
1672  flux2 = this%uzflst(this%nwavst(icell) + 1, icell)
1673  flux1 = this%uzflst(this%nwavst(icell), icell)
1674  theta1 = this%uzthst(this%nwavst(icell), icell)
1675  this%uzspst(this%nwavst(icell) + 1, icell) = &
1676  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1677  this%thtr(icell), this%eps(icell), this%vks(icell))
1678  this%uzdpst(this%nwavst(icell) + 1, icell) = this%extdpuz(icell)
1679  this%nwavst(icell) = this%nwavst(icell) + 1
1680  if (this%nwavst(icell) > this%nwav(icell)) then
1681  !
1682  ! -- too many waves error, deallocate temp arrays and return
1683  ierr = 1
1684  goto 500
1685  end if
1686  else
1687  numadd = 0
1688  end if
1689  !
1690  ! -- one wave below extinction depth
1691  else if (this%nwavst(icell) == 1) then
1692  if (this%nwavst(icell) + 1 > this%nwav(icell)) then
1693  !
1694  ! -- too many waves error
1695  ierr = 1
1696  goto 500
1697  end if
1698  if (ietflag == 2) then
1699  tho = this%uzthst(1, icell)
1700  fktho = this%uzflst(1, icell)
1701  hcap = this%caph(icell, tho)
1702  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1703  end if
1704  if ((this%uzthst(1, icell) - thetaout) > this%thtr(icell) + extwc1) then
1705  if (thetaout > dem30) then
1706  this%uzthst(2, icell) = this%uzthst(1, icell) - thetaout
1707  this%uzflst(2, icell) = &
1708  this%vks(icell) * (((this%uzthst(2, icell) - this%thtr(icell)) * &
1709  thtsrinv)**this%eps(icell))
1710  this%uzdpst(2, icell) = this%extdpuz(icell)
1711  theta2 = this%uzthst(2, icell)
1712  flux2 = this%uzflst(2, icell)
1713  flux1 = this%uzflst(1, icell)
1714  theta1 = this%uzthst(1, icell)
1715  this%uzspst(2, icell) = &
1716  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1717  this%thtr(icell), this%eps(icell), this%vks(icell))
1718  this%nwavst(icell) = this%nwavst(icell) + 1
1719  if (this%nwavst(icell) > this%nwav(icell)) then
1720  !
1721  ! -- too many waves error
1722  ierr = 1
1723  goto 500
1724  end if
1725  end if
1726  else if (this%uzthst(1, icell) > this%thtr(icell) + extwc1) then
1727  if (thetaout > dem30) then
1728  this%uzthst(2, icell) = this%thtr(icell) + extwc1
1729  this%uzflst(2, icell) = &
1730  this%vks(icell) * (((this%uzthst(2, icell) - &
1731  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1732  this%uzdpst(2, icell) = this%extdpuz(icell)
1733  theta2 = this%uzthst(2, icell)
1734  flux2 = this%uzflst(2, icell)
1735  flux1 = this%uzflst(1, icell)
1736  theta1 = this%uzthst(1, icell)
1737  this%uzspst(2, icell) = &
1738  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1739  this%thtr(icell), this%eps(icell), this%vks(icell))
1740  this%nwavst(icell) = this%nwavst(icell) + 1
1741  if (this%nwavst(icell) > this%nwav(icell)) then
1742  !
1743  ! -- too many waves error
1744  ierr = 1
1745  goto 500
1746  end if
1747  end if
1748  end if
1749  else
1750  !
1751  ! -- extinction depth splits waves
1752  if (this%uzdpst(1, icell) - this%extdpuz(icell) > dem7) then
1753  j = 2
1754  jk = 0
1755  !
1756  ! -- locate extinction depth between waves
1757  do while (jk == 0)
1758  diff = this%uzdpst(j, icell) - this%extdpuz(icell)
1759  if (diff > dzero) then
1760  j = j + 1
1761  else
1762  jk = 1
1763  end if
1764  end do
1765  kk = j
1766  if (this%uzthst(j, icell) > this%thtr(icell) + extwc1) then
1767  !
1768  ! -- create a wave at extinction depth
1769  if (abs(diff) > dem5) then
1770  if (this%nwavst(icell) + 1 > this%nwav(icell)) then
1771  !
1772  ! -- too many waves error
1773  ierr = 1
1774  goto 500
1775  end if
1776  call this%wave_shift(this, icell, icell, -1, &
1777  this%nwavst(icell) + 1, j, -1)
1778  this%uzdpst(j, icell) = this%extdpuz(icell)
1779  this%nwavst(icell) = this%nwavst(icell) + 1
1780  if (this%nwavst(icell) > this%nwav(icell)) then
1781  !
1782  ! -- too many waves error
1783  ierr = 1
1784  goto 500
1785  end if
1786  end if
1787  kk = j
1788  else
1789  jhold = this%nwavst(icell)
1790  i = j + 1
1791  do while (i < this%nwavst(icell))
1792  if (this%uzthst(i, icell) > this%thtr(icell) + extwc1) then
1793  jhold = i
1794  i = this%nwavst(icell) + 1
1795  end if
1796  i = i + 1
1797  end do
1798  j = jhold
1799  kk = jhold
1800  end if
1801  else
1802  kk = 1
1803  end if
1804  !
1805  ! -- all waves above extinction depth
1806  do while (kk <= this%nwavst(icell))
1807  if (ietflag == 2) then
1808  tho = this%uzthst(kk, icell)
1809  fktho = this%uzflst(kk, icell)
1810  hcap = this%caph(icell, tho)
1811  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1812  end if
1813  if (this%uzthst(kk, icell) > this%thtr(icell) + extwc1) then
1814  if (this%uzthst(kk, icell) - thetaout > &
1815  this%thtr(icell) + extwc1) then
1816  this%uzthst(kk, icell) = this%uzthst(kk, icell) - thetaout
1817  else if (this%uzthst(kk, icell) > this%thtr(icell) + extwc1) then
1818  this%uzthst(kk, icell) = this%thtr(icell) + extwc1
1819  end if
1820  if (kk == 1) then
1821  this%uzflst(kk, icell) = &
1822  this%vks(icell) * &
1823  (((this%uzthst(kk, icell) - &
1824  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1825  end if
1826  if (kk > 1) then
1827  flux1 = &
1828  this%vks(icell) * ((this%uzthst(kk - 1, icell) - &
1829  this%thtr(icell)) * thtsrinv)**this%eps(icell)
1830  flux2 = &
1831  this%vks(icell) * ((this%uzthst(kk, icell) - &
1832  this%thtr(icell)) * thtsrinv)**this%eps(icell)
1833  this%uzflst(kk, icell) = flux2
1834  theta2 = this%uzthst(kk, icell)
1835  theta1 = this%uzthst(kk - 1, icell)
1836  this%uzspst(kk, icell) = leadspeed(theta1, theta2, flux1, flux2, &
1837  this%thts(icell), &
1838  this%thtr(icell), &
1839  this%eps(icell), this%vks(icell))
1840  end if
1841  end if
1842  kk = kk + 1
1843  end do
1844  end if
1845  !
1846  ! -- calculate aet
1847  kj = 1
1848  do while (kj <= this%nwavst(icell) - 1)
1849  if (abs(this%uzthst(kj, icell) - this%uzthst(kj + 1, icell)) < dem6) then
1850  call this%wave_shift(this, icell, icell, 1, kj + 1, &
1851  this%nwavst(icell) - 1, 1)
1852  kj = kj - 1
1853  this%nwavst(icell) = this%nwavst(icell) - 1
1854  end if
1855  kj = kj + 1
1856  end do
1857  depth = this%uzdpst(1, icell)
1858  fm = this%unsat_stor(icell, depth)
1859  this%etact(icell) = st - fm
1860  fm = this%etact(icell) / delt
1861  if (this%etact(icell) < dzero) then
1862  call this%wave_shift(uzfktemp, icell, 1, 0, 1, nwv, 1)
1863  this%nwavst(icell) = nwv
1864  this%etact(icell) = dzero
1865  elseif (petsub - fm < -dem15 .AND. ietflag == 2) then
1866  !
1867  ! -- aet greater than pet, reset and try again
1868  call this%wave_shift(uzfktemp, icell, 1, 0, 1, nwv, 1)
1869  this%nwavst(icell) = nwv
1870  this%etact(icell) = dzero
1871  else
1872  itest = 1
1873  end if
1874  !
1875  ! -- end aet-pet loop for head dependent et
1876  fmp = fm
1877  if (k > 100) then
1878  itest = 1
1879  elseif (ietflag < 2) then
1880  fmp = petsub
1881  itest = 1
1882  end if
1883  end do
1884 500 continue
1885  !
1886  ! -- deallocate temporary worker
1887  call uzfktemp%dealloc()
1888  end subroutine uzet
1889 
1890  !> @brief Calculate capillary pressure head from B-C equation
1891  !<
1892  function caph(this, icell, tho)
1893  ! -- dummy
1894  class(uzfcellgrouptype) :: this
1895  integer(I4B), intent(in) :: icell
1896  real(dp), intent(in) :: tho
1897  ! -- local
1898  real(dp) :: caph, lambda, star
1899  !
1900  caph = -dem6
1901  star = (tho - this%thtr(icell)) / (this%thts(icell) - this%thtr(icell))
1902  if (star < dem15) star = dem15
1903  lambda = dtwo / (this%eps(icell) - dthree)
1904  if (star > dem15) then
1905  if (tho - this%thts(icell) < dem15) then
1906  caph = this%ha(icell) * star**(-done / lambda)
1907  else
1908  caph = dzero
1909  end if
1910  end if
1911  end function caph
1912 
1913  !> @brief Calculate capillary pressure-based uz et
1914  function rate_et_z(this, icell, factor, fktho, h)
1915  ! -- Return
1916  real(dp) :: rate_et_z
1917  ! -- dummy
1918  class(uzfcellgrouptype) :: this
1919  integer(I4B), intent(in) :: icell
1920  real(dp), intent(in) :: factor, fktho, h
1921  !
1922  rate_et_z = factor * fktho * (h - this%hroot(icell))
1923  if (rate_et_z < dzero) rate_et_z = dzero
1924  end function rate_et_z
1925 
1926  !> @brief Determine the water content at a specific depth
1927  !!
1928  !! Because UZF-calculated waves are internal to UZF objects, different water
1929  !! contents exists at different depths.
1930  !<
1931  function get_water_content_at_depth(this, icell, depth) result(theta_at_depth)
1932  ! -- dummy
1933  class(uzfcellgrouptype) :: this
1934  integer(I4B), intent(in) :: icell !< uzf cell containing depth
1935  real(dp), intent(in) :: depth !< depth within the cell
1936  ! -- return
1937  real(dp) :: theta_at_depth
1938  ! -- local
1939  real(dp) :: d1
1940  real(dp) :: d2
1941  real(dp) :: f1
1942  real(dp) :: f2
1943  !
1944  if (this%watab(icell) < this%celtop(icell)) then
1945  if (this%celtop(icell) - depth > this%watab(icell)) then
1946  d1 = depth - dem3
1947  d2 = depth + dem3
1948  f1 = this%unsat_stor(icell, d1)
1949  f2 = this%unsat_stor(icell, d2)
1950  theta_at_depth = this%thtr(icell) + (f2 - f1) / (d2 - d1)
1951  else
1952  theta_at_depth = this%thts(icell)
1953  end if
1954  else
1955  theta_at_depth = this%thts(icell)
1956  end if
1957  end function get_water_content_at_depth
1958 
1959  !> @brief Calculate and return the cell-based water content value
1960  !<
1961  function get_wcnew(this, icell) result(watercontent)
1962  ! -- dummy
1963  class(uzfcellgrouptype) :: this
1964  integer(I4B), intent(in) :: icell !< uzf cell containing depth
1965  ! -- return
1966  real(dp) :: watercontent
1967  ! -- local
1968  real(dp) :: top
1969  real(dp) :: bot
1970  real(dp) :: theta_r
1971  real(dp) :: thk
1972  real(dp) :: hgwf
1973  real(dp) :: fm
1974  real(dp) :: d
1975  !
1976  hgwf = this%watab(icell)
1977  top = this%celtop(icell)
1978  bot = this%celbot(icell)
1979  thk = top - max(bot, hgwf)
1980  if (thk > dzero) then
1981  theta_r = this%thtr(icell)
1982  d = thk
1983  fm = this%unsat_stor(icell, d)
1984  watercontent = fm / thk
1985  watercontent = watercontent + theta_r
1986  else
1987  watercontent = dzero
1988  end if
1989  end function get_wcnew
1990 
1991 end module uzfcellgroupmodule
init
subroutine init()
Definition: GridSorting.f90:25
memorymanagermodule::mem_allocate
Definition: MemoryManager.f90:55
memorymanagermodule::mem_deallocate
Definition: MemoryManager.f90:129
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:145
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:25
tdismodule::kper
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:26
tdismodule::delt
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:32
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:1465
uzfcellgroupmodule::factors
subroutine factors(feps1, feps2)
Calculate unit specific tolerances.
Definition: UzfCellGroup.f90:1070
uzfcellgroupmodule::uzflow
subroutine uzflow(this, thick, thickold, delt, ietflag, icell, ierr)
Method of Characteristics solution for kinematic wave equation.
Definition: UzfCellGroup.f90:981
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:772
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:1185
uzfcellgroupmodule::caph
real(dp) function caph(this, icell, tho)
Calculate capillary pressure head from B-C equation.
Definition: UzfCellGroup.f90:1893
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:705
uzfcellgroupmodule::setwaves
subroutine setwaves(this, icell)
Reset waves to default values at start of simulation.
Definition: UzfCellGroup.f90:869
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:1385
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:1932
uzfcellgroupmodule::wave_shift
subroutine wave_shift(this, this2, icell, icell2, shft, strt, stp, cntr)
Copy waves or shift waves in arrays.
Definition: UzfCellGroup.f90:956
uzfcellgroupmodule::rejfinf
subroutine rejfinf(this, icell, deriv, hgwf, trhs, thcof, finfact)
Reject applied infiltration due to low vks.
Definition: UzfCellGroup.f90:744
uzfcellgroupmodule::rate_et_z
real(dp) function rate_et_z(this, icell, factor, fktho, h)
Calculate capillary pressure-based uz et.
Definition: UzfCellGroup.f90:1915
uzfcellgroupmodule::simgwet
subroutine simgwet(this, igwetflag, icell, hgwf, trhs, thcof, det)
Calculate gwf et using residual uzf pet.
Definition: UzfCellGroup.f90:812
uzfcellgroupmodule::get_wcnew
real(dp) function get_wcnew(this, icell)
Calculate and return the cell-based water content value.
Definition: UzfCellGroup.f90:1962
uzfcellgroupmodule::uz_rise
subroutine uz_rise(this, icell, totfluxtot)
Calculate recharge due to a rise in the gwf head.
Definition: UzfCellGroup.f90:849
uzfcellgroupmodule::uzet
subroutine uzet(this, icell, delt, ietflag, ierr)
Remove water from uz due to et.
Definition: UzfCellGroup.f90:1541
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:918
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:1421
uzfcellgroupmodule::trailwav
subroutine trailwav(this, icell, ierr)
Create and set trail waves.
Definition: UzfCellGroup.f90:1098
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