15 character(len=LENMEMPATH) :: memorypath
16 integer(I4B),
pointer :: inunit => null()
17 integer(I4B),
pointer :: iout => null()
18 integer(I4B),
pointer :: inewton => null()
19 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
20 integer(I4B),
dimension(:),
pointer,
contiguous :: iax => null()
21 integer(I4B),
dimension(:),
pointer,
contiguous :: jax => null()
22 integer(I4B),
dimension(:),
pointer,
contiguous :: idxglox => null()
23 integer(I4B),
dimension(:),
pointer,
contiguous :: ia_xt3d => null()
24 integer(I4B),
dimension(:),
pointer,
contiguous :: ja_xt3d => null()
25 integer(I4B),
pointer :: numextnbrs => null()
26 integer(I4B),
pointer :: ixt3d => null()
27 logical,
pointer :: nozee => null()
28 real(dp),
pointer :: vcthresh => null()
29 real(dp),
dimension(:, :),
pointer,
contiguous :: rmatck => null()
30 real(dp),
dimension(:),
pointer,
contiguous :: qsat => null()
31 integer(I4B),
pointer :: nbrmax => null()
32 real(dp),
dimension(:),
pointer,
contiguous :: amatpc => null()
33 real(dp),
dimension(:),
pointer,
contiguous :: amatpcx => null()
34 integer(I4B),
dimension(:),
pointer,
contiguous :: iallpc => null()
35 logical,
pointer :: lamatsaved => null()
38 real(dp),
dimension(:),
pointer,
contiguous :: k11 => null()
39 real(dp),
dimension(:),
pointer,
contiguous :: k22 => null()
40 real(dp),
dimension(:),
pointer,
contiguous :: k33 => null()
41 integer(I4B),
pointer :: ik22 => null()
42 integer(I4B),
pointer :: ik33 => null()
43 real(dp),
dimension(:),
pointer,
contiguous :: sat => null()
44 integer(I4B),
dimension(:),
pointer,
contiguous :: icelltype => null()
45 integer(I4B),
pointer :: iangle1 => null()
46 integer(I4B),
pointer :: iangle2 => null()
47 integer(I4B),
pointer :: iangle3 => null()
48 real(dp),
dimension(:),
pointer,
contiguous :: angle1 => null()
49 real(dp),
dimension(:),
pointer,
contiguous :: angle2 => null()
50 real(dp),
dimension(:),
pointer,
contiguous :: angle3 => null()
51 logical,
pointer :: ldispersion => null()
89 subroutine xt3d_cr(xt3dobj, name_model, inunit, iout, ldispopt)
92 character(len=*),
intent(in) :: name_model
93 integer(I4B),
intent(in) :: inunit
94 integer(I4B),
intent(in) :: iout
95 logical,
optional,
intent(in) :: ldispopt
104 call xt3dobj%allocate_scalars()
107 xt3dobj%inunit = inunit
109 if (
present(ldispopt)) xt3dobj%ldispersion = ldispopt
136 integer(I4B),
intent(in) :: moffset
140 integer(I4B) :: i, j, k, jj, kk, iglo, kglo, iadded
142 integer(I4B) :: ierror
145 if (this%ixt3d == 1)
then
149 call sparse_xt3d%init(this%dis%nodes, this%dis%nodes, nnz)
154 do i = 1, this%dis%nodes
157 do jj = this%dis%con%ia(i) + 1, this%dis%con%ia(i + 1) - 1
158 j = this%dis%con%ja(jj)
160 do kk = this%dis%con%ia(j) + 1, this%dis%con%ia(j + 1) - 1
161 k = this%dis%con%ja(kk)
163 call sparse_xt3d%addconnection(i, k, 1)
170 call mem_allocate(this%ia_xt3d, this%dis%nodes + 1,
'IA_XT3D', &
171 trim(this%memoryPath))
172 call mem_allocate(this%ja_xt3d, sparse_xt3d%nnz,
'JA_XT3D', &
173 trim(this%memoryPath))
174 call sparse_xt3d%filliaja(this%ia_xt3d, this%ja_xt3d, ierror)
175 call sparse_xt3d%destroy()
179 do i = 1, this%dis%nodes
181 do kk = this%ia_xt3d(i), this%ia_xt3d(i + 1) - 1
184 call sparse%addconnection(iglo, kglo, 1, iadded)
185 this%numextnbrs = this%numextnbrs + 1
190 call mem_allocate(this%ia_xt3d, 0,
'IA_XT3D', trim(this%memoryPath))
191 call mem_allocate(this%ja_xt3d, 0,
'JA_XT3D', trim(this%memoryPath))
205 integer(I4B),
intent(in) :: moffset
208 integer(I4B) :: i, j, iglo, jglo, niax, njax, ipos
209 integer(I4B) :: ipos_sln, icol_first, icol_last
210 integer(I4B) :: jj_xt3d
211 integer(I4B) :: igfirstnod, iglastnod
216 if (this%ixt3d == 1)
then
220 igfirstnod = moffset + 1
221 iglastnod = moffset + this%dis%nodes
224 niax = this%dis%nodes + 1
225 njax = this%numextnbrs
226 call mem_allocate(this%iax, niax,
'IAX', trim(this%memoryPath))
227 call mem_allocate(this%jax, njax,
'JAX', trim(this%memoryPath))
228 call mem_allocate(this%idxglox, njax,
'IDXGLOX', trim(this%memoryPath))
235 do i = 1, this%dis%nodes
239 icol_first = matrix_sln%get_first_col_pos(iglo)
240 icol_last = matrix_sln%get_last_col_pos(iglo)
241 do ipos_sln = icol_first, icol_last
245 jglo = matrix_sln%get_column(ipos_sln)
246 if (jglo < igfirstnod .or. jglo > iglastnod)
then
255 do jj_xt3d = this%ia_xt3d(i), this%ia_xt3d(i + 1) - 1
256 if (j == this%ja_xt3d(jj_xt3d))
then
264 this%jax(ipos) = matrix_sln%get_column(ipos_sln) - moffset
265 this%idxglox(ipos) = ipos_sln
270 this%iax(i + 1) = ipos
275 call mem_allocate(this%iax, 0,
'IAX', trim(this%memoryPath))
276 call mem_allocate(this%jax, 0,
'JAX', trim(this%memoryPath))
277 call mem_allocate(this%idxglox, 0,
'IDXGLOX', trim(this%memoryPath))
287 subroutine xt3d_ar(this, ibound, k11, ik33, k33, sat, ik22, k22, iangle1, &
288 iangle2, iangle3, angle1, angle2, angle3, inewton, icelltype)
293 integer(I4B),
dimension(:),
pointer,
contiguous,
intent(inout) :: ibound
294 real(DP),
dimension(:),
intent(in),
pointer,
contiguous :: k11
295 integer(I4B),
intent(in),
pointer :: ik33
296 real(DP),
dimension(:),
intent(in),
pointer,
contiguous :: k33
297 real(DP),
dimension(:),
intent(in),
pointer,
contiguous :: sat
298 integer(I4B),
intent(in),
pointer :: ik22
299 real(DP),
dimension(:),
intent(in),
pointer,
contiguous :: k22
300 integer(I4B),
intent(in),
pointer :: iangle1
301 integer(I4B),
intent(in),
pointer :: iangle2
302 integer(I4B),
intent(in),
pointer :: iangle3
303 real(DP),
dimension(:),
intent(in),
pointer,
contiguous :: angle1
304 real(DP),
dimension(:),
intent(in),
pointer,
contiguous :: angle2
305 real(DP),
dimension(:),
intent(in),
pointer,
contiguous :: angle3
306 integer(I4B),
intent(in),
pointer,
optional :: inewton
307 integer(I4B),
dimension(:),
intent(in),
pointer, &
308 contiguous,
optional :: icelltype
310 integer(I4B) :: n, nnbrs
312 character(len=*),
parameter :: fmtheader = &
313 "(1x, /1x, 'XT3D is active.'//)"
316 if (this%iout > 0)
then
317 write (this%iout, fmtheader)
321 this%ibound => ibound
328 this%iangle1 => iangle1
329 this%iangle2 => iangle2
330 this%iangle3 => iangle3
331 this%angle1 => angle1
332 this%angle2 => angle2
333 this%angle3 => angle3
335 if (
present(inewton))
then
337 this%inewton = inewton
339 if (
present(icelltype))
then
343 this%icelltype => icelltype
348 if (this%iangle2 == 0) this%nozee = .true.
352 do n = 1, this%dis%nodes
353 nnbrs = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
354 this%nbrmax = max(nnbrs, this%nbrmax)
358 if (this%dis%icondir == 0)
then
359 call store_error(
'Vertices not specified for discretization package, '// &
360 'but XT3D is active: '//trim(adjustl(this%memoryPath))// &
361 '. Vertices must be specified in discretization '// &
362 'package in order to use XT3D.', terminate=.true.)
366 if (this%dis%con%ianglex == 0)
then
367 call store_error(
'ANGLDEGX is not specified in the DIS package, '// &
368 'but XT3D is active: '//trim(adjustl(this%memoryPath))// &
369 '. ANGLDEGX must be provided in discretization '// &
370 'package in order to use XT3D.', terminate=.true.)
374 call this%allocate_arrays()
378 if (this%lamatsaved .and. .not. this%ldispersion) &
379 call this%xt3d_fcpc(this%dis%nodes, .true.)
386 subroutine xt3d_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
392 integer(I4B) :: kiter
394 integer(I4B),
intent(in),
dimension(:) :: idxglo
395 real(DP),
intent(inout),
dimension(:) :: rhs
396 real(DP),
intent(inout),
dimension(:) :: hnew
398 integer(I4B) :: nodes, nja
399 integer(I4B) :: n, m, ipos
400 logical :: allhc0, allhc1
401 integer(I4B) :: nnbr0, nnbr1
402 integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
404 integer(I4B),
dimension(this%nbrmax) :: inbr0, inbr1
405 real(DP) :: ar01, ar10
406 real(DP),
dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
407 real(DP),
dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
408 real(DP),
dimension(3, 3) :: ck0, ck1
410 real(DP),
dimension(this%nbrmax) :: chati0, chat1j
411 real(DP) :: qnm, qnbrs
416 nodes = this%dis%nodes
417 nja = this%dis%con%nja
418 if (this%lamatsaved)
then
419 do i = 1, this%dis%con%nja
420 call matrix_sln%add_value_pos(idxglo(i), this%amatpc(i))
422 do i = 1, this%numextnbrs
423 call matrix_sln%add_value_pos(this%idxglox(i), this%amatpcx(i))
429 if (this%ibound(n) .eq. 0) cycle
431 if (this%lamatsaved)
then
432 if (this%iallpc(n) == 1) cycle
434 nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
436 call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
441 ipos = this%dis%con%ia(n) + il0
442 if (this%dis%con%mask(ipos) == 0) cycle
446 if ((m .eq. 0) .or. (m .lt. n)) cycle
447 nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
449 call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
452 call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
455 if (this%inewton /= 0)
then
459 call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
463 call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, ck0, &
464 nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
465 this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
469 if (this%inewton /= 0)
then
471 qnm = chat01 * (hnew(m) - hnew(n))
473 call this%xt3d_qnbrs(nodes, n, m, nnbr0, inbr0, chati0, hnew, qnbrs)
476 call this%xt3d_qnbrs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, qnbrs)
479 call this%xt3d_areas(nodes, n, m, jjs01, .true., ar01, ar10, hnew)
480 this%qsat(ii01) = qnm * ar01
482 call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
483 chat01 = chat01 * ar01
484 chati0 = chati0 * ar01
485 chat1j = chat1j * ar01
489 call matrix_sln%add_value_pos(idxglo(ii00), -chat01)
490 call matrix_sln%add_value_pos(idxglo(ii01), chat01)
491 call matrix_sln%add_value_pos(idxglo(ii11), -chat01)
492 call matrix_sln%add_value_pos(idxglo(ii10), chat01)
493 if (this%ixt3d == 1)
then
494 call this%xt3d_amat_nbrs(nodes, n, ii00, nnbr0, nja, matrix_sln, &
495 inbr0, idxglo, chati0)
496 call this%xt3d_amat_nbrnbrs(nodes, n, m, ii01, nnbr1, nja, matrix_sln, &
497 inbr1, idxglo, chat1j)
498 call this%xt3d_amat_nbrs(nodes, m, ii11, nnbr1, nja, matrix_sln, &
499 inbr1, idxglo, chat1j)
500 call this%xt3d_amat_nbrnbrs(nodes, m, n, ii10, nnbr0, nja, matrix_sln, &
501 inbr0, idxglo, chati0)
503 call this%xt3d_rhs(nodes, n, m, nnbr0, inbr0, chati0, hnew, rhs)
504 call this%xt3d_rhs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, rhs)
522 integer(I4B),
intent(in) :: nodes
523 logical,
intent(in) :: lsat
525 integer(I4B) :: n, m, ipos
527 logical :: allhc0, allhc1
528 integer(I4B) :: nnbr0, nnbr1
529 integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
530 integer(I4B),
dimension(this%nbrmax) :: inbr0, inbr1
531 real(DP) :: ar01, ar10
532 real(DP),
dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
533 real(DP),
dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
534 real(DP),
dimension(3, 3) :: ck0, ck1
536 real(DP),
dimension(this%nbrmax) :: chati0, chat1j
539 do n = 1,
size(this%amatpc)
540 this%amatpc(n) =
dzero
542 do n = 1,
size(this%amatpcx)
543 this%amatpcx(n) =
dzero
550 if (this%iallpc(n) == 0) cycle
551 if (this%ibound(n) == 0) cycle
552 nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
554 call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
559 ipos = this%dis%con%ia(n) + il0
560 if (this%dis%con%mask(ipos) == 0) cycle
565 nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
567 call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
570 call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
573 call this%xt3d_areas(nodes, n, m, jjs01, lsat, ar01, ar10)
576 call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, ck0, &
577 nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
578 this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
580 this%amatpc(ii00) = this%amatpc(ii00) - chat01
581 this%amatpc(ii01) = this%amatpc(ii01) + chat01
582 this%amatpc(ii11) = this%amatpc(ii11) - chat01
583 this%amatpc(ii10) = this%amatpc(ii10) + chat01
584 call this%xt3d_amatpc_nbrs(nodes, n, ii00, nnbr0, inbr0, chati0)
585 call this%xt3d_amatpcx_nbrnbrs(nodes, n, m, ii01, nnbr1, inbr1, chat1j)
586 call this%xt3d_amatpc_nbrs(nodes, m, ii11, nnbr1, inbr1, chat1j)
587 call this%xt3d_amatpcx_nbrnbrs(nodes, m, n, ii10, nnbr0, inbr0, chati0)
597 subroutine xt3d_fhfb(this, kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew, &
604 integer(I4B) :: kiter
605 integer(I4B),
intent(in) :: nodes
606 integer(I4B),
intent(in) :: nja
609 integer(I4B),
intent(in),
dimension(nja) :: idxglo
610 real(DP),
intent(inout),
dimension(nodes) :: rhs
611 real(DP),
intent(inout),
dimension(nodes) :: hnew
614 logical :: allhc0, allhc1
615 integer(I4B) :: nnbr0, nnbr1
616 integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10, il
617 integer(I4B),
dimension(this%nbrmax) :: inbr0, inbr1
618 real(DP) :: ar01, ar10
619 real(DP),
dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
620 real(DP),
dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
621 real(DP),
dimension(3, 3) :: ck0, ck1
623 real(DP),
dimension(this%nbrmax) :: chati0, chat1j
624 real(DP) :: qnm, qnbrs
630 nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
632 call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
636 if (inbr0(il) .eq. m)
then
641 nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
643 call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
646 call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
649 if (this%inewton /= 0)
then
653 call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
658 call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, &
659 ck0, nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
660 this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
663 if (condhfb > dzero)
then
664 term = chat01 / (chat01 + condhfb)
668 chat01 = -chat01 * term
669 chati0 = -chati0 * term
670 chat1j = -chat1j * term
674 if (this%inewton /= 0)
then
676 qnm = chat01 * (hnew(m) - hnew(n))
678 call this%xt3d_qnbrs(nodes, n, m, nnbr0, inbr0, chati0, hnew, qnbrs)
681 call this%xt3d_qnbrs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, qnbrs)
684 call this%xt3d_areas(nodes, n, m, jjs01, .true., ar01, ar10, hnew)
685 this%qsat(ii01) = this%qsat(ii01) + qnm * ar01
687 call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
688 chat01 = chat01 * ar01
689 chati0 = chati0 * ar01
690 chat1j = chat1j * ar01
694 call matrix_sln%add_value_pos(idxglo(ii00), -chat01)
695 call matrix_sln%add_value_pos(idxglo(ii01), chat01)
696 call matrix_sln%add_value_pos(idxglo(ii11), -chat01)
697 call matrix_sln%add_value_pos(idxglo(ii10), chat01)
698 if (this%ixt3d == 1)
then
699 call this%xt3d_amat_nbrs(nodes, n, ii00, nnbr0, nja, matrix_sln, &
700 inbr0, idxglo, chati0)
701 call this%xt3d_amat_nbrnbrs(nodes, n, m, ii01, nnbr1, nja, matrix_sln, &
702 inbr1, idxglo, chat1j)
703 call this%xt3d_amat_nbrs(nodes, m, ii11, nnbr1, nja, matrix_sln, &
704 inbr1, idxglo, chat1j)
705 call this%xt3d_amat_nbrnbrs(nodes, m, n, ii10, nnbr0, nja, matrix_sln, &
706 inbr0, idxglo, chati0)
708 call this%xt3d_rhs(nodes, n, m, nnbr0, inbr0, chati0, hnew, rhs)
709 call this%xt3d_rhs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, rhs)
717 subroutine xt3d_fn(this, kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew)
723 integer(I4B) :: kiter
724 integer(I4B),
intent(in) :: nodes
725 integer(I4B),
intent(in) :: nja
727 integer(I4B),
intent(in),
dimension(nja) :: idxglo
728 real(DP),
intent(inout),
dimension(nodes) :: rhs
729 real(DP),
intent(inout),
dimension(nodes) :: hnew
731 integer(I4B) :: n, m, ipos
732 integer(I4B) :: nnbr0
733 integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
734 integer(I4B),
dimension(this%nbrmax) :: inbr0
735 integer(I4B) :: iups, idn
736 real(DP) :: topup, botup, derv, term
742 if (this%ibound(n) .eq. 0) cycle
746 if (this%lamatsaved)
then
747 if (this%iallpc(n) == 1) cycle
749 nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
753 call this%xt3d_load_inbr(n, nnbr0, inbr0)
758 ipos = this%dis%con%ia(n) + il0
759 if (this%dis%con%mask(ipos) == 0) cycle
764 if ((inbr0(il0) .eq. 0) .or. (m .lt. n)) cycle
767 call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
772 if (hnew(m) < hnew(n)) iups = n
774 if (iups == n) idn = m
777 if ((this%icelltype(iups) == 0) .and. (this%ixt3d .eq. 1)) cycle
781 topup = this%dis%top(iups)
782 botup = this%dis%bot(iups)
783 if (this%dis%con%ihc(jjs01) == 2)
then
784 topup = min(this%dis%top(n), this%dis%top(m))
785 botup = max(this%dis%bot(n), this%dis%bot(m))
790 term = this%qsat(ii01) * derv
796 call matrix_sln%add_value_pos(idxglo(ii00), term)
797 rhs(n) = rhs(n) + term * hnew(n)
800 call matrix_sln%add_value_pos(idxglo(ii10), -term)
801 rhs(m) = rhs(m) - term * hnew(n)
807 call matrix_sln%add_value_pos(idxglo(ii01), term)
808 rhs(n) = rhs(n) + term * hnew(m)
811 call matrix_sln%add_value_pos(idxglo(ii11), -term)
812 rhs(m) = rhs(m) - term * hnew(m)
829 real(DP),
intent(inout),
dimension(:) :: hnew
830 real(DP),
intent(inout),
dimension(:) :: flowja
832 integer(I4B) :: n, ipos, m, nodes
833 real(DP) :: qnm, qnbrs
834 logical :: allhc0, allhc1
835 integer(I4B) :: nnbr0, nnbr1
836 integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
837 integer(I4B),
dimension(this%nbrmax) :: inbr0, inbr1
838 real(DP) :: ar01, ar10
839 real(DP),
dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
840 real(DP),
dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
841 real(DP),
dimension(3, 3) :: ck0, ck1
843 real(DP),
dimension(this%nbrmax) :: chati0, chat1j
846 nodes = this%dis%nodes
850 if (this%ibound(n) .eq. 0) cycle
851 nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
854 call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
863 if ((inbr0(il0) .eq. 0) .or. (m .lt. n)) cycle
864 nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
867 call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
871 call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
875 if (this%inewton /= 0) &
876 call this%xt3d_areas(nodes, n, m, jjs01, .true., ar01, ar10, hnew)
877 call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
881 call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, ck0, &
882 nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
883 this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
886 qnm = chat01 * (hnew(m) - hnew(n))
889 call this%xt3d_qnbrs(nodes, n, m, nnbr0, inbr0, chati0, hnew, qnbrs)
893 call this%xt3d_qnbrs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, qnbrs)
896 flowja(ipos) = flowja(ipos) + qnm
897 flowja(this%dis%con%isym(ipos)) = flowja(this%dis%con%isym(ipos)) - qnm
914 real(DP),
intent(inout),
dimension(:) :: hnew
915 real(DP),
intent(inout),
dimension(:) :: flowja
918 integer(I4B) :: nodes
919 logical :: allhc0, allhc1
921 integer(I4B) :: nnbr0, nnbr1
922 integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10, il
923 integer(I4B),
dimension(this%nbrmax) :: inbr0, inbr1
925 real(DP) :: ar01, ar10
926 real(DP),
dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
927 real(DP),
dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
928 real(DP),
dimension(3, 3) :: ck0, ck1
930 real(DP),
dimension(this%nbrmax) :: chati0, chat1j
931 real(DP) :: qnm, qnbrs
936 nodes = this%dis%nodes
937 nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
940 call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
945 if (inbr0(il) .eq. m)
then
951 nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
954 call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
958 call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
962 if (this%inewton /= 0)
then
966 call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
971 call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, &
972 ck0, nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
973 this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
976 if (condhfb > dzero)
then
977 term = chat01 / (chat01 + condhfb)
981 chat01 = -chat01 * term
982 chati0 = -chati0 * term
983 chat1j = -chat1j * term
986 qnm = chat01 * (hnew(m) - hnew(n))
989 call this%xt3d_qnbrs(nodes, n, m, nnbr0, inbr0, chati0, hnew, qnbrs)
993 call this%xt3d_qnbrs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, qnbrs)
998 if (this%inewton /= 0)
then
999 call this%xt3d_areas(nodes, n, m, jjs01, .true., ar01, ar10, hnew)
1000 call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
1005 flowja(ipos) = flowja(ipos) + qnm
1006 flowja(this%dis%con%isym(ipos)) = flowja(this%dis%con%isym(ipos)) - qnm
1021 if (this%ixt3d /= 0)
then
1059 call mem_allocate(this%ixt3d,
'IXT3D', this%memoryPath)
1060 call mem_allocate(this%nbrmax,
'NBRMAX', this%memoryPath)
1061 call mem_allocate(this%inunit,
'INUNIT', this%memoryPath)
1063 call mem_allocate(this%inewton,
'INEWTON', this%memoryPath)
1064 call mem_allocate(this%numextnbrs,
'NUMEXTNBRS', this%memoryPath)
1065 call mem_allocate(this%nozee,
'NOZEE', this%memoryPath)
1066 call mem_allocate(this%vcthresh,
'VCTHRESH', this%memoryPath)
1067 call mem_allocate(this%lamatsaved,
'LAMATSAVED', this%memoryPath)
1068 call mem_allocate(this%ldispersion,
'LDISPERSION', this%memoryPath)
1077 this%nozee = .false.
1078 this%vcthresh = 1.d-10
1079 this%lamatsaved = .false.
1080 this%ldispersion = .false.
1094 integer(I4B) :: njax
1098 if (this%inewton /= 0)
then
1099 call mem_allocate(this%qsat, this%dis%nja,
'QSAT', this%memoryPath)
1101 call mem_allocate(this%qsat, 0,
'QSAT', this%memoryPath)
1107 if (this%ldispersion)
then
1111 this%lamatsaved = .true.
1112 call mem_allocate(this%iallpc, this%dis%nodes,
'IALLPC', this%memoryPath)
1113 do n = 1, this%dis%nodes
1121 call this%xt3d_iallpc()
1125 if (this%lamatsaved)
then
1126 call mem_allocate(this%amatpc, this%dis%nja,
'AMATPC', this%memoryPath)
1127 njax = this%numextnbrs
1128 call mem_allocate(this%amatpcx, njax,
'AMATPCX', this%memoryPath)
1130 call mem_allocate(this%amatpc, 0,
'AMATPC', this%memoryPath)
1131 call mem_allocate(this%amatpcx, 0,
'AMATPCX', this%memoryPath)
1135 call mem_allocate(this%rmatck, 3, 3,
'RMATCK', this%memoryPath)
1139 if (this%inewton /= 0)
then
1141 else if (this%lamatsaved)
then
1143 this%amatpcx =
dzero
1157 integer(I4B) :: n, m, mm, il0, il1
1158 integer(I4B) :: nnbr0, nnbr1
1159 integer(I4B),
dimension(this%nbrmax) :: inbr0, inbr1
1161 if (this%ixt3d == 2)
then
1162 this%lamatsaved = .false.
1163 call mem_allocate(this%iallpc, 0,
'IALLPC', this%memoryPath)
1167 call mem_allocate(this%iallpc, this%dis%nodes,
'IALLPC', this%memoryPath)
1168 do n = 1, this%dis%nodes
1174 do n = 1, this%dis%nodes
1175 if (this%icelltype(n) /= 0)
then
1179 nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
1180 call this%xt3d_load_inbr(n, nnbr0, inbr0)
1184 if (this%icelltype(m) /= 0)
then
1189 nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
1190 call this%xt3d_load_inbr(m, nnbr1, inbr1)
1193 if (this%icelltype(mm) /= 0)
then
1205 this%lamatsaved = .false.
1206 do n = 1, this%dis%nodes
1207 if (this%iallpc(n) == 1)
then
1208 this%lamatsaved = .true.
1214 if (.not. this%lamatsaved)
then
1216 call mem_allocate(this%iallpc, 0,
'IALLPC', this%memoryPath)
1229 integer(I4B) :: n, m, il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
1231 integer(I4B) :: iinm
1238 call this%xt3d_get_iinm(m, n, iinm)
1239 il10 = iinm - this%dis%con%ia(m)
1241 ii00 = this%dis%con%ia(n)
1245 jjs01 = this%dis%con%jas(ii01)
1247 ii11 = this%dis%con%ia(m)
1258 subroutine xt3d_load(this, nodes, n, nnbr, inbr, vc, vn, dl, dln, ck, allhc)
1264 integer(I4B),
intent(in) :: nodes
1265 integer(I4B) :: n, nnbr
1266 integer(I4B),
dimension(this%nbrmax) :: inbr
1267 real(DP),
dimension(this%nbrmax, 3) :: vc, vn
1268 real(DP),
dimension(this%nbrmax) :: dl, dln
1269 real(DP),
dimension(3, 3) :: ck
1271 integer(I4B) :: il, ii, jj, jjs
1272 integer(I4B) :: ihcnjj
1273 real(DP) :: satn, satjj
1274 real(DP) :: cl1njj, cl2njj, dltot, ooclsum
1278 ck(1, 1) = this%k11(n)
1279 ck(2, 2) = this%k22(n)
1280 ck(3, 3) = this%k33(n)
1281 call this%xt3d_fillrmatck(n)
1282 ck = matmul(this%rmatck, ck)
1283 ck = matmul(ck, transpose(this%rmatck))
1291 ii = il + this%dis%con%ia(n)
1292 jj = this%dis%con%ja(ii)
1293 jjs = this%dis%con%jas(ii)
1294 if (this%ibound(jj) .ne. 0)
then
1297 satjj = this%sat(jj)
1300 ihcnjj = this%dis%con%ihc(jjs)
1301 call this%dis%connection_normal(n, jj, ihcnjj, vn(il, 1), vn(il, 2), &
1303 call this%dis%connection_vector(n, jj, this%nozee, satn, satjj, ihcnjj, &
1304 vc(il, 1), vc(il, 2), vc(il, 3), dltot)
1306 cl1njj = this%dis%con%cl1(jjs)
1307 cl2njj = this%dis%con%cl2(jjs)
1309 cl1njj = this%dis%con%cl2(jjs)
1310 cl2njj = this%dis%con%cl1(jjs)
1312 ooclsum = 1d0 / (cl1njj + cl2njj)
1313 dl(il) = dltot * cl1njj * ooclsum
1314 dln(il) = dltot * cl2njj * ooclsum
1315 if (this%dis%con%ihc(jjs) .eq. 0) allhc = .false.
1330 integer(I4B) :: n, nnbr
1331 integer(I4B),
dimension(this%nbrmax) :: inbr
1333 integer(I4B) :: il, ii, jj
1338 ii = il + this%dis%con%ia(n)
1339 jj = this%dis%con%ja(ii)
1340 if (this%ibound(jj) .ne. 0)
then
1353 subroutine xt3d_areas(this, nodes, n, m, jjs01, lsat, ar01, ar10, hnew)
1357 integer(I4B) :: nodes, n, m, jjs01
1358 real(DP) :: ar01, ar10
1359 real(DP),
intent(inout),
dimension(:),
optional :: hnew
1361 real(DP) :: topn, botn, topm, botm, thksatn, thksatm
1362 real(DP) :: sill_top, sill_bot, tpn, tpm
1366 if (this%dis%con%ihc(jjs01) == 0)
then
1369 ar01 = this%dis%con%hwva(jjs01)
1371 else if (this%inewton /= 0)
then
1377 topn = this%dis%top(n)
1378 botn = this%dis%bot(n)
1379 topm = this%dis%top(m)
1380 botm = this%dis%bot(m)
1381 thksatn = topn - botn
1382 thksatm = topm - botm
1383 if (this%dis%con%ihc(jjs01) .eq. 2)
then
1385 sill_top = min(topn, topm)
1386 sill_bot = max(botn, botm)
1387 tpn = botn + thksatn
1388 tpm = botm + thksatm
1389 thksatn = max(min(tpn, sill_top) - sill_bot,
dzero)
1390 thksatm = max(min(tpm, sill_top) - sill_bot,
dzero)
1392 ar01 = this%dis%con%hwva(jjs01) *
dhalf * (thksatn + thksatm)
1398 if (hnew(m) < hnew(n))
then
1399 satups = this%sat(n)
1401 satups = this%sat(m)
1403 ar01 = ar01 * satups
1409 topn = this%dis%top(n)
1410 botn = this%dis%bot(n)
1411 topm = this%dis%top(m)
1412 botm = this%dis%bot(m)
1413 thksatn = topn - botn
1414 thksatm = topm - botm
1415 if (.not. lsat)
then
1416 thksatn = this%sat(n) * thksatn
1417 thksatm = this%sat(m) * thksatm
1419 if (this%dis%con%ihc(jjs01) == 2)
then
1421 sill_top = min(topn, topm)
1422 sill_bot = max(botn, botm)
1423 tpn = botn + thksatn
1424 tpm = botm + thksatm
1425 thksatn = max(min(tpn, sill_top) - sill_bot,
dzero)
1426 thksatm = max(min(tpm, sill_top) - sill_bot,
dzero)
1428 ar01 = this%dis%con%hwva(jjs01) * thksatn
1429 ar10 = this%dis%con%hwva(jjs01) * thksatm
1439 matrix_sln, inbr, idxglo, chat)
1442 integer(I4B),
intent(in) :: nodes
1443 integer(I4B) :: n, idiag, nnbr, nja
1445 integer(I4B),
dimension(this%nbrmax) :: inbr
1446 integer(I4B),
intent(in),
dimension(nja) :: idxglo
1447 real(DP),
dimension(this%nbrmax) :: chat
1449 integer(I4B) :: iil, iii
1452 if (inbr(iil) .ne. 0)
then
1453 iii = this%dis%con%ia(n) + iil
1454 call matrix_sln%add_value_pos(idxglo(idiag), -chat(iil))
1455 call matrix_sln%add_value_pos(idxglo(iii), chat(iil))
1466 matrix_sln, inbr, idxglo, chat)
1469 integer(I4B),
intent(in) :: nodes
1470 integer(I4B) :: n, m, ii01, nnbr, nja
1472 integer(I4B),
dimension(this%nbrmax) :: inbr
1473 integer(I4B),
intent(in),
dimension(nja) :: idxglo
1474 real(DP),
dimension(this%nbrmax) :: chat
1476 integer(I4B) :: iil, iii, jjj, iixjjj, iijjj
1479 if (inbr(iil) .ne. 0)
then
1480 call matrix_sln%add_value_pos(idxglo(ii01), chat(iil))
1481 iii = this%dis%con%ia(m) + iil
1482 jjj = this%dis%con%ja(iii)
1483 call this%xt3d_get_iinmx(n, jjj, iixjjj)
1484 if (iixjjj .ne. 0)
then
1485 call matrix_sln%add_value_pos(this%idxglox(iixjjj), -chat(iil))
1487 call this%xt3d_get_iinm(n, jjj, iijjj)
1488 call matrix_sln%add_value_pos(idxglo(iijjj), -chat(iil))
1502 integer(I4B),
intent(in) :: nodes
1503 integer(I4B) :: n, idiag, nnbr
1504 integer(I4B),
dimension(this%nbrmax) :: inbr
1505 real(DP),
dimension(this%nbrmax) :: chat
1507 integer(I4B) :: iil, iii
1510 iii = this%dis%con%ia(n) + iil
1511 this%amatpc(idiag) = this%amatpc(idiag) - chat(iil)
1512 this%amatpc(iii) = this%amatpc(iii) + chat(iil)
1524 integer(I4B),
intent(in) :: nodes
1525 integer(I4B) :: n, m, ii01, nnbr
1526 integer(I4B),
dimension(this%nbrmax) :: inbr
1527 real(DP),
dimension(this%nbrmax) :: chat
1529 integer(I4B) :: iil, iii, jjj, iixjjj, iijjj
1532 this%amatpc(ii01) = this%amatpc(ii01) + chat(iil)
1533 iii = this%dis%con%ia(m) + iil
1534 jjj = this%dis%con%ja(iii)
1535 call this%xt3d_get_iinmx(n, jjj, iixjjj)
1536 if (iixjjj .ne. 0)
then
1537 this%amatpcx(iixjjj) = this%amatpcx(iixjjj) - chat(iil)
1539 call this%xt3d_get_iinm(n, jjj, iijjj)
1540 this%amatpc(iijjj) = this%amatpc(iijjj) - chat(iil)
1554 integer(I4B) :: n, m, iinm
1556 integer(I4B) :: ii, jj
1559 do ii = this%dis%con%ia(n), this%dis%con%ia(n + 1) - 1
1560 jj = this%dis%con%ja(ii)
1577 integer(I4B) :: n, m, iinmx
1579 integer(I4B) :: iix, jjx
1582 do iix = this%iax(n), this%iax(n + 1) - 1
1584 if (jjx .eq. m)
then
1596 subroutine xt3d_rhs(this, nodes, n, m, nnbr, inbr, chat, hnew, &
1600 integer(I4B),
intent(in) :: nodes
1601 integer(I4B) :: n, m, nnbr
1602 integer(I4B),
dimension(this%nbrmax) :: inbr
1603 real(DP),
dimension(this%nbrmax) :: chat
1604 real(DP),
intent(inout),
dimension(nodes) :: hnew, rhs
1606 integer(I4B) :: iil, iii, jjj
1610 if (inbr(iil) .ne. 0)
then
1611 iii = iil + this%dis%con%ia(n)
1612 jjj = this%dis%con%ja(iii)
1613 term = chat(iil) * (hnew(jjj) - hnew(n))
1614 rhs(n) = rhs(n) - term
1615 rhs(m) = rhs(m) + term
1625 subroutine xt3d_qnbrs(this, nodes, n, m, nnbr, inbr, chat, hnew, &
1629 integer(I4B),
intent(in) :: nodes
1630 integer(I4B) :: n, m, nnbr
1631 integer(I4B),
dimension(this%nbrmax) :: inbr
1633 real(DP),
dimension(this%nbrmax) :: chat
1634 real(DP),
intent(inout),
dimension(nodes) :: hnew
1636 integer(I4B) :: iil, iii, jjj
1641 if (inbr(iil) .ne. 0)
then
1642 iii = iil + this%dis%con%ia(n)
1643 jjj = this%dis%con%ja(iii)
1644 term = chat(iil) * (hnew(jjj) - hnew(n))
1645 qnbrs = qnbrs + term
1660 integer(I4B),
intent(in) :: n
1662 real(DP) :: ang1, ang2, ang3, ang2d, ang3d
1663 real(DP) :: s1, c1, s2, c2, s3, c3
1665 if (this%nozee)
then
1668 ang1 = this%angle1(n)
1672 ang1 = this%angle1(n)
1673 ang2 = this%angle2(n)
1674 ang3 = this%angle3(n)
1682 this%rmatck(1, 1) = c1 * c2
1683 this%rmatck(1, 2) = c1 * s2 * s3 - s1 * c3
1684 this%rmatck(1, 3) = -c1 * s2 * c3 - s1 * s3
1685 this%rmatck(2, 1) = s1 * c2
1686 this%rmatck(2, 2) = s1 * s2 * s3 + c1 * c3
1687 this%rmatck(2, 3) = -s1 * s2 * c3 + c1 * s3
1688 this%rmatck(3, 1) = s2
1689 this%rmatck(3, 2) = -c2 * s3
1690 this%rmatck(3, 3) = c2 * c3
This module contains simulation constants.
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dzero
real constant zero
integer(i4b), parameter lenmempath
maximum length of the memory path
real(dp), parameter done
real constant 1
This module defines variable data types.
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
This module contains simulation methods.
subroutine, public store_error(msg, terminate)
Store an error message.
real(dp) function squadraticsaturationderivative(top, bot, x, eps)
@ brief Derivative of the quadratic saturation function
subroutine qconds(nnbrmx, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, ck0, nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
Compute the "conductances" in the normal-flux expression for an interface (modflow-usg version)....
subroutine xt3d_load(this, nodes, n, nnbr, inbr, vc, vn, dl, dln, ck, allhc)
Load conductivity and connection info for a cell into arrays used by XT3D.
subroutine xt3d_ac(this, moffset, sparse)
Add connections for extended neighbors to the sparse matrix.
subroutine xt3d_get_iinmx(this, n, m, iinmx)
Get position of n-m "extended connection" in jax array (return 0 if not connected)
subroutine xt3d_amatpcx_nbrnbrs(this, nodes, n, m, ii01, nnbr, inbr, chat)
Add contributions from neighbors of neighbor to amatpc and amatpcx.
subroutine allocate_arrays(this)
Allocate xt3d arrays.
subroutine xt3d_amat_nbrnbrs(this, nodes, n, m, ii01, nnbr, nja, matrix_sln, inbr, idxglo, chat)
Add contributions from neighbors of neighbor to amat.
subroutine xt3d_da(this)
Deallocate variables.
subroutine allocate_scalars(this)
Allocate scalar pointer variables.
subroutine xt3d_amat_nbrs(this, nodes, n, idiag, nnbr, nja, matrix_sln, inbr, idxglo, chat)
Add contributions from neighbors to amat.
subroutine xt3d_areas(this, nodes, n, m, jjs01, lsat, ar01, ar10, hnew)
Compute interfacial areas.
subroutine xt3d_flowja(this, hnew, flowja)
Budget.
subroutine xt3d_get_iinm(this, n, m, iinm)
Get position of n-m connection in ja array (return 0 if not connected)
subroutine xt3d_fhfb(this, kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew, n, m, condhfb)
Formulate HFB correction.
subroutine xt3d_flowjahfb(this, n, m, hnew, flowja, condhfb)
hfb contribution to flowja when xt3d is used
subroutine xt3d_indices(this, n, m, il0, ii01, jjs01, il01, il10, ii00, ii11, ii10)
Set various indices for XT3D.
subroutine xt3d_load_inbr(this, n, nnbr, inbr)
Load neighbor list for a cell.
subroutine xt3d_fcpc(this, nodes, lsat)
Formulate for permanently confined connections and save in amatpc and amatpcx.
subroutine xt3d_fn(this, kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew)
Fill Newton terms for xt3d.
subroutine xt3d_fillrmatck(this, n)
Fill rmat array for cell n.
subroutine xt3d_mc(this, moffset, matrix_sln)
Map connections and construct iax, jax, and idxglox.
subroutine, public xt3d_cr(xt3dobj, name_model, inunit, iout, ldispopt)
Create a new xt3d object.
subroutine xt3d_qnbrs(this, nodes, n, m, nnbr, inbr, chat, hnew, qnbrs)
Add contributions to flow from neighbors.
subroutine xt3d_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
Formulate.
subroutine xt3d_rhs(this, nodes, n, m, nnbr, inbr, chat, hnew, rhs)
Add contributions to rhs.
subroutine xt3d_df(this, dis)
Define the xt3d object.
subroutine xt3d_iallpc(this)
Allocate and populate iallpc array. Set lamatsaved.
subroutine xt3d_ar(this, ibound, k11, ik33, k33, sat, ik22, k22, iangle1, iangle2, iangle3, angle1, angle2, angle3, inewton, icelltype)
Allocate and Read.
subroutine xt3d_amatpc_nbrs(this, nodes, n, idiag, nnbr, inbr, chat)
Add contributions from neighbors to amatpc.