MODFLOW 6  version 6.9.0.dev0
USGS Modular Hydrologic Model
Data Types | Functions/Subroutines
uzfcellgroupmodule Module Reference

Data Types

type  uzfwavestoretype
 Snapshot of one cell's wave train. More...
 
type  uzfcellgrouptype
 

Functions/Subroutines

subroutine init (this, ncells, nwav, memory_path)
 Allocate and set uzf object variables. More...
 
subroutine dealloc (this)
 Deallocate uzf object variables. More...
 
subroutine setdata (this, icell, area, top, bot, surfdep, vks, thtr, thts, thti, eps, ntrail, landflag, ivertcon)
 Set uzf object material properties. More...
 
subroutine sethead (this, icell, hgwf)
 Set initial head for uzf object. More...
 
subroutine setdatafinf (this, icell, finf)
 Set infiltration. More...
 
subroutine setdatauzfarea (this, icell, areamult)
 Set uzfarea using cellarea and areamult. More...
 
subroutine setdataet (this, icell, jbelow, pet, extdp)
 Set unsaturated ET-related variables. More...
 
subroutine setgwpet (this, icell)
 Subtract aet from pet to calculate residual et for gw. More...
 
subroutine setbelowpet (this, icell, jbelow)
 Subtract aet from pet to calculate residual et for deeper cells. More...
 
subroutine setdataetwc (this, icell, jbelow, extwc)
 Set extinction water content. More...
 
subroutine setdataetha (this, icell, jbelow, ha, hroot, rootact)
 Set variables for head-based unsaturated flow. More...
 
subroutine advance (this, icell)
 Set variables to advance to new time step. nothing yet. More...
 
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. More...
 
subroutine addrech (this, icell, jbelow, hgwf, trhs, thcof, deriv, delt)
 Add recharge or infiltration to cells. More...
 
subroutine rejfinf (this, icell, deriv, hgwf, trhs, thcof, finfact)
 Reject applied infiltration due to low vks. More...
 
subroutine gwseep (this, icell, deriv, scale, hgwf, trhs, thcof, seep)
 Calculate groundwater discharge to land surface. More...
 
subroutine simgwet (this, igwetflag, icell, hgwf, trhs, thcof, det)
 Calculate gwf et using residual uzf pet. More...
 
subroutine uz_rise (this, icell, totfluxtot)
 Calculate recharge due to a rise in the gwf head. More...
 
subroutine setwaves (this, icell)
 Reset waves to default values at start of simulation. More...
 
subroutine routewaves (this, totfluxtot, delt, ietflag, icell, ierr)
 Prepare and route waves over time step. More...
 
subroutine shift_waves (this, icell, shft, strt, stp, cntr)
 Move waves within a cell by shft positions. More...
 
subroutine store_waves (this, icell, store)
 Save the wave train of one cell into a snapshot. More...
 
subroutine load_waves (this, icell, store)
 Restore the wave train of one cell from a snapshot. More...
 
subroutine uzflow (this, thick, thickold, delt, ietflag, icell, ierr)
 Method of Characteristics solution for kinematic wave equation. More...
 
subroutine factors (feps1, feps2)
 Calculate unit specific tolerances. More...
 
subroutine trailwav (this, icell, ierr)
 Create and set trail waves. More...
 
subroutine leadwav (this, time, itester, itrailflg, thetab, fluxb, ffcheck, feps2, delt, icell)
 Create a lead wave and route over time step. More...
 
real(dp) function leadspeed (theta1, theta2, flux1, flux2, thts, thtr, eps, vks)
 Calculates waves speed from dflux/dtheta. More...
 
real(dp) function unsat_stor (this, icell, d1)
 Sums up mobile water over depth interval. More...
 
subroutine update_wav (this, icell, delt, iss, itest)
 Update to new state of uz at end of time step. More...
 
subroutine uzet (this, icell, delt, ietflag, ierr)
 Remove water from uz due to et. More...
 
real(dp) function caph (this, icell, tho)
 Calculate capillary pressure head from B-C equation. More...
 
real(dp) function rate_et_z (this, icell, factor, fktho, h)
 Calculate capillary pressure-based uz et. More...
 
real(dp) function get_water_content_at_depth (this, icell, depth)
 Determine the water content at a specific depth. More...
 
real(dp) function get_wcnew (this, icell)
 Calculate and return the cell-based water content value. More...
 

Function/Subroutine Documentation

◆ addrech()

subroutine uzfcellgroupmodule::addrech ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
integer(i4b), intent(in)  jbelow,
real(dp), intent(in)  hgwf,
real(dp), intent(inout)  trhs,
real(dp), intent(inout)  thcof,
real(dp), intent(inout)  deriv,
real(dp), intent(in)  delt 
)

Definition at line 684 of file UzfCellGroup.f90.

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
Here is the call graph for this function:

◆ advance()

subroutine uzfcellgroupmodule::advance ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell 
)
private

Definition at line 532 of file UzfCellGroup.f90.

533  ! -- dummy
534  class(UzfCellGroupType) :: this
535  integer(I4B), intent(in) :: icell
536  !
537  ! -- set variables
538  this%surfseep(icell) = dzero

◆ caph()

real(dp) function uzfcellgroupmodule::caph ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  tho 
)
private

Definition at line 1897 of file UzfCellGroup.f90.

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

◆ dealloc()

subroutine uzfcellgroupmodule::dealloc ( class(uzfcellgrouptype)  this)

Definition at line 254 of file UzfCellGroup.f90.

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)
memorymanagermodule::mem_deallocate
Definition: MemoryManager.f90:129

◆ factors()

subroutine uzfcellgroupmodule::factors ( real(dp), intent(out)  feps1,
real(dp), intent(out)  feps2 
)
private

Definition at line 1087 of file UzfCellGroup.f90.

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
Here is the caller graph for this function:

◆ get_water_content_at_depth()

real(dp) function uzfcellgroupmodule::get_water_content_at_depth ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  depth 
)
private

Because UZF-calculated waves are internal to UZF objects, different water contents exists at different depths.

Parameters
[in]icelluzf cell containing depth
[in]depthdepth within the cell

Definition at line 1936 of file UzfCellGroup.f90.

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

◆ get_wcnew()

real(dp) function uzfcellgroupmodule::get_wcnew ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell 
)
private
Parameters
[in]icelluzf cell containing depth

Definition at line 1966 of file UzfCellGroup.f90.

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

◆ gwseep()

subroutine uzfcellgroupmodule::gwseep ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(inout)  deriv,
real(dp), intent(out)  scale,
real(dp), intent(in)  hgwf,
real(dp), intent(inout)  trhs,
real(dp), intent(inout)  thcof,
real(dp), intent(inout)  seep 
)
private

Definition at line 751 of file UzfCellGroup.f90.

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
Here is the call graph for this function:

◆ init()

subroutine uzfcellgroupmodule::init ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  ncells,
integer(i4b), intent(in)  nwav,
character(len=*), intent(in)  memory_path 
)
private

Definition at line 121 of file UzfCellGroup.f90.

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
memorymanagermodule::mem_allocate
Definition: MemoryManager.f90:55

◆ leadspeed()

real(dp) function uzfcellgroupmodule::leadspeed ( real(dp), intent(in)  theta1,
real(dp), intent(in)  theta2,
real(dp), intent(in)  flux1,
real(dp), intent(inout)  flux2,
real(dp), intent(in)  thts,
real(dp), intent(in)  thtr,
real(dp), intent(in)  eps,
real(dp), intent(in)  vks 
)
private

Definition at line 1396 of file UzfCellGroup.f90.

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
Here is the caller graph for this function:

◆ leadwav()

subroutine uzfcellgroupmodule::leadwav ( class(uzfcellgrouptype)  this,
real(dp), intent(inout)  time,
integer(i4b), intent(inout)  itester,
integer(i4b), intent(inout)  itrailflg,
real(dp), intent(inout)  thetab,
real(dp), intent(inout)  fluxb,
real(dp), intent(inout)  ffcheck,
real(dp), intent(in)  feps2,
real(dp), intent(in)  delt,
integer(i4b), intent(in)  icell 
)
private

Definition at line 1201 of file UzfCellGroup.f90.

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
Here is the call graph for this function:

◆ load_waves()

subroutine uzfcellgroupmodule::load_waves ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
type(uzfwavestoretype), intent(in)  store 
)
private

Definition at line 979 of file UzfCellGroup.f90.

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

◆ rate_et_z()

real(dp) function uzfcellgroupmodule::rate_et_z ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  factor,
real(dp), intent(in)  fktho,
real(dp), intent(in)  h 
)
private

Definition at line 1919 of file UzfCellGroup.f90.

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

◆ rejfinf()

subroutine uzfcellgroupmodule::rejfinf ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(inout)  deriv,
real(dp), intent(in)  hgwf,
real(dp), intent(inout)  trhs,
real(dp), intent(inout)  thcof,
real(dp), intent(inout)  finfact 
)
private

Definition at line 723 of file UzfCellGroup.f90.

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
Here is the call graph for this function:

◆ routewaves()

subroutine uzfcellgroupmodule::routewaves ( class(uzfcellgrouptype)  this,
real(dp), intent(inout)  totfluxtot,
real(dp), intent(in)  delt,
integer(i4b), intent(in)  ietflag,
integer(i4b), intent(in)  icell,
integer(i4b), intent(inout)  ierr 
)
private

Definition at line 897 of file UzfCellGroup.f90.

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

◆ setbelowpet()

subroutine uzfcellgroupmodule::setbelowpet ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
integer(i4b), intent(in)  jbelow 
)

Definition at line 472 of file UzfCellGroup.f90.

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
tdismodule::delt
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:32

◆ setdata()

subroutine uzfcellgroupmodule::setdata ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  area,
real(dp), intent(in)  top,
real(dp), intent(in)  bot,
real(dp), intent(in)  surfdep,
real(dp), intent(in)  vks,
real(dp), intent(in)  thtr,
real(dp), intent(in)  thts,
real(dp), intent(in)  thti,
real(dp), intent(in)  eps,
integer(i4b), intent(in)  ntrail,
integer(i4b), intent(in)  landflag,
integer(i4b), intent(in)  ivertcon 
)

Definition at line 316 of file UzfCellGroup.f90.

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

◆ setdataet()

subroutine uzfcellgroupmodule::setdataet ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
integer(i4b), intent(in)  jbelow,
real(dp), intent(in)  pet,
real(dp), intent(in)  extdp 
)
private

Definition at line 409 of file UzfCellGroup.f90.

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

◆ setdataetha()

subroutine uzfcellgroupmodule::setdataetha ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
integer(i4b), intent(in)  jbelow,
real(dp), intent(in)  ha,
real(dp), intent(in)  hroot,
real(dp), intent(in)  rootact 
)
private

Definition at line 510 of file UzfCellGroup.f90.

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

◆ setdataetwc()

subroutine uzfcellgroupmodule::setdataetwc ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
integer(i4b), intent(in)  jbelow,
real(dp), intent(in)  extwc 
)

Definition at line 496 of file UzfCellGroup.f90.

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

◆ setdatafinf()

subroutine uzfcellgroupmodule::setdatafinf ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  finf 
)
private

Definition at line 377 of file UzfCellGroup.f90.

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

◆ setdatauzfarea()

subroutine uzfcellgroupmodule::setdatauzfarea ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  areamult 
)
private

Definition at line 397 of file UzfCellGroup.f90.

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

◆ setgwpet()

subroutine uzfcellgroupmodule::setgwpet ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell 
)
private

Definition at line 453 of file UzfCellGroup.f90.

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

◆ sethead()

subroutine uzfcellgroupmodule::sethead ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  hgwf 
)
private

Definition at line 361 of file UzfCellGroup.f90.

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)

◆ setwaves()

subroutine uzfcellgroupmodule::setwaves ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell 
)
private

Definition at line 848 of file UzfCellGroup.f90.

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

◆ shift_waves()

subroutine uzfcellgroupmodule::shift_waves ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
integer(i4b), intent(in)  shft,
integer(i4b), intent(in)  strt,
integer(i4b), intent(in)  stp,
integer(i4b), intent(in)  cntr 
)
private

Waves strt to stp are overwritten by the waves shft positions away. cntr is the loop increment; it must run the loop in the direction that keeps the source ahead of the destination.

Definition at line 939 of file UzfCellGroup.f90.

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

◆ simgwet()

subroutine uzfcellgroupmodule::simgwet ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  igwetflag,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  hgwf,
real(dp), intent(inout)  trhs,
real(dp), intent(inout)  thcof,
real(dp), intent(inout)  det 
)
private

Definition at line 791 of file UzfCellGroup.f90.

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)
Here is the call graph for this function:

◆ solve()

subroutine uzfcellgroupmodule::solve ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  jbelow,
integer(i4b), intent(in)  icell,
real(dp), intent(inout)  totfluxtot,
integer(i4b), intent(in)  ietflag,
integer(i4b), intent(in)  issflag,
integer(i4b), intent(in)  iseepflag,
real(dp), intent(in)  hgwf,
real(dp), intent(in)  qfrommvr,
integer(i4b), intent(inout)  ierr,
logical, intent(in)  reset_state,
real(dp), intent(inout), optional  trhs,
real(dp), intent(inout), optional  thcof,
real(dp), intent(inout), optional  deriv,
real(dp), intent(inout), optional  watercontent 
)
private
Parameters
[in]jbelownumber of underlying uzf object or 0 if none
[in]icellnumber of this uzf object
[in]ietflaget is off (0) or based one water content (1) or pressure (2)
[in]issflagsteady state flag
[in]iseepflagdischarge to land is active (1) or not (0)
[in]hgwfhead for cell icell
[in]qfrommvrwater inflow from mover
[in,out]ierrflag indicating not enough waves
[in]reset_stateflag indicating that waves should be reset after solution
[in,out]trhstotal uzf rhs contribution to GWF model
[in,out]thcoftotal uzf hcof contribution to GWF model
[in,out]derivderivate term for contribution to GWF model
[in,out]watercontentcalculated water content

Definition at line 544 of file UzfCellGroup.f90.

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

◆ store_waves()

subroutine uzfcellgroupmodule::store_waves ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
type(uzfwavestoretype), intent(inout)  store 
)
private

Definition at line 960 of file UzfCellGroup.f90.

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)

◆ trailwav()

subroutine uzfcellgroupmodule::trailwav ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
integer(i4b), intent(inout)  ierr 
)
private

Definition at line 1115 of file UzfCellGroup.f90.

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
Here is the call graph for this function:

◆ unsat_stor()

real(dp) function uzfcellgroupmodule::unsat_stor ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(inout)  d1 
)
private

Definition at line 1432 of file UzfCellGroup.f90.

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

◆ update_wav()

subroutine uzfcellgroupmodule::update_wav ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  delt,
integer(i4b), intent(in)  iss,
integer(i4b), intent(in)  itest 
)
private

Definition at line 1476 of file UzfCellGroup.f90.

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

◆ uz_rise()

subroutine uzfcellgroupmodule::uz_rise ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(inout)  totfluxtot 
)
private

Definition at line 828 of file UzfCellGroup.f90.

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

◆ uzet()

subroutine uzfcellgroupmodule::uzet ( class(uzfcellgrouptype)  this,
integer(i4b), intent(in)  icell,
real(dp), intent(in)  delt,
integer(i4b), intent(in)  ietflag,
integer(i4b), intent(inout)  ierr 
)
private

Definition at line 1552 of file UzfCellGroup.f90.

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
Here is the call graph for this function:

◆ uzflow()

subroutine uzfcellgroupmodule::uzflow ( class(uzfcellgrouptype)  this,
real(dp), intent(inout)  thick,
real(dp), intent(inout)  thickold,
real(dp), intent(in)  delt,
integer(i4b), intent(in)  ietflag,
integer(i4b), intent(in)  icell,
integer(i4b), intent(inout)  ierr 
)
private

Definition at line 998 of file UzfCellGroup.f90.

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
Here is the call graph for this function: