215 class(GwfHfbType) :: this
216 integer(I4B) :: kiter
217 class(MatrixBaseType),
pointer :: matrix_sln
218 integer(I4B),
intent(in),
dimension(:) :: idxglo
219 real(DP),
intent(inout),
dimension(:) :: rhs
220 real(DP),
intent(inout),
dimension(:) :: hnew
222 integer(I4B) :: nodes, nja
223 integer(I4B) :: ihfb, n, m
225 integer(I4B) :: idiag, isymcon
226 integer(I4B) :: ixt3d
227 real(DP) :: cond, condhfb, aterm
228 real(DP) :: fawidth, faheight
229 real(DP) :: topn, topm, botn, botm
230 real(DP) :: viscratio
234 nodes = this%dis%nodes
235 nja = this%dis%con%nja
236 if (
associated(this%xt3d%ixt3d))
then
237 ixt3d = this%xt3d%ixt3d
244 do ihfb = 1, this%nhfb
245 n = min(this%noden(ihfb), this%nodem(ihfb))
246 m = max(this%noden(ihfb), this%nodem(ihfb))
248 if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle
250 if (this%ivsc /= 0)
then
251 call this%vsc%get_visc_ratio(n, m, hnew(n), hnew(m), viscratio)
254 if (this%hydchr(ihfb) >
dzero)
then
255 if (this%inewton == 0)
then
256 ipos = this%idxloc(ihfb)
261 if (this%icelltype(n) == 1)
then
262 if (hnew(n) < topn) topn = hnew(n)
264 if (this%icelltype(m) == 1)
then
265 if (hnew(m) < topm) topm = hnew(m)
267 if (this%ihc(this%jas(ipos)) == 2)
then
268 faheight = min(topn, topm) - max(botn, botm)
270 faheight =
dhalf * ((topn - botn) + (topm - botm))
272 fawidth = this%hwva(this%jas(ipos))
273 condhfb = this%hydchr(ihfb) * viscratio * &
276 condhfb = this%hydchr(ihfb) * viscratio
279 condhfb = this%hydchr(ihfb)
282 call this%xt3d%xt3d_fhfb(kiter, nodes, nja, matrix_sln, idxglo, &
283 rhs, hnew, n, m, condhfb)
289 if (this%inewton == 0)
then
290 do ihfb = 1, this%nhfb
291 ipos = this%idxloc(ihfb)
292 aterm = matrix_sln%get_value_pos(idxglo(ipos))
295 if (this%ibound(n) == 0 .or. this%ibound(m) == 0) cycle
297 if (this%ivsc /= 0)
then
298 call this%vsc%get_visc_ratio(n, m, hnew(n), hnew(m), viscratio)
301 if (this%icelltype(n) == 1 .or. this%icelltype(m) == 1 .or. &
309 if (this%icelltype(n) == 1)
then
310 if (hnew(n) < topn) topn = hnew(n)
312 if (this%icelltype(m) == 1)
then
313 if (hnew(m) < topm) topm = hnew(m)
315 if (this%ihc(this%jas(ipos)) == 2)
then
316 faheight = min(topn, topm) - max(botn, botm)
318 faheight =
dhalf * ((topn - botn) + (topm - botm))
320 if (this%hydchr(ihfb) >
dzero)
then
321 fawidth = this%hwva(this%jas(ipos))
322 condhfb = this%hydchr(ihfb) * viscratio * &
324 cond = aterm * condhfb / (aterm + condhfb)
326 cond = -aterm * this%hydchr(ihfb)
330 this%condsav(ihfb) = cond
334 call matrix_sln%add_value_pos(idxglo(idiag), aterm - cond)
335 call matrix_sln%set_value_pos(idxglo(ipos), cond)
338 isymcon = this%isym(ipos)
340 call matrix_sln%add_value_pos(idxglo(idiag), aterm - cond)
341 call matrix_sln%set_value_pos(idxglo(isymcon), cond)