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