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

Data Types

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, 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. 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 wave_shift (this, this2, icell, icell2, shft, strt, stp, cntr)
 Copy waves or shift waves in arrays. 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 704 of file UzfCellGroup.f90.

705  ! -- dummy
706  class(UzfCellGroupType) :: this
707  integer(I4B), intent(in) :: icell
708  integer(I4B), intent(in) :: jbelow
709  real(DP), intent(inout) :: trhs
710  real(DP), intent(inout) :: thcof
711  real(DP), intent(inout) :: deriv
712  real(DP), intent(in) :: delt
713  real(DP), intent(in) :: hgwf
714  ! -- local
715  real(DP) :: fcheck
716  real(DP) :: x, scale, range
717  !
718  ! -- initialize
719  range = dem5
720  deriv = dzero
721  thcof = dzero
722  trhs = this%uzfarea(icell) * this%totflux(icell) / delt
723  if (this%totflux(icell) < dem14) return
724  scale = done
725  !
726  ! -- smoothly reduce flow between cells when head close to cell top
727  x = hgwf - (this%celbot(icell) - range)
728  call sscurve(x, range, deriv, scale)
729  deriv = this%uzfarea(icell) * deriv * this%totflux(icell) / delt
730  this%finf(jbelow) = (done - scale) * this%totflux(icell) / delt
731  fcheck = this%finf(jbelow) - this%vks(jbelow)
732  !
733  ! -- reduce flow between cells when vks is too small
734  if (fcheck < dem14) fcheck = dzero
735  this%finf(jbelow) = this%finf(jbelow) - fcheck
736  this%surfluxbelow(icell) = this%finf(jbelow)
737  this%totflux(icell) = scale * this%totflux(icell) + fcheck * delt
738  trhs = this%uzfarea(icell) * this%totflux(icell) / delt
Here is the call graph for this function:

◆ advance()

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

Definition at line 551 of file UzfCellGroup.f90.

552  ! -- dummy
553  class(UzfCellGroupType) :: this
554  integer(I4B), intent(in) :: icell
555  !
556  ! -- set variables
557  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 1892 of file UzfCellGroup.f90.

1893  ! -- dummy
1894  class(UzfCellGroupType) :: this
1895  integer(I4B), intent(in) :: icell
1896  real(DP), intent(in) :: tho
1897  ! -- local
1898  real(DP) :: caph, lambda, star
1899  !
1900  caph = -dem6
1901  star = (tho - this%thtr(icell)) / (this%thts(icell) - this%thtr(icell))
1902  if (star < dem15) star = dem15
1903  lambda = dtwo / (this%eps(icell) - dthree)
1904  if (star > dem15) then
1905  if (tho - this%thts(icell) < dem15) then
1906  caph = this%ha(icell) * star**(-done / lambda)
1907  else
1908  caph = dzero
1909  end if
1910  end if

◆ dealloc()

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

Definition at line 241 of file UzfCellGroup.f90.

242  ! -- modules
243  use memorymanagermodule, only: mem_deallocate
244  ! -- dummy
245  class(UzfCellGroupType) :: this
246  !
247  ! -- deallocate based on whether or not memory manager was used
248  if (this%imem_manager == 0) then
249  deallocate (this%uzdpst)
250  deallocate (this%uzthst)
251  deallocate (this%uzflst)
252  deallocate (this%uzspst)
253  deallocate (this%nwavst)
254  deallocate (this%thtr)
255  deallocate (this%thts)
256  deallocate (this%thti)
257  deallocate (this%eps)
258  deallocate (this%ha)
259  deallocate (this%hroot)
260  deallocate (this%rootact)
261  deallocate (this%extwc)
262  deallocate (this%etact)
263  deallocate (this%nwav)
264  deallocate (this%ntrail)
265  deallocate (this%totflux)
266  deallocate (this%sinf)
267  deallocate (this%finf)
268  deallocate (this%finf_rej)
269  deallocate (this%gwet)
270  deallocate (this%uzfarea)
271  deallocate (this%cellarea)
272  deallocate (this%celtop)
273  deallocate (this%celbot)
274  deallocate (this%landtop)
275  deallocate (this%watab)
276  deallocate (this%watabold)
277  deallocate (this%surfdep)
278  deallocate (this%vks)
279  deallocate (this%surflux)
280  deallocate (this%surfluxbelow)
281  deallocate (this%surfseep)
282  deallocate (this%gwpet)
283  deallocate (this%pet)
284  deallocate (this%petmax)
285  deallocate (this%extdp)
286  deallocate (this%extdpuz)
287  deallocate (this%landflag)
288  deallocate (this%ivertcon)
289  else
290  call mem_deallocate(this%uzdpst)
291  call mem_deallocate(this%uzthst)
292  call mem_deallocate(this%uzflst)
293  call mem_deallocate(this%uzspst)
294  call mem_deallocate(this%nwavst)
295  call mem_deallocate(this%thtr)
296  call mem_deallocate(this%thts)
297  call mem_deallocate(this%thti)
298  call mem_deallocate(this%eps)
299  call mem_deallocate(this%ha)
300  call mem_deallocate(this%hroot)
301  call mem_deallocate(this%rootact)
302  call mem_deallocate(this%extwc)
303  call mem_deallocate(this%etact)
304  call mem_deallocate(this%nwav)
305  call mem_deallocate(this%ntrail)
306  call mem_deallocate(this%totflux)
307  call mem_deallocate(this%sinf)
308  call mem_deallocate(this%finf)
309  call mem_deallocate(this%finf_rej)
310  call mem_deallocate(this%gwet)
311  call mem_deallocate(this%uzfarea)
312  call mem_deallocate(this%cellarea)
313  call mem_deallocate(this%celtop)
314  call mem_deallocate(this%celbot)
315  call mem_deallocate(this%landtop)
316  call mem_deallocate(this%watab)
317  call mem_deallocate(this%watabold)
318  call mem_deallocate(this%surfdep)
319  call mem_deallocate(this%vks)
320  call mem_deallocate(this%surflux)
321  call mem_deallocate(this%surfluxbelow)
322  call mem_deallocate(this%surfseep)
323  call mem_deallocate(this%gwpet)
324  call mem_deallocate(this%pet)
325  call mem_deallocate(this%petmax)
326  call mem_deallocate(this%extdp)
327  call mem_deallocate(this%extdpuz)
328  call mem_deallocate(this%landflag)
329  call mem_deallocate(this%ivertcon)
330  end if
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 1069 of file UzfCellGroup.f90.

1070  ! -- dummy
1071  real(DP), intent(out) :: feps1
1072  real(DP), intent(out) :: feps2
1073  real(DP) :: factor1
1074  real(DP) :: factor2
1075  !
1076  ! calculate constants for uzflow
1077  factor1 = done
1078  factor2 = done
1079  feps1 = dem9
1080  feps2 = dem9
1081  if (itmuni == 1) then
1082  factor1 = done / 86400.d0
1083  else if (itmuni == 2) then
1084  factor1 = done / 1440.d0
1085  else if (itmuni == 3) then
1086  factor1 = done / 24.0d0
1087  else if (itmuni == 5) then
1088  factor1 = 365.0d0
1089  end if
1090  factor2 = done / 0.3048
1091  feps1 = feps1 * factor1 * factor2
1092  feps2 = feps2 * factor1 * factor2
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 1931 of file UzfCellGroup.f90.

1932  ! -- dummy
1933  class(UzfCellGroupType) :: this
1934  integer(I4B), intent(in) :: icell !< uzf cell containing depth
1935  real(DP), intent(in) :: depth !< depth within the cell
1936  ! -- return
1937  real(DP) :: theta_at_depth
1938  ! -- local
1939  real(DP) :: d1
1940  real(DP) :: d2
1941  real(DP) :: f1
1942  real(DP) :: f2
1943  !
1944  if (this%watab(icell) < this%celtop(icell)) then
1945  if (this%celtop(icell) - depth > this%watab(icell)) then
1946  d1 = depth - dem3
1947  d2 = depth + dem3
1948  f1 = this%unsat_stor(icell, d1)
1949  f2 = this%unsat_stor(icell, d2)
1950  theta_at_depth = this%thtr(icell) + (f2 - f1) / (d2 - d1)
1951  else
1952  theta_at_depth = this%thts(icell)
1953  end if
1954  else
1955  theta_at_depth = this%thts(icell)
1956  end if

◆ 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 1961 of file UzfCellGroup.f90.

1962  ! -- dummy
1963  class(UzfCellGroupType) :: this
1964  integer(I4B), intent(in) :: icell !< uzf cell containing depth
1965  ! -- return
1966  real(DP) :: watercontent
1967  ! -- local
1968  real(DP) :: top
1969  real(DP) :: bot
1970  real(DP) :: theta_r
1971  real(DP) :: thk
1972  real(DP) :: hgwf
1973  real(DP) :: fm
1974  real(DP) :: d
1975  !
1976  hgwf = this%watab(icell)
1977  top = this%celtop(icell)
1978  bot = this%celbot(icell)
1979  thk = top - max(bot, hgwf)
1980  if (thk > dzero) then
1981  theta_r = this%thtr(icell)
1982  d = thk
1983  fm = this%unsat_stor(icell, d)
1984  watercontent = fm / thk
1985  watercontent = watercontent + theta_r
1986  else
1987  watercontent = dzero
1988  end if

◆ 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 771 of file UzfCellGroup.f90.

772  ! -- dummy
773  class(UzfCellGroupType) :: this
774  integer(I4B), intent(in) :: icell
775  real(DP), intent(inout) :: deriv
776  real(DP), intent(inout) :: trhs
777  real(DP), intent(inout) :: thcof
778  real(DP), intent(inout) :: seep
779  real(DP), intent(out) :: scale
780  real(DP), intent(in) :: hgwf
781  ! -- local
782  real(DP) :: x, range, y, deriv1, d1, d2, Q
783  !
784  seep = dzero
785  deriv = dzero
786  deriv1 = dzero
787  d1 = dzero
788  d2 = dzero
789  scale = dzero
790  q = this%uzfarea(icell) * this%vks(icell)
791  range = this%surfdep(icell)
792  x = hgwf - this%celtop(icell)
793  call scubiclinear(x, range, deriv1, y)
794  scale = y
795  seep = scale * q * (hgwf - this%celtop(icell)) / range
796  trhs = scale * q * this%celtop(icell) / range
797  thcof = -scale * q / range
798  d1 = -deriv1 * q * x / range
799  d2 = -scale * q / range
800  deriv = d1 + d2
801  if (seep < dzero) then
802  seep = dzero
803  deriv = dzero
804  trhs = dzero
805  thcof = dzero
806  end if
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), optional  memory_path 
)
private

Definition at line 98 of file UzfCellGroup.f90.

99  ! -- modules
100  use memorymanagermodule, only: mem_allocate
101  ! -- dummy
102  class(UzfCellGroupType) :: this
103  integer(I4B), intent(in) :: nwav
104  integer(I4B), intent(in) :: ncells
105  character(len=*), intent(in), optional :: memory_path
106  ! -- local
107  integer(I4B) :: icell
108  !
109  ! -- Use mem_allocate if memory path is passed in, otherwise it's a temp object
110  if (present(memory_path)) then
111  this%imem_manager = 1
112  call mem_allocate(this%uzdpst, nwav, ncells, 'UZDPST', memory_path)
113  call mem_allocate(this%uzthst, nwav, ncells, 'UZTHST', memory_path)
114  call mem_allocate(this%uzflst, nwav, ncells, 'UZFLST', memory_path)
115  call mem_allocate(this%uzspst, nwav, ncells, 'UZSPST', memory_path)
116  call mem_allocate(this%nwavst, ncells, 'NWAVST', memory_path)
117  call mem_allocate(this%thtr, ncells, 'THTR', memory_path)
118  call mem_allocate(this%thts, ncells, 'THTS', memory_path)
119  call mem_allocate(this%thti, ncells, 'THTI', memory_path)
120  call mem_allocate(this%eps, ncells, 'EPS', memory_path)
121  call mem_allocate(this%ha, ncells, 'HA', memory_path)
122  call mem_allocate(this%hroot, ncells, 'HROOT', memory_path)
123  call mem_allocate(this%rootact, ncells, 'ROOTACT', memory_path)
124  call mem_allocate(this%extwc, ncells, 'EXTWC', memory_path)
125  call mem_allocate(this%etact, ncells, 'ETACT', memory_path)
126  call mem_allocate(this%nwav, ncells, 'NWAV', memory_path)
127  call mem_allocate(this%ntrail, ncells, 'NTRAIL', memory_path)
128  call mem_allocate(this%totflux, ncells, 'TOTFLUX', memory_path)
129  call mem_allocate(this%sinf, ncells, 'SINF', memory_path)
130  call mem_allocate(this%finf, ncells, 'FINF', memory_path)
131  call mem_allocate(this%finf_rej, ncells, 'FINF_REJ', memory_path)
132  call mem_allocate(this%gwet, ncells, 'GWET', memory_path)
133  call mem_allocate(this%uzfarea, ncells, 'UZFAREA', memory_path)
134  call mem_allocate(this%cellarea, ncells, 'CELLAREA', memory_path)
135  call mem_allocate(this%celtop, ncells, 'CELTOP', memory_path)
136  call mem_allocate(this%celbot, ncells, 'CELBOT', memory_path)
137  call mem_allocate(this%landtop, ncells, 'LANDTOP', memory_path)
138  call mem_allocate(this%watab, ncells, 'WATAB', memory_path)
139  call mem_allocate(this%watabold, ncells, 'WATABOLD', memory_path)
140  call mem_allocate(this%surfdep, ncells, 'SURFDEP', memory_path)
141  call mem_allocate(this%vks, ncells, 'VKS', memory_path)
142  call mem_allocate(this%surflux, ncells, 'SURFLUX', memory_path)
143  call mem_allocate(this%surfluxbelow, ncells, 'SURFLUXBELOW', memory_path)
144  call mem_allocate(this%surfseep, ncells, 'SURFSEEP', memory_path)
145  call mem_allocate(this%gwpet, ncells, 'GWPET', memory_path)
146  call mem_allocate(this%pet, ncells, 'PET', memory_path)
147  call mem_allocate(this%petmax, ncells, 'PETMAX', memory_path)
148  call mem_allocate(this%extdp, ncells, 'EXTDP', memory_path)
149  call mem_allocate(this%extdpuz, ncells, 'EXTDPUZ', memory_path)
150  call mem_allocate(this%landflag, ncells, 'LANDFLAG', memory_path)
151  call mem_allocate(this%ivertcon, ncells, 'IVERTCON', memory_path)
152  else
153  this%imem_manager = 0
154  allocate (this%uzdpst(nwav, ncells))
155  allocate (this%uzthst(nwav, ncells))
156  allocate (this%uzflst(nwav, ncells))
157  allocate (this%uzspst(nwav, ncells))
158  allocate (this%nwavst(ncells))
159  allocate (this%thtr(ncells))
160  allocate (this%thts(ncells))
161  allocate (this%thti(ncells))
162  allocate (this%eps(ncells))
163  allocate (this%ha(ncells))
164  allocate (this%hroot(ncells))
165  allocate (this%rootact(ncells))
166  allocate (this%extwc(ncells))
167  allocate (this%etact(ncells))
168  allocate (this%nwav(ncells))
169  allocate (this%ntrail(ncells))
170  allocate (this%totflux(ncells))
171  allocate (this%sinf(ncells))
172  allocate (this%finf(ncells))
173  allocate (this%finf_rej(ncells))
174  allocate (this%gwet(ncells))
175  allocate (this%uzfarea(ncells))
176  allocate (this%cellarea(ncells))
177  allocate (this%celtop(ncells))
178  allocate (this%celbot(ncells))
179  allocate (this%landtop(ncells))
180  allocate (this%watab(ncells))
181  allocate (this%watabold(ncells))
182  allocate (this%surfdep(ncells))
183  allocate (this%vks(ncells))
184  allocate (this%surflux(ncells))
185  allocate (this%surfluxbelow(ncells))
186  allocate (this%surfseep(ncells))
187  allocate (this%gwpet(ncells))
188  allocate (this%pet(ncells))
189  allocate (this%petmax(ncells))
190  allocate (this%extdp(ncells))
191  allocate (this%extdpuz(ncells))
192  allocate (this%landflag(ncells))
193  allocate (this%ivertcon(ncells))
194  end if
195  do icell = 1, ncells
196  this%uzdpst(:, icell) = dzero
197  this%uzthst(:, icell) = dzero
198  this%uzflst(:, icell) = dzero
199  this%uzspst(:, icell) = dzero
200  this%nwavst(icell) = 1
201  this%thtr(icell) = dzero
202  this%thts(icell) = dzero
203  this%thti(icell) = dzero
204  this%eps(icell) = dzero
205  this%ha(icell) = dzero
206  this%hroot(icell) = dzero
207  this%rootact(icell) = dzero
208  this%extwc(icell) = dzero
209  this%etact(icell) = dzero
210  this%nwav(icell) = nwav
211  this%ntrail(icell) = 0
212  this%totflux(icell) = dzero
213  this%sinf(icell) = dzero
214  this%finf(icell) = dzero
215  this%finf_rej(icell) = dzero
216  this%gwet(icell) = dzero
217  this%uzfarea(icell) = dzero
218  this%cellarea(icell) = dzero
219  this%celtop(icell) = dzero
220  this%celbot(icell) = dzero
221  this%landtop(icell) = dzero
222  this%watab(icell) = dzero
223  this%watabold(icell) = dzero
224  this%surfdep(icell) = dzero
225  this%vks(icell) = dzero
226  this%surflux(icell) = dzero
227  this%surfluxbelow(icell) = dzero
228  this%surfseep(icell) = dzero
229  this%gwpet(icell) = dzero
230  this%pet(icell) = dzero
231  this%petmax(icell) = dzero
232  this%extdp(icell) = dzero
233  this%extdpuz(icell) = dzero
234  this%landflag(icell) = 0
235  this%ivertcon(icell) = 0
236  end do
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 1384 of file UzfCellGroup.f90.

1385  ! -- Return
1386  real(DP) :: leadspeed
1387  ! -- dummy
1388  real(DP), intent(in) :: theta1
1389  real(DP), intent(in) :: theta2
1390  real(DP), intent(in) :: flux1
1391  real(DP), intent(inout) :: flux2
1392  real(DP), intent(in) :: thts
1393  real(DP), intent(in) :: thtr
1394  real(DP), intent(in) :: eps
1395  real(DP), intent(in) :: vks
1396  ! -- local
1397  real(DP) :: comp1, comp2, thsrinv, epsfksths
1398  real(DP) :: eps_m1, fhold, comp3
1399  !
1400  eps_m1 = eps - done
1401  thsrinv = done / (thts - thtr)
1402  epsfksths = eps * vks * thsrinv
1403  comp1 = theta2 - theta1
1404  comp2 = abs(flux2 - flux1)
1405  comp3 = theta1 - thtr
1406  if (comp2 < dem15) flux2 = flux1 + dem15
1407  if (abs(comp1) < dem30) then
1408  fhold = dem30
1409  if (comp3 > dem30) fhold = (comp3 * thsrinv)**eps
1410  if (fhold < dem30) fhold = dem30
1411  leadspeed = epsfksths * (fhold**eps_m1)
1412  else
1413  leadspeed = (flux2 - flux1) / (theta2 - theta1)
1414  end if
1415  if (leadspeed < dem30) leadspeed = dem30
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 1183 of file UzfCellGroup.f90.

1185  ! -- dummy
1186  class(UzfCellGroupType) :: this
1187  real(DP), intent(inout) :: thetab
1188  real(DP), intent(inout) :: fluxb
1189  real(DP), intent(in) :: feps2
1190  real(DP), intent(inout) :: time
1191  integer(I4B), intent(inout) :: itester
1192  integer(I4B), intent(inout) :: itrailflg
1193  real(DP), intent(inout) :: ffcheck
1194  real(DP), intent(in) :: delt
1195  integer(I4B), intent(in) :: icell
1196  ! -- local
1197  real(DP) :: bottomtime, shortest, fcheck
1198  real(DP) :: eps_m1, timenew, bottom, timedt
1199  real(DP) :: thtsrinv, diff, fluxhld2
1200  real(DP) :: flux1, flux2, theta1, theta2, ftest
1201  real(DP), allocatable, dimension(:) :: checktime
1202  integer(I4B) :: iflx, iremove, j, l
1203  integer(I4B) :: nwavp1, jshort
1204  integer(I4B), allocatable, dimension(:) :: more
1205  !
1206  allocate (checktime(this%nwavst(icell)))
1207  allocate (more(this%nwavst(icell)))
1208  ftest = dzero
1209  eps_m1 = dble(this%eps(icell)) - done
1210  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1211  !
1212  ! -- initialize new wave
1213  if (itrailflg == 0) then
1214  if (ffcheck > feps2) then
1215  this%uzflst(this%nwavst(icell), icell) = this%surflux(icell)
1216  if (this%uzflst(this%nwavst(icell), icell) < dem30) &
1217  this%uzflst(this%nwavst(icell), icell) = dzero
1218  this%uzthst(this%nwavst(icell), icell) = &
1219  (((this%uzflst(this%nwavst(icell), icell) / this%vks(icell))** &
1220  (done / this%eps(icell))) * (this%thts(icell) - this%thtr(icell))) &
1221  + this%thtr(icell)
1222  theta2 = this%uzthst(this%nwavst(icell), icell)
1223  flux2 = this%uzflst(this%nwavst(icell), icell)
1224  flux1 = this%uzflst(this%nwavst(icell) - 1, icell)
1225  theta1 = this%uzthst(this%nwavst(icell) - 1, icell)
1226  this%uzspst(this%nwavst(icell), icell) = &
1227  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1228  this%thtr(icell), this%eps(icell), this%vks(icell))
1229  this%uzdpst(this%nwavst(icell), icell) = dzero
1230  end if
1231  end if
1232  !
1233  ! -- route all waves and interception of waves over times step
1234  diff = done
1235  timedt = dzero
1236  iflx = 0
1237  fluxhld2 = this%uzflst(1, icell)
1238  if (this%nwavst(icell) == 0) itester = 1
1239  if (itester /= 1) then
1240  do while (diff > dem6)
1241  nwavp1 = this%nwavst(icell) + 1
1242  timedt = delt - time
1243  do j = 1, this%nwavst(icell)
1244  checktime(j) = dep20
1245  more(j) = 0
1246  end do
1247  shortest = timedt
1248  if (this%nwavst(icell) > 2) then
1249  j = 2
1250  !
1251  ! -- calculate time until wave overtakes wave ahead
1252  nwavp1 = this%nwavst(icell) + 1
1253  do while (j < nwavp1)
1254  ftest = this%uzspst(j - 1, icell) - this%uzspst(j, icell)
1255  if (abs(ftest) > dem30) then
1256  checktime(j) = (this%uzdpst(j, icell) - &
1257  this%uzdpst(j - 1, icell)) / ftest
1258  if (checktime(j) < dem30) checktime(j) = dep20
1259  end if
1260  j = j + 1
1261  end do
1262  end if
1263  !
1264  ! - calc time until wave reaches bottom of cell
1265  bottomtime = dep20
1266  if (this%nwavst(icell) > 1) then
1267  if (this%uzspst(2, icell) > dzero) then
1268  bottom = this%uzspst(2, icell)
1269  if (bottom < dem15) bottom = dem15
1270  bottomtime = (this%uzdpst(1, icell) - this%uzdpst(2, icell)) / bottom
1271  if (bottomtime < dzero) bottomtime = dem12
1272  end if
1273  end if
1274  !
1275  ! -- calc time for wave interception
1276  jshort = 0
1277  do j = this%nwavst(icell), 3, -1
1278  if (shortest - checktime(j) > -dem9) then
1279  more(j) = 1
1280  jshort = j
1281  shortest = checktime(j)
1282  end if
1283  end do
1284  do j = 3, this%nwavst(icell)
1285  if (shortest - checktime(j) < dem9) then
1286  if (j /= jshort) more(j) = 0
1287  end if
1288  end do
1289  !
1290  ! -- what happens first, waves hits bottom or interception
1291  iremove = 0
1292  timenew = time
1293  fcheck = (time + shortest) - delt
1294  if (shortest < dem7) fcheck = -done
1295  if (bottomtime < shortest .AND. time + bottomtime < delt) then
1296  j = 2
1297  do while (j < nwavp1)
1298  !
1299  ! -- route waves
1300  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1301  this%uzspst(j, icell) * bottomtime
1302  j = j + 1
1303  end do
1304  fluxb = this%uzflst(2, icell)
1305  thetab = this%uzthst(2, icell)
1306  iflx = 1
1307  call this%wave_shift(this, icell, icell, 1, 1, &
1308  this%nwavst(icell) - 1, 1)
1309  iremove = 1
1310  timenew = time + bottomtime
1311  this%uzspst(1, icell) = dzero
1312  !
1313  ! -- do waves intercept before end of time step
1314  else if (fcheck < dzero .AND. this%nwavst(icell) > 2) then
1315  j = 2
1316  do while (j < nwavp1)
1317  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1318  this%uzspst(j, icell) * shortest
1319  j = j + 1
1320  end do
1321  !
1322  ! -- combine waves that intercept, remove a wave
1323  j = 3
1324  l = j
1325  do while (j < this%nwavst(icell) + 1)
1326  if (more(j) == 1) then
1327  l = j
1328  theta2 = this%uzthst(j, icell)
1329  flux2 = this%uzflst(j, icell)
1330  if (j == 3) then
1331  flux1 = fluxb
1332  theta1 = thetab
1333  else
1334  flux1 = this%uzflst(j - 2, icell)
1335  theta1 = this%uzthst(j - 2, icell)
1336  end if
1337  this%uzspst(j, icell) = leadspeed(theta1, theta2, flux1, flux2, &
1338  this%thts(icell), &
1339  this%thtr(icell), &
1340  this%eps(icell), this%vks(icell))
1341  !
1342  ! -- update waves
1343  call this%wave_shift(this, icell, icell, 1, l - 1, &
1344  this%nwavst(icell) - 1, 1)
1345  l = this%nwavst(icell) + 1
1346  iremove = iremove + 1
1347  end if
1348  j = j + 1
1349  end do
1350  timenew = timenew + shortest
1351  !
1352  ! -- calc. total flux to bottom during remaining time in step
1353  else
1354  j = 2
1355  do while (j < nwavp1)
1356  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1357  this%uzspst(j, icell) * timedt
1358  j = j + 1
1359  end do
1360  timenew = delt
1361  end if
1362  this%totflux(icell) = this%totflux(icell) + fluxhld2 * (timenew - time)
1363  if (iflx == 1) then
1364  fluxhld2 = this%uzflst(1, icell)
1365  iflx = 0
1366  end if
1367  !
1368  ! -- remove dead waves
1369  this%nwavst(icell) = this%nwavst(icell) - iremove
1370  time = timenew
1371  diff = delt - time
1372  if (this%nwavst(icell) == 1) then
1373  itester = 1
1374  exit
1375  end if
1376  end do
1377  end if
1378  deallocate (checktime)
1379  deallocate (more)
Here is the call graph for this function:

◆ 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 1914 of file UzfCellGroup.f90.

1915  ! -- Return
1916  real(DP) :: rate_et_z
1917  ! -- dummy
1918  class(UzfCellGroupType) :: this
1919  integer(I4B), intent(in) :: icell
1920  real(DP), intent(in) :: factor, fktho, h
1921  !
1922  rate_et_z = factor * fktho * (h - this%hroot(icell))
1923  if (rate_et_z < dzero) rate_et_z = dzero

◆ 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 743 of file UzfCellGroup.f90.

744  ! -- dummy
745  class(UzfCellGroupType) :: this
746  integer(I4B), intent(in) :: icell
747  real(DP), intent(inout) :: deriv
748  real(DP), intent(inout) :: finfact
749  real(DP), intent(inout) :: thcof
750  real(DP), intent(inout) :: trhs
751  real(DP), intent(in) :: hgwf
752  ! -- local
753  real(DP) :: x, range, scale, q
754  !
755  range = this%surfdep(icell)
756  q = this%surflux(icell)
757  finfact = q
758  trhs = finfact * this%uzfarea(icell)
759  x = this%celtop(icell) - hgwf
760  call slinear(x, range, deriv, scale)
761  deriv = -q * deriv * this%uzfarea(icell) * scale
762  if (scale < done) then
763  finfact = q * scale
764  trhs = finfact * this%uzfarea(icell) * this%celtop(icell) / range
765  thcof = finfact * this%uzfarea(icell) / range
766  end if
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 917 of file UzfCellGroup.f90.

918  ! -- dummy
919  class(UzfCellGroupType) :: this
920  real(DP), intent(inout) :: totfluxtot
921  real(DP), intent(in) :: delt
922  integer(I4B), intent(in) :: ietflag
923  integer(I4B), intent(in) :: icell
924  integer(I4B), intent(inout) :: ierr
925  ! -- local
926  real(DP) :: thick, thickold
927  integer(I4B) :: idelt, iwav, ik
928  !
929  ! -- initialize
930  this%totflux(icell) = dzero
931  this%etact(icell) = dzero
932  thick = this%celtop(icell) - this%watab(icell)
933  thickold = this%celtop(icell) - this%watabold(icell)
934  !
935  ! -- no uz, clear waves
936  if (thickold < dzero) then
937  do iwav = 1, this%nwavst(icell)
938  this%uzthst(iwav, icell) = this%thtr(icell)
939  this%uzdpst(iwav, icell) = dzero
940  this%uzspst(iwav, icell) = dzero
941  this%uzflst(iwav, icell) = dzero
942  end do
943  this%nwavst(icell) = 1
944  end if
945  idelt = 1
946  do ik = 1, idelt
947  call this%uzflow(thick, thickold, delt, ietflag, icell, ierr)
948  if (ierr > 0) return
949  totfluxtot = totfluxtot + this%totflux(icell)
950  end do

◆ setbelowpet()

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

Definition at line 491 of file UzfCellGroup.f90.

492  ! -- modules
493  use tdismodule, only: delt
494  ! -- dummy
495  class(UzfCellGroupType) :: this
496  integer(I4B), intent(in) :: icell
497  integer(I4B), intent(in) :: jbelow
498  ! -- dummy
499  real(DP) :: pet
500  !
501  pet = dzero
502  !
503  ! -- transfer unmet pet to lower cell
504  !
505  if (this%extdpuz(jbelow) > dem3) then
506  pet = this%pet(icell) - this%etact(icell) / delt - &
507  this%gwet(icell) / this%uzfarea(icell)
508  if (pet < dzero) pet = dzero
509  end if
510  this%pet(jbelow) = pet
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 335 of file UzfCellGroup.f90.

337  ! -- dummy
338  class(UzfCellGroupType) :: this
339  integer(I4B), intent(in) :: icell
340  real(DP), intent(in) :: area
341  real(DP), intent(in) :: top
342  real(DP), intent(in) :: bot
343  real(DP), intent(in) :: surfdep
344  real(DP), intent(in) :: vks
345  real(DP), intent(in) :: thtr
346  real(DP), intent(in) :: thts
347  real(DP), intent(in) :: thti
348  real(DP), intent(in) :: eps
349  integer(I4B), intent(in) :: ntrail
350  integer(I4B), intent(in) :: landflag
351  integer(I4B), intent(in) :: ivertcon
352  !
353  ! -- set the values for uzf cell icell
354  this%landflag(icell) = landflag
355  this%ivertcon(icell) = ivertcon
356  this%surfdep(icell) = surfdep
357  this%uzfarea(icell) = area
358  this%cellarea(icell) = area
359  if (this%landflag(icell) == 1) then
360  this%celtop(icell) = top - dhalf * this%surfdep(icell)
361  else
362  this%celtop(icell) = top
363  end if
364  this%celbot(icell) = bot
365  this%vks(icell) = vks
366  this%thtr(icell) = thtr
367  this%thts(icell) = thts
368  this%thti(icell) = thti
369  this%eps(icell) = eps
370  this%ntrail(icell) = ntrail
371  this%pet(icell) = dzero
372  this%extdp(icell) = dzero
373  this%extwc(icell) = dzero
374  this%ha(icell) = dzero
375  this%hroot(icell) = dzero

◆ 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 428 of file UzfCellGroup.f90.

429  ! -- dummy
430  class(UzfCellGroupType) :: this
431  integer(I4B), intent(in) :: icell
432  integer(I4B), intent(in) :: jbelow
433  real(DP), intent(in) :: pet
434  real(DP), intent(in) :: extdp
435  ! -- local
436  real(DP) :: thick
437  !
438  if (this%landflag(icell) == 1) then
439  this%pet(icell) = pet
440  this%gwpet(icell) = pet
441  else
442  this%pet(icell) = dzero
443  this%gwpet(icell) = dzero
444  end if
445  thick = this%celtop(icell) - this%celbot(icell)
446  this%extdp(icell) = extdp
447  if (this%landflag(icell) > 0) then
448  this%landtop(icell) = this%celtop(icell)
449  this%petmax(icell) = this%pet(icell)
450  end if
451  !
452  ! -- set uz extinction depth
453  if (this%landtop(icell) - this%extdp(icell) < this%celbot(icell)) then
454  this%extdpuz(icell) = thick
455  else
456  this%extdpuz(icell) = this%celtop(icell) - &
457  (this%landtop(icell) - this%extdp(icell))
458  end if
459  if (this%extdpuz(icell) < dzero) this%extdpuz(icell) = dzero
460  if (this%extdpuz(icell) > dem7 .and. this%extdp(icell) < dem7) &
461  this%extdp(icell) = this%extdpuz(icell)
462  !
463  ! -- set pet for underlying cell
464  if (jbelow > 0) then
465  this%landtop(jbelow) = this%landtop(icell)
466  this%petmax(jbelow) = this%petmax(icell)
467  end if

◆ 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 529 of file UzfCellGroup.f90.

530  ! -- dummy
531  class(UzfCellGroupType) :: this
532  integer(I4B), intent(in) :: icell
533  integer(I4B), intent(in) :: jbelow
534  real(DP), intent(in) :: ha
535  real(DP), intent(in) :: hroot
536  real(DP), intent(in) :: rootact
537  !
538  ! -- set variables
539  this%ha(icell) = ha
540  this%hroot(icell) = hroot
541  this%rootact(icell) = rootact
542  if (jbelow > 0) then
543  this%ha(jbelow) = ha
544  this%hroot(jbelow) = hroot
545  this%rootact(jbelow) = rootact
546  end if

◆ 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 515 of file UzfCellGroup.f90.

516  ! -- dummy
517  class(UzfCellGroupType) :: this
518  integer(I4B), intent(in) :: icell
519  integer(I4B), intent(in) :: jbelow
520  real(DP), intent(in) :: extwc
521  !
522  ! -- set extinction water content
523  this%extwc(icell) = extwc
524  if (jbelow > 0) this%extwc(jbelow) = extwc

◆ setdatafinf()

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

Definition at line 396 of file UzfCellGroup.f90.

397  ! -- dummy
398  class(UzfCellGroupType) :: this
399  integer(I4B), intent(in) :: icell
400  real(DP), intent(in) :: finf
401  !
402  if (this%landflag(icell) == 1) then
403  this%sinf(icell) = finf
404  this%finf(icell) = finf
405  else
406  this%sinf(icell) = dzero
407  this%finf(icell) = dzero
408  end if
409  this%finf_rej(icell) = dzero
410  this%surflux(icell) = dzero
411  this%surfluxbelow(icell) = dzero

◆ setdatauzfarea()

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

Definition at line 416 of file UzfCellGroup.f90.

417  ! -- dummy
418  class(UzfCellGroupType) :: this
419  integer(I4B), intent(in) :: icell
420  real(DP), intent(in) :: areamult
421  !
422  ! -- set uzf area
423  this%uzfarea(icell) = this%cellarea(icell) * areamult

◆ setgwpet()

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

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  ! -- dummy
479  real(DP) :: pet
480  !
481  pet = dzero
482  !
483  ! -- reduce pet for gw by uzet
484  pet = this%pet(icell) - this%etact(icell) / delt
485  if (pet < dzero) pet = dzero
486  this%gwpet(icell) = pet

◆ sethead()

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

Definition at line 380 of file UzfCellGroup.f90.

381  ! -- dummy
382  class(UzfCellGroupType) :: this
383  integer(I4B), intent(in) :: icell
384  real(DP), intent(in) :: hgwf
385  !
386  ! -- set initial head
387  this%watab(icell) = this%celbot(icell)
388  if (hgwf > this%celbot(icell)) this%watab(icell) = hgwf
389  if (this%watab(icell) > this%celtop(icell)) &
390  this%watab(icell) = this%celtop(icell)
391  this%watabold(icell) = this%watab(icell)

◆ setwaves()

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

Definition at line 868 of file UzfCellGroup.f90.

869  ! -- dummy
870  class(UzfCellGroupType) :: this
871  ! -- local
872  integer(I4B), intent(in) :: icell
873  real(DP) :: bottom, top
874  integer(I4B) :: jk
875  real(DP) :: thick
876  !
877  ! -- initialize
878  this%totflux(icell) = dzero
879  this%nwavst(icell) = 1
880  this%uzdpst(:, icell) = dzero
881  thick = this%celtop(icell) - this%watab(icell)
882  do jk = 1, this%nwav(icell)
883  this%uzthst(jk, icell) = this%thtr(icell)
884  end do
885  !
886  ! -- initialize waves for first stress period
887  if (thick > dzero) then
888  this%uzdpst(1, icell) = thick
889  this%uzthst(1, icell) = this%thti(icell)
890  top = this%uzthst(1, icell) - this%thtr(icell)
891  if (top < dzero) top = dzero
892  bottom = this%thts(icell) - this%thtr(icell)
893  if (bottom < dzero) bottom = dzero
894  this%uzflst(1, icell) = this%vks(icell) * (top / bottom)**this%eps(icell)
895  if (this%uzthst(1, icell) < this%thtr(icell)) &
896  this%uzthst(1, icell) = this%thtr(icell)
897  !
898  ! -- calculate water stored in the unsaturated zone
899  if (top > dzero) then
900  this%uzspst(1, icell) = dzero
901  else
902  this%uzflst(1, icell) = dzero
903  this%uzspst(1, icell) = dzero
904  end if
905  !
906  ! no unsaturated zone
907  else
908  this%uzflst(1, icell) = dzero
909  this%uzdpst(1, icell) = dzero
910  this%uzspst(1, icell) = dzero
911  this%uzthst(1, icell) = this%thtr(icell)
912  end if

◆ 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 811 of file UzfCellGroup.f90.

812  ! -- dummy
813  class(UzfCellGroupType) :: this
814  integer(I4B), intent(in) :: igwetflag
815  integer(I4B), intent(in) :: icell
816  real(DP), intent(in) :: hgwf
817  real(DP), intent(inout) :: trhs
818  real(DP), intent(inout) :: thcof
819  real(DP), intent(inout) :: det
820  ! -- local
821  real(DP) :: s, x, c, b, et
822  !
823  this%gwet(icell) = dzero
824  trhs = dzero
825  thcof = dzero
826  det = dzero
827  s = this%landtop(icell)
828  x = this%extdp(icell)
829  c = this%gwpet(icell)
830  b = this%celbot(icell)
831  if (b > hgwf) return
832  if (x < dem6) return
833  if (igwetflag == 1) then
834  et = etfunc_lin(s, x, c, det, trhs, thcof, hgwf, &
835  this%celtop(icell), this%celbot(icell))
836  else if (igwetflag == 2) then
837  et = etfunc_nlin(s, x, c, det, trhs, thcof, hgwf)
838  end if
839  ! this%gwet(icell) = et * this%uzfarea(icell)
840  trhs = trhs * this%uzfarea(icell)
841  thcof = thcof * this%uzfarea(icell)
842  this%gwet(icell) = trhs - (thcof * hgwf)
843  ! write(99,*)'in group', icell, this%gwet(icell)
Here is the call graph for this function:

◆ solve()

subroutine uzfcellgroupmodule::solve ( class(uzfcellgrouptype)  this,
type(uzfcellgrouptype)  thiswork,
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
thisworkwork object for resetting wave state
[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 563 of file UzfCellGroup.f90.

566  ! -- modules
567  use tdismodule, only: delt
568  ! -- dummy
569  class(UzfCellGroupType) :: this
570  type(UzfCellGroupType) :: thiswork !< work object for resetting wave state
571  integer(I4B), intent(in) :: jbelow !< number of underlying uzf object or 0 if none
572  integer(I4B), intent(in) :: icell !< number of this uzf object
573  real(DP), intent(inout) :: totfluxtot !<
574  integer(I4B), intent(in) :: ietflag !< et is off (0) or based one water content (1) or pressure (2)
575  integer(I4B), intent(in) :: issflag !< steady state flag
576  integer(I4B), intent(in) :: iseepflag !< discharge to land is active (1) or not (0)
577  real(DP), intent(in) :: hgwf !< head for cell icell
578  real(DP), intent(in) :: qfrommvr !< water inflow from mover
579  integer(I4B), intent(inout) :: ierr !< flag indicating not enough waves
580  logical, intent(in) :: reset_state !< flag indicating that waves should be reset after solution
581  real(DP), intent(inout), optional :: trhs !< total uzf rhs contribution to GWF model
582  real(DP), intent(inout), optional :: thcof !< total uzf hcof contribution to GWF model
583  real(DP), intent(inout), optional :: deriv !< derivate term for contribution to GWF model
584  real(DP), intent(inout), optional :: watercontent !< calculated water content
585  ! -- local
586  real(DP) :: test
587  real(DP) :: scale
588  real(DP) :: seep
589  real(DP) :: finfact
590  real(DP) :: derivfinf
591  real(DP) :: trhsfinf
592  real(DP) :: thcoffinf
593  real(DP) :: trhsseep
594  real(DP) :: thcofseep
595  real(DP) :: deriv1
596  real(DP) :: deriv2
597  !
598  ! -- initialize variables
599  totfluxtot = dzero
600  trhsfinf = dzero
601  thcoffinf = dzero
602  trhsseep = dzero
603  thcofseep = dzero
604  this%finf_rej(icell) = dzero
605  this%surflux(icell) = this%finf(icell) + qfrommvr / this%uzfarea(icell)
606  this%watab(icell) = hgwf
607  this%surfseep(icell) = dzero
608  seep = dzero
609  finfact = dzero
610  this%etact(icell) = dzero
611  this%surfluxbelow(icell) = dzero
612  if (this%ivertcon(icell) > 0) then
613  this%finf(jbelow) = dzero
614  end if
615  if (this%watab(icell) < this%celbot(icell)) &
616  this%watab(icell) = this%celbot(icell)
617  !
618  ! -- initialize derivative variables
619  deriv1 = dzero
620  deriv2 = dzero
621  derivfinf = dzero
622  !
623  ! -- save wave states for resetting after iteration.
624  if (reset_state) then
625  call thiswork%wave_shift(this, 1, icell, 0, 1, this%nwavst(icell), 1)
626  end if
627  !
628  if (this%watab(icell) > this%celtop(icell)) &
629  this%watab(icell) = this%celtop(icell)
630  !
631  ! -- add water from mover to applied infiltration.
632  if (this%surflux(icell) > this%vks(icell)) then
633  this%surflux(icell) = this%vks(icell)
634  end if
635  !
636  ! -- saturation excess rejected infiltration
637  if (this%landflag(icell) == 1) then
638  call this%rejfinf(icell, deriv1, hgwf, trhsfinf, thcoffinf, finfact)
639  this%surflux(icell) = finfact
640  end if
641  !
642  ! -- calculate rejected infiltration
643  this%finf_rej(icell) = this%finf(icell) + &
644  (qfrommvr / this%uzfarea(icell)) - this%surflux(icell)
645  !
646  ! -- calculate groundwater discharge
647  if (iseepflag > 0 .and. this%landflag(icell) == 1) then
648  call this%gwseep(icell, deriv2, scale, hgwf, trhsseep, thcofseep, seep)
649  this%surfseep(icell) = seep
650  end if
651  !
652  ! -- route water through unsat zone, calc. storage change and recharge
653  test = this%watab(icell)
654  if (this%watabold(icell) - test < -dem15) test = this%watabold(icell)
655  if (this%celtop(icell) - test > dem15) then
656  if (issflag == 0) then
657  call this%routewaves(totfluxtot, delt, ietflag, icell, ierr)
658  if (ierr > 0) then
659  if (reset_state) &
660  call this%wave_shift(thiswork, icell, 1, 0, 1, thiswork%nwavst(1), 1)
661  return
662  end if
663  call this%uz_rise(icell, totfluxtot)
664  this%totflux(icell) = totfluxtot
665  if (this%ivertcon(icell) > 0) then
666  call this%addrech(icell, jbelow, hgwf, trhsfinf, thcoffinf, &
667  derivfinf, delt)
668  end if
669  else
670  this%totflux(icell) = this%surflux(icell) * delt
671  totfluxtot = this%surflux(icell) * delt
672  end if
673  thcoffinf = dzero
674  trhsfinf = this%totflux(icell) * this%uzfarea(icell) / delt
675  if (.not. reset_state) then
676  call this%update_wav(icell, delt, issflag, 0)
677  end if
678  else
679  this%totflux(icell) = this%surflux(icell) * delt
680  totfluxtot = this%surflux(icell) * delt
681  if (.not. reset_state) then
682  call this%update_wav(icell, delt, issflag, 1)
683  end if
684  end if
685  !
686  ! -- If formulating, then these variables will be present
687  if (present(deriv)) deriv = deriv1 + deriv2 + derivfinf
688  if (present(trhs)) trhs = trhsfinf + trhsseep
689  if (present(thcof)) thcof = thcoffinf + thcofseep
690  !
691  ! -- Assign water content prior to resetting waves
692  if (present(watercontent)) then
693  watercontent = this%get_wcnew(icell)
694  end if
695  !
696  ! -- reset waves to previous state for next iteration
697  if (reset_state) then
698  call this%wave_shift(thiswork, icell, 1, 0, 1, thiswork%nwavst(1), 1)
699  end if

◆ trailwav()

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

Definition at line 1097 of file UzfCellGroup.f90.

1098  ! -- dummy
1099  class(UzfCellGroupType) :: this
1100  integer(I4B), intent(in) :: icell
1101  integer(I4B), intent(inout) :: ierr
1102  ! -- local
1103  real(DP) :: smoist, smoistinc, ftrail, eps_m1
1104  real(DP) :: thtsrinv
1105  real(DP) :: flux1, flux2, theta1, theta2
1106  real(DP) :: fnuminc
1107  integer(I4B) :: j, jj, jk, nwavstm1
1108  !
1109  ! -- initialize
1110  eps_m1 = dble(this%eps(icell)) - done
1111  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1112  nwavstm1 = this%nwavst(icell) - 1
1113  !
1114  ! -- initialize trailwaves
1115  smoist = (((this%surflux(icell) / this%vks(icell))** &
1116  (done / this%eps(icell))) * &
1117  (this%thts(icell) - this%thtr(icell))) + this%thtr(icell)
1118  if (this%uzthst(nwavstm1, icell) - smoist > dem9) then
1119  fnuminc = dzero
1120  do jk = 1, this%ntrail(icell)
1121  fnuminc = fnuminc + float(jk)
1122  end do
1123  smoistinc = (this%uzthst(nwavstm1, icell) - smoist) / (fnuminc - done)
1124  jj = this%ntrail(icell)
1125  ftrail = dble(this%ntrail(icell)) + done
1126  do j = this%nwavst(icell), this%nwavst(icell) + this%ntrail(icell) - 1
1127  if (j > this%nwav(icell)) then
1128  ! -- too many waves error
1129  ierr = 1
1130  return
1131  end if
1132  if (j > this%nwavst(icell)) then
1133  this%uzthst(j, icell) = this%uzthst(j - 1, icell) &
1134  - ((ftrail - float(jj)) * smoistinc)
1135  else
1136  this%uzthst(j, icell) = this%uzthst(j - 1, icell) - dem9
1137  end if
1138  jj = jj - 1
1139  if (this%uzthst(j, icell) <= this%thtr(icell) + dem9) &
1140  this%uzthst(j, icell) = this%thtr(icell) + dem9
1141  this%uzflst(j, icell) = &
1142  this%vks(icell) * (((this%uzthst(j, icell) - this%thtr(icell)) * &
1143  thtsrinv)**this%eps(icell))
1144  theta2 = this%uzthst(j - 1, icell)
1145  flux2 = this%uzflst(j - 1, icell)
1146  flux1 = this%uzflst(j, icell)
1147  theta1 = this%uzthst(j, icell)
1148  this%uzspst(j, icell) = leadspeed(theta1, theta2, flux1, flux2, &
1149  this%thts(icell), this%thtr(icell), &
1150  this%eps(icell), this%vks(icell))
1151  this%uzdpst(j, icell) = dzero
1152  if (j == this%nwavst(icell)) then
1153  this%uzdpst(j, icell) = this%uzdpst(j, icell) + &
1154  (this%ntrail(icell) + 1) * dem9
1155  else
1156  this%uzdpst(j, icell) = this%uzdpst(j - 1, icell) - dem9
1157  end if
1158  end do
1159  this%nwavst(icell) = this%nwavst(icell) + this%ntrail(icell) - 1
1160  if (this%nwavst(icell) >= this%nwav(icell)) then
1161  ! -- too many waves error
1162  ierr = 1
1163  return
1164  end if
1165  else
1166  this%uzdpst(this%nwavst(icell), icell) = dzero
1167  this%uzflst(this%nwavst(icell), icell) = &
1168  this%vks(icell) * (((this%uzthst(this%nwavst(icell), icell) - &
1169  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1170  this%uzthst(this%nwavst(icell), icell) = smoist
1171  theta2 = this%uzthst(this%nwavst(icell) - 1, icell)
1172  flux2 = this%uzflst(this%nwavst(icell) - 1, icell)
1173  flux1 = this%uzflst(this%nwavst(icell), icell)
1174  theta1 = this%uzthst(this%nwavst(icell), icell)
1175  this%uzspst(this%nwavst(icell), icell) = &
1176  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1177  this%thtr(icell), this%eps(icell), this%vks(icell))
1178  end if
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 1420 of file UzfCellGroup.f90.

1421  ! -- Return
1422  real(DP) :: unsat_stor
1423  ! -- dummy
1424  class(UzfCellGroupType) :: this
1425  integer(I4B), intent(in) :: icell
1426  real(DP), intent(inout) :: d1
1427  ! -- local
1428  real(DP) :: fm
1429  integer(I4B) :: j, k, nwavm1, jj
1430  !
1431  fm = dzero
1432  j = this%nwavst(icell) + 1
1433  k = this%nwavst(icell)
1434  nwavm1 = k - 1
1435  if (d1 > this%uzdpst(1, icell)) d1 = this%uzdpst(1, icell)
1436  !
1437  ! -- find deepest wave above depth d1, counter held as j
1438  do while (k > 0)
1439  if (this%uzdpst(k, icell) - d1 < -dem30) j = k
1440  k = k - 1
1441  end do
1442  if (j > this%nwavst(icell)) then
1443  fm = fm + (this%uzthst(this%nwavst(icell), icell) - this%thtr(icell)) * d1
1444  elseif (this%nwavst(icell) > 1) then
1445  if (j > 1) then
1446  fm = fm + (this%uzthst(j - 1, icell) - this%thtr(icell)) &
1447  * (d1 - this%uzdpst(j, icell))
1448  end if
1449  do jj = j, nwavm1
1450  fm = fm + (this%uzthst(jj, icell) - this%thtr(icell)) &
1451  * (this%uzdpst(jj, icell) &
1452  - this%uzdpst(jj + 1, icell))
1453  end do
1454  fm = fm + (this%uzthst(this%nwavst(icell), icell) - this%thtr(icell)) &
1455  * (this%uzdpst(this%nwavst(icell), icell))
1456  else
1457  fm = fm + (this%uzthst(1, icell) - this%thtr(icell)) * d1
1458  end if
1459  unsat_stor = fm

◆ 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 1464 of file UzfCellGroup.f90.

1465  ! -- dummy
1466  class(UzfCellGroupType) :: this
1467  integer(I4B), intent(in) :: icell
1468  integer(I4B), intent(in) :: itest
1469  integer(I4B), intent(in) :: iss
1470  real(DP), intent(in) :: delt
1471  ! -- local
1472  real(DP) :: bot, depthsave, top
1473  real(DP) :: thick, thtsrinv
1474  integer(I4B) :: nwavhld, k, j
1475  !
1476  bot = this%watab(icell)
1477  top = this%celtop(icell)
1478  thick = top - bot
1479  nwavhld = this%nwavst(icell)
1480  if (itest == 1) then
1481  this%uzflst(1, icell) = dzero
1482  this%uzthst(1, icell) = this%thtr(icell)
1483  return
1484  end if
1485  if (iss == 1) then
1486  if (this%thts(icell) - this%thtr(icell) < dem7) then
1487  thtsrinv = done / dem7
1488  else
1489  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1490  end if
1491  this%totflux(icell) = this%surflux(icell) * delt
1492  this%watabold(icell) = this%watab(icell)
1493  this%uzthst(1, icell) = this%thti(icell)
1494  this%uzflst(1, icell) = &
1495  this%vks(icell) * (((this%uzthst(1, icell) - this%thtr(icell)) &
1496  * thtsrinv)**this%eps(icell))
1497  this%uzdpst(1, icell) = thick
1498  this%uzspst(1, icell) = thick
1499  this%nwavst(icell) = 1
1500  else
1501  !
1502  ! -- water table rises through waves
1503  if (this%watab(icell) - this%watabold(icell) > dem30) then
1504  depthsave = this%uzdpst(1, icell)
1505  j = 0
1506  k = this%nwavst(icell)
1507  do while (k > 0)
1508  if (this%uzdpst(k, icell) - thick < -dem30) j = k
1509  k = k - 1
1510  end do
1511  this%uzdpst(1, icell) = thick
1512  if (j > 1) then
1513  this%uzspst(1, icell) = dzero
1514  this%nwavst(icell) = this%nwavst(icell) - j + 2
1515  this%uzthst(1, icell) = this%uzthst(j - 1, icell)
1516  this%uzflst(1, icell) = this%uzflst(j - 1, icell)
1517  if (j > 2) call this%wave_shift(this, icell, icell, j - 2, 2, &
1518  nwavhld - (j - 2), 1)
1519  elseif (j == 0) then
1520  this%uzspst(1, icell) = dzero
1521  this%uzthst(1, icell) = this%uzthst(this%nwavst(icell), icell)
1522  this%uzflst(1, icell) = this%uzflst(this%nwavst(icell), icell)
1523  this%nwavst(icell) = 1
1524  end if
1525  end if
1526  !
1527  ! -- calculate new unsat. storage
1528  if (thick <= dzero) then
1529  this%uzspst(1, icell) = dzero
1530  this%nwavst(icell) = 1
1531  this%uzthst(1, icell) = this%thtr(icell)
1532  this%uzflst(1, icell) = dzero
1533  end if
1534  this%watabold(icell) = this%watab(icell)
1535  end if

◆ uz_rise()

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

Definition at line 848 of file UzfCellGroup.f90.

849  ! -- dummy
850  class(UzfCellGroupType) :: this
851  integer(I4B), intent(in) :: icell
852  real(DP), intent(inout) :: totfluxtot
853  ! -- local
854  real(DP) :: fm1, fm2, d1
855  !
856  ! -- additional recharge from a rising water table
857  if (this%watab(icell) - this%watabold(icell) > dem30) then
858  d1 = this%celtop(icell) - this%watabold(icell)
859  fm1 = this%unsat_stor(icell, d1)
860  d1 = this%celtop(icell) - this%watab(icell)
861  fm2 = this%unsat_stor(icell, d1)
862  totfluxtot = totfluxtot + (fm1 - fm2)
863  end if

◆ 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 1540 of file UzfCellGroup.f90.

1541  ! -- dummy
1542  class(UzfCellGroupType) :: this
1543  integer(I4B), intent(in) :: icell
1544  real(DP), intent(in) :: delt
1545  integer(I4B), intent(in) :: ietflag
1546  integer(I4B), intent(inout) :: ierr
1547  ! -- local
1548  type(UzfCellGroupType) :: uzfktemp
1549  real(DP) :: diff
1550  real(DP) :: thetaout
1551  real(DP) :: fm
1552  real(DP) :: st
1553  real(DP) :: thtsrinv
1554  real(DP) :: epsfksthts
1555  real(DP) :: fmp
1556  real(DP) :: fktho
1557  real(DP) :: theta1
1558  real(DP) :: theta2
1559  real(DP) :: flux1
1560  real(DP) :: flux2
1561  real(DP) :: hcap
1562  real(DP) :: factor
1563  real(DP) :: tho
1564  real(DP) :: depth
1565  real(DP) :: extwc1
1566  real(DP) :: petsub
1567  integer(I4B) :: i
1568  integer(I4B) :: j
1569  integer(I4B) :: jhold
1570  integer(I4B) :: jk
1571  integer(I4B) :: kj
1572  integer(I4B) :: kk
1573  integer(I4B) :: numadd
1574  integer(I4B) :: k
1575  integer(I4B) :: nwv
1576  integer(I4B) :: itest
1577  !
1578  ! -- initialize
1579  this%etact(icell) = dzero
1580  if (this%extdpuz(icell) < dem7) return
1581  petsub = this%rootact(icell) * this%pet(icell) * &
1582  this%extdpuz(icell) / this%extdp(icell)
1583  thetaout = delt * petsub / this%extdp(icell)
1584  if (ietflag == 1) thetaout = delt * this%pet(icell) / this%extdp(icell)
1585  if (thetaout < dem10) return
1586  depth = this%uzdpst(1, icell)
1587  st = this%unsat_stor(icell, depth)
1588  if (st < dem4) return
1589  !
1590  ! -- allocate temporary wave storage.
1591  nwv = this%nwavst(icell)
1592  itest = 0
1593  call uzfktemp%init(1, nwv)
1594  !
1595  ! store original wave characteristics
1596  call uzfktemp%wave_shift(this, 1, icell, 0, 1, nwv, 1)
1597  factor = done
1598  this%etact(icell) = dzero
1599  if (this%thts(icell) - this%thtr(icell) < dem7) then
1600  thtsrinv = 1.0 / dem7
1601  else
1602  thtsrinv = done / (this%thts(icell) - this%thtr(icell))
1603  end if
1604  epsfksthts = this%eps(icell) * this%vks(icell) * thtsrinv
1605  this%etact(icell) = dzero
1606  fmp = dzero
1607  extwc1 = this%extwc(icell) - this%thtr(icell)
1608  if (extwc1 < dem6) extwc1 = dem7
1609  numadd = 0
1610  fm = st
1611  k = 0
1612  !
1613  ! -- loop for reducing aet to pet when et is head dependent
1614  do while (itest == 0)
1615  k = k + 1
1616  if (k > 1 .AND. abs(fmp - petsub) > dem5 * petsub) then
1617  factor = factor / (fm / petsub)
1618  end if
1619  !
1620  ! -- one wave shallower than extdp
1621  if (this%nwavst(icell) == 1 .AND. &
1622  this%uzdpst(1, icell) <= this%extdpuz(icell)) then
1623  if (ietflag == 2) then
1624  tho = this%uzthst(1, icell)
1625  fktho = this%uzflst(1, icell)
1626  hcap = this%caph(icell, tho)
1627  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1628  end if
1629  if ((this%uzthst(1, icell) - thetaout) > this%thtr(icell) + extwc1) then
1630  this%uzthst(1, icell) = this%uzthst(1, icell) - thetaout
1631  this%uzflst(1, icell) = &
1632  this%vks(icell) * (((this%uzthst(1, icell) - &
1633  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1634  else if (this%uzthst(1, icell) > this%thtr(icell) + extwc1) then
1635  this%uzthst(1, icell) = this%thtr(icell) + extwc1
1636  this%uzflst(1, icell) = &
1637  this%vks(icell) * (((this%uzthst(1, icell) - &
1638  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1639  end if
1640  !
1641  ! -- all waves shallower than extinction depth
1642  else if (this%nwavst(icell) > 1 .AND. &
1643  this%uzdpst(this%nwavst(icell), icell) > this%extdpuz(icell)) then
1644  if (ietflag == 2) then
1645  tho = this%uzthst(this%nwavst(icell), icell)
1646  fktho = this%uzflst(this%nwavst(icell), icell)
1647  hcap = this%caph(icell, tho)
1648  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1649  end if
1650  if (this%nwavst(icell) + 1 > this%nwav(icell)) then
1651  !
1652  ! -- too many waves error
1653  ierr = 1
1654  goto 500
1655  end if
1656  if (this%uzthst(this%nwavst(icell), icell) - thetaout > &
1657  this%thtr(icell) + extwc1) then
1658  this%uzthst(this%nwavst(icell) + 1, icell) = &
1659  this%uzthst(this%nwavst(icell), icell) - thetaout
1660  numadd = 1
1661  else if (this%uzthst(this%nwavst(icell), icell) > &
1662  this%thtr(icell) + extwc1) then
1663  this%uzthst(this%nwavst(icell) + 1, icell) = this%thtr(icell) + extwc1
1664  numadd = 1
1665  end if
1666  if (numadd == 1) then
1667  this%uzflst(this%nwavst(icell) + 1, icell) = &
1668  this%vks(icell) * &
1669  (((this%uzthst(this%nwavst(icell) + 1, icell) - &
1670  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1671  theta2 = this%uzthst(this%nwavst(icell) + 1, icell)
1672  flux2 = this%uzflst(this%nwavst(icell) + 1, icell)
1673  flux1 = this%uzflst(this%nwavst(icell), icell)
1674  theta1 = this%uzthst(this%nwavst(icell), icell)
1675  this%uzspst(this%nwavst(icell) + 1, icell) = &
1676  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1677  this%thtr(icell), this%eps(icell), this%vks(icell))
1678  this%uzdpst(this%nwavst(icell) + 1, icell) = this%extdpuz(icell)
1679  this%nwavst(icell) = this%nwavst(icell) + 1
1680  if (this%nwavst(icell) > this%nwav(icell)) then
1681  !
1682  ! -- too many waves error, deallocate temp arrays and return
1683  ierr = 1
1684  goto 500
1685  end if
1686  else
1687  numadd = 0
1688  end if
1689  !
1690  ! -- one wave below extinction depth
1691  else if (this%nwavst(icell) == 1) then
1692  if (this%nwavst(icell) + 1 > this%nwav(icell)) then
1693  !
1694  ! -- too many waves error
1695  ierr = 1
1696  goto 500
1697  end if
1698  if (ietflag == 2) then
1699  tho = this%uzthst(1, icell)
1700  fktho = this%uzflst(1, icell)
1701  hcap = this%caph(icell, tho)
1702  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1703  end if
1704  if ((this%uzthst(1, icell) - thetaout) > this%thtr(icell) + extwc1) then
1705  if (thetaout > dem30) then
1706  this%uzthst(2, icell) = this%uzthst(1, icell) - thetaout
1707  this%uzflst(2, icell) = &
1708  this%vks(icell) * (((this%uzthst(2, icell) - this%thtr(icell)) * &
1709  thtsrinv)**this%eps(icell))
1710  this%uzdpst(2, icell) = this%extdpuz(icell)
1711  theta2 = this%uzthst(2, icell)
1712  flux2 = this%uzflst(2, icell)
1713  flux1 = this%uzflst(1, icell)
1714  theta1 = this%uzthst(1, icell)
1715  this%uzspst(2, icell) = &
1716  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1717  this%thtr(icell), this%eps(icell), this%vks(icell))
1718  this%nwavst(icell) = this%nwavst(icell) + 1
1719  if (this%nwavst(icell) > this%nwav(icell)) then
1720  !
1721  ! -- too many waves error
1722  ierr = 1
1723  goto 500
1724  end if
1725  end if
1726  else if (this%uzthst(1, icell) > this%thtr(icell) + extwc1) then
1727  if (thetaout > dem30) then
1728  this%uzthst(2, icell) = this%thtr(icell) + extwc1
1729  this%uzflst(2, icell) = &
1730  this%vks(icell) * (((this%uzthst(2, icell) - &
1731  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1732  this%uzdpst(2, icell) = this%extdpuz(icell)
1733  theta2 = this%uzthst(2, icell)
1734  flux2 = this%uzflst(2, icell)
1735  flux1 = this%uzflst(1, icell)
1736  theta1 = this%uzthst(1, icell)
1737  this%uzspst(2, icell) = &
1738  leadspeed(theta1, theta2, flux1, flux2, this%thts(icell), &
1739  this%thtr(icell), this%eps(icell), this%vks(icell))
1740  this%nwavst(icell) = this%nwavst(icell) + 1
1741  if (this%nwavst(icell) > this%nwav(icell)) then
1742  !
1743  ! -- too many waves error
1744  ierr = 1
1745  goto 500
1746  end if
1747  end if
1748  end if
1749  else
1750  !
1751  ! -- extinction depth splits waves
1752  if (this%uzdpst(1, icell) - this%extdpuz(icell) > dem7) then
1753  j = 2
1754  jk = 0
1755  !
1756  ! -- locate extinction depth between waves
1757  do while (jk == 0)
1758  diff = this%uzdpst(j, icell) - this%extdpuz(icell)
1759  if (diff > dzero) then
1760  j = j + 1
1761  else
1762  jk = 1
1763  end if
1764  end do
1765  kk = j
1766  if (this%uzthst(j, icell) > this%thtr(icell) + extwc1) then
1767  !
1768  ! -- create a wave at extinction depth
1769  if (abs(diff) > dem5) then
1770  if (this%nwavst(icell) + 1 > this%nwav(icell)) then
1771  !
1772  ! -- too many waves error
1773  ierr = 1
1774  goto 500
1775  end if
1776  call this%wave_shift(this, icell, icell, -1, &
1777  this%nwavst(icell) + 1, j, -1)
1778  this%uzdpst(j, icell) = this%extdpuz(icell)
1779  this%nwavst(icell) = this%nwavst(icell) + 1
1780  if (this%nwavst(icell) > this%nwav(icell)) then
1781  !
1782  ! -- too many waves error
1783  ierr = 1
1784  goto 500
1785  end if
1786  end if
1787  kk = j
1788  else
1789  jhold = this%nwavst(icell)
1790  i = j + 1
1791  do while (i < this%nwavst(icell))
1792  if (this%uzthst(i, icell) > this%thtr(icell) + extwc1) then
1793  jhold = i
1794  i = this%nwavst(icell) + 1
1795  end if
1796  i = i + 1
1797  end do
1798  j = jhold
1799  kk = jhold
1800  end if
1801  else
1802  kk = 1
1803  end if
1804  !
1805  ! -- all waves above extinction depth
1806  do while (kk <= this%nwavst(icell))
1807  if (ietflag == 2) then
1808  tho = this%uzthst(kk, icell)
1809  fktho = this%uzflst(kk, icell)
1810  hcap = this%caph(icell, tho)
1811  thetaout = this%rate_et_z(icell, factor, fktho, hcap)
1812  end if
1813  if (this%uzthst(kk, icell) > this%thtr(icell) + extwc1) then
1814  if (this%uzthst(kk, icell) - thetaout > &
1815  this%thtr(icell) + extwc1) then
1816  this%uzthst(kk, icell) = this%uzthst(kk, icell) - thetaout
1817  else if (this%uzthst(kk, icell) > this%thtr(icell) + extwc1) then
1818  this%uzthst(kk, icell) = this%thtr(icell) + extwc1
1819  end if
1820  if (kk == 1) then
1821  this%uzflst(kk, icell) = &
1822  this%vks(icell) * &
1823  (((this%uzthst(kk, icell) - &
1824  this%thtr(icell)) * thtsrinv)**this%eps(icell))
1825  end if
1826  if (kk > 1) then
1827  flux1 = &
1828  this%vks(icell) * ((this%uzthst(kk - 1, icell) - &
1829  this%thtr(icell)) * thtsrinv)**this%eps(icell)
1830  flux2 = &
1831  this%vks(icell) * ((this%uzthst(kk, icell) - &
1832  this%thtr(icell)) * thtsrinv)**this%eps(icell)
1833  this%uzflst(kk, icell) = flux2
1834  theta2 = this%uzthst(kk, icell)
1835  theta1 = this%uzthst(kk - 1, icell)
1836  this%uzspst(kk, icell) = leadspeed(theta1, theta2, flux1, flux2, &
1837  this%thts(icell), &
1838  this%thtr(icell), &
1839  this%eps(icell), this%vks(icell))
1840  end if
1841  end if
1842  kk = kk + 1
1843  end do
1844  end if
1845  !
1846  ! -- calculate aet
1847  kj = 1
1848  do while (kj <= this%nwavst(icell) - 1)
1849  if (abs(this%uzthst(kj, icell) - this%uzthst(kj + 1, icell)) < dem6) then
1850  call this%wave_shift(this, icell, icell, 1, kj + 1, &
1851  this%nwavst(icell) - 1, 1)
1852  kj = kj - 1
1853  this%nwavst(icell) = this%nwavst(icell) - 1
1854  end if
1855  kj = kj + 1
1856  end do
1857  depth = this%uzdpst(1, icell)
1858  fm = this%unsat_stor(icell, depth)
1859  this%etact(icell) = st - fm
1860  fm = this%etact(icell) / delt
1861  if (this%etact(icell) < dzero) then
1862  call this%wave_shift(uzfktemp, icell, 1, 0, 1, nwv, 1)
1863  this%nwavst(icell) = nwv
1864  this%etact(icell) = dzero
1865  elseif (petsub - fm < -dem15 .AND. ietflag == 2) then
1866  !
1867  ! -- aet greater than pet, reset and try again
1868  call this%wave_shift(uzfktemp, icell, 1, 0, 1, nwv, 1)
1869  this%nwavst(icell) = nwv
1870  this%etact(icell) = dzero
1871  else
1872  itest = 1
1873  end if
1874  !
1875  ! -- end aet-pet loop for head dependent et
1876  fmp = fm
1877  if (k > 100) then
1878  itest = 1
1879  elseif (ietflag < 2) then
1880  fmp = petsub
1881  itest = 1
1882  end if
1883  end do
1884 500 continue
1885  !
1886  ! -- deallocate temporary worker
1887  call uzfktemp%dealloc()
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 980 of file UzfCellGroup.f90.

981  ! -- dummy
982  class(UzfCellGroupType) :: this
983  real(DP), intent(inout) :: thickold
984  real(DP), intent(inout) :: thick
985  real(DP), intent(in) :: delt
986  integer(I4B), intent(in) :: ietflag
987  integer(I4B), intent(in) :: icell
988  integer(I4B), intent(inout) :: ierr
989  ! -- local
990  real(DP) :: ffcheck, time, feps1, feps2
991  real(DP) :: thetadif, thetab, fluxb, oldsflx
992  integer(I4B) :: itrailflg, itester
993  !
994  time = dzero
995  this%totflux(icell) = dzero
996  itrailflg = 0
997  oldsflx = this%uzflst(this%nwavst(icell), icell)
998  call factors(feps1, feps2)
999  !
1000  ! -- check for falling or rising water table
1001  if ((thick - thickold) > feps1) then
1002  thetadif = abs(this%uzthst(1, icell) - this%thtr(icell))
1003  if (thetadif > dem6) then
1004  call this%wave_shift(this, icell, icell, -1, &
1005  this%nwavst(icell) + 1, 2, -1)
1006  if (this%uzdpst(2, icell) < dem30) &
1007  this%uzdpst(2, icell) = (this%ntrail(icell) + dtwo) * dem6
1008  if (this%uzthst(2, icell) > this%thtr(icell)) then
1009  this%uzspst(2, icell) = this%uzflst(2, icell) / &
1010  (this%uzthst(2, icell) - this%thtr(icell))
1011  else
1012  this%uzspst(2, icell) = dzero
1013  end if
1014  this%uzthst(1, icell) = this%thtr(icell)
1015  this%uzflst(1, icell) = dzero
1016  this%uzspst(1, icell) = dzero
1017  this%uzdpst(1, icell) = thick
1018  this%nwavst(icell) = this%nwavst(icell) + 1
1019  if (this%nwavst(icell) >= this%nwav(icell)) then
1020  ! -- too many waves error
1021  ierr = 1
1022  return
1023  end if
1024  else
1025  this%uzdpst(1, icell) = thick
1026  end if
1027  end if
1028  thetab = this%uzthst(1, icell)
1029  fluxb = this%uzflst(1, icell)
1030  this%totflux(icell) = dzero
1031  itester = 0
1032  ffcheck = (this%surflux(icell) - this%uzflst(this%nwavst(icell), icell))
1033  !
1034  ! -- increase new waves in infiltration changes
1035  if (ffcheck > feps2 .OR. ffcheck < -feps2) then
1036  this%nwavst(icell) = this%nwavst(icell) + 1
1037  if (this%nwavst(icell) >= this%nwav(icell)) then
1038  !
1039  ! -- too many waves error
1040  ierr = 1
1041  return
1042  end if
1043  else if (this%nwavst(icell) == 1) then
1044  itester = 1
1045  end if
1046  if (this%nwavst(icell) > 1) then
1047  if (ffcheck < -feps2) then
1048  call this%trailwav(icell, ierr)
1049  if (ierr > 0) return
1050  itrailflg = 1
1051  end if
1052  call this%leadwav(time, itester, itrailflg, thetab, fluxb, ffcheck, &
1053  feps2, delt, icell)
1054  end if
1055  if (itester == 1) then
1056  this%totflux(icell) = this%totflux(icell) + &
1057  (delt - time) * this%uzflst(1, icell)
1058  time = dzero
1059  itester = 0
1060  end if
1061  !
1062  ! -- simulate et
1063  if (ietflag > 0) call this%uzet(icell, delt, ietflag, ierr)
1064  if (ierr > 0) return
Here is the call graph for this function:

◆ wave_shift()

subroutine uzfcellgroupmodule::wave_shift ( class(uzfcellgrouptype)  this,
type(uzfcellgrouptype)  this2,
integer(i4b), intent(in)  icell,
integer(i4b), intent(in)  icell2,
integer(i4b), intent(in)  shft,
integer(i4b), intent(in)  strt,
integer(i4b), intent(in)  stp,
integer(i4b), intent(in)  cntr 
)
private

Definition at line 955 of file UzfCellGroup.f90.

956  ! -- dummy
957  class(UzfCellGroupType) :: this
958  type(UzfCellGroupType) :: this2
959  integer(I4B), intent(in) :: icell
960  integer(I4B), intent(in) :: icell2
961  integer(I4B), intent(in) :: shft
962  integer(I4B), intent(in) :: strt
963  integer(I4B), intent(in) :: stp
964  integer(I4B), intent(in) :: cntr
965  ! -- local
966  integer(I4B) :: j
967  !
968  ! -- copy waves from one uzf cell group to another
969  do j = strt, stp, cntr
970  this%uzthst(j, icell) = this2%uzthst(j + shft, icell2)
971  this%uzdpst(j, icell) = this2%uzdpst(j + shft, icell2)
972  this%uzflst(j, icell) = this2%uzflst(j + shft, icell2)
973  this%uzspst(j, icell) = this2%uzspst(j + shft, icell2)
974  end do
975  this%nwavst(icell) = this2%nwavst(icell2)