18 integer(I4B),
pointer :: iadvwt => null()
19 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
21 real(dp),
pointer :: eqnsclfac => null()
46 subroutine adv_cr(advobj, name_model, inunit, iout, fmi, eqnsclfac)
49 character(len=*),
intent(in) :: name_model
50 integer(I4B),
intent(in) :: inunit
51 integer(I4B),
intent(in) :: iout
53 real(dp),
intent(in),
pointer :: eqnsclfac
59 call advobj%set_names(1, name_model,
'ADV',
'ADV')
62 call advobj%allocate_scalars()
65 advobj%inunit = inunit
68 advobj%eqnsclfac => eqnsclfac
83 character(len=*),
parameter :: fmtadv = &
84 "(1x,/1x,'ADV-- ADVECTION PACKAGE, VERSION 1, 8/25/2017', &
85 &' INPUT READ FROM UNIT ', i0, //)"
88 if (.not.
present(adv_options))
then
92 call this%parser%Initialize(this%inunit, this%iout)
95 write (this%iout, fmtadv) this%inunit
98 call this%read_options()
102 this%iadvwt = adv_options%iAdvScheme
118 integer(I4B),
dimension(:),
pointer,
contiguous,
intent(in) :: ibound
124 this%ibound => ibound
134 subroutine adv_fc(this, nodes, matrix_sln, idxglo, cnew, rhs)
138 integer,
intent(in) :: nodes
140 integer(I4B),
intent(in),
dimension(:) :: idxglo
141 real(DP),
intent(in),
dimension(:) :: cnew
142 real(DP),
dimension(:),
intent(inout) :: rhs
144 integer(I4B) :: n, m, idiag, ipos
145 real(DP) :: omega, qnm
150 if (this%ibound(n) == 0) cycle
151 idiag = this%dis%con%ia(n)
152 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
153 if (this%dis%con%mask(ipos) == 0) cycle
154 m = this%dis%con%ja(ipos)
155 if (this%ibound(m) == 0) cycle
156 qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
157 omega = this%adv_weight(this%iadvwt, ipos, n, m, qnm)
158 call matrix_sln%add_value_pos(idxglo(ipos), qnm * (
done - omega))
159 call matrix_sln%add_value_pos(idxglo(idiag), qnm * omega)
164 if (this%iadvwt == 2)
then
166 if (this%ibound(n) == 0) cycle
167 call this%advtvd(n, cnew, rhs)
184 integer(I4B),
intent(in) :: n
185 real(DP),
dimension(:),
intent(in) :: cnew
186 real(DP),
dimension(:),
intent(inout) :: rhs
189 integer(I4B) :: m, ipos
192 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
193 if (this%dis%con%mask(ipos) == 0) cycle
194 m = this%dis%con%ja(ipos)
195 if (m > n .and. this%ibound(m) /= 0)
then
196 qtvd = this%advqtvd(n, m, ipos, cnew)
197 rhs(n) = rhs(n) - qtvd
198 rhs(m) = rhs(m) + qtvd
211 function advqtvd(this, n, m, iposnm, cnew)
result(qtvd)
218 integer(I4B),
intent(in) :: n
219 integer(I4B),
intent(in) :: m
220 integer(I4B),
intent(in) :: iposnm
221 real(dp),
dimension(:),
intent(in) :: cnew
223 integer(I4B) :: ipos, isympos, iup, idn, i2up, j
224 real(dp) :: qnm, qmax, qupj, elupdn, elup2up
225 real(dp) :: smooth, cdiff, alimiter
231 isympos = this%dis%con%jas(iposnm)
232 qnm = this%fmi%gwfflowja(iposnm)
233 if (qnm > dzero)
then
241 elupdn = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos)
246 do ipos = this%dis%con%ia(iup) + 1, this%dis%con%ia(iup + 1) - 1
247 j = this%dis%con%ja(ipos)
248 if (this%ibound(j) == 0) cycle
249 qupj = this%fmi%gwfflowja(ipos)
250 isympos = this%dis%con%jas(ipos)
251 if (qupj > qmax)
then
254 elup2up = this%dis%con%cl1(isympos) + this%dis%con%cl2(isympos)
261 cdiff = abs(cnew(idn) - cnew(iup))
262 if (cdiff >
dprec)
then
263 smooth = (cnew(iup) - cnew(i2up)) / elup2up * &
264 elupdn / (cnew(idn) - cnew(iup))
266 if (smooth > dzero)
then
267 alimiter = dtwo * smooth / (done + smooth)
268 qtvd = dhalf * alimiter * qnm * (cnew(idn) - cnew(iup))
269 qtvd = qtvd * this%eqnsclfac
283 real(DP),
intent(in),
dimension(:) :: cnew
284 real(DP),
intent(inout),
dimension(:) :: flowja
286 integer(I4B) :: nodes
287 integer(I4B) :: n, m, idiag, ipos
288 real(DP) :: omega, qnm
292 nodes = this%dis%nodes
294 if (this%ibound(n) == 0) cycle
295 idiag = this%dis%con%ia(n)
296 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
297 m = this%dis%con%ja(ipos)
298 if (this%ibound(m) == 0) cycle
299 qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
300 omega = this%adv_weight(this%iadvwt, ipos, n, m, qnm)
301 flowja(ipos) = flowja(ipos) + qnm * omega * cnew(n) + &
302 qnm * (
done - omega) * cnew(m)
307 if (this%iadvwt == 2)
call this%advtvd_bd(cnew, flowja)
318 real(DP),
dimension(:),
intent(in) :: cnew
319 real(DP),
dimension(:),
intent(inout) :: flowja
321 real(DP) :: qtvd, qnm
322 integer(I4B) :: nodes, n, m, ipos
324 nodes = this%dis%nodes
326 if (this%ibound(n) == 0) cycle
327 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
328 m = this%dis%con%ja(ipos)
329 if (this%ibound(m) /= 0)
then
330 qnm = this%fmi%gwfflowja(ipos)
331 qtvd = this%advqtvd(n, m, ipos, cnew)
332 flowja(ipos) = flowja(ipos) + qtvd
350 if (this%inunit > 0)
then
354 this%ibound => null()
360 call this%NumericalPackageType%da()
378 call this%NumericalPackageType%allocate_scalars()
381 call mem_allocate(this%iadvwt,
'IADVWT', this%memoryPath)
404 character(len=LINELENGTH) :: errmsg, keyword
406 logical :: isfound, endOfBlock
408 character(len=*),
parameter :: fmtiadvwt = &
409 &
"(4x,'ADVECTION WEIGHTING SCHEME HAS BEEN SET TO: ', a)"
412 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, blockrequired=.false., &
413 supportopenclose=.true.)
417 write (this%iout,
'(1x,a)')
'PROCESSING ADVECTION OPTIONS'
419 call this%parser%GetNextLine(endofblock)
421 call this%parser%GetStringCaps(keyword)
422 select case (keyword)
424 call this%parser%GetStringCaps(keyword)
425 select case (keyword)
428 write (this%iout, fmtiadvwt)
'UPSTREAM'
431 write (this%iout, fmtiadvwt)
'CENTRAL'
434 write (this%iout, fmtiadvwt)
'TVD'
436 write (errmsg,
'(a, a)') &
437 'Unknown scheme: ', trim(keyword)
439 write (errmsg,
'(a, a)') &
440 'Scheme must be "UPSTREAM", "CENTRAL" or "TVD"'
442 call this%parser%StoreErrorUnit()
445 write (errmsg,
'(a,a)')
'Unknown ADVECTION option: ', &
450 write (this%iout,
'(1x,a)')
'END OF ADVECTION OPTIONS'
461 function adv_weight(this, iadvwt, ipos, n, m, qnm)
result(omega)
466 integer,
intent(in) :: iadvwt
467 integer,
intent(in) :: ipos
468 integer,
intent(in) :: n
469 integer,
intent(in) :: m
470 real(dp),
intent(in) :: qnm
478 if (this%dis%con%ihc(this%dis%con%jas(ipos)) == 0)
then
480 lnm =
dhalf * (this%dis%top(n) - this%dis%bot(n))
481 lmn =
dhalf * (this%dis%top(m) - this%dis%bot(m))
484 lnm = this%dis%con%cl1(this%dis%con%jas(ipos))
485 lmn = this%dis%con%cl2(this%dis%con%jas(ipos))
487 omega = lmn / (lnm + lmn)
490 if (qnm >
dzero)
then
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dzero
real constant zero
real(dp), parameter dprec
real constant machine precision
real(dp), parameter dtwo
real constant 2
real(dp), parameter done
real constant 1
This module defines variable data types.
This module contains the base numerical package type.
This module contains simulation methods.
subroutine, public store_error(msg, terminate)
Store an error message.
subroutine advtvd_bd(this, cnew, flowja)
Add TVD contribution to flowja.
real(dp) function advqtvd(this, n, m, iposnm, cnew)
Calculate TVD.
subroutine adv_df(this, adv_options)
Define ADV object.
subroutine adv_fc(this, nodes, matrix_sln, idxglo, cnew, rhs)
Fill coefficient method for ADV package.
subroutine adv_cq(this, cnew, flowja)
Calculate advection contribution to flowja.
subroutine adv_da(this)
Deallocate memory.
subroutine advtvd(this, n, cnew, rhs)
Calculate TVD.
subroutine adv_ar(this, dis, ibound)
Allocate and read method for package.
real(dp) function adv_weight(this, iadvwt, ipos, n, m, qnm)
@ brief Advection weight
subroutine read_options(this)
Read options.
subroutine, public adv_cr(advobj, name_model, inunit, iout, fmi, eqnsclfac)
@ brief Create a new ADV object
subroutine allocate_scalars(this)
Allocate scalars specific to the streamflow energy transport (SFE) package.