MODFLOW 6  version 6.7.0.dev0
USGS Modular Hydrologic Model
xt3dmodule Module Reference

Data Types

type  xt3dtype
 

Functions/Subroutines

subroutine, public xt3d_cr (xt3dobj, name_model, inunit, iout, ldispopt)
 Create a new xt3d object. More...
 
subroutine xt3d_df (this, dis)
 Define the xt3d object. More...
 
subroutine xt3d_ac (this, moffset, sparse)
 Add connections for extended neighbors to the sparse matrix. More...
 
subroutine xt3d_mc (this, moffset, matrix_sln)
 Map connections and construct iax, jax, and idxglox. More...
 
subroutine xt3d_ar (this, ibound, k11, ik33, k33, sat, ik22, k22, iangle1, iangle2, iangle3, angle1, angle2, angle3, inewton, icelltype)
 Allocate and Read. More...
 
subroutine xt3d_fc (this, kiter, matrix_sln, idxglo, rhs, hnew)
 Formulate. More...
 
subroutine xt3d_fcpc (this, nodes, lsat)
 Formulate for permanently confined connections and save in amatpc and amatpcx. More...
 
subroutine xt3d_fhfb (this, kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew, n, m, condhfb)
 Formulate HFB correction. More...
 
subroutine xt3d_fn (this, kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew)
 Fill Newton terms for xt3d. More...
 
subroutine xt3d_flowja (this, hnew, flowja)
 Budget. More...
 
subroutine xt3d_flowjahfb (this, n, m, hnew, flowja, condhfb)
 hfb contribution to flowja when xt3d is used More...
 
subroutine xt3d_da (this)
 Deallocate variables. More...
 
subroutine allocate_scalars (this)
 Allocate scalar pointer variables. More...
 
subroutine allocate_arrays (this)
 Allocate xt3d arrays. More...
 
subroutine xt3d_iallpc (this)
 Allocate and populate iallpc array. Set lamatsaved. More...
 
subroutine xt3d_indices (this, n, m, il0, ii01, jjs01, il01, il10, ii00, ii11, ii10)
 Set various indices for XT3D. More...
 
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. More...
 
subroutine xt3d_load_inbr (this, n, nnbr, inbr)
 Load neighbor list for a cell. More...
 
subroutine xt3d_areas (this, nodes, n, m, jjs01, lsat, ar01, ar10, hnew)
 Compute interfacial areas. More...
 
subroutine xt3d_amat_nbrs (this, nodes, n, idiag, nnbr, nja, matrix_sln, inbr, idxglo, chat)
 Add contributions from neighbors to amat. More...
 
subroutine xt3d_amat_nbrnbrs (this, nodes, n, m, ii01, nnbr, nja, matrix_sln, inbr, idxglo, chat)
 Add contributions from neighbors of neighbor to amat. More...
 
subroutine xt3d_amatpc_nbrs (this, nodes, n, idiag, nnbr, inbr, chat)
 Add contributions from neighbors to amatpc. More...
 
subroutine xt3d_amatpcx_nbrnbrs (this, nodes, n, m, ii01, nnbr, inbr, chat)
 Add contributions from neighbors of neighbor to amatpc and amatpcx. More...
 
subroutine xt3d_get_iinm (this, n, m, iinm)
 Get position of n-m connection in ja array (return 0 if not connected) More...
 
subroutine xt3d_get_iinmx (this, n, m, iinmx)
 Get position of n-m "extended connection" in jax array (return 0 if not connected) More...
 
subroutine xt3d_rhs (this, nodes, n, m, nnbr, inbr, chat, hnew, rhs)
 Add contributions to rhs. More...
 
subroutine xt3d_qnbrs (this, nodes, n, m, nnbr, inbr, chat, hnew, qnbrs)
 Add contributions to flow from neighbors. More...
 
subroutine xt3d_fillrmatck (this, n)
 Fill rmat array for cell n. More...
 

Function/Subroutine Documentation

◆ allocate_arrays()

subroutine xt3dmodule::allocate_arrays ( class(xt3dtype this)

Definition at line 1049 of file Xt3dInterface.f90.

1050  ! -- modules
1052  ! -- dummy
1053  class(Xt3dType) :: this
1054  ! -- local
1055  integer(I4B) :: njax
1056  integer(I4B) :: n
1057  !
1058  ! -- Allocate Newton-dependent arrays
1059  if (this%inewton /= 0) then
1060  call mem_allocate(this%qsat, this%dis%nja, 'QSAT', this%memoryPath)
1061  else
1062  call mem_allocate(this%qsat, 0, 'QSAT', this%memoryPath)
1063  end if
1064  !
1065  ! -- If dispersion, set iallpc to 1 otherwise call xt3d_iallpc to go through
1066  ! each connection and mark cells that are permanenntly confined and can
1067  ! have their coefficients precalculated
1068  if (this%ldispersion) then
1069  !
1070  ! -- xt3d is being used for dispersion; all matrix terms are precalculated
1071  ! and used repeatedly until flows change
1072  this%lamatsaved = .true.
1073  call mem_allocate(this%iallpc, this%dis%nodes, 'IALLPC', this%memoryPath)
1074  do n = 1, this%dis%nodes
1075  this%iallpc(n) = 1
1076  end do
1077  else
1078  !
1079  ! -- xt3d is being used for flow so find where connections are
1080  ! permanently confined and precalculate matrix terms case where
1081  ! conductances do not depend on head
1082  call this%xt3d_iallpc()
1083  end if
1084  !
1085  ! -- Allocate space for precalculated matrix terms
1086  if (this%lamatsaved) then
1087  call mem_allocate(this%amatpc, this%dis%nja, 'AMATPC', this%memoryPath)
1088  njax = this%numextnbrs ! + 1
1089  call mem_allocate(this%amatpcx, njax, 'AMATPCX', this%memoryPath)
1090  else
1091  call mem_allocate(this%amatpc, 0, 'AMATPC', this%memoryPath)
1092  call mem_allocate(this%amatpcx, 0, 'AMATPCX', this%memoryPath)
1093  end if
1094  !
1095  ! -- Allocate work arrays
1096  call mem_allocate(this%rmatck, 3, 3, 'RMATCK', this%memoryPath)
1097  !
1098  ! -- Initialize arrays to zero
1099  this%rmatck = dzero
1100  if (this%inewton /= 0) then
1101  this%qsat = dzero
1102  else if (this%lamatsaved) then
1103  this%amatpc = dzero
1104  this%amatpcx = dzero
1105  end if

◆ allocate_scalars()

subroutine xt3dmodule::allocate_scalars ( class(xt3dtype this)

Definition at line 1016 of file Xt3dInterface.f90.

1017  ! -- modules
1019  ! -- dummy
1020  class(Xt3dType) :: this
1021  !
1022  ! -- Allocate scalars
1023  call mem_allocate(this%ixt3d, 'IXT3D', this%memoryPath)
1024  call mem_allocate(this%nbrmax, 'NBRMAX', this%memoryPath)
1025  call mem_allocate(this%inunit, 'INUNIT', this%memoryPath)
1026  call mem_allocate(this%iout, 'IOUT', this%memoryPath)
1027  call mem_allocate(this%inewton, 'INEWTON', this%memoryPath)
1028  call mem_allocate(this%numextnbrs, 'NUMEXTNBRS', this%memoryPath)
1029  call mem_allocate(this%nozee, 'NOZEE', this%memoryPath)
1030  call mem_allocate(this%vcthresh, 'VCTHRESH', this%memoryPath)
1031  call mem_allocate(this%lamatsaved, 'LAMATSAVED', this%memoryPath)
1032  call mem_allocate(this%ldispersion, 'LDISPERSION', this%memoryPath)
1033  !
1034  ! -- Initialize value
1035  this%ixt3d = 0
1036  this%nbrmax = 0
1037  this%inunit = 0
1038  this%iout = 0
1039  this%inewton = 0
1040  this%numextnbrs = 0
1041  this%nozee = .false.
1042  this%vcthresh = 1.d-10
1043  this%lamatsaved = .false.
1044  this%ldispersion = .false.

◆ xt3d_ac()

subroutine xt3dmodule::xt3d_ac ( class(xt3dtype this,
integer(i4b), intent(in)  moffset,
type(sparsematrix), intent(inout)  sparse 
)

Definition at line 124 of file Xt3dInterface.f90.

125  ! -- modules
126  use sparsemodule, only: sparsematrix
128  ! -- dummy
129  class(Xt3dType) :: this
130  integer(I4B), intent(in) :: moffset
131  type(sparsematrix), intent(inout) :: sparse
132  ! -- local
133  type(sparsematrix) :: sparse_xt3d
134  integer(I4B) :: i, j, k, jj, kk, iglo, kglo, iadded
135  integer(I4B) :: nnz
136  integer(I4B) :: ierror
137  !
138  ! -- If not rhs, add connections
139  if (this%ixt3d == 1) then
140  ! -- Assume nnz is 19, which is an approximate value
141  ! based on a 3d structured grid
142  nnz = 19
143  call sparse_xt3d%init(this%dis%nodes, this%dis%nodes, nnz)
144  !
145  ! -- Loop over nodes and store extended xt3d neighbors
146  ! temporarily in sparse_xt3d; this will be converted to
147  ! ia_xt3d and ja_xt3d next
148  do i = 1, this%dis%nodes
149  iglo = i + moffset
150  ! -- loop over neighbors
151  do jj = this%dis%con%ia(i) + 1, this%dis%con%ia(i + 1) - 1
152  j = this%dis%con%ja(jj)
153  ! -- loop over neighbors of neighbors
154  do kk = this%dis%con%ia(j) + 1, this%dis%con%ia(j + 1) - 1
155  k = this%dis%con%ja(kk)
156  kglo = k + moffset
157  call sparse_xt3d%addconnection(i, k, 1)
158  end do
159  end do
160  end do
161  !
162  ! -- calculate ia_xt3d and ja_xt3d from sparse_xt3d and
163  ! then destroy sparse
164  call mem_allocate(this%ia_xt3d, this%dis%nodes + 1, 'IA_XT3D', &
165  trim(this%memoryPath))
166  call mem_allocate(this%ja_xt3d, sparse_xt3d%nnz, 'JA_XT3D', &
167  trim(this%memoryPath))
168  call sparse_xt3d%filliaja(this%ia_xt3d, this%ja_xt3d, ierror)
169  call sparse_xt3d%destroy()
170  !
171  ! -- add extended neighbors to sparse and count number of
172  ! extended neighbors
173  do i = 1, this%dis%nodes
174  iglo = i + moffset
175  do kk = this%ia_xt3d(i), this%ia_xt3d(i + 1) - 1
176  k = this%ja_xt3d(kk)
177  kglo = k + moffset
178  call sparse%addconnection(iglo, kglo, 1, iadded)
179  this%numextnbrs = this%numextnbrs + 1
180  end do
181  end do
182  else
183  ! -- Arrays not needed; set to size zero
184  call mem_allocate(this%ia_xt3d, 0, 'IA_XT3D', trim(this%memoryPath))
185  call mem_allocate(this%ja_xt3d, 0, 'JA_XT3D', trim(this%memoryPath))
186  end if

◆ xt3d_amat_nbrnbrs()

subroutine xt3dmodule::xt3d_amat_nbrnbrs ( class(xt3dtype this,
integer(i4b), intent(in)  nodes,
integer(i4b)  n,
integer(i4b)  m,
integer(i4b)  ii01,
integer(i4b)  nnbr,
integer(i4b)  nja,
class(matrixbasetype), pointer  matrix_sln,
integer(i4b), dimension(this%nbrmax)  inbr,
integer(i4b), dimension(nja), intent(in)  idxglo,
real(dp), dimension(this%nbrmax)  chat 
)

Definition at line 1406 of file Xt3dInterface.f90.

1408  ! -- dummy
1409  class(Xt3dType) :: this
1410  integer(I4B), intent(in) :: nodes
1411  integer(I4B) :: n, m, ii01, nnbr, nja
1412  class(MatrixBaseType), pointer :: matrix_sln
1413  integer(I4B), dimension(this%nbrmax) :: inbr
1414  integer(I4B), intent(in), dimension(nja) :: idxglo
1415  real(DP), dimension(this%nbrmax) :: chat
1416  ! -- local
1417  integer(I4B) :: iil, iii, jjj, iixjjj, iijjj
1418  !
1419  do iil = 1, nnbr
1420  if (inbr(iil) .ne. 0) then
1421  call matrix_sln%add_value_pos(idxglo(ii01), chat(iil))
1422  iii = this%dis%con%ia(m) + iil
1423  jjj = this%dis%con%ja(iii)
1424  call this%xt3d_get_iinmx(n, jjj, iixjjj)
1425  if (iixjjj .ne. 0) then
1426  call matrix_sln%add_value_pos(this%idxglox(iixjjj), -chat(iil))
1427  else
1428  call this%xt3d_get_iinm(n, jjj, iijjj)
1429  call matrix_sln%add_value_pos(idxglo(iijjj), -chat(iil))
1430  end if
1431  end if
1432  end do

◆ xt3d_amat_nbrs()

subroutine xt3dmodule::xt3d_amat_nbrs ( class(xt3dtype this,
integer(i4b), intent(in)  nodes,
integer(i4b)  n,
integer(i4b)  idiag,
integer(i4b)  nnbr,
integer(i4b)  nja,
class(matrixbasetype), pointer  matrix_sln,
integer(i4b), dimension(this%nbrmax)  inbr,
integer(i4b), dimension(nja), intent(in)  idxglo,
real(dp), dimension(this%nbrmax)  chat 
)

Definition at line 1382 of file Xt3dInterface.f90.

1384  ! -- dummy
1385  class(Xt3dType) :: this
1386  integer(I4B), intent(in) :: nodes
1387  integer(I4B) :: n, idiag, nnbr, nja
1388  class(MatrixBaseType), pointer :: matrix_sln
1389  integer(I4B), dimension(this%nbrmax) :: inbr
1390  integer(I4B), intent(in), dimension(nja) :: idxglo
1391  real(DP), dimension(this%nbrmax) :: chat
1392  ! -- local
1393  integer(I4B) :: iil, iii
1394  !
1395  do iil = 1, nnbr
1396  if (inbr(iil) .ne. 0) then
1397  iii = this%dis%con%ia(n) + iil
1398  call matrix_sln%add_value_pos(idxglo(idiag), -chat(iil))
1399  call matrix_sln%add_value_pos(idxglo(iii), chat(iil))
1400  end if
1401  end do

◆ xt3d_amatpc_nbrs()

subroutine xt3dmodule::xt3d_amatpc_nbrs ( class(xt3dtype this,
integer(i4b), intent(in)  nodes,
integer(i4b)  n,
integer(i4b)  idiag,
integer(i4b)  nnbr,
integer(i4b), dimension(this%nbrmax)  inbr,
real(dp), dimension(this%nbrmax)  chat 
)

Definition at line 1437 of file Xt3dInterface.f90.

1438  ! -- dummy
1439  class(Xt3dType) :: this
1440  integer(I4B), intent(in) :: nodes
1441  integer(I4B) :: n, idiag, nnbr
1442  integer(I4B), dimension(this%nbrmax) :: inbr
1443  real(DP), dimension(this%nbrmax) :: chat
1444  ! -- local
1445  integer(I4B) :: iil, iii
1446  !
1447  do iil = 1, nnbr
1448  iii = this%dis%con%ia(n) + iil
1449  this%amatpc(idiag) = this%amatpc(idiag) - chat(iil)
1450  this%amatpc(iii) = this%amatpc(iii) + chat(iil)
1451  end do

◆ xt3d_amatpcx_nbrnbrs()

subroutine xt3dmodule::xt3d_amatpcx_nbrnbrs ( class(xt3dtype this,
integer(i4b), intent(in)  nodes,
integer(i4b)  n,
integer(i4b)  m,
integer(i4b)  ii01,
integer(i4b)  nnbr,
integer(i4b), dimension(this%nbrmax)  inbr,
real(dp), dimension(this%nbrmax)  chat 
)

Definition at line 1456 of file Xt3dInterface.f90.

1457  ! -- dummy
1458  class(Xt3dType) :: this
1459  integer(I4B), intent(in) :: nodes
1460  integer(I4B) :: n, m, ii01, nnbr
1461  integer(I4B), dimension(this%nbrmax) :: inbr
1462  real(DP), dimension(this%nbrmax) :: chat
1463  ! -- local
1464  integer(I4B) :: iil, iii, jjj, iixjjj, iijjj
1465  !
1466  do iil = 1, nnbr
1467  this%amatpc(ii01) = this%amatpc(ii01) + chat(iil)
1468  iii = this%dis%con%ia(m) + iil
1469  jjj = this%dis%con%ja(iii)
1470  call this%xt3d_get_iinmx(n, jjj, iixjjj)
1471  if (iixjjj .ne. 0) then
1472  this%amatpcx(iixjjj) = this%amatpcx(iixjjj) - chat(iil)
1473  else
1474  call this%xt3d_get_iinm(n, jjj, iijjj)
1475  this%amatpc(iijjj) = this%amatpc(iijjj) - chat(iil)
1476  end if
1477  end do

◆ xt3d_ar()

subroutine xt3dmodule::xt3d_ar ( class(xt3dtype this,
integer(i4b), dimension(:), intent(inout), pointer, contiguous  ibound,
real(dp), dimension(:), intent(in), pointer, contiguous  k11,
integer(i4b), intent(in), pointer  ik33,
real(dp), dimension(:), intent(in), pointer, contiguous  k33,
real(dp), dimension(:), intent(in), pointer, contiguous  sat,
integer(i4b), intent(in), pointer  ik22,
real(dp), dimension(:), intent(in), pointer, contiguous  k22,
integer(i4b), intent(in), pointer  iangle1,
integer(i4b), intent(in), pointer  iangle2,
integer(i4b), intent(in), pointer  iangle3,
real(dp), dimension(:), intent(in), pointer, contiguous  angle1,
real(dp), dimension(:), intent(in), pointer, contiguous  angle2,
real(dp), dimension(:), intent(in), pointer, contiguous  angle3,
integer(i4b), intent(in), optional, pointer  inewton,
integer(i4b), dimension(:), intent(in), optional, pointer, contiguous  icelltype 
)

Definition at line 275 of file Xt3dInterface.f90.

277  ! -- modules
278  use simmodule, only: store_error
279  ! -- dummy
280  class(Xt3dType) :: this
281  integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: ibound
282  real(DP), dimension(:), intent(in), pointer, contiguous :: k11
283  integer(I4B), intent(in), pointer :: ik33
284  real(DP), dimension(:), intent(in), pointer, contiguous :: k33
285  real(DP), dimension(:), intent(in), pointer, contiguous :: sat
286  integer(I4B), intent(in), pointer :: ik22
287  real(DP), dimension(:), intent(in), pointer, contiguous :: k22
288  integer(I4B), intent(in), pointer :: iangle1
289  integer(I4B), intent(in), pointer :: iangle2
290  integer(I4B), intent(in), pointer :: iangle3
291  real(DP), dimension(:), intent(in), pointer, contiguous :: angle1
292  real(DP), dimension(:), intent(in), pointer, contiguous :: angle2
293  real(DP), dimension(:), intent(in), pointer, contiguous :: angle3
294  integer(I4B), intent(in), pointer, optional :: inewton
295  integer(I4B), dimension(:), intent(in), pointer, &
296  contiguous, optional :: icelltype
297  ! -- local
298  integer(I4B) :: n, nnbrs
299  ! -- formats
300  character(len=*), parameter :: fmtheader = &
301  "(1x, /1x, 'XT3D is active.'//)"
302  !
303  ! -- Print a message identifying the xt3d module.
304  if (this%iout > 0) then
305  write (this%iout, fmtheader)
306  end if
307  !
308  ! -- Store pointers to arguments that were passed in
309  this%ibound => ibound
310  this%k11 => k11
311  this%ik33 => ik33
312  this%k33 => k33
313  this%sat => sat
314  this%ik22 => ik22
315  this%k22 => k22
316  this%iangle1 => iangle1
317  this%iangle2 => iangle2
318  this%iangle3 => iangle3
319  this%angle1 => angle1
320  this%angle2 => angle2
321  this%angle3 => angle3
322  !
323  if (present(inewton)) then
324  ! -- inewton is not needed for transport so it's optional.
325  this%inewton = inewton
326  end if
327  if (present(icelltype)) then
328  ! -- icelltype is not needed for transport, so it's optional.
329  ! It is only needed to determine if cell connections are permanently
330  ! confined, which means that some matrix terms can be precalculated
331  this%icelltype => icelltype
332  end if
333  !
334  ! -- If angle1 and angle2 were not specified, then there is no z
335  ! component in the xt3d formulation for horizontal connections.
336  if (this%iangle2 == 0) this%nozee = .true.
337  !
338  ! -- Determine the maximum number of neighbors for any cell.
339  this%nbrmax = 0
340  do n = 1, this%dis%nodes
341  nnbrs = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
342  this%nbrmax = max(nnbrs, this%nbrmax)
343  end do
344  !
345  ! -- Check to make sure dis package can calculate connection direction info
346  if (this%dis%icondir == 0) then
347  call store_error('Vertices not specified for discretization package, '// &
348  'but XT3D is active: '//trim(adjustl(this%memoryPath))// &
349  '. Vertices must be specified in discretization '// &
350  'package in order to use XT3D.', terminate=.true.)
351  end if
352  !
353  ! -- Check to make sure ANGLEDEGX is available for interface normals
354  if (this%dis%con%ianglex == 0) then
355  call store_error('ANGLDEGX is not specified in the DIS package, '// &
356  'but XT3D is active: '//trim(adjustl(this%memoryPath))// &
357  '. ANGLDEGX must be provided in discretization '// &
358  'package in order to use XT3D.', terminate=.true.)
359  end if
360  !
361  ! -- allocate arrays
362  call this%allocate_arrays()
363  !
364  ! -- If not Newton and not rhs, precalculate amatpc and amatpcx for
365  ! -- permanently confined connections
366  if (this%lamatsaved .and. .not. this%ldispersion) &
367  call this%xt3d_fcpc(this%dis%nodes, .true.)
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
Here is the call graph for this function:

◆ xt3d_areas()

subroutine xt3dmodule::xt3d_areas ( class(xt3dtype this,
integer(i4b)  nodes,
integer(i4b)  n,
integer(i4b)  m,
integer(i4b)  jjs01,
logical  lsat,
real(dp)  ar01,
real(dp)  ar10,
real(dp), dimension(:), intent(inout), optional  hnew 
)

Definition at line 1300 of file Xt3dInterface.f90.

1301  ! -- dummy
1302  class(Xt3dType) :: this
1303  logical :: lsat
1304  integer(I4B) :: nodes, n, m, jjs01
1305  real(DP) :: ar01, ar10
1306  real(DP), intent(inout), dimension(:), optional :: hnew
1307  ! -- local
1308  real(DP) :: topn, botn, topm, botm, thksatn, thksatm
1309  real(DP) :: sill_top, sill_bot, tpn, tpm
1310  real(DP) :: satups
1311  !
1312  ! -- Compute area depending on connection type
1313  if (this%dis%con%ihc(jjs01) == 0) then
1314  !
1315  ! -- Vertical connection
1316  ar01 = this%dis%con%hwva(jjs01)
1317  ar10 = ar01
1318  else if (this%inewton /= 0) then
1319  !
1320  ! -- If Newton horizontal connection
1321  if (lsat) then
1322  !
1323  ! -- lsat true means use full saturation
1324  topn = this%dis%top(n)
1325  botn = this%dis%bot(n)
1326  topm = this%dis%top(m)
1327  botm = this%dis%bot(m)
1328  thksatn = topn - botn
1329  thksatm = topm - botm
1330  if (this%dis%con%ihc(jjs01) .eq. 2) then
1331  ! -- Vertically staggered
1332  sill_top = min(topn, topm)
1333  sill_bot = max(botn, botm)
1334  tpn = botn + thksatn
1335  tpm = botm + thksatm
1336  thksatn = max(min(tpn, sill_top) - sill_bot, dzero)
1337  thksatm = max(min(tpm, sill_top) - sill_bot, dzero)
1338  end if
1339  ar01 = this%dis%con%hwva(jjs01) * dhalf * (thksatn + thksatm)
1340  else
1341  ! -- If Newton and lsat=.false., it is assumed that the fully saturated
1342  ! areas have already been calculated and are being passed in through
1343  ! ar01 and ar10. The actual areas are obtained simply by scaling by
1344  ! the upstream saturation.
1345  if (hnew(m) < hnew(n)) then
1346  satups = this%sat(n)
1347  else
1348  satups = this%sat(m)
1349  end if
1350  ar01 = ar01 * satups
1351  end if
1352  ar10 = ar01
1353  else
1354  !
1355  ! -- Non-Newton horizontal connection
1356  topn = this%dis%top(n)
1357  botn = this%dis%bot(n)
1358  topm = this%dis%top(m)
1359  botm = this%dis%bot(m)
1360  thksatn = topn - botn
1361  thksatm = topm - botm
1362  if (.not. lsat) then
1363  thksatn = this%sat(n) * thksatn
1364  thksatm = this%sat(m) * thksatm
1365  end if
1366  if (this%dis%con%ihc(jjs01) == 2) then
1367  ! -- Vertically staggered
1368  sill_top = min(topn, topm)
1369  sill_bot = max(botn, botm)
1370  tpn = botn + thksatn
1371  tpm = botm + thksatm
1372  thksatn = max(min(tpn, sill_top) - sill_bot, dzero)
1373  thksatm = max(min(tpm, sill_top) - sill_bot, dzero)
1374  end if
1375  ar01 = this%dis%con%hwva(jjs01) * thksatn
1376  ar10 = this%dis%con%hwva(jjs01) * thksatm
1377  end if

◆ xt3d_cr()

subroutine, public xt3dmodule::xt3d_cr ( type(xt3dtype), pointer  xt3dobj,
character(len=*), intent(in)  name_model,
integer(i4b), intent(in)  inunit,
integer(i4b), intent(in)  iout,
logical, intent(in), optional  ldispopt 
)

Definition at line 89 of file Xt3dInterface.f90.

90  ! -- dummy
91  type(Xt3dType), pointer :: xt3dobj
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
96  !
97  ! -- Create the object
98  allocate (xt3dobj)
99  !
100  ! -- Assign the memory path
101  xt3dobj%memoryPath = create_mem_path(name_model, 'XT3D')
102  !
103  ! -- Allocate scalars
104  call xt3dobj%allocate_scalars()
105  !
106  ! -- Set variables
107  xt3dobj%inunit = inunit
108  xt3dobj%iout = iout
109  if (present(ldispopt)) xt3dobj%ldispersion = ldispopt
Here is the call graph for this function:
Here is the caller graph for this function:

◆ xt3d_da()

subroutine xt3dmodule::xt3d_da ( class(xt3dtype this)

Definition at line 981 of file Xt3dInterface.f90.

982  ! -- modules
984  ! -- dummy
985  class(Xt3dType) :: this
986  !
987  ! -- Deallocate arrays
988  if (this%ixt3d /= 0) then
989  call mem_deallocate(this%iax)
990  call mem_deallocate(this%jax)
991  call mem_deallocate(this%idxglox)
992  call mem_deallocate(this%ia_xt3d)
993  call mem_deallocate(this%ja_xt3d)
994  call mem_deallocate(this%rmatck)
995  call mem_deallocate(this%qsat)
996  call mem_deallocate(this%amatpc)
997  call mem_deallocate(this%amatpcx)
998  call mem_deallocate(this%iallpc)
999  end if
1000  !
1001  ! -- Scalars
1002  call mem_deallocate(this%ixt3d)
1003  call mem_deallocate(this%inunit)
1004  call mem_deallocate(this%iout)
1005  call mem_deallocate(this%inewton)
1006  call mem_deallocate(this%numextnbrs)
1007  call mem_deallocate(this%nozee)
1008  call mem_deallocate(this%vcthresh)
1009  call mem_deallocate(this%lamatsaved)
1010  call mem_deallocate(this%nbrmax)
1011  call mem_deallocate(this%ldispersion)

◆ xt3d_df()

subroutine xt3dmodule::xt3d_df ( class(xt3dtype this,
class(disbasetype), intent(inout), pointer  dis 
)

Definition at line 114 of file Xt3dInterface.f90.

115  ! -- dummy
116  class(Xt3dType) :: this
117  class(DisBaseType), pointer, intent(inout) :: dis
118  !
119  this%dis => dis

◆ xt3d_fc()

subroutine xt3dmodule::xt3d_fc ( class(xt3dtype this,
integer(i4b)  kiter,
class(matrixbasetype), pointer  matrix_sln,
integer(i4b), dimension(:), intent(in)  idxglo,
real(dp), dimension(:), intent(inout)  rhs,
real(dp), dimension(:), intent(inout)  hnew 
)

Definition at line 371 of file Xt3dInterface.f90.

372  ! -- modules
373  use constantsmodule, only: done
374  use xt3dalgorithmmodule, only: qconds
375  ! -- dummy
376  class(Xt3dType) :: this
377  integer(I4B) :: kiter
378  class(MatrixBaseType), pointer :: matrix_sln
379  integer(I4B), intent(in), dimension(:) :: idxglo
380  real(DP), intent(inout), dimension(:) :: rhs
381  real(DP), intent(inout), dimension(:) :: hnew
382  ! -- local
383  integer(I4B) :: nodes, nja
384  integer(I4B) :: n, m, ipos
385  logical :: allhc0, allhc1
386  integer(I4B) :: nnbr0, nnbr1
387  integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
388  integer(I4B) :: i
389  integer(I4B), dimension(this%nbrmax) :: inbr0, inbr1
390  real(DP) :: ar01, ar10
391  real(DP), dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
392  real(DP), dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
393  real(DP), dimension(3, 3) :: ck0, ck1
394  real(DP) :: chat01
395  real(DP), dimension(this%nbrmax) :: chati0, chat1j
396  real(DP) :: qnm, qnbrs
397  !
398  ! -- Calculate xt3d conductance-like coefficients and put into amat and rhs
399  ! -- as appropriate
400  !
401  nodes = this%dis%nodes
402  nja = this%dis%con%nja
403  if (this%lamatsaved) then
404  do i = 1, this%dis%con%nja
405  call matrix_sln%add_value_pos(idxglo(i), this%amatpc(i))
406  end do
407  do i = 1, this%numextnbrs
408  call matrix_sln%add_value_pos(this%idxglox(i), this%amatpcx(i))
409  end do
410  end if
411  !
412  do n = 1, nodes
413  ! -- Skip if inactive.
414  if (this%ibound(n) .eq. 0) cycle
415  ! -- Skip if all connections are permanently confined
416  if (this%lamatsaved) then
417  if (this%iallpc(n) == 1) cycle
418  end if
419  nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
420  ! -- Load conductivity and connection info for cell 0.
421  call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
422  ck0, allhc0)
423  ! -- Loop over active neighbors of cell 0 that have a higher
424  ! cell number (taking advantage of reciprocity).
425  do il0 = 1, nnbr0
426  ipos = this%dis%con%ia(n) + il0
427  if (this%dis%con%mask(ipos) == 0) cycle
428 
429  m = inbr0(il0)
430  ! -- Skip if neighbor is inactive or has lower cell number.
431  if ((m .eq. 0) .or. (m .lt. n)) cycle
432  nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
433  ! -- Load conductivity and connection info for cell 1.
434  call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
435  ck1, allhc1)
436  ! -- Set various indices.
437  call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
438  ii00, ii11, ii10)
439  ! -- Compute areas.
440  if (this%inewton /= 0) then
441  ar01 = done
442  ar10 = done
443  else
444  call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
445  end if
446  ! -- Compute "conductances" for interface between
447  ! cells 0 and 1.
448  call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, ck0, &
449  nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
450  this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
451  ! -- If Newton, compute and save saturated flow, then scale
452  ! conductance-like coefficients by the actual area for
453  ! subsequent amat and rhs assembly.
454  if (this%inewton /= 0) then
455  ! -- Contribution to flow from primary connection.
456  qnm = chat01 * (hnew(m) - hnew(n))
457  ! -- Contribution from immediate neighbors of node 0.
458  call this%xt3d_qnbrs(nodes, n, m, nnbr0, inbr0, chati0, hnew, qnbrs)
459  qnm = qnm + qnbrs
460  ! -- Contribution from immediate neighbors of node 1.
461  call this%xt3d_qnbrs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, qnbrs)
462  qnm = qnm - qnbrs
463  ! -- Multiply by saturated area and save in qsat.
464  call this%xt3d_areas(nodes, n, m, jjs01, .true., ar01, ar10, hnew)
465  this%qsat(ii01) = qnm * ar01
466  ! -- Scale coefficients by actual area.
467  call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
468  chat01 = chat01 * ar01
469  chati0 = chati0 * ar01
470  chat1j = chat1j * ar01
471  end if
472  !
473  ! -- Contribute to rows for cells 0 and 1.
474  call matrix_sln%add_value_pos(idxglo(ii00), -chat01)
475  call matrix_sln%add_value_pos(idxglo(ii01), chat01)
476  call matrix_sln%add_value_pos(idxglo(ii11), -chat01)
477  call matrix_sln%add_value_pos(idxglo(ii10), chat01)
478  if (this%ixt3d == 1) then
479  call this%xt3d_amat_nbrs(nodes, n, ii00, nnbr0, nja, matrix_sln, &
480  inbr0, idxglo, chati0)
481  call this%xt3d_amat_nbrnbrs(nodes, n, m, ii01, nnbr1, nja, matrix_sln, &
482  inbr1, idxglo, chat1j)
483  call this%xt3d_amat_nbrs(nodes, m, ii11, nnbr1, nja, matrix_sln, &
484  inbr1, idxglo, chat1j)
485  call this%xt3d_amat_nbrnbrs(nodes, m, n, ii10, nnbr0, nja, matrix_sln, &
486  inbr0, idxglo, chati0)
487  else
488  call this%xt3d_rhs(nodes, n, m, nnbr0, inbr0, chati0, hnew, rhs)
489  call this%xt3d_rhs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, rhs)
490  end if
491  !
492  end do
493  end do
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
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)....
Here is the call graph for this function:

◆ xt3d_fcpc()

subroutine xt3dmodule::xt3d_fcpc ( class(xt3dtype this,
integer(i4b), intent(in)  nodes,
logical, intent(in)  lsat 
)
Parameters
[in]lsatif true, then calculations made with saturated areas (should be false for solute dispersion; should be true for heat)

Definition at line 499 of file Xt3dInterface.f90.

500  ! -- modules
501  use xt3dalgorithmmodule, only: qconds
502  ! -- dummy
503  class(Xt3dType) :: this
504  integer(I4B), intent(in) :: nodes
505  logical, intent(in) :: lsat !< if true, then calculations made with saturated areas (should be false for solute dispersion; should be true for heat)
506  ! -- local
507  integer(I4B) :: n, m, ipos
508  !
509  logical :: allhc0, allhc1
510  integer(I4B) :: nnbr0, nnbr1
511  integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
512  integer(I4B), dimension(this%nbrmax) :: inbr0, inbr1
513  real(DP) :: ar01, ar10
514  real(DP), dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
515  real(DP), dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
516  real(DP), dimension(3, 3) :: ck0, ck1
517  real(DP) :: chat01
518  real(DP), dimension(this%nbrmax) :: chati0, chat1j
519  !
520  ! -- Initialize amatpc and amatpcx to zero
521  do n = 1, size(this%amatpc)
522  this%amatpc(n) = dzero
523  end do
524  do n = 1, size(this%amatpcx)
525  this%amatpcx(n) = dzero
526  end do
527  !
528  ! -- Calculate xt3d conductance-like coefficients for permanently confined
529  ! -- connections and put into amatpc and amatpcx as appropriate
530  do n = 1, nodes
531  ! -- Skip if not iallpc.
532  if (this%iallpc(n) == 0) cycle
533  if (this%ibound(n) == 0) cycle
534  nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
535  ! -- Load conductivity and connection info for cell 0.
536  call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
537  ck0, allhc0)
538  ! -- Loop over active neighbors of cell 0 that have a higher
539  ! -- cell number (taking advantage of reciprocity).
540  do il0 = 1, nnbr0
541  ipos = this%dis%con%ia(n) + il0
542  if (this%dis%con%mask(ipos) == 0) cycle
543 
544  m = inbr0(il0)
545  ! -- Skip if neighbor has lower cell number.
546  if (m .lt. n) cycle
547  nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
548  ! -- Load conductivity and connection info for cell 1.
549  call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
550  ck1, allhc1)
551  ! -- Set various indices.
552  call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
553  ii00, ii11, ii10)
554  ! -- Compute confined areas.
555  call this%xt3d_areas(nodes, n, m, jjs01, lsat, ar01, ar10)
556  ! -- Compute "conductances" for interface between
557  ! -- cells 0 and 1.
558  call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, ck0, &
559  nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
560  this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
561  ! -- Contribute to rows for cells 0 and 1.
562  this%amatpc(ii00) = this%amatpc(ii00) - chat01
563  this%amatpc(ii01) = this%amatpc(ii01) + chat01
564  this%amatpc(ii11) = this%amatpc(ii11) - chat01
565  this%amatpc(ii10) = this%amatpc(ii10) + chat01
566  call this%xt3d_amatpc_nbrs(nodes, n, ii00, nnbr0, inbr0, chati0)
567  call this%xt3d_amatpcx_nbrnbrs(nodes, n, m, ii01, nnbr1, inbr1, chat1j)
568  call this%xt3d_amatpc_nbrs(nodes, m, ii11, nnbr1, inbr1, chat1j)
569  call this%xt3d_amatpcx_nbrnbrs(nodes, m, n, ii10, nnbr0, inbr0, chati0)
570  end do
571  end do
Here is the call graph for this function:

◆ xt3d_fhfb()

subroutine xt3dmodule::xt3d_fhfb ( class(xt3dtype this,
integer(i4b)  kiter,
integer(i4b), intent(in)  nodes,
integer(i4b), intent(in)  nja,
class(matrixbasetype), pointer  matrix_sln,
integer(i4b), dimension(nja), intent(in)  idxglo,
real(dp), dimension(nodes), intent(inout)  rhs,
real(dp), dimension(nodes), intent(inout)  hnew,
integer(i4b)  n,
integer(i4b)  m,
real(dp)  condhfb 
)

Definition at line 576 of file Xt3dInterface.f90.

578  ! -- modules
579  use constantsmodule, only: done
580  use xt3dalgorithmmodule, only: qconds
581  ! -- dummy
582  class(Xt3dType) :: this
583  integer(I4B) :: kiter
584  integer(I4B), intent(in) :: nodes
585  integer(I4B), intent(in) :: nja
586  class(MatrixBaseType), pointer :: matrix_sln
587  integer(I4B) :: n, m
588  integer(I4B), intent(in), dimension(nja) :: idxglo
589  real(DP), intent(inout), dimension(nodes) :: rhs
590  real(DP), intent(inout), dimension(nodes) :: hnew
591  real(DP) :: condhfb
592  ! -- local
593  logical :: allhc0, allhc1
594  integer(I4B) :: nnbr0, nnbr1
595  integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10, il
596  integer(I4B), dimension(this%nbrmax) :: inbr0, inbr1
597  real(DP) :: ar01, ar10
598  real(DP), dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
599  real(DP), dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
600  real(DP), dimension(3, 3) :: ck0, ck1
601  real(DP) :: chat01
602  real(DP), dimension(this%nbrmax) :: chati0, chat1j
603  real(DP) :: qnm, qnbrs
604  real(DP) :: term
605  !
606  ! -- Calculate hfb corrections to xt3d conductance-like coefficients and
607  ! put into amat and rhs as appropriate
608  !
609  nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
610  ! -- Load conductivity and connection info for cell 0.
611  call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
612  ck0, allhc0)
613  ! -- Find local neighbor number of cell 1.
614  do il = 1, nnbr0
615  if (inbr0(il) .eq. m) then
616  il0 = il
617  exit
618  end if
619  end do
620  nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
621  ! -- Load conductivity and connection info for cell 1.
622  call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
623  ck1, allhc1)
624  ! -- Set various indices.
625  call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
626  ii00, ii11, ii10)
627  ! -- Compute areas.
628  if (this%inewton /= 0) then
629  ar01 = done
630  ar10 = done
631  else
632  call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
633  end if
634  !
635  ! -- Compute "conductances" for interface between
636  ! cells 0 and 1.
637  call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, &
638  ck0, nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
639  this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
640  !
641  ! -- Apply scale factor to compute "conductances" for hfb correction
642  if (condhfb > dzero) then
643  term = chat01 / (chat01 + condhfb)
644  else
645  term = -condhfb
646  end if
647  chat01 = -chat01 * term
648  chati0 = -chati0 * term
649  chat1j = -chat1j * term
650  !
651  ! -- If Newton, compute and save saturated flow, then scale conductance-like
652  ! coefficients by the actual area for subsequent amat and rhs assembly.
653  if (this%inewton /= 0) then
654  ! -- Contribution to flow from primary connection.
655  qnm = chat01 * (hnew(m) - hnew(n))
656  ! -- Contribution from immediate neighbors of node 0.
657  call this%xt3d_qnbrs(nodes, n, m, nnbr0, inbr0, chati0, hnew, qnbrs)
658  qnm = qnm + qnbrs
659  ! -- Contribution from immediate neighbors of node 1.
660  call this%xt3d_qnbrs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, qnbrs)
661  qnm = qnm - qnbrs
662  ! -- Multiply by saturated area and add correction to qsat.
663  call this%xt3d_areas(nodes, n, m, jjs01, .true., ar01, ar10, hnew)
664  this%qsat(ii01) = this%qsat(ii01) + qnm * ar01
665  ! -- Scale coefficients by actual area.
666  call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
667  chat01 = chat01 * ar01
668  chati0 = chati0 * ar01
669  chat1j = chat1j * ar01
670  end if
671  !
672  ! -- Contribute to rows for cells 0 and 1.
673  call matrix_sln%add_value_pos(idxglo(ii00), -chat01)
674  call matrix_sln%add_value_pos(idxglo(ii01), chat01)
675  call matrix_sln%add_value_pos(idxglo(ii11), -chat01)
676  call matrix_sln%add_value_pos(idxglo(ii10), chat01)
677  if (this%ixt3d == 1) then
678  call this%xt3d_amat_nbrs(nodes, n, ii00, nnbr0, nja, matrix_sln, &
679  inbr0, idxglo, chati0)
680  call this%xt3d_amat_nbrnbrs(nodes, n, m, ii01, nnbr1, nja, matrix_sln, &
681  inbr1, idxglo, chat1j)
682  call this%xt3d_amat_nbrs(nodes, m, ii11, nnbr1, nja, matrix_sln, &
683  inbr1, idxglo, chat1j)
684  call this%xt3d_amat_nbrnbrs(nodes, m, n, ii10, nnbr0, nja, matrix_sln, &
685  inbr0, idxglo, chati0)
686  else
687  call this%xt3d_rhs(nodes, n, m, nnbr0, inbr0, chati0, hnew, rhs)
688  call this%xt3d_rhs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, rhs)
689  end if
Here is the call graph for this function:

◆ xt3d_fillrmatck()

subroutine xt3dmodule::xt3d_fillrmatck ( class(xt3dtype this,
integer(i4b), intent(in)  n 
)

angle1, 2, and 3 must be in radians.

Definition at line 1577 of file Xt3dInterface.f90.

1578  ! -- dummy
1579  class(Xt3dType) :: this
1580  integer(I4B), intent(in) :: n
1581  ! -- local
1582  real(DP) :: ang1, ang2, ang3, ang2d, ang3d
1583  real(DP) :: s1, c1, s2, c2, s3, c3
1584  !
1585  if (this%nozee) then
1586  ang2d = 0d0
1587  ang3d = 0d0
1588  ang1 = this%angle1(n)
1589  ang2 = 0d0
1590  ang3 = 0d0
1591  else
1592  ang1 = this%angle1(n)
1593  ang2 = this%angle2(n)
1594  ang3 = this%angle3(n)
1595  end if
1596  s1 = sin(ang1)
1597  c1 = cos(ang1)
1598  s2 = sin(ang2)
1599  c2 = cos(ang2)
1600  s3 = sin(ang3)
1601  c3 = cos(ang3)
1602  this%rmatck(1, 1) = c1 * c2
1603  this%rmatck(1, 2) = c1 * s2 * s3 - s1 * c3
1604  this%rmatck(1, 3) = -c1 * s2 * c3 - s1 * s3
1605  this%rmatck(2, 1) = s1 * c2
1606  this%rmatck(2, 2) = s1 * s2 * s3 + c1 * c3
1607  this%rmatck(2, 3) = -s1 * s2 * c3 + c1 * s3
1608  this%rmatck(3, 1) = s2
1609  this%rmatck(3, 2) = -c2 * s3
1610  this%rmatck(3, 3) = c2 * c3

◆ xt3d_flowja()

subroutine xt3dmodule::xt3d_flowja ( class(xt3dtype this,
real(dp), dimension(:), intent(inout)  hnew,
real(dp), dimension(:), intent(inout)  flowja 
)

Definition at line 797 of file Xt3dInterface.f90.

798  ! -- modules
799  use xt3dalgorithmmodule, only: qconds
800  ! -- dummy
801  class(Xt3dType) :: this
802  real(DP), intent(inout), dimension(:) :: hnew
803  real(DP), intent(inout), dimension(:) :: flowja
804  ! -- local
805  integer(I4B) :: n, ipos, m, nodes
806  real(DP) :: qnm, qnbrs
807  logical :: allhc0, allhc1
808  integer(I4B) :: nnbr0, nnbr1
809  integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
810  integer(I4B), dimension(this%nbrmax) :: inbr0, inbr1
811  real(DP) :: ar01, ar10
812  real(DP), dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
813  real(DP), dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
814  real(DP), dimension(3, 3) :: ck0, ck1
815  real(DP) :: chat01
816  real(DP), dimension(this%nbrmax) :: chati0, chat1j
817  !
818  ! -- Calculate the flow across each cell face and store in flowja
819  nodes = this%dis%nodes
820  do n = 1, nodes
821  !
822  ! -- Skip if inactive.
823  if (this%ibound(n) .eq. 0) cycle
824  nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
825  !
826  ! -- Load conductivity and connection info for cell 0.
827  call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
828  ck0, allhc0)
829  !
830  ! -- Loop over active neighbors of cell 0 that have a higher
831  ! cell number (taking advantage of reciprocity).
832  do il0 = 1, nnbr0
833  m = inbr0(il0)
834  !
835  ! -- Skip if neighbor is inactive or has lower cell number.
836  if ((inbr0(il0) .eq. 0) .or. (m .lt. n)) cycle
837  nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
838  !
839  ! -- Load conductivity and connection info for cell 1.
840  call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
841  ck1, allhc1)
842  !
843  ! -- Set various indices.
844  call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
845  ii00, ii11, ii10)
846  !
847  ! -- Compute areas.
848  if (this%inewton /= 0) &
849  call this%xt3d_areas(nodes, n, m, jjs01, .true., ar01, ar10, hnew)
850  call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
851  !
852  ! -- Compute "conductances" for interface between
853  ! cells 0 and 1.
854  call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, ck0, &
855  nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
856  this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
857  !
858  ! -- Contribution to flow from primary connection.
859  qnm = chat01 * (hnew(m) - hnew(n))
860  !
861  ! -- Contribution from immediate neighbors of node 0.
862  call this%xt3d_qnbrs(nodes, n, m, nnbr0, inbr0, chati0, hnew, qnbrs)
863  qnm = qnm + qnbrs
864  !
865  ! -- Contribution from immediate neighbors of node 1.
866  call this%xt3d_qnbrs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, qnbrs)
867  qnm = qnm - qnbrs
868  ipos = ii01
869  flowja(ipos) = flowja(ipos) + qnm
870  flowja(this%dis%con%isym(ipos)) = flowja(this%dis%con%isym(ipos)) - qnm
871  end do
872  end do
Here is the call graph for this function:

◆ xt3d_flowjahfb()

subroutine xt3dmodule::xt3d_flowjahfb ( class(xt3dtype this,
integer(i4b)  n,
integer(i4b)  m,
real(dp), dimension(:), intent(inout)  hnew,
real(dp), dimension(:), intent(inout)  flowja,
real(dp)  condhfb 
)

Definition at line 877 of file Xt3dInterface.f90.

878  ! -- modules
879  use constantsmodule, only: done
880  use xt3dalgorithmmodule, only: qconds
881  ! -- dummy
882  class(Xt3dType) :: this
883  integer(I4B) :: n, m
884  real(DP), intent(inout), dimension(:) :: hnew
885  real(DP), intent(inout), dimension(:) :: flowja
886  real(DP) :: condhfb
887  ! -- local
888  integer(I4B) :: nodes
889  logical :: allhc0, allhc1
890  ! integer(I4B), parameter :: nbrmax = 10
891  integer(I4B) :: nnbr0, nnbr1
892  integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10, il
893  integer(I4B), dimension(this%nbrmax) :: inbr0, inbr1
894  integer(I4B) :: ipos
895  real(DP) :: ar01, ar10
896  real(DP), dimension(this%nbrmax, 3) :: vc0, vn0, vc1, vn1
897  real(DP), dimension(this%nbrmax) :: dl0, dl0n, dl1, dl1n
898  real(DP), dimension(3, 3) :: ck0, ck1
899  real(DP) :: chat01
900  real(DP), dimension(this%nbrmax) :: chati0, chat1j
901  real(DP) :: qnm, qnbrs
902  real(DP) :: term
903  !
904  ! -- Calculate hfb corrections to xt3d conductance-like coefficients and
905  ! put into amat and rhs as appropriate
906  nodes = this%dis%nodes
907  nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
908  !
909  ! -- Load conductivity and connection info for cell 0.
910  call this%xt3d_load(nodes, n, nnbr0, inbr0, vc0, vn0, dl0, dl0n, &
911  ck0, allhc0)
912  !
913  ! -- Find local neighbor number of cell 1.
914  do il = 1, nnbr0
915  if (inbr0(il) .eq. m) then
916  il0 = il
917  exit
918  end if
919  end do
920  !
921  nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
922  !
923  ! -- Load conductivity and connection info for cell 1.
924  call this%xt3d_load(nodes, m, nnbr1, inbr1, vc1, vn1, dl1, dl1n, &
925  ck1, allhc1)
926  !
927  ! -- Set various indices.
928  call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
929  ii00, ii11, ii10)
930  !
931  ! -- Compute areas.
932  if (this%inewton /= 0) then
933  ar01 = done
934  ar10 = done
935  else
936  call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
937  end if
938  !
939  ! -- Compute "conductances" for interface between
940  ! cells 0 and 1.
941  call qconds(this%nbrmax, nnbr0, inbr0, il01, vc0, vn0, dl0, dl0n, &
942  ck0, nnbr1, inbr1, il10, vc1, vn1, dl1, dl1n, ck1, ar01, ar10, &
943  this%vcthresh, allhc0, allhc1, chat01, chati0, chat1j)
944  !
945  ! -- Apply scale factor to compute "conductances" for hfb correction
946  if (condhfb > dzero) then
947  term = chat01 / (chat01 + condhfb)
948  else
949  term = -condhfb
950  end if
951  chat01 = -chat01 * term
952  chati0 = -chati0 * term
953  chat1j = -chat1j * term
954  !
955  ! -- Contribution to flow from primary connection.
956  qnm = chat01 * (hnew(m) - hnew(n))
957  !
958  ! -- Contribution from immediate neighbors of node 0.
959  call this%xt3d_qnbrs(nodes, n, m, nnbr0, inbr0, chati0, hnew, qnbrs)
960  qnm = qnm + qnbrs
961  !
962  ! -- Contribution from immediate neighbors of node 1.
963  call this%xt3d_qnbrs(nodes, m, n, nnbr1, inbr1, chat1j, hnew, qnbrs)
964  qnm = qnm - qnbrs
965  !
966  ! -- If Newton, scale conductance-like coefficients by the
967  ! actual area.
968  if (this%inewton /= 0) then
969  call this%xt3d_areas(nodes, n, m, jjs01, .true., ar01, ar10, hnew)
970  call this%xt3d_areas(nodes, n, m, jjs01, .false., ar01, ar10, hnew)
971  qnm = qnm * ar01
972  end if
973  !
974  ipos = ii01
975  flowja(ipos) = flowja(ipos) + qnm
976  flowja(this%dis%con%isym(ipos)) = flowja(this%dis%con%isym(ipos)) - qnm
Here is the call graph for this function:

◆ xt3d_fn()

subroutine xt3dmodule::xt3d_fn ( class(xt3dtype this,
integer(i4b)  kiter,
integer(i4b), intent(in)  nodes,
integer(i4b), intent(in)  nja,
class(matrixbasetype), pointer  matrix_sln,
integer(i4b), dimension(nja), intent(in)  idxglo,
real(dp), dimension(nodes), intent(inout)  rhs,
real(dp), dimension(nodes), intent(inout)  hnew 
)

Definition at line 693 of file Xt3dInterface.f90.

694  ! -- modules
695  use constantsmodule, only: done
697  ! -- dummy
698  class(Xt3dType) :: this
699  integer(I4B) :: kiter
700  integer(I4B), intent(in) :: nodes
701  integer(I4B), intent(in) :: nja
702  class(MatrixBaseType), pointer :: matrix_sln
703  integer(I4B), intent(in), dimension(nja) :: idxglo
704  real(DP), intent(inout), dimension(nodes) :: rhs
705  real(DP), intent(inout), dimension(nodes) :: hnew
706  ! -- local
707  integer(I4B) :: n, m, ipos
708  integer(I4B) :: nnbr0
709  integer(I4B) :: il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
710  integer(I4B), dimension(this%nbrmax) :: inbr0
711  integer(I4B) :: iups, idn
712  real(DP) :: topup, botup, derv, term
713  !
714  ! -- Update amat and rhs with Newton terms
715  do n = 1, nodes
716  !
717  ! -- Skip if inactive.
718  if (this%ibound(n) .eq. 0) cycle
719  !
720  ! -- No Newton correction if amat saved (which implies no rhs option)
721  ! and all connections for the cell are permanently confined.
722  if (this%lamatsaved) then
723  if (this%iallpc(n) == 1) cycle
724  end if
725  nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
726  !
727  ! -- Load neighbors of cell. Set cell numbers for inactive
728  ! neighbors to zero.
729  call this%xt3d_load_inbr(n, nnbr0, inbr0)
730  !
731  ! -- Loop over active neighbors of cell 0 that have a higher
732  ! cell number (taking advantage of reciprocity).
733  do il0 = 1, nnbr0
734  ipos = this%dis%con%ia(n) + il0
735  if (this%dis%con%mask(ipos) == 0) cycle
736  !
737  m = inbr0(il0)
738  !
739  ! -- Skip if neighbor is inactive or has lower cell number.
740  if ((inbr0(il0) .eq. 0) .or. (m .lt. n)) cycle
741  !
742  ! -- Set various indices.
743  call this%xt3d_indices(n, m, il0, ii01, jjs01, il01, il10, &
744  ii00, ii11, ii10)
745  !
746  ! -- Determine upstream node
747  iups = m
748  if (hnew(m) < hnew(n)) iups = n
749  idn = n
750  if (iups == n) idn = m
751  !
752  ! -- No Newton terms if upstream cell is confined and no rhs option
753  if ((this%icelltype(iups) == 0) .and. (this%ixt3d .eq. 1)) cycle
754  !
755  ! -- Set the upstream top and bot, and then recalculate for a
756  ! vertically staggered horizontal connection
757  topup = this%dis%top(iups)
758  botup = this%dis%bot(iups)
759  if (this%dis%con%ihc(jjs01) == 2) then
760  topup = min(this%dis%top(n), this%dis%top(m))
761  botup = max(this%dis%bot(n), this%dis%bot(m))
762  end if
763  !
764  ! -- Derivative term
765  derv = squadraticsaturationderivative(topup, botup, hnew(iups))
766  term = this%qsat(ii01) * derv
767  !
768  ! -- Fill Jacobian for n being the upstream node
769  if (iups == n) then
770  !
771  ! -- Fill in row of n
772  call matrix_sln%add_value_pos(idxglo(ii00), term)
773  rhs(n) = rhs(n) + term * hnew(n)
774  !
775  ! -- Fill in row of m
776  call matrix_sln%add_value_pos(idxglo(ii10), -term)
777  rhs(m) = rhs(m) - term * hnew(n)
778  !
779  ! -- Fill Jacobian for m being the upstream node
780  else
781  !
782  ! -- Fill in row of n
783  call matrix_sln%add_value_pos(idxglo(ii01), term)
784  rhs(n) = rhs(n) + term * hnew(m)
785  !
786  ! -- Fill in row of m
787  call matrix_sln%add_value_pos(idxglo(ii11), -term)
788  rhs(m) = rhs(m) - term * hnew(m)
789  !
790  end if
791  end do
792  end do
real(dp) function squadraticsaturationderivative(top, bot, x, eps)
@ brief Derivative of the quadratic saturation function
Here is the call graph for this function:

◆ xt3d_get_iinm()

subroutine xt3dmodule::xt3d_get_iinm ( class(xt3dtype this,
integer(i4b)  n,
integer(i4b)  m,
integer(i4b)  iinm 
)

Definition at line 1483 of file Xt3dInterface.f90.

1484  ! -- dummy
1485  class(Xt3dType) :: this
1486  integer(I4B) :: n, m, iinm
1487  ! -- local
1488  integer(I4B) :: ii, jj
1489  !
1490  iinm = 0
1491  do ii = this%dis%con%ia(n), this%dis%con%ia(n + 1) - 1
1492  jj = this%dis%con%ja(ii)
1493  if (jj .eq. m) then
1494  iinm = ii
1495  exit
1496  end if
1497  end do

◆ xt3d_get_iinmx()

subroutine xt3dmodule::xt3d_get_iinmx ( class(xt3dtype this,
integer(i4b)  n,
integer(i4b)  m,
integer(i4b)  iinmx 
)

Definition at line 1503 of file Xt3dInterface.f90.

1504  ! -- dummy
1505  class(Xt3dType) :: this
1506  integer(I4B) :: n, m, iinmx
1507  ! -- local
1508  integer(I4B) :: iix, jjx
1509  !
1510  iinmx = 0
1511  do iix = this%iax(n), this%iax(n + 1) - 1
1512  jjx = this%jax(iix)
1513  if (jjx .eq. m) then
1514  iinmx = iix
1515  exit
1516  end if
1517  end do

◆ xt3d_iallpc()

subroutine xt3dmodule::xt3d_iallpc ( class(xt3dtype this)

Definition at line 1109 of file Xt3dInterface.f90.

1110  ! -- modules
1112  ! -- dummy
1113  class(Xt3dType) :: this
1114  ! -- local
1115  integer(I4B) :: n, m, mm, il0, il1
1116  integer(I4B) :: nnbr0, nnbr1
1117  integer(I4B), dimension(this%nbrmax) :: inbr0, inbr1
1118  !
1119  if (this%ixt3d == 2) then
1120  this%lamatsaved = .false.
1121  call mem_allocate(this%iallpc, 0, 'IALLPC', this%memoryPath)
1122  else
1123  !
1124  ! -- allocate memory for iallpc and initialize to 1
1125  call mem_allocate(this%iallpc, this%dis%nodes, 'IALLPC', this%memoryPath)
1126  do n = 1, this%dis%nodes
1127  this%iallpc(n) = 1
1128  end do
1129  !
1130  ! -- Go through cells and connections and set iallpc to 0 if any
1131  ! connected cell is not confined
1132  do n = 1, this%dis%nodes
1133  if (this%icelltype(n) /= 0) then
1134  this%iallpc(n) = 0
1135  cycle
1136  end if
1137  nnbr0 = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
1138  call this%xt3d_load_inbr(n, nnbr0, inbr0)
1139  do il0 = 1, nnbr0
1140  m = inbr0(il0)
1141  if (m .lt. n) cycle
1142  if (this%icelltype(m) /= 0) then
1143  this%iallpc(n) = 0
1144  this%iallpc(m) = 0
1145  cycle
1146  end if
1147  nnbr1 = this%dis%con%ia(m + 1) - this%dis%con%ia(m) - 1
1148  call this%xt3d_load_inbr(m, nnbr1, inbr1)
1149  do il1 = 1, nnbr1
1150  mm = inbr1(il1)
1151  if (this%icelltype(mm) /= 0) then
1152  this%iallpc(n) = 0
1153  this%iallpc(m) = 0
1154  this%iallpc(mm) = 0
1155  end if
1156  end do
1157  end do
1158  end do
1159  !
1160  ! -- If any cells have all permanently confined connections then
1161  ! performance can be improved by precalculating coefficients
1162  ! so set lamatsaved to true.
1163  this%lamatsaved = .false.
1164  do n = 1, this%dis%nodes
1165  if (this%iallpc(n) == 1) then
1166  this%lamatsaved = .true.
1167  exit
1168  end if
1169  end do
1170  end if
1171  !
1172  if (.not. this%lamatsaved) then
1173  ! there are no permanently confined connections so deallocate iallpc
1174  ! in order to save memory
1175  call mem_reallocate(this%iallpc, 0, 'IALLPC', this%memoryPath)
1176  end if

◆ xt3d_indices()

subroutine xt3dmodule::xt3d_indices ( class(xt3dtype this,
integer(i4b)  n,
integer(i4b)  m,
integer(i4b)  il0,
integer(i4b)  ii01,
integer(i4b)  jjs01,
integer(i4b)  il01,
integer(i4b)  il10,
integer(i4b)  ii00,
integer(i4b)  ii11,
integer(i4b)  ii10 
)

Definition at line 1181 of file Xt3dInterface.f90.

1183  ! -- dummy
1184  class(Xt3dType) :: this
1185  integer(I4B) :: n, m, il0, ii01, jjs01, il01, il10, ii00, ii11, ii10
1186  ! -- local
1187  integer(I4B) :: iinm
1188  !
1189  ! -- Set local number of node 0-1 connection (local cell number of cell 1
1190  ! in cell 0's neighbor list).
1191  il01 = il0
1192  ! -- Set local number of node 1-0 connection (local cell number of cell 0
1193  ! in cell 1's neighbor list).
1194  call this%xt3d_get_iinm(m, n, iinm)
1195  il10 = iinm - this%dis%con%ia(m)
1196  ! -- Set index of node 0 diagonal in the ja array.
1197  ii00 = this%dis%con%ia(n)
1198  ! -- Set index of node 0-1 connection in the ja array.
1199  ii01 = ii00 + il01
1200  ! -- Set symmetric index of node 0-1 connection.
1201  jjs01 = this%dis%con%jas(ii01)
1202  ! -- Set index of node 1 diagonal in the ja array.
1203  ii11 = this%dis%con%ia(m)
1204  ! -- Set index of node 1-0 connection in the ja array.
1205  ii10 = ii11 + il10

◆ xt3d_load()

subroutine xt3dmodule::xt3d_load ( class(xt3dtype this,
integer(i4b), intent(in)  nodes,
integer(i4b)  n,
integer(i4b)  nnbr,
integer(i4b), dimension(this%nbrmax)  inbr,
real(dp), dimension(this%nbrmax, 3)  vc,
real(dp), dimension(this%nbrmax, 3)  vn,
real(dp), dimension(this%nbrmax)  dl,
real(dp), dimension(this%nbrmax)  dln,
real(dp), dimension(3, 3)  ck,
logical  allhc 
)

Definition at line 1211 of file Xt3dInterface.f90.

1212  ! -- module
1213  use constantsmodule, only: dzero, dhalf, done
1214  ! -- dummy
1215  class(Xt3dType) :: this
1216  logical :: allhc
1217  integer(I4B), intent(in) :: nodes
1218  integer(I4B) :: n, nnbr
1219  integer(I4B), dimension(this%nbrmax) :: inbr
1220  real(DP), dimension(this%nbrmax, 3) :: vc, vn
1221  real(DP), dimension(this%nbrmax) :: dl, dln
1222  real(DP), dimension(3, 3) :: ck
1223  ! -- local
1224  integer(I4B) :: il, ii, jj, jjs
1225  integer(I4B) :: ihcnjj
1226  real(DP) :: satn, satjj
1227  real(DP) :: cl1njj, cl2njj, dltot, ooclsum
1228  !
1229  ! -- Set conductivity tensor for cell.
1230  ck = dzero
1231  ck(1, 1) = this%k11(n)
1232  ck(2, 2) = this%k22(n)
1233  ck(3, 3) = this%k33(n)
1234  call this%xt3d_fillrmatck(n)
1235  ck = matmul(this%rmatck, ck)
1236  ck = matmul(ck, transpose(this%rmatck))
1237  !
1238  ! -- Load neighbors of cell. Set cell numbers for inactive neighbors to
1239  ! zero so xt3d knows to ignore them. Compute direct connection lengths
1240  ! from perpendicular connection lengths. Also determine if all active
1241  ! connections are horizontal.
1242  allhc = .true.
1243  do il = 1, nnbr
1244  ii = il + this%dis%con%ia(n)
1245  jj = this%dis%con%ja(ii)
1246  jjs = this%dis%con%jas(ii)
1247  if (this%ibound(jj) .ne. 0) then
1248  inbr(il) = jj
1249  satn = this%sat(n)
1250  satjj = this%sat(jj)
1251  !
1252  ! -- DISV and DIS
1253  ihcnjj = this%dis%con%ihc(jjs)
1254  call this%dis%connection_normal(n, jj, ihcnjj, vn(il, 1), vn(il, 2), &
1255  vn(il, 3), ii)
1256  call this%dis%connection_vector(n, jj, this%nozee, satn, satjj, ihcnjj, &
1257  vc(il, 1), vc(il, 2), vc(il, 3), dltot)
1258  if (jj > n) then
1259  cl1njj = this%dis%con%cl1(jjs)
1260  cl2njj = this%dis%con%cl2(jjs)
1261  else
1262  cl1njj = this%dis%con%cl2(jjs)
1263  cl2njj = this%dis%con%cl1(jjs)
1264  end if
1265  ooclsum = 1d0 / (cl1njj + cl2njj)
1266  dl(il) = dltot * cl1njj * ooclsum
1267  dln(il) = dltot * cl2njj * ooclsum
1268  if (this%dis%con%ihc(jjs) .eq. 0) allhc = .false.
1269  else
1270  inbr(il) = 0
1271  end if
1272  end do
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65

◆ xt3d_load_inbr()

subroutine xt3dmodule::xt3d_load_inbr ( class(xt3dtype this,
integer(i4b)  n,
integer(i4b)  nnbr,
integer(i4b), dimension(this%nbrmax)  inbr 
)

Definition at line 1277 of file Xt3dInterface.f90.

1278  ! -- dummy
1279  class(Xt3dType) :: this
1280  integer(I4B) :: n, nnbr
1281  integer(I4B), dimension(this%nbrmax) :: inbr
1282  ! -- local
1283  integer(I4B) :: il, ii, jj
1284  !
1285  ! -- Load neighbors of cell. Set cell numbers for inactive
1286  ! neighbors to zero so xt3d knows to ignore them.
1287  do il = 1, nnbr
1288  ii = il + this%dis%con%ia(n)
1289  jj = this%dis%con%ja(ii)
1290  if (this%ibound(jj) .ne. 0) then
1291  inbr(il) = jj
1292  else
1293  inbr(il) = 0
1294  end if
1295  end do

◆ xt3d_mc()

subroutine xt3dmodule::xt3d_mc ( class(xt3dtype this,
integer(i4b), intent(in)  moffset,
class(matrixbasetype), pointer  matrix_sln 
)

Definition at line 191 of file Xt3dInterface.f90.

192  ! -- modules
194  ! -- dummy
195  class(Xt3dType) :: this
196  integer(I4B), intent(in) :: moffset
197  class(MatrixBaseType), pointer :: matrix_sln
198  ! -- local
199  integer(I4B) :: i, j, iglo, jglo, niax, njax, ipos
200  integer(I4B) :: ipos_sln, icol_first, icol_last
201  integer(I4B) :: jj_xt3d
202  integer(I4B) :: igfirstnod, iglastnod
203  logical :: isextnbr
204  !
205  ! -- If not rhs, map connections for extended neighbors and construct iax,
206  ! -- jax, and idxglox
207  if (this%ixt3d == 1) then
208  !
209  ! -- calculate the first node for the model and the last node in global
210  ! numbers
211  igfirstnod = moffset + 1
212  iglastnod = moffset + this%dis%nodes
213  !
214  ! -- allocate iax, jax, and idxglox
215  niax = this%dis%nodes + 1
216  njax = this%numextnbrs ! + 1
217  call mem_allocate(this%iax, niax, 'IAX', trim(this%memoryPath))
218  call mem_allocate(this%jax, njax, 'JAX', trim(this%memoryPath))
219  call mem_allocate(this%idxglox, njax, 'IDXGLOX', trim(this%memoryPath))
220  !
221  ! -- load first iax entry
222  ipos = 1
223  this%iax(1) = ipos
224  !
225  ! -- loop over nodes
226  do i = 1, this%dis%nodes
227  !
228  ! -- calculate global node number
229  iglo = i + moffset
230  icol_first = matrix_sln%get_first_col_pos(iglo)
231  icol_last = matrix_sln%get_last_col_pos(iglo)
232  do ipos_sln = icol_first, icol_last
233  !
234  ! -- if jglo is in a different model, then it cannot be an extended
235  ! neighbor, so skip over it
236  jglo = matrix_sln%get_column(ipos_sln)
237  if (jglo < igfirstnod .or. jglo > iglastnod) then
238  cycle
239  end if
240  !
241  ! -- Check to see if this local connection was added by
242  ! xt3d. If not, then this connection was added by
243  ! something else, such as an interface model.
244  j = jglo - moffset
245  isextnbr = .false.
246  do jj_xt3d = this%ia_xt3d(i), this%ia_xt3d(i + 1) - 1
247  if (j == this%ja_xt3d(jj_xt3d)) then
248  isextnbr = .true.
249  exit
250  end if
251  end do
252  !
253  ! -- if an extended neighbor, add it to jax and idxglox
254  if (isextnbr) then
255  this%jax(ipos) = matrix_sln%get_column(ipos_sln) - moffset
256  this%idxglox(ipos) = ipos_sln
257  ipos = ipos + 1
258  end if
259  end do
260  ! -- load next iax entry
261  this%iax(i + 1) = ipos
262  end do
263  !
264  else
265  !
266  call mem_allocate(this%iax, 0, 'IAX', trim(this%memoryPath))
267  call mem_allocate(this%jax, 0, 'JAX', trim(this%memoryPath))
268  call mem_allocate(this%idxglox, 0, 'IDXGLOX', trim(this%memoryPath))
269  !
270  end if

◆ xt3d_qnbrs()

subroutine xt3dmodule::xt3d_qnbrs ( class(xt3dtype this,
integer(i4b), intent(in)  nodes,
integer(i4b)  n,
integer(i4b)  m,
integer(i4b)  nnbr,
integer(i4b), dimension(this%nbrmax)  inbr,
real(dp), dimension(this%nbrmax)  chat,
real(dp), dimension(nodes), intent(inout)  hnew,
real(dp)  qnbrs 
)

Definition at line 1548 of file Xt3dInterface.f90.

1550  ! -- dummy
1551  class(Xt3dType) :: this
1552  integer(I4B), intent(in) :: nodes
1553  integer(I4B) :: n, m, nnbr
1554  integer(I4B), dimension(this%nbrmax) :: inbr
1555  real(DP) :: qnbrs
1556  real(DP), dimension(this%nbrmax) :: chat
1557  real(DP), intent(inout), dimension(nodes) :: hnew
1558  ! -- local
1559  integer(I4B) :: iil, iii, jjj
1560  real(DP) :: term
1561  !
1562  qnbrs = 0d0
1563  do iil = 1, nnbr
1564  if (inbr(iil) .ne. 0) then
1565  iii = iil + this%dis%con%ia(n)
1566  jjj = this%dis%con%ja(iii)
1567  term = chat(iil) * (hnew(jjj) - hnew(n))
1568  qnbrs = qnbrs + term
1569  end if
1570  end do

◆ xt3d_rhs()

subroutine xt3dmodule::xt3d_rhs ( class(xt3dtype this,
integer(i4b), intent(in)  nodes,
integer(i4b)  n,
integer(i4b)  m,
integer(i4b)  nnbr,
integer(i4b), dimension(this%nbrmax)  inbr,
real(dp), dimension(this%nbrmax)  chat,
real(dp), dimension(nodes), intent(inout)  hnew,
real(dp), dimension(nodes), intent(inout)  rhs 
)

Definition at line 1522 of file Xt3dInterface.f90.

1524  ! -- dummy
1525  class(Xt3dType) :: this
1526  integer(I4B), intent(in) :: nodes
1527  integer(I4B) :: n, m, nnbr
1528  integer(I4B), dimension(this%nbrmax) :: inbr
1529  real(DP), dimension(this%nbrmax) :: chat
1530  real(DP), intent(inout), dimension(nodes) :: hnew, rhs
1531  ! -- local
1532  integer(I4B) :: iil, iii, jjj
1533  real(DP) :: term
1534  !
1535  do iil = 1, nnbr
1536  if (inbr(iil) .ne. 0) then
1537  iii = iil + this%dis%con%ia(n)
1538  jjj = this%dis%con%ja(iii)
1539  term = chat(iil) * (hnew(jjj) - hnew(n))
1540  rhs(n) = rhs(n) - term
1541  rhs(m) = rhs(m) + term
1542  end if
1543  end do