MODFLOW 6  version 6.7.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 700 of file UzfCellGroup.f90.

701  ! -- dummy
702  class(UzfCellGroupType) :: this
703  integer(I4B), intent(in) :: icell
704  integer(I4B), intent(in) :: jbelow
705  real(DP), intent(inout) :: trhs
706  real(DP), intent(inout) :: thcof
707  real(DP), intent(inout) :: deriv
708  real(DP), intent(in) :: delt
709  real(DP), intent(in) :: hgwf
710  ! -- local
711  real(DP) :: fcheck
712  real(DP) :: x, scale, range
713  !
714  ! -- initialize
715  range = dem5
716  deriv = dzero
717  thcof = dzero
718  trhs = this%uzfarea(icell) * this%totflux(icell) / delt
719  if (this%totflux(icell) < dem14) return
720  scale = done
721  !
722  ! -- smoothly reduce flow between cells when head close to cell top
723  x = hgwf - (this%celbot(icell) - range)
724  call sscurve(x, range, deriv, scale)
725  deriv = this%uzfarea(icell) * deriv * this%totflux(icell) / delt
726  this%finf(jbelow) = (done - scale) * this%totflux(icell) / delt
727  fcheck = this%finf(jbelow) - this%vks(jbelow)
728  !
729  ! -- reduce flow between cells when vks is too small
730  if (fcheck < dem14) fcheck = dzero
731  this%finf(jbelow) = this%finf(jbelow) - fcheck
732  this%surfluxbelow(icell) = this%finf(jbelow)
733  this%totflux(icell) = scale * this%totflux(icell) + fcheck * delt
734  trhs = this%uzfarea(icell) * this%totflux(icell) / delt
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 1869 of file UzfCellGroup.f90.

1870  ! -- dummy
1871  class(UzfCellGroupType) :: this
1872  integer(I4B), intent(in) :: icell
1873  real(DP), intent(in) :: tho
1874  ! -- local
1875  real(DP) :: caph, lambda, star
1876  !
1877  caph = -dem6
1878  star = (tho - this%thtr(icell)) / (this%thts(icell) - this%thtr(icell))
1879  if (star < dem15) star = dem15
1880  lambda = dtwo / (this%eps(icell) - dthree)
1881  if (star > dem15) then
1882  if (tho - this%thts(icell) < dem15) then
1883  caph = this%ha(icell) * star**(-done / lambda)
1884  else
1885  caph = dzero
1886  end if
1887  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:122

◆ factors()

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

Definition at line 1065 of file UzfCellGroup.f90.

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

1909  ! -- dummy
1910  class(UzfCellGroupType) :: this
1911  integer(I4B), intent(in) :: icell !< uzf cell containing depth
1912  real(DP), intent(in) :: depth !< depth within the cell
1913  ! -- return
1914  real(DP) :: theta_at_depth
1915  ! -- local
1916  real(DP) :: d1
1917  real(DP) :: d2
1918  real(DP) :: f1
1919  real(DP) :: f2
1920  !
1921  if (this%watab(icell) < this%celtop(icell)) then
1922  if (this%celtop(icell) - depth > this%watab(icell)) then
1923  d1 = depth - dem3
1924  d2 = depth + dem3
1925  f1 = this%unsat_stor(icell, d1)
1926  f2 = this%unsat_stor(icell, d2)
1927  theta_at_depth = this%thtr(icell) + (f2 - f1) / (d2 - d1)
1928  else
1929  theta_at_depth = this%thts(icell)
1930  end if
1931  else
1932  theta_at_depth = this%thts(icell)
1933  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 1938 of file UzfCellGroup.f90.

1939  ! -- dummy
1940  class(UzfCellGroupType) :: this
1941  integer(I4B), intent(in) :: icell !< uzf cell containing depth
1942  ! -- return
1943  real(DP) :: watercontent
1944  ! -- local
1945  real(DP) :: top
1946  real(DP) :: bot
1947  real(DP) :: theta_r
1948  real(DP) :: thk
1949  real(DP) :: hgwf
1950  real(DP) :: fm
1951  real(DP) :: d
1952  !
1953  hgwf = this%watab(icell)
1954  top = this%celtop(icell)
1955  bot = this%celbot(icell)
1956  thk = top - max(bot, hgwf)
1957  if (thk > dzero) then
1958  theta_r = this%thtr(icell)
1959  d = thk
1960  fm = this%unsat_stor(icell, d)
1961  watercontent = fm / thk
1962  watercontent = watercontent + theta_r
1963  else
1964  watercontent = dzero
1965  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 767 of file UzfCellGroup.f90.

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

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

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

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

1892  ! -- Return
1893  real(DP) :: rate_et_z
1894  ! -- dummy
1895  class(UzfCellGroupType) :: this
1896  integer(I4B), intent(in) :: icell
1897  real(DP), intent(in) :: factor, fktho, h
1898  !
1899  rate_et_z = factor * fktho * (h - this%hroot(icell))
1900  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 739 of file UzfCellGroup.f90.

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

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

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

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

808  ! -- dummy
809  class(UzfCellGroupType) :: this
810  integer(I4B), intent(in) :: igwetflag
811  integer(I4B), intent(in) :: icell
812  real(DP), intent(in) :: hgwf
813  real(DP), intent(inout) :: trhs
814  real(DP), intent(inout) :: thcof
815  real(DP), intent(inout) :: det
816  ! -- local
817  real(DP) :: s, x, c, b, et
818  !
819  this%gwet(icell) = dzero
820  trhs = dzero
821  thcof = dzero
822  det = dzero
823  s = this%landtop(icell)
824  x = this%extdp(icell)
825  c = this%gwpet(icell)
826  b = this%celbot(icell)
827  if (b > hgwf) return
828  if (x < dem6) return
829  if (igwetflag == 1) then
830  et = etfunc_lin(s, x, c, det, trhs, thcof, hgwf, &
831  this%celtop(icell), this%celbot(icell))
832  else if (igwetflag == 2) then
833  et = etfunc_nlin(s, x, c, det, trhs, thcof, hgwf)
834  end if
835  ! this%gwet(icell) = et * this%uzfarea(icell)
836  trhs = trhs * this%uzfarea(icell)
837  thcof = thcof * this%uzfarea(icell)
838  this%gwet(icell) = trhs - (thcof * hgwf)
839  ! write(99,*)'in group', icell, this%gwet(icell)
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) return
659  call this%uz_rise(icell, totfluxtot)
660  this%totflux(icell) = totfluxtot
661  if (this%ivertcon(icell) > 0) then
662  call this%addrech(icell, jbelow, hgwf, trhsfinf, thcoffinf, &
663  derivfinf, delt)
664  end if
665  else
666  this%totflux(icell) = this%surflux(icell) * delt
667  totfluxtot = this%surflux(icell) * delt
668  end if
669  thcoffinf = dzero
670  trhsfinf = this%totflux(icell) * this%uzfarea(icell) / delt
671  if (.not. reset_state) then
672  call this%update_wav(icell, delt, issflag, 0)
673  end if
674  else
675  this%totflux(icell) = this%surflux(icell) * delt
676  totfluxtot = this%surflux(icell) * delt
677  if (.not. reset_state) then
678  call this%update_wav(icell, delt, issflag, 1)
679  end if
680  end if
681  !
682  ! -- If formulating, then these variables will be present
683  if (present(deriv)) deriv = deriv1 + deriv2 + derivfinf
684  if (present(trhs)) trhs = trhsfinf + trhsseep
685  if (present(thcof)) thcof = thcoffinf + thcofseep
686  !
687  ! -- Assign water content prior to resetting waves
688  if (present(watercontent)) then
689  watercontent = this%get_wcnew(icell)
690  end if
691  !
692  ! -- reset waves to previous state for next iteration
693  if (reset_state) then
694  call this%wave_shift(thiswork, icell, 1, 0, 1, thiswork%nwavst(1), 1)
695  end if

◆ trailwav()

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

Definition at line 1093 of file UzfCellGroup.f90.

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

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

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

845  ! -- dummy
846  class(UzfCellGroupType) :: this
847  integer(I4B), intent(in) :: icell
848  real(DP), intent(inout) :: totfluxtot
849  ! -- local
850  real(DP) :: fm1, fm2, d1
851  !
852  ! -- additional recharge from a rising water table
853  if (this%watab(icell) - this%watabold(icell) > dem30) then
854  d1 = this%celtop(icell) - this%watabold(icell)
855  fm1 = this%unsat_stor(icell, d1)
856  d1 = this%celtop(icell) - this%watab(icell)
857  fm2 = this%unsat_stor(icell, d1)
858  totfluxtot = totfluxtot + (fm1 - fm2)
859  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 1535 of file UzfCellGroup.f90.

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

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

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