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