16 integer(I4B) :: imem_manager
17 real(dp),
pointer,
dimension(:),
contiguous :: thtr => null()
18 real(dp),
pointer,
dimension(:),
contiguous :: thts => null()
19 real(dp),
pointer,
dimension(:),
contiguous :: thti => null()
20 real(dp),
pointer,
dimension(:),
contiguous :: eps => null()
21 real(dp),
pointer,
dimension(:),
contiguous :: extwc => null()
22 real(dp),
pointer,
dimension(:),
contiguous :: ha => null()
23 real(dp),
pointer,
dimension(:),
contiguous :: hroot => null()
24 real(dp),
pointer,
dimension(:),
contiguous :: rootact => null()
25 real(dp),
pointer,
dimension(:),
contiguous :: etact => null()
26 real(dp),
dimension(:, :),
pointer,
contiguous :: uzspst => null()
27 real(dp),
dimension(:, :),
pointer,
contiguous :: uzthst => null()
28 real(dp),
dimension(:, :),
pointer,
contiguous :: uzflst => null()
29 real(dp),
dimension(:, :),
pointer,
contiguous :: uzdpst => null()
30 integer(I4B),
pointer,
dimension(:),
contiguous :: nwavst => null()
31 real(dp),
pointer,
dimension(:),
contiguous :: totflux => null()
32 integer(I4B),
pointer,
dimension(:),
contiguous :: nwav => null()
33 integer(I4B),
pointer,
dimension(:),
contiguous :: ntrail => null()
34 real(dp),
pointer,
dimension(:),
contiguous :: sinf => null()
35 real(dp),
pointer,
dimension(:),
contiguous :: finf => null()
36 real(dp),
pointer,
dimension(:),
contiguous :: pet => null()
37 real(dp),
pointer,
dimension(:),
contiguous :: petmax => null()
38 real(dp),
pointer,
dimension(:),
contiguous :: extdp => null()
39 real(dp),
pointer,
dimension(:),
contiguous :: extdpuz => null()
40 real(dp),
pointer,
dimension(:),
contiguous :: finf_rej => null()
41 real(dp),
pointer,
dimension(:),
contiguous :: gwet => null()
42 real(dp),
pointer,
dimension(:),
contiguous :: uzfarea => null()
43 real(dp),
pointer,
dimension(:),
contiguous :: cellarea => null()
44 real(dp),
pointer,
dimension(:),
contiguous :: celtop => null()
45 real(dp),
pointer,
dimension(:),
contiguous :: celbot => null()
46 real(dp),
pointer,
dimension(:),
contiguous :: landtop => null()
47 real(dp),
pointer,
dimension(:),
contiguous :: watab => null()
48 real(dp),
pointer,
dimension(:),
contiguous :: watabold => null()
49 real(dp),
pointer,
dimension(:),
contiguous :: vks => null()
50 real(dp),
pointer,
dimension(:),
contiguous :: surfdep => null()
51 real(dp),
pointer,
dimension(:),
contiguous :: surflux => null()
52 real(dp),
pointer,
dimension(:),
contiguous :: surfluxbelow => null()
53 real(dp),
pointer,
dimension(:),
contiguous :: surfseep => null()
54 real(dp),
pointer,
dimension(:),
contiguous :: gwpet => null()
55 integer(I4B),
pointer,
dimension(:),
contiguous :: landflag => null()
56 integer(I4B),
pointer,
dimension(:),
contiguous :: ivertcon => null()
97 subroutine init(this, ncells, nwav, memory_path)
102 integer(I4B),
intent(in) :: nwav
103 integer(I4B),
intent(in) :: ncells
104 character(len=*),
intent(in),
optional :: memory_path
106 integer(I4B) :: icell
109 if (
present(memory_path))
then
110 this%imem_manager = 1
111 call mem_allocate(this%uzdpst, nwav, ncells,
'UZDPST', memory_path)
112 call mem_allocate(this%uzthst, nwav, ncells,
'UZTHST', memory_path)
113 call mem_allocate(this%uzflst, nwav, ncells,
'UZFLST', memory_path)
114 call mem_allocate(this%uzspst, nwav, ncells,
'UZSPST', memory_path)
115 call mem_allocate(this%nwavst, ncells,
'NWAVST', memory_path)
116 call mem_allocate(this%thtr, ncells,
'THTR', memory_path)
117 call mem_allocate(this%thts, ncells,
'THTS', memory_path)
118 call mem_allocate(this%thti, ncells,
'THTI', memory_path)
121 call mem_allocate(this%hroot, ncells,
'HROOT', memory_path)
122 call mem_allocate(this%rootact, ncells,
'ROOTACT', memory_path)
123 call mem_allocate(this%extwc, ncells,
'EXTWC', memory_path)
124 call mem_allocate(this%etact, ncells,
'ETACT', memory_path)
125 call mem_allocate(this%nwav, ncells,
'NWAV', memory_path)
126 call mem_allocate(this%ntrail, ncells,
'NTRAIL', memory_path)
127 call mem_allocate(this%totflux, ncells,
'TOTFLUX', memory_path)
128 call mem_allocate(this%sinf, ncells,
'SINF', memory_path)
129 call mem_allocate(this%finf, ncells,
'FINF', memory_path)
130 call mem_allocate(this%finf_rej, ncells,
'FINF_REJ', memory_path)
131 call mem_allocate(this%gwet, ncells,
'GWET', memory_path)
132 call mem_allocate(this%uzfarea, ncells,
'UZFAREA', memory_path)
133 call mem_allocate(this%cellarea, ncells,
'CELLAREA', memory_path)
134 call mem_allocate(this%celtop, ncells,
'CELTOP', memory_path)
135 call mem_allocate(this%celbot, ncells,
'CELBOT', memory_path)
136 call mem_allocate(this%landtop, ncells,
'LANDTOP', memory_path)
137 call mem_allocate(this%watab, ncells,
'WATAB', memory_path)
138 call mem_allocate(this%watabold, ncells,
'WATABOLD', memory_path)
139 call mem_allocate(this%surfdep, ncells,
'SURFDEP', memory_path)
141 call mem_allocate(this%surflux, ncells,
'SURFLUX', memory_path)
142 call mem_allocate(this%surfluxbelow, ncells,
'SURFLUXBELOW', memory_path)
143 call mem_allocate(this%surfseep, ncells,
'SURFSEEP', memory_path)
144 call mem_allocate(this%gwpet, ncells,
'GWPET', memory_path)
146 call mem_allocate(this%petmax, ncells,
'PETMAX', memory_path)
147 call mem_allocate(this%extdp, ncells,
'EXTDP', memory_path)
148 call mem_allocate(this%extdpuz, ncells,
'EXTDPUZ', memory_path)
149 call mem_allocate(this%landflag, ncells,
'LANDFLAG', memory_path)
150 call mem_allocate(this%ivertcon, ncells,
'IVERTCON', memory_path)
152 this%imem_manager = 0
153 allocate (this%uzdpst(nwav, ncells))
154 allocate (this%uzthst(nwav, ncells))
155 allocate (this%uzflst(nwav, ncells))
156 allocate (this%uzspst(nwav, ncells))
157 allocate (this%nwavst(ncells))
158 allocate (this%thtr(ncells))
159 allocate (this%thts(ncells))
160 allocate (this%thti(ncells))
161 allocate (this%eps(ncells))
162 allocate (this%ha(ncells))
163 allocate (this%hroot(ncells))
164 allocate (this%rootact(ncells))
165 allocate (this%extwc(ncells))
166 allocate (this%etact(ncells))
167 allocate (this%nwav(ncells))
168 allocate (this%ntrail(ncells))
169 allocate (this%totflux(ncells))
170 allocate (this%sinf(ncells))
171 allocate (this%finf(ncells))
172 allocate (this%finf_rej(ncells))
173 allocate (this%gwet(ncells))
174 allocate (this%uzfarea(ncells))
175 allocate (this%cellarea(ncells))
176 allocate (this%celtop(ncells))
177 allocate (this%celbot(ncells))
178 allocate (this%landtop(ncells))
179 allocate (this%watab(ncells))
180 allocate (this%watabold(ncells))
181 allocate (this%surfdep(ncells))
182 allocate (this%vks(ncells))
183 allocate (this%surflux(ncells))
184 allocate (this%surfluxbelow(ncells))
185 allocate (this%surfseep(ncells))
186 allocate (this%gwpet(ncells))
187 allocate (this%pet(ncells))
188 allocate (this%petmax(ncells))
189 allocate (this%extdp(ncells))
190 allocate (this%extdpuz(ncells))
191 allocate (this%landflag(ncells))
192 allocate (this%ivertcon(ncells))
195 this%uzdpst(:, icell) =
dzero
196 this%uzthst(:, icell) =
dzero
197 this%uzflst(:, icell) =
dzero
198 this%uzspst(:, icell) =
dzero
199 this%nwavst(icell) = 1
200 this%thtr(icell) =
dzero
201 this%thts(icell) =
dzero
202 this%thti(icell) =
dzero
203 this%eps(icell) =
dzero
204 this%ha(icell) =
dzero
205 this%hroot(icell) =
dzero
206 this%rootact(icell) =
dzero
207 this%extwc(icell) =
dzero
208 this%etact(icell) =
dzero
209 this%nwav(icell) = nwav
210 this%ntrail(icell) = 0
211 this%totflux(icell) =
dzero
212 this%sinf(icell) =
dzero
213 this%finf(icell) =
dzero
214 this%finf_rej(icell) =
dzero
215 this%gwet(icell) =
dzero
216 this%uzfarea(icell) =
dzero
217 this%cellarea(icell) =
dzero
218 this%celtop(icell) =
dzero
219 this%celbot(icell) =
dzero
220 this%landtop(icell) =
dzero
221 this%watab(icell) =
dzero
222 this%watabold(icell) =
dzero
223 this%surfdep(icell) =
dzero
224 this%vks(icell) =
dzero
225 this%surflux(icell) =
dzero
226 this%surfluxbelow(icell) =
dzero
227 this%surfseep(icell) =
dzero
228 this%gwpet(icell) =
dzero
229 this%pet(icell) =
dzero
230 this%petmax(icell) =
dzero
231 this%extdp(icell) =
dzero
232 this%extdpuz(icell) =
dzero
233 this%landflag(icell) = 0
234 this%ivertcon(icell) = 0
250 if (this%imem_manager == 0)
then
251 deallocate (this%uzdpst)
252 deallocate (this%uzthst)
253 deallocate (this%uzflst)
254 deallocate (this%uzspst)
255 deallocate (this%nwavst)
256 deallocate (this%thtr)
257 deallocate (this%thts)
258 deallocate (this%thti)
259 deallocate (this%eps)
261 deallocate (this%hroot)
262 deallocate (this%rootact)
263 deallocate (this%extwc)
264 deallocate (this%etact)
265 deallocate (this%nwav)
266 deallocate (this%ntrail)
267 deallocate (this%totflux)
268 deallocate (this%sinf)
269 deallocate (this%finf)
270 deallocate (this%finf_rej)
271 deallocate (this%gwet)
272 deallocate (this%uzfarea)
273 deallocate (this%cellarea)
274 deallocate (this%celtop)
275 deallocate (this%celbot)
276 deallocate (this%landtop)
277 deallocate (this%watab)
278 deallocate (this%watabold)
279 deallocate (this%surfdep)
280 deallocate (this%vks)
281 deallocate (this%surflux)
282 deallocate (this%surfluxbelow)
283 deallocate (this%surfseep)
284 deallocate (this%gwpet)
285 deallocate (this%pet)
286 deallocate (this%petmax)
287 deallocate (this%extdp)
288 deallocate (this%extdpuz)
289 deallocate (this%landflag)
290 deallocate (this%ivertcon)
340 subroutine setdata(this, icell, area, top, bot, surfdep, vks, thtr, thts, &
341 thti, eps, ntrail, landflag, ivertcon)
344 integer(I4B),
intent(in) :: icell
345 real(DP),
intent(in) :: area
346 real(DP),
intent(in) :: top
347 real(DP),
intent(in) :: bot
348 real(DP),
intent(in) :: surfdep
349 real(DP),
intent(in) :: vks
350 real(DP),
intent(in) :: thtr
351 real(DP),
intent(in) :: thts
352 real(DP),
intent(in) :: thti
353 real(DP),
intent(in) :: eps
354 integer(I4B),
intent(in) :: ntrail
355 integer(I4B),
intent(in) :: landflag
356 integer(I4B),
intent(in) :: ivertcon
359 this%landflag(icell) = landflag
360 this%ivertcon(icell) = ivertcon
361 this%surfdep(icell) = surfdep
362 this%uzfarea(icell) = area
363 this%cellarea(icell) = area
364 if (this%landflag(icell) == 1)
then
365 this%celtop(icell) = top -
dhalf * this%surfdep(icell)
367 this%celtop(icell) = top
369 this%celbot(icell) = bot
370 this%vks(icell) = vks
371 this%thtr(icell) = thtr
372 this%thts(icell) = thts
373 this%thti(icell) = thti
374 this%eps(icell) = eps
375 this%ntrail(icell) = ntrail
376 this%pet(icell) =
dzero
377 this%extdp(icell) =
dzero
378 this%extwc(icell) =
dzero
379 this%ha(icell) =
dzero
380 this%hroot(icell) =
dzero
388 integer(I4B),
intent(in) :: icell
389 real(DP),
intent(in) :: hgwf
392 this%watab(icell) = this%celbot(icell)
393 if (hgwf > this%celbot(icell)) this%watab(icell) = hgwf
394 if (this%watab(icell) > this%celtop(icell)) &
395 this%watab(icell) = this%celtop(icell)
396 this%watabold(icell) = this%watab(icell)
407 integer(I4B),
intent(in) :: icell
408 real(DP),
intent(in) :: finf
410 if (this%landflag(icell) == 1)
then
411 this%sinf(icell) = finf
412 this%finf(icell) = finf
414 this%sinf(icell) =
dzero
415 this%finf(icell) =
dzero
417 this%finf_rej(icell) =
dzero
418 this%surflux(icell) =
dzero
419 this%surfluxbelow(icell) =
dzero
430 integer(I4B),
intent(in) :: icell
431 real(DP),
intent(in) :: areamult
434 this%uzfarea(icell) = this%cellarea(icell) * areamult
445 integer(I4B),
intent(in) :: icell
446 integer(I4B),
intent(in) :: jbelow
447 real(DP),
intent(in) :: pet
448 real(DP),
intent(in) :: extdp
452 if (this%landflag(icell) == 1)
then
453 this%pet(icell) = pet
454 this%gwpet(icell) = pet
456 this%pet(icell) =
dzero
457 this%gwpet(icell) =
dzero
459 thick = this%celtop(icell) - this%celbot(icell)
460 this%extdp(icell) = extdp
461 if (this%landflag(icell) > 0)
then
462 this%landtop(icell) = this%celtop(icell)
463 this%petmax(icell) = this%pet(icell)
467 if (this%landtop(icell) - this%extdp(icell) < this%celbot(icell))
then
468 this%extdpuz(icell) = thick
470 this%extdpuz(icell) = this%celtop(icell) - &
471 (this%landtop(icell) - this%extdp(icell))
473 if (this%extdpuz(icell) <
dzero) this%extdpuz(icell) =
dzero
474 if (this%extdpuz(icell) >
dem7 .and. this%extdp(icell) <
dem7) &
475 this%extdp(icell) = this%extdpuz(icell)
479 this%landtop(jbelow) = this%landtop(icell)
480 this%petmax(jbelow) = this%petmax(icell)
494 integer(I4B),
intent(in) :: icell
501 pet = this%pet(icell) - this%etact(icell) /
delt
503 this%gwpet(icell) = pet
516 integer(I4B),
intent(in) :: icell
517 integer(I4B),
intent(in) :: jbelow
525 if (this%extdpuz(jbelow) >
dem3)
then
526 pet = this%pet(icell) - this%etact(icell) /
delt - &
527 this%gwet(icell) / this%uzfarea(icell)
530 this%pet(jbelow) = pet
541 integer(I4B),
intent(in) :: icell
542 integer(I4B),
intent(in) :: jbelow
543 real(DP),
intent(in) :: extwc
546 this%extwc(icell) = extwc
547 if (jbelow > 0) this%extwc(jbelow) = extwc
558 integer(I4B),
intent(in) :: icell
559 integer(I4B),
intent(in) :: jbelow
560 real(DP),
intent(in) :: ha
561 real(DP),
intent(in) :: hroot
562 real(DP),
intent(in) :: rootact
566 this%hroot(icell) = hroot
567 this%rootact(icell) = rootact
570 this%hroot(jbelow) = hroot
571 this%rootact(jbelow) = rootact
583 integer(I4B),
intent(in) :: icell
586 this%surfseep(icell) =
dzero
595 subroutine solve(this, thiswork, jbelow, icell, totfluxtot, ietflag, &
596 issflag, iseepflag, hgwf, qfrommvr, ierr, &
597 reset_state, trhs, thcof, deriv, watercontent)
603 integer(I4B),
intent(in) :: jbelow
604 integer(I4B),
intent(in) :: icell
605 real(DP),
intent(inout) :: totfluxtot
606 integer(I4B),
intent(in) :: ietflag
607 integer(I4B),
intent(in) :: issflag
608 integer(I4B),
intent(in) :: iseepflag
609 real(DP),
intent(in) :: hgwf
610 real(DP),
intent(in) :: qfrommvr
611 integer(I4B),
intent(inout) :: ierr
612 logical,
intent(in) :: reset_state
613 real(DP),
intent(inout),
optional :: trhs
614 real(DP),
intent(inout),
optional :: thcof
615 real(DP),
intent(inout),
optional :: deriv
616 real(DP),
intent(inout),
optional :: watercontent
622 real(DP) :: derivfinf
624 real(DP) :: thcoffinf
626 real(DP) :: thcofseep
636 this%finf_rej(icell) =
dzero
637 this%surflux(icell) = this%finf(icell) + qfrommvr / this%uzfarea(icell)
638 this%watab(icell) = hgwf
639 this%surfseep(icell) =
dzero
642 this%etact(icell) =
dzero
643 this%surfluxbelow(icell) =
dzero
644 if (this%ivertcon(icell) > 0)
then
645 this%finf(jbelow) =
dzero
647 if (this%watab(icell) < this%celbot(icell)) &
648 this%watab(icell) = this%celbot(icell)
656 if (reset_state)
then
657 call thiswork%wave_shift(this, 1, icell, 0, 1, this%nwavst(icell), 1)
660 if (this%watab(icell) > this%celtop(icell)) &
661 this%watab(icell) = this%celtop(icell)
664 if (this%surflux(icell) > this%vks(icell))
then
665 this%surflux(icell) = this%vks(icell)
669 if (this%landflag(icell) == 1)
then
670 call this%rejfinf(icell, deriv1, hgwf, trhsfinf, thcoffinf, finfact)
671 this%surflux(icell) = finfact
675 this%finf_rej(icell) = this%finf(icell) + &
676 (qfrommvr / this%uzfarea(icell)) - this%surflux(icell)
679 if (iseepflag > 0 .and. this%landflag(icell) == 1)
then
680 call this%gwseep(icell, deriv2, scale, hgwf, trhsseep, thcofseep, seep)
681 this%surfseep(icell) = seep
685 test = this%watab(icell)
686 if (this%watabold(icell) - test < -
dem15) test = this%watabold(icell)
687 if (this%celtop(icell) - test >
dem15)
then
688 if (issflag == 0)
then
689 call this%routewaves(totfluxtot,
delt, ietflag, icell, ierr)
691 call this%uz_rise(icell, totfluxtot)
692 this%totflux(icell) = totfluxtot
693 if (this%ivertcon(icell) > 0)
then
694 call this%addrech(icell, jbelow, hgwf, trhsfinf, thcoffinf, &
698 this%totflux(icell) = this%surflux(icell) *
delt
699 totfluxtot = this%surflux(icell) *
delt
702 trhsfinf = this%totflux(icell) * this%uzfarea(icell) /
delt
703 if (.not. reset_state)
then
704 call this%update_wav(icell,
delt, issflag, 0)
707 this%totflux(icell) = this%surflux(icell) *
delt
708 totfluxtot = this%surflux(icell) *
delt
709 if (.not. reset_state)
then
710 call this%update_wav(icell,
delt, issflag, 1)
715 if (
present(deriv)) deriv = deriv1 + deriv2 + derivfinf
716 if (
present(trhs)) trhs = trhsfinf + trhsseep
717 if (
present(thcof)) thcof = thcoffinf + thcofseep
720 if (
present(watercontent))
then
721 watercontent = this%get_wcnew(icell)
725 if (reset_state)
then
726 call this%wave_shift(thiswork, icell, 1, 0, 1, thiswork%nwavst(1), 1)
735 subroutine addrech(this, icell, jbelow, hgwf, trhs, thcof, deriv, delt)
738 integer(I4B),
intent(in) :: icell
739 integer(I4B),
intent(in) :: jbelow
740 real(DP),
intent(inout) :: trhs
741 real(DP),
intent(inout) :: thcof
742 real(DP),
intent(inout) :: deriv
743 real(DP),
intent(in) :: delt
744 real(DP),
intent(in) :: hgwf
747 real(DP) :: x, scale, range
753 trhs = this%uzfarea(icell) * this%totflux(icell) / delt
754 if (this%totflux(icell) <
dem14)
return
758 x = hgwf - (this%celbot(icell) - range)
759 call sscurve(x, range, deriv, scale)
760 deriv = this%uzfarea(icell) * deriv * this%totflux(icell) / delt
761 this%finf(jbelow) = (
done - scale) * this%totflux(icell) / delt
762 fcheck = this%finf(jbelow) - this%vks(jbelow)
766 this%finf(jbelow) = this%finf(jbelow) - fcheck
767 this%surfluxbelow(icell) = this%finf(jbelow)
768 this%totflux(icell) = scale * this%totflux(icell) + fcheck * delt
769 trhs = this%uzfarea(icell) * this%totflux(icell) / delt
777 subroutine rejfinf(this, icell, deriv, hgwf, trhs, thcof, finfact)
780 integer(I4B),
intent(in) :: icell
781 real(DP),
intent(inout) :: deriv
782 real(DP),
intent(inout) :: finfact
783 real(DP),
intent(inout) :: thcof
784 real(DP),
intent(inout) :: trhs
785 real(DP),
intent(in) :: hgwf
787 real(DP) :: x, range, scale, q
789 range = this%surfdep(icell)
790 q = this%surflux(icell)
792 trhs = finfact * this%uzfarea(icell)
793 x = this%celtop(icell) - hgwf
794 call slinear(x, range, deriv, scale)
795 deriv = -q * deriv * this%uzfarea(icell) * scale
796 if (scale <
done)
then
798 trhs = finfact * this%uzfarea(icell) * this%celtop(icell) / range
799 thcof = finfact * this%uzfarea(icell) / range
808 subroutine gwseep(this, icell, deriv, scale, hgwf, trhs, thcof, seep)
811 integer(I4B),
intent(in) :: icell
812 real(DP),
intent(inout) :: deriv
813 real(DP),
intent(inout) :: trhs
814 real(DP),
intent(inout) :: thcof
815 real(DP),
intent(inout) :: seep
816 real(DP),
intent(out) :: scale
817 real(DP),
intent(in) :: hgwf
819 real(DP) :: x, range, y, deriv1, d1, d2, Q
827 q = this%uzfarea(icell) * this%vks(icell)
828 range = this%surfdep(icell)
829 x = hgwf - this%celtop(icell)
832 seep = scale * q * (hgwf - this%celtop(icell)) / range
833 trhs = scale * q * this%celtop(icell) / range
834 thcof = -scale * q / range
835 d1 = -deriv1 * q * x / range
836 d2 = -scale * q / range
838 if (seep <
dzero)
then
851 subroutine simgwet(this, igwetflag, icell, hgwf, trhs, thcof, det)
854 integer(I4B),
intent(in) :: igwetflag
855 integer(I4B),
intent(in) :: icell
856 real(DP),
intent(in) :: hgwf
857 real(DP),
intent(inout) :: trhs
858 real(DP),
intent(inout) :: thcof
859 real(DP),
intent(inout) :: det
861 real(DP) :: s, x, c, b, et
863 this%gwet(icell) =
dzero
867 s = this%landtop(icell)
868 x = this%extdp(icell)
869 c = this%gwpet(icell)
870 b = this%celbot(icell)
873 if (igwetflag == 1)
then
874 et =
etfunc_lin(s, x, c, det, trhs, thcof, hgwf, &
875 this%celtop(icell), this%celbot(icell))
876 else if (igwetflag == 2)
then
880 trhs = trhs * this%uzfarea(icell)
881 thcof = thcof * this%uzfarea(icell)
882 this%gwet(icell) = trhs - (thcof * hgwf)
891 function etfunc_lin(s, x, c, det, trhs, thcof, hgwf, celtop, celbot)
895 real(dp),
intent(in) :: s
896 real(dp),
intent(in) :: x
897 real(dp),
intent(in) :: c
898 real(dp),
intent(inout) :: det
899 real(dp),
intent(inout) :: trhs
900 real(dp),
intent(inout) :: thcof
901 real(dp),
intent(in) :: hgwf
902 real(dp),
intent(in) :: celtop
903 real(dp),
intent(in) :: celbot
907 real(dp) :: depth, scale, thick
910 if (hgwf > (s - x) .and. hgwf < s)
THEN
911 etgw = (c * (hgwf - (s - x)) / x)
917 etgw = trhs - (thcof * hgwf)
921 else if (hgwf >= s)
then
931 depth = hgwf - (s - x)
932 thick = celtop - celbot
933 if (depth > thick) depth = thick
936 call scubic(depth, range, det, scale)
938 thcof = scale * thcof
939 etgw = trhs - (thcof * hgwf)
953 real(dp),
intent(in) :: s
954 real(dp),
intent(in) :: x
955 real(dp),
intent(in) :: c
956 real(dp),
intent(inout) :: det
957 real(dp),
intent(inout) :: trhs
958 real(dp),
intent(inout) :: thcof
959 real(dp),
intent(in) :: hgwf
963 real(dp) :: depth, scale
965 depth = hgwf - (s - x)
969 call scubic(depth, range, det, scale)
984 integer(I4B),
intent(in) :: icell
985 real(DP),
intent(inout) :: totfluxtot
987 real(DP) :: fm1, fm2, d1
990 if (this%watab(icell) - this%watabold(icell) >
dem30)
then
991 d1 = this%celtop(icell) - this%watabold(icell)
992 fm1 = this%unsat_stor(icell, d1)
993 d1 = this%celtop(icell) - this%watab(icell)
994 fm2 = this%unsat_stor(icell, d1)
995 totfluxtot = totfluxtot + (fm1 - fm2)
1008 integer(I4B),
intent(in) :: icell
1009 real(DP) :: bottom, top
1014 this%totflux(icell) =
dzero
1015 this%nwavst(icell) = 1
1016 this%uzdpst(:, icell) =
dzero
1017 thick = this%celtop(icell) - this%watab(icell)
1018 do jk = 1, this%nwav(icell)
1019 this%uzthst(jk, icell) = this%thtr(icell)
1023 if (thick >
dzero)
then
1024 this%uzdpst(1, icell) = thick
1025 this%uzthst(1, icell) = this%thti(icell)
1026 top = this%uzthst(1, icell) - this%thtr(icell)
1028 bottom = this%thts(icell) - this%thtr(icell)
1030 this%uzflst(1, icell) = this%vks(icell) * (top / bottom)**this%eps(icell)
1031 if (this%uzthst(1, icell) < this%thtr(icell)) &
1032 this%uzthst(1, icell) = this%thtr(icell)
1035 if (top >
dzero)
then
1036 this%uzspst(1, icell) =
dzero
1038 this%uzflst(1, icell) =
dzero
1039 this%uzspst(1, icell) =
dzero
1044 this%uzflst(1, icell) =
dzero
1045 this%uzdpst(1, icell) =
dzero
1046 this%uzspst(1, icell) =
dzero
1047 this%uzthst(1, icell) = this%thtr(icell)
1056 subroutine routewaves(this, totfluxtot, delt, ietflag, icell, ierr)
1059 real(DP),
intent(inout) :: totfluxtot
1060 real(DP),
intent(in) :: delt
1061 integer(I4B),
intent(in) :: ietflag
1062 integer(I4B),
intent(in) :: icell
1063 integer(I4B),
intent(inout) :: ierr
1065 real(DP) :: thick, thickold
1066 integer(I4B) :: idelt, iwav, ik
1069 this%totflux(icell) =
dzero
1070 this%etact(icell) =
dzero
1071 thick = this%celtop(icell) - this%watab(icell)
1072 thickold = this%celtop(icell) - this%watabold(icell)
1075 if (thickold <
dzero)
then
1077 this%uzthst(iwav, icell) = this%thtr(icell)
1078 this%uzdpst(iwav, icell) =
dzero
1079 this%uzspst(iwav, icell) =
dzero
1080 this%uzflst(iwav, icell) =
dzero
1081 this%nwavst(icell) = 1
1086 call this%uzflow(thick, thickold, delt, ietflag, icell, ierr)
1087 if (ierr > 0)
return
1088 totfluxtot = totfluxtot + this%totflux(icell)
1097 subroutine wave_shift(this, this2, icell, icell2, shft, strt, stp, cntr)
1101 integer(I4B),
intent(in) :: icell
1102 integer(I4B),
intent(in) :: icell2
1103 integer(I4B),
intent(in) :: shft
1104 integer(I4B),
intent(in) :: strt
1105 integer(I4B),
intent(in) :: stp
1106 integer(I4B),
intent(in) :: cntr
1111 do j = strt, stp, cntr
1112 this%uzthst(j, icell) = this2%uzthst(j + shft, icell2)
1113 this%uzdpst(j, icell) = this2%uzdpst(j + shft, icell2)
1114 this%uzflst(j, icell) = this2%uzflst(j + shft, icell2)
1115 this%uzspst(j, icell) = this2%uzspst(j + shft, icell2)
1117 this%nwavst(icell) = this2%nwavst(icell2)
1125 subroutine uzflow(this, thick, thickold, delt, ietflag, icell, ierr)
1128 real(DP),
intent(inout) :: thickold
1129 real(DP),
intent(inout) :: thick
1130 real(DP),
intent(in) :: delt
1131 integer(I4B),
intent(in) :: ietflag
1132 integer(I4B),
intent(in) :: icell
1133 integer(I4B),
intent(inout) :: ierr
1135 real(DP) :: ffcheck, time, feps1, feps2
1136 real(DP) :: thetadif, thetab, fluxb, oldsflx
1137 integer(I4B) :: itrailflg, itester
1140 this%totflux(icell) =
dzero
1142 oldsflx = this%uzflst(this%nwavst(icell), icell)
1146 if ((thick - thickold) > feps1)
then
1147 thetadif = abs(this%uzthst(1, icell) - this%thtr(icell))
1148 if (thetadif >
dem6)
then
1149 call this%wave_shift(this, icell, icell, -1, &
1150 this%nwavst(icell) + 1, 2, -1)
1151 if (this%uzdpst(2, icell) <
dem30) &
1152 this%uzdpst(2, icell) = (this%ntrail(icell) +
dtwo) *
dem6
1153 if (this%uzthst(2, icell) > this%thtr(icell))
then
1154 this%uzspst(2, icell) = this%uzflst(2, icell) / &
1155 (this%uzthst(2, icell) - this%thtr(icell))
1157 this%uzspst(2, icell) =
dzero
1159 this%uzthst(1, icell) = this%thtr(icell)
1160 this%uzflst(1, icell) =
dzero
1161 this%uzspst(1, icell) =
dzero
1162 this%uzdpst(1, icell) = thick
1163 this%nwavst(icell) = this%nwavst(icell) + 1
1164 if (this%nwavst(icell) >= this%nwav(icell))
then
1170 this%uzdpst(1, icell) = thick
1173 thetab = this%uzthst(1, icell)
1174 fluxb = this%uzflst(1, icell)
1175 this%totflux(icell) =
dzero
1177 ffcheck = (this%surflux(icell) - this%uzflst(this%nwavst(icell), icell))
1180 if (ffcheck > feps2 .OR. ffcheck < -feps2)
then
1181 this%nwavst(icell) = this%nwavst(icell) + 1
1182 if (this%nwavst(icell) >= this%nwav(icell))
then
1188 else if (this%nwavst(icell) == 1)
then
1191 if (this%nwavst(icell) > 1)
then
1192 if (ffcheck < -feps2)
then
1193 call this%trailwav(icell, ierr)
1194 if (ierr > 0)
return
1197 call this%leadwav(time, itester, itrailflg, thetab, fluxb, ffcheck, &
1200 if (itester == 1)
then
1201 this%totflux(icell) = this%totflux(icell) + &
1202 (delt - time) * this%uzflst(1, icell)
1208 if (ietflag > 0)
call this%uzet(icell, delt, ietflag, ierr)
1209 if (ierr > 0)
return
1219 real(DP),
intent(out) :: feps1
1220 real(DP),
intent(out) :: feps2
1230 factor1 =
done / 86400.d0
1231 else if (
itmuni == 2)
then
1232 factor1 =
done / 1440.d0
1233 else if (
itmuni == 3)
then
1234 factor1 =
done / 24.0d0
1235 else if (
itmuni == 5)
then
1238 factor2 =
done / 0.3048
1239 feps1 = feps1 * factor1 * factor2
1240 feps2 = feps2 * factor1 * factor2
1251 integer(I4B),
intent(in) :: icell
1252 integer(I4B),
intent(inout) :: ierr
1254 real(DP) :: smoist, smoistinc, ftrail, eps_m1
1255 real(DP) :: thtsrinv
1256 real(DP) :: flux1, flux2, theta1, theta2
1258 integer(I4B) :: j, jj, jk, nwavstm1
1261 eps_m1 = dble(this%eps(icell)) -
done
1262 thtsrinv =
done / (this%thts(icell) - this%thtr(icell))
1263 nwavstm1 = this%nwavst(icell) - 1
1266 smoist = (((this%surflux(icell) / this%vks(icell))** &
1267 (
done / this%eps(icell))) * &
1268 (this%thts(icell) - this%thtr(icell))) + this%thtr(icell)
1269 if (this%uzthst(nwavstm1, icell) - smoist >
dem9)
then
1271 do jk = 1, this%ntrail(icell)
1272 fnuminc = fnuminc + float(jk)
1274 smoistinc = (this%uzthst(nwavstm1, icell) - smoist) / (fnuminc -
done)
1275 jj = this%ntrail(icell)
1276 ftrail = dble(this%ntrail(icell)) +
done
1277 do j = this%nwavst(icell), this%nwavst(icell) + this%ntrail(icell) - 1
1278 if (j > this%nwav(icell))
then
1283 if (j > this%nwavst(icell))
then
1284 this%uzthst(j, icell) = this%uzthst(j - 1, icell) &
1285 - ((ftrail - float(jj)) * smoistinc)
1287 this%uzthst(j, icell) = this%uzthst(j - 1, icell) -
dem9
1290 if (this%uzthst(j, icell) <= this%thtr(icell) +
dem9) &
1291 this%uzthst(j, icell) = this%thtr(icell) +
dem9
1292 this%uzflst(j, icell) = &
1293 this%vks(icell) * (((this%uzthst(j, icell) - this%thtr(icell)) * &
1294 thtsrinv)**this%eps(icell))
1295 theta2 = this%uzthst(j - 1, icell)
1296 flux2 = this%uzflst(j - 1, icell)
1297 flux1 = this%uzflst(j, icell)
1298 theta1 = this%uzthst(j, icell)
1299 this%uzspst(j, icell) =
leadspeed(theta1, theta2, flux1, flux2, &
1300 this%thts(icell), this%thtr(icell), &
1301 this%eps(icell), this%vks(icell))
1302 this%uzdpst(j, icell) =
dzero
1303 if (j == this%nwavst(icell))
then
1304 this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1305 (this%ntrail(icell) + 1) *
dem9
1307 this%uzdpst(j, icell) = this%uzdpst(j - 1, icell) -
dem9
1310 this%nwavst(icell) = this%nwavst(icell) + this%ntrail(icell) - 1
1311 if (this%nwavst(icell) >= this%nwav(icell))
then
1317 this%uzdpst(this%nwavst, icell) =
dzero
1318 this%uzflst(this%nwavst, icell) = &
1319 this%vks(icell) * (((this%uzthst(this%nwavst, icell) - &
1320 this%thtr(icell)) * thtsrinv)**this%eps(icell))
1321 this%uzthst(this%nwavst, icell) = smoist
1322 theta2 = this%uzthst(this%nwavst(icell) - 1, icell)
1323 flux2 = this%uzflst(this%nwavst(icell) - 1, icell)
1324 flux1 = this%uzflst(this%nwavst(icell), icell)
1325 theta1 = this%uzthst(this%nwavst(icell), icell)
1326 this%uzspst(this%nwavst(icell), icell) = &
1327 leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1328 this%thtr(icell), this%eps(icell), this%vks(icell))
1337 subroutine leadwav(this, time, itester, itrailflg, thetab, fluxb, &
1338 ffcheck, feps2, delt, icell)
1341 real(DP),
intent(inout) :: thetab
1342 real(DP),
intent(inout) :: fluxb
1343 real(DP),
intent(in) :: feps2
1344 real(DP),
intent(inout) :: time
1345 integer(I4B),
intent(inout) :: itester
1346 integer(I4B),
intent(inout) :: itrailflg
1347 real(DP),
intent(inout) :: ffcheck
1348 real(DP),
intent(in) :: delt
1349 integer(I4B),
intent(in) :: icell
1351 real(DP) :: bottomtime, shortest, fcheck
1352 real(DP) :: eps_m1, timenew, bottom, timedt
1353 real(DP) :: thtsrinv, diff, fluxhld2
1354 real(DP) :: flux1, flux2, theta1, theta2, ftest
1355 real(DP),
allocatable,
dimension(:) :: checktime
1356 integer(I4B) :: iflx, iremove, j, l
1357 integer(I4B) :: nwavp1, jshort
1358 integer(I4B),
allocatable,
dimension(:) :: more
1360 allocate (checktime(this%nwavst(icell)))
1361 allocate (more(this%nwavst(icell)))
1363 eps_m1 = dble(this%eps(icell)) -
done
1364 thtsrinv =
done / (this%thts(icell) - this%thtr(icell))
1367 if (itrailflg == 0)
then
1368 if (ffcheck > feps2)
then
1369 this%uzflst(this%nwavst(icell), icell) = this%surflux(icell)
1370 if (this%uzflst(this%nwavst(icell), icell) <
dem30) &
1371 this%uzflst(this%nwavst(icell), icell) =
dzero
1372 this%uzthst(this%nwavst(icell), icell) = &
1373 (((this%uzflst(this%nwavst(icell), icell) / this%vks(icell))** &
1374 (
done / this%eps(icell))) * (this%thts(icell) - this%thtr(icell))) &
1376 theta2 = this%uzthst(this%nwavst(icell), icell)
1377 flux2 = this%uzflst(this%nwavst(icell), icell)
1378 flux1 = this%uzflst(this%nwavst(icell) - 1, icell)
1379 theta1 = this%uzthst(this%nwavst(icell) - 1, icell)
1380 this%uzspst(this%nwavst(icell), icell) = &
1381 leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1382 this%thtr(icell), this%eps(icell), this%vks(icell))
1383 this%uzdpst(this%nwavst(icell), icell) =
dzero
1391 fluxhld2 = this%uzflst(1, icell)
1392 if (this%nwavst(icell) == 0) itester = 1
1393 if (itester /= 1)
then
1394 do while (diff >
dem6)
1395 nwavp1 = this%nwavst(icell) + 1
1396 timedt = delt - time
1397 do j = 1, this%nwavst(icell)
1398 checktime(j) =
dep20
1402 if (this%nwavst(icell) > 2)
then
1406 nwavp1 = this%nwavst(icell) + 1
1407 do while (j < nwavp1)
1408 ftest = this%uzspst(j - 1, icell) - this%uzspst(j, icell)
1409 if (abs(ftest) >
dem30)
then
1410 checktime(j) = (this%uzdpst(j, icell) - &
1411 this%uzdpst(j - 1, icell)) / ftest
1412 if (checktime(j) <
dem30) checktime(j) =
dep20
1420 if (this%nwavst(icell) > 1)
then
1421 if (this%uzspst(2, icell) >
dzero)
then
1422 bottom = this%uzspst(2, icell)
1424 bottomtime = (this%uzdpst(1, icell) - this%uzdpst(2, icell)) / bottom
1431 do j = this%nwavst(icell), 3, -1
1432 if (shortest - checktime(j) > -
dem9)
then
1435 shortest = checktime(j)
1438 do j = 3, this%nwavst(icell)
1439 if (shortest - checktime(j) <
dem9)
then
1440 if (j /= jshort) more(j) = 0
1447 fcheck = (time + shortest) - delt
1448 if (shortest <
dem7) fcheck = -
done
1449 if (bottomtime < shortest .AND. time + bottomtime < delt)
then
1451 do while (j < nwavp1)
1454 this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1455 this%uzspst(j, icell) * bottomtime
1458 fluxb = this%uzflst(2, icell)
1459 thetab = this%uzthst(2, icell)
1461 call this%wave_shift(this, icell, icell, 1, 1, &
1462 this%nwavst(icell) - 1, 1)
1464 timenew = time + bottomtime
1465 this%uzspst(1, icell) =
dzero
1468 else if (fcheck <
dzero .AND. this%nwavst(icell) > 2)
then
1470 do while (j < nwavp1)
1471 this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1472 this%uzspst(j, icell) * shortest
1479 do while (j < this%nwavst(icell) + 1)
1480 if (more(j) == 1)
then
1482 theta2 = this%uzthst(j, icell)
1483 flux2 = this%uzflst(j, icell)
1488 flux1 = this%uzflst(j - 2, icell)
1489 theta1 = this%uzthst(j - 2, icell)
1491 this%uzspst(j, icell) =
leadspeed(theta1, theta2, flux1, flux2, &
1494 this%eps(icell), this%vks(icell))
1497 call this%wave_shift(this, icell, icell, 1, l - 1, &
1498 this%nwavst(icell) - 1, 1)
1499 l = this%nwavst(icell) + 1
1500 iremove = iremove + 1
1504 timenew = timenew + shortest
1509 do while (j < nwavp1)
1510 this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1511 this%uzspst(j, icell) * timedt
1516 this%totflux(icell) = this%totflux(icell) + fluxhld2 * (timenew - time)
1518 fluxhld2 = this%uzflst(1, icell)
1523 this%nwavst(icell) = this%nwavst(icell) - iremove
1526 if (this%nwavst(icell) == 1)
then
1532 deallocate (checktime)
1541 function leadspeed(theta1, theta2, flux1, flux2, thts, thtr, eps, vks)
1545 real(dp),
intent(in) :: theta1
1546 real(dp),
intent(in) :: theta2
1547 real(dp),
intent(in) :: flux1
1548 real(dp),
intent(inout) :: flux2
1549 real(dp),
intent(in) :: thts
1550 real(dp),
intent(in) :: thtr
1551 real(dp),
intent(in) :: eps
1552 real(dp),
intent(in) :: vks
1554 real(dp) :: comp1, comp2, thsrinv, epsfksths
1555 real(dp) :: eps_m1, fhold, comp3
1558 thsrinv =
done / (thts - thtr)
1559 epsfksths = eps * vks * thsrinv
1560 comp1 = theta2 - theta1
1561 comp2 = abs(flux2 - flux1)
1562 comp3 = theta1 - thtr
1564 if (abs(comp1) <
dem30)
then
1565 if (comp3 >
dem30) fhold = (comp3 * thsrinv)**eps
1569 leadspeed = (flux2 - flux1) / (theta2 - theta1)
1584 integer(I4B),
intent(in) :: icell
1585 real(dp),
intent(inout) :: d1
1588 integer(I4B) :: j, k, nwavm1, jj
1591 j = this%nwavst(icell) + 1
1592 k = this%nwavst(icell)
1594 if (d1 > this%uzdpst(1, icell)) d1 = this%uzdpst(1, icell)
1598 if (this%uzdpst(k, icell) - d1 < -
dem30) j = k
1601 if (j > this%nwavst(icell))
then
1602 fm = fm + (this%uzthst(this%nwavst(icell), icell) - this%thtr(icell)) * d1
1603 elseif (this%nwavst(icell) > 1)
then
1605 fm = fm + (this%uzthst(j - 1, icell) - this%thtr(icell)) &
1606 * (d1 - this%uzdpst(j, icell))
1609 fm = fm + (this%uzthst(jj, icell) - this%thtr(icell)) &
1610 * (this%uzdpst(jj, icell) &
1611 - this%uzdpst(jj + 1, icell))
1613 fm = fm + (this%uzthst(this%nwavst(icell), icell) - this%thtr(icell)) &
1614 * (this%uzdpst(this%nwavst(icell), icell))
1616 fm = fm + (this%uzthst(1, icell) - this%thtr(icell)) * d1
1629 integer(I4B),
intent(in) :: icell
1630 integer(I4B),
intent(in) :: itest
1631 integer(I4B),
intent(in) :: iss
1632 real(DP),
intent(in) :: delt
1634 real(DP) :: bot, depthsave, top
1635 real(DP) :: thick, thtsrinv
1636 integer(I4B) :: nwavhld, k, j
1638 bot = this%watab(icell)
1639 top = this%celtop(icell)
1641 nwavhld = this%nwavst(icell)
1642 if (itest == 1)
then
1643 this%uzflst(1, icell) =
dzero
1644 this%uzthst(1, icell) = this%thtr(icell)
1648 if (this%thts(icell) - this%thtr(icell) <
dem7)
then
1651 thtsrinv =
done / (this%thts(icell) - this%thtr(icell))
1653 this%totflux(icell) = this%surflux(icell) * delt
1654 this%watabold(icell) = this%watab(icell)
1655 this%uzthst(1, icell) = this%thti(icell)
1656 this%uzflst(1, icell) = &
1657 this%vks(icell) * (((this%uzthst(1, icell) - this%thtr(icell)) &
1658 * thtsrinv)**this%eps(icell))
1659 this%uzdpst(1, icell) = thick
1660 this%uzspst(1, icell) = thick
1661 this%nwavst(icell) = 1
1665 if (this%watab(icell) - this%watabold(icell) >
dem30)
then
1666 depthsave = this%uzdpst(1, icell)
1668 k = this%nwavst(icell)
1670 if (this%uzdpst(k, icell) - thick < -
dem30) j = k
1673 this%uzdpst(1, icell) = thick
1675 this%uzspst(1, icell) =
dzero
1676 this%nwavst(icell) = this%nwavst(icell) - j + 2
1677 this%uzthst(1, icell) = this%uzthst(j - 1, icell)
1678 this%uzflst(1, icell) = this%uzflst(j - 1, icell)
1679 if (j > 2)
call this%wave_shift(this, icell, icell, j - 2, 2, &
1680 nwavhld - (j - 2), 1)
1681 elseif (j == 0)
then
1682 this%uzspst(1, icell) =
dzero
1683 this%uzthst(1, icell) = this%uzthst(this%nwavst(icell), icell)
1684 this%uzflst(1, icell) = this%uzflst(this%nwavst(icell), icell)
1685 this%nwavst(icell) = 1
1690 if (thick <=
dzero)
then
1691 this%uzspst(1, icell) =
dzero
1692 this%nwavst(icell) = 1
1693 this%uzthst(1, icell) = this%thtr(icell)
1694 this%uzflst(1, icell) =
dzero
1696 this%watabold(icell) = this%watab(icell)
1705 subroutine uzet(this, icell, delt, ietflag, ierr)
1708 integer(I4B),
intent(in) :: icell
1709 real(DP),
intent(in) :: delt
1710 integer(I4B),
intent(in) :: ietflag
1711 integer(I4B),
intent(inout) :: ierr
1715 real(DP) :: thetaout
1718 real(DP) :: thtsrinv
1719 real(DP) :: epsfksthts
1734 integer(I4B) :: jhold
1738 integer(I4B) :: numadd
1741 integer(I4B) :: itest
1744 this%etact(icell) =
dzero
1745 if (this%extdpuz(icell) <
dem7)
return
1746 petsub = this%rootact(icell) * this%pet(icell) * &
1747 this%extdpuz(icell) / this%extdp(icell)
1748 thetaout = delt * petsub / this%extdp(icell)
1749 if (ietflag == 1) thetaout = delt * this%pet(icell) / this%extdp(icell)
1750 if (thetaout <
dem10)
return
1751 depth = this%uzdpst(1, icell)
1752 st = this%unsat_stor(icell, depth)
1753 if (st <
dem4)
return
1756 nwv = this%nwavst(icell)
1758 call uzfktemp%init(1, nwv)
1761 call uzfktemp%wave_shift(this, 1, icell, 0, 1, nwv, 1)
1763 this%etact(icell) =
dzero
1764 if (this%thts(icell) - this%thtr(icell) <
dem7)
then
1765 thtsrinv = 1.0 /
dem7
1767 thtsrinv =
done / (this%thts(icell) - this%thtr(icell))
1769 epsfksthts = this%eps(icell) * this%vks(icell) * thtsrinv
1770 this%etact(icell) =
dzero
1772 extwc1 = this%extwc(icell) - this%thtr(icell)
1779 do while (itest == 0)
1781 if (k > 1 .AND. abs(fmp - petsub) >
dem5 * petsub)
then
1782 factor = factor / (fm / petsub)
1786 if (this%nwavst(icell) == 1 .AND. &
1787 this%uzdpst(1, icell) <= this%extdpuz(icell))
then
1788 if (ietflag == 2)
then
1789 tho = this%uzthst(1, icell)
1790 fktho = this%uzflst(1, icell)
1791 hcap = this%caph(icell, tho)
1792 thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1794 if ((this%uzthst(1, icell) - thetaout) > this%thtr(icell) + extwc1)
then
1795 this%uzthst(1, icell) = this%uzthst(1, icell) - thetaout
1796 this%uzflst(1, icell) = &
1797 this%vks(icell) * (((this%uzthst(1, icell) - &
1798 this%thtr(icell)) * thtsrinv)**this%eps(icell))
1799 else if (this%uzthst(1, icell) > this%thtr(icell) + extwc1)
then
1800 this%uzthst(1, icell) = this%thtr(icell) + extwc1
1801 this%uzflst(1, icell) = &
1802 this%vks(icell) * (((this%uzthst(1, icell) - &
1803 this%thtr(icell)) * thtsrinv)**this%eps(icell))
1807 else if (this%nwavst(icell) > 1 .AND. &
1808 this%uzdpst(this%nwavst(icell), icell) > this%extdpuz(icell))
then
1809 if (ietflag == 2)
then
1810 tho = this%uzthst(this%nwavst(icell), icell)
1811 fktho = this%uzflst(this%nwavst(icell), icell)
1812 hcap = this%caph(icell, tho)
1813 thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1815 if (this%uzthst(this%nwavst(icell), icell) - thetaout > &
1816 this%thtr(icell) + extwc1)
then
1817 this%uzthst(this%nwavst(icell) + 1, icell) = &
1818 this%uzthst(this%nwavst(icell), icell) - thetaout
1820 else if (this%uzthst(this%nwavst(icell), icell) > &
1821 this%thtr(icell) + extwc1)
then
1822 this%uzthst(this%nwavst(icell) + 1, icell) = this%thtr(icell) + extwc1
1825 if (numadd == 1)
then
1826 this%uzflst(this%nwavst(icell) + 1, icell) = &
1828 (((this%uzthst(this%nwavst(icell) + 1, icell) - &
1829 this%thtr(icell)) * thtsrinv)**this%eps(icell))
1830 theta2 = this%uzthst(this%nwavst(icell) + 1, icell)
1831 flux2 = this%uzflst(this%nwavst(icell) + 1, icell)
1832 flux1 = this%uzflst(this%nwavst(icell), icell)
1833 theta1 = this%uzthst(this%nwavst(icell), icell)
1834 this%uzspst(this%nwavst(icell) + 1, icell) = &
1835 leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1836 this%thtr(icell), this%eps(icell), this%vks(icell))
1837 this%uzdpst(this%nwavst(icell) + 1, icell) = this%extdpuz(icell)
1838 this%nwavst(icell) = this%nwavst(icell) + 1
1839 if (this%nwavst(icell) > this%nwav(icell))
then
1850 else if (this%nwavst(icell) == 1)
then
1851 if (ietflag == 2)
then
1852 tho = this%uzthst(1, icell)
1853 fktho = this%uzflst(1, icell)
1854 hcap = this%caph(icell, tho)
1855 thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1857 if ((this%uzthst(1, icell) - thetaout) > this%thtr(icell) + extwc1)
then
1858 if (thetaout >
dem30)
then
1859 this%uzthst(2, icell) = this%uzthst(1, icell) - thetaout
1860 this%uzflst(2, icell) = &
1861 this%vks(icell) * (((this%uzthst(2, icell) - this%thtr(icell)) * &
1862 thtsrinv)**this%eps(icell))
1863 this%uzdpst(2, icell) = this%extdpuz(icell)
1864 theta2 = this%uzthst(2, icell)
1865 flux2 = this%uzflst(2, icell)
1866 flux1 = this%uzflst(1, icell)
1867 theta1 = this%uzthst(1, icell)
1868 this%uzspst(2, icell) = &
1869 leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1870 this%thtr(icell), this%eps(icell), this%vks(icell))
1871 this%nwavst(icell) = this%nwavst(icell) + 1
1872 if (this%nwavst(icell) > this%nwav(icell))
then
1879 else if (this%uzthst(1, icell) > this%thtr(icell) + extwc1)
then
1880 if (thetaout >
dem30)
then
1881 this%uzthst(2, icell) = this%thtr(icell) + extwc1
1882 this%uzflst(2, icell) = &
1883 this%vks(icell) * (((this%uzthst(2, icell) - &
1884 this%thtr(icell)) * thtsrinv)**this%eps(icell))
1885 this%uzdpst(2, icell) = this%extdpuz(icell)
1886 theta2 = this%uzthst(2, icell)
1887 flux2 = this%uzflst(2, icell)
1888 flux1 = this%uzflst(1, icell)
1889 theta1 = this%uzthst(1, icell)
1890 this%uzspst(2, icell) = &
1891 leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1892 this%thtr(icell), this%eps(icell), this%vks(icell))
1893 this%nwavst(icell) = this%nwavst(icell) + 1
1894 if (this%nwavst(icell) > this%nwav(icell))
then
1905 if (this%uzdpst(1, icell) - this%extdpuz(icell) >
dem7)
then
1911 diff = this%uzdpst(j, icell) - this%extdpuz(icell)
1912 if (diff >
dzero)
then
1919 if (this%uzthst(j, icell) > this%thtr(icell) + extwc1)
then
1922 if (abs(diff) >
dem5)
then
1923 call this%wave_shift(this, icell, icell, -1, &
1924 this%nwavst(icell) + 1, j, -1)
1925 this%uzdpst(j, icell) = this%extdpuz(icell)
1926 this%nwavst(icell) = this%nwavst(icell) + 1
1927 if (this%nwavst(icell) > this%nwav(icell))
then
1936 jhold = this%nwavst(icell)
1938 do while (i < this%nwavst(icell))
1939 if (this%uzthst(i, icell) > this%thtr(icell) + extwc1)
then
1941 i = this%nwavst(icell) + 1
1953 do while (kk <= this%nwavst(icell))
1954 if (ietflag == 2)
then
1955 tho = this%uzthst(kk, icell)
1956 fktho = this%uzflst(kk, icell)
1957 hcap = this%caph(icell, tho)
1958 thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1960 if (this%uzthst(kk, icell) > this%thtr(icell) + extwc1)
then
1961 if (this%uzthst(kk, icell) - thetaout > &
1962 this%thtr(icell) + extwc1)
then
1963 this%uzthst(kk, icell) = this%uzthst(kk, icell) - thetaout
1964 else if (this%uzthst(kk, icell) > this%thtr(icell) + extwc1)
then
1965 this%uzthst(kk, icell) = this%thtr(icell) + extwc1
1968 this%uzflst(kk, icell) = &
1970 (((this%uzthst(kk, icell) - &
1971 this%thtr(icell)) * thtsrinv)**this%eps(icell))
1975 this%vks(icell) * ((this%uzthst(kk - 1, icell) - &
1976 this%thtr(icell)) * thtsrinv)**this%eps(icell)
1978 this%vks(icell) * ((this%uzthst(kk, icell) - &
1979 this%thtr(icell)) * thtsrinv)**this%eps(icell)
1980 this%uzflst(kk, icell) = flux2
1981 theta2 = this%uzthst(kk, icell)
1982 theta1 = this%uzthst(kk - 1, icell)
1983 this%uzspst(kk, icell) =
leadspeed(theta1, theta2, flux1, flux2, &
1986 this%eps(icell), this%vks(icell))
1995 do while (kj <= this%nwavst(icell) - 1)
1996 if (abs(this%uzthst(kj, icell) - this%uzthst(kj + 1, icell)) <
dem6)
then
1997 call this%wave_shift(this, icell, icell, 1, kj + 1, &
1998 this%nwavst(icell) - 1, 1)
2000 this%nwavst(icell) = this%nwavst(icell) - 1
2004 depth = this%uzdpst(1, icell)
2005 fm = this%unsat_stor(icell, depth)
2006 this%etact(icell) = st - fm
2007 fm = this%etact(icell) / delt
2008 if (this%etact(icell) <
dzero)
then
2009 call this%wave_shift(uzfktemp, icell, 1, 0, 1, nwv, 1)
2010 this%nwavst(icell) = nwv
2011 this%etact(icell) =
dzero
2012 elseif (petsub - fm < -
dem15 .AND. ietflag == 2)
then
2015 call this%wave_shift(uzfktemp, icell, 1, 0, 1, nwv, 1)
2016 this%nwavst(icell) = nwv
2017 this%etact(icell) =
dzero
2026 elseif (ietflag < 2)
then
2034 call uzfktemp%dealloc()
2045 integer(I4B),
intent(in) :: icell
2046 real(dp),
intent(in) :: tho
2048 real(dp) ::
caph, lambda, star
2051 star = (tho - this%thtr(icell)) / (this%thts(icell) - this%thtr(icell))
2054 if (star >
dem15)
then
2055 if (tho - this%thts(icell) <
dem15)
then
2056 caph = this%ha(icell) * star**(-
done / lambda)
2072 integer(I4B),
intent(in) :: icell
2073 real(dp),
intent(in) :: factor, fktho, h
2075 rate_et_z = factor * fktho * (h - this%hroot(icell))
2090 integer(I4B),
intent(in) :: icell
2091 real(dp),
intent(in) :: depth
2093 real(dp) :: theta_at_depth
2100 if (this%watab(icell) < this%celtop(icell))
then
2101 if (this%celtop(icell) - depth > this%watab(icell))
then
2104 f1 = this%unsat_stor(icell, d1)
2105 f2 = this%unsat_stor(icell, d2)
2106 theta_at_depth = this%thtr(icell) + (f2 - f1) / (d2 - d1)
2108 theta_at_depth = this%thts(icell)
2111 theta_at_depth = this%thts(icell)
2123 integer(I4B),
intent(in) :: icell
2125 real(dp) :: watercontent
2135 hgwf = this%watab(icell)
2136 top = this%celtop(icell)
2137 bot = this%celbot(icell)
2138 thk = top - max(bot, hgwf)
2139 if (thk >
dzero)
then
2140 theta_r = this%thtr(icell)
2142 fm = this%unsat_stor(icell, d)
2143 watercontent = fm / thk
2144 watercontent = watercontent + theta_r
2146 watercontent =
dzero
This module contains simulation constants.
real(dp), parameter dem12
real constant 1e-12
real(dp), parameter dem20
real constant 1e-20
real(dp), parameter dem10
real constant 1e-10
real(dp), parameter dem7
real constant 1e-7
real(dp), parameter dep20
real constant 1e20
real(dp), parameter dem14
real constant 1e-14
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dem3
real constant 1e-3
real(dp), parameter dem4
real constant 1e-4
real(dp), parameter dem30
real constant 1e-30
real(dp), parameter dem6
real constant 1e-6
real(dp), parameter dzero
real constant zero
real(dp), parameter dem5
real constant 1e-5
real(dp), parameter dem9
real constant 1e-9
real(dp), parameter dem15
real constant 1e-15
real(dp), parameter dtwo
real constant 2
real(dp), parameter dthree
real constant 3
real(dp), parameter done
real constant 1
This module defines variable data types.
subroutine slinear(x, range, dydx, y)
@ brief sLinear
subroutine scubiclinear(x, range, dydx, y)
@ brief sCubicLinear
subroutine sscurve(x, range, dydx, y)
@ brief SCurve
subroutine scubic(x, range, dydx, y)
@ brief sCubic
integer(i4b), pointer, public itmuni
flag indicating time units
integer(i4b), pointer, public kper
current stress period number
real(dp), pointer, public delt
length of the current time step
subroutine update_wav(this, icell, delt, iss, itest)
Update to new state of uz at end of time step.
subroutine factors(feps1, feps2)
Calculate unit specific tolerances.
subroutine uzflow(this, thick, thickold, delt, ietflag, icell, ierr)
Method of Characteristics solution for kinematic wave equation.
subroutine advance(this, icell)
Set variables to advance to new time step. nothing yet.
subroutine setbelowpet(this, icell, jbelow)
Subtract aet from pet to calculate residual et for deeper cells.
subroutine gwseep(this, icell, deriv, scale, hgwf, trhs, thcof, seep)
Calculate groudwater discharge to land surface.
subroutine sethead(this, icell, hgwf)
Set initial head for uzf object.
subroutine leadwav(this, time, itester, itrailflg, thetab, fluxb, ffcheck, feps2, delt, icell)
Create a lead wave and route over time step.
real(dp) function caph(this, icell, tho)
Calculate capillary pressure head from B-C equation.
subroutine setdatauzfarea(this, icell, areamult)
Set uzfarea using cellarea and areamult.
subroutine setdataetha(this, icell, jbelow, ha, hroot, rootact)
Set variables for head-based unsaturated flow.
subroutine addrech(this, icell, jbelow, hgwf, trhs, thcof, deriv, delt)
Add recharge or infiltration to cells.
subroutine setwaves(this, icell)
Reset waves to default values at start of simulation.
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.
subroutine setdata(this, icell, area, top, bot, surfdep, vks, thtr, thts, thti, eps, ntrail, landflag, ivertcon)
Set uzf object material properties.
real(dp) function leadspeed(theta1, theta2, flux1, flux2, thts, thtr, eps, vks)
Calculates waves speed from dflux/dtheta.
real(dp) function get_water_content_at_depth(this, icell, depth)
Determine the water content at a specific depth.
real(dp) function etfunc_nlin(s, x, c, det, trhs, thcof, hgwf)
Square-wave ET function with smoothing at extinction depth.
subroutine wave_shift(this, this2, icell, icell2, shft, strt, stp, cntr)
Copy waves or shift waves in arrays.
subroutine rejfinf(this, icell, deriv, hgwf, trhs, thcof, finfact)
Reject applied infiltration due to low vks.
real(dp) function rate_et_z(this, icell, factor, fktho, h)
Calculate capillary pressure-based uz et.
subroutine simgwet(this, igwetflag, icell, hgwf, trhs, thcof, det)
Calculate gwf et using residual uzf pet.
real(dp) function get_wcnew(this, icell)
Calculate and return the cell-based water content value.
subroutine uz_rise(this, icell, totfluxtot)
Calculate recharge due to a rise in the gwf head.
subroutine uzet(this, icell, delt, ietflag, ierr)
Remove water from uz due to et.
subroutine setgwpet(this, icell)
Subtract aet from pet to calculate residual et for gw.
subroutine routewaves(this, totfluxtot, delt, ietflag, icell, ierr)
Prepare and route waves over time step.
subroutine setdataet(this, icell, jbelow, pet, extdp)
Set unsaturated ET-related variables.
real(dp) function etfunc_lin(s, x, c, det, trhs, thcof, hgwf, celtop, celbot)
Calculate gwf et using linear ET function from mf-2005.
subroutine dealloc(this)
Deallocate uzf object variables.
real(dp) function unsat_stor(this, icell, d1)
Sums up mobile water over depth interval.
subroutine trailwav(this, icell, ierr)
Create and set trail waves.
subroutine setdataetwc(this, icell, jbelow, extwc)
Set extinction water content.
subroutine setdatafinf(this, icell, finf)
Set infiltration.