47 subroutine qconds(nnbrmx, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, ck0, &
48 nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
49 vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
51 integer(I4B) :: nnbrmx
53 integer(I4B),
dimension(nnbrmx) :: inbr0
55 real(DP),
dimension(nnbrmx, 3) :: vc0
56 real(DP),
dimension(nnbrmx, 3) :: vn0
57 real(DP),
dimension(nnbrmx) :: dl0
58 real(DP),
dimension(nnbrmx) :: dl0n
59 real(DP),
dimension(3, 3) :: ck0
61 integer(I4B),
dimension(nnbrmx) :: inbr1
63 real(DP),
dimension(nnbrmx) :: vc1
64 real(DP),
dimension(nnbrmx) :: vn1
65 real(DP),
dimension(nnbrmx) :: dl1
66 real(DP),
dimension(nnbrmx) :: dl1n
67 real(DP),
dimension(3, 3) :: ck1
74 real(DP),
dimension(nnbrmx) :: chati0
75 real(DP),
dimension(nnbrmx) :: chat1j
83 real(DP),
dimension(nnbrmx) :: bhat0
84 real(DP),
dimension(nnbrmx) :: bhat1
95 if (ar01 .eq. 0d0)
then
104 call abhats(nnbrmx, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, ck0, &
105 vcthresh, allhc0, ar01, ahat0, bhat0)
107 call abhats(nnbrmx, nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, &
108 vcthresh, allhc1, ar10, ahat1, bhat1)
110 denom = (ahat0 + ahat1)
111 if (abs(denom) >
dprec)
then
112 wght1 = ahat0 / (ahat0 + ahat1)
117 chat01 = wght1 * ahat1
119 chati0(i) = wght0 * bhat0(i)
120 chat1j(i) = wght1 * bhat1(i)
130 subroutine abhats(nnbrmx, nnbr, inbr, il01, vc, vn, dl0, dln, ck, &
131 vcthresh, allhc, ar01, ahat, bhat)
133 integer(I4B) :: nnbrmx
135 integer(I4B),
dimension(nnbrmx) :: inbr
137 real(DP),
dimension(nnbrmx, 3) :: vc
138 real(DP),
dimension(nnbrmx, 3) :: vn
139 real(DP),
dimension(nnbrmx) :: dl0
140 real(DP),
dimension(nnbrmx) :: dln
141 real(DP),
dimension(3, 3) :: ck
146 real(DP),
dimension(nnbrmx) :: bhat
149 real(DP),
dimension(nnbrmx, 3) :: vccde
150 real(DP),
dimension(3, 3) :: rmat
151 real(DP),
dimension(3) :: sigma
152 real(DP),
dimension(nnbrmx) :: bd
153 real(DP),
dimension(nnbrmx) :: be
154 real(DP),
dimension(nnbrmx) :: betad
155 real(DP),
dimension(nnbrmx) :: betae
176 call getrot(nnbrmx, nnbr, inbr, vc, il01, rmat, iml0)
182 if (iml0 .eq. 0)
then
186 sigma(1) = dot_product(vn(il01, :), matmul(ck, rmat(:, 1)))
187 ahat = sigma(1) / dl0(il01)
194 call tranvc(nnbrmx, nnbr, rmat, vc, vccde)
197 call abwts(nnbrmx, nnbr, inbr, il01, 2, vccde, &
198 vcthresh, dl0, dln, acd, add, aed, bd)
213 if ((il == il01) .or. (inbr(il) == 0))
then
215 else if (dabs(vccde(il, 3)) > 1d-10)
then
221 call abwts(nnbrmx, nnbr, inbr, il01, 3, vccde, &
222 vcthresh, dl0, dln, ace, aee, ade, be)
232 determ = add * aee - ade * aed
234 alphad = (acd * aee - ace * aed) * oodet
235 alphae = (ace * add - acd * ade) * oodet
240 if ((il == il01) .or. (inbr(il) == 0)) cycle
241 betad(il) = (bd(il) * aee - be(il) * aed) * oodet
242 betae(il) = (be(il) * add - bd(il) * ade) * oodet
246 sigma = matmul(vn(il01, :), matmul(ck, rmat))
249 ahat = (sigma(1) - sigma(2) * alphad - sigma(3) * alphae) / dl0(il01)
253 if ((il == il01) .or. (inbr(il) == 0)) cycle
254 dl0il = dl0(il) + dln(il)
255 bhat(il) = (sigma(2) * betad(il) + sigma(3) * betae(il)) / dl0il
289 subroutine getrot(nnbrmx, nnbr, inbr, vc, il01, rmat, iml0)
291 integer(I4B) :: nnbrmx
293 integer(I4B),
dimension(nnbrmx) :: inbr
294 real(DP),
dimension(nnbrmx, 3) :: vc
296 real(DP),
dimension(3, 3) :: rmat
299 real(DP),
dimension(3) :: vcc
300 real(DP),
dimension(3) :: vcd
301 real(DP),
dimension(3) :: vce
302 real(DP),
dimension(3) :: vcmax
318 if ((il .eq. il01) .or. (inbr(il) .eq. 0))
then
321 cmp = dot_product(vc(il, :), vcc)
323 if (acmp .lt. acmpmn)
then
334 vcmax(:) = vc(iml0, :)
339 vcd = vcmax - cmpmn * vcc
340 vcd = vcd / dsqrt(1d0 - cmpmn * cmpmn)
344 vce(1) = vcc(2) * vcd(3) - vcc(3) * vcd(2)
345 vce(2) = vcc(3) * vcd(1) - vcc(1) * vcd(3)
346 vce(3) = vcc(1) * vcd(2) - vcc(2) * vcd(1)
366 subroutine tranvc(nnbrmx, nnbrs, rmat, vc, vccde)
368 integer(I4B) :: nnbrmx
369 integer(I4B) :: nnbrs
370 real(DP),
dimension(3, 3) :: rmat
371 real(DP),
dimension(nnbrmx, 3) :: vc
372 real(DP),
dimension(nnbrmx, 3) :: vccde
380 vccde(il, :) = matmul(transpose(rmat), vc(il, :))
398 subroutine abwts(nnbrmx, nnbr, inbr, il01, nde1, vccde, &
399 vcthresh, dl0, dln, acd, add, aed, bd)
401 integer(I4B) :: nnbrmx
403 integer(I4B),
dimension(nnbrmx) :: inbr
406 real(DP),
dimension(nnbrmx, 3) :: vccde
408 real(DP),
dimension(nnbrmx) :: dl0
409 real(DP),
dimension(nnbrmx) :: dln
413 real(DP),
dimension(nnbrmx) :: bd
425 real(DP),
dimension(nnbrmx) :: omwt
436 if ((il .eq. il01) .or. (inbr(il) .eq. 0)) cycle
437 vcmx = max(dabs(vccde(il, nde1)), vcmx)
438 dlm = 5d-1 * (dl0(il) + dln(il))
444 cosang = vccde(il, 1)
445 dl4wt = dsqrt(dlm * dlm + dl0(il01) * dl0(il01) &
446 - 2d0 * dlm * dl0(il01) * cosang)
447 omwt(il) = dabs(vccde(il, nde1)) * dl4wt
448 dsum = dsum + omwt(il)
455 dsum = dsum + 1d-10 * dsum
458 if ((il .eq. il01) .or. (inbr(il) .eq. 0)) cycle
459 fact = dsum - omwt(il)
460 omwt(il) = fact * dabs(vccde(il, nde1))
468 if ((il .eq. il01) .or. (inbr(il) .eq. 0)) cycle
469 bd(il) = omwt(il) * sign(1d0, vccde(il, nde1))
470 dsum = dsum + omwt(il) * dabs(vccde(il, nde1))
476 if ((il .eq. il01) .or. (inbr(il) .eq. 0)) cycle
477 bd(il) = bd(il) * oodsum
486 if ((il .eq. il01) .or. (inbr(il) .eq. 0)) cycle
487 acd = acd + bd(il) * vccde(il, 1)
488 aed = aed + bd(il) * vccde(il, nde2)
492 if (vcmx .lt. vcthresh)
then
493 fatten = vcmx / vcthresh
This module contains simulation constants.
real(dp), parameter dprec
real constant machine precision
real(dp), parameter done
real constant 1
This module defines variable data types.
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 abwts(nnbrmx, nnbr, inbr, il01, nde1, vccde, vcthresh, dl0, dln, acd, add, aed, bd)
Compute "a" and "b" weights for the local connections with respect to the perpendicular direction of ...
subroutine getrot(nnbrmx, nnbr, inbr, vc, il01, rmat, iml0)
Compute the matrix that rotates the model-coordinate axes to the "(c, d, e)-coordinate" axes associat...
subroutine tranvc(nnbrmx, nnbrs, rmat, vc, vccde)
Transform local connection unit-vectors from model coordinates to "(c, d, e)" coordinates associated ...
subroutine abhats(nnbrmx, nnbr, inbr, il01, vc, vn, dl0, dln, ck, vcthresh, allhc, ar01, ahat, bhat)
Compute "ahat" and "bhat" coefficients for one side of an interface.