MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
gwf-npf.f90
Go to the documentation of this file.
2  use kindmodule, only: dp, i4b
3  use simvariablesmodule, only: errmsg
4  use constantsmodule, only: dzero, dem9, dem8, dem7, dem6, dem2, &
5  dhalf, dp9, done, dtwo, &
6  dhnoflo, dhdry, dem10, &
13  use basedismodule, only: disbasetype
14  use gwficmodule, only: gwfictype
15  use gwfvscmodule, only: gwfvsctype
16  use xt3dmodule, only: xt3dtype
19  use tvkmodule, only: tvktype, tvk_cr
24  use hgeoutilmodule, only: hyeff
26  condmean, thksatnm, &
28 
29  implicit none
30 
31  private
32  public :: gwfnpftype
33  public :: npf_cr
34 
35  type, extends(numericalpackagetype) :: gwfnpftype
36 
37  type(gwfictype), pointer :: ic => null() !< initial conditions object
38  type(gwfvsctype), pointer :: vsc => null() !< viscosity object
39  type(xt3dtype), pointer :: xt3d => null() !< xt3d pointer
40  integer(I4B), pointer :: iname => null() !< length of variable names
41  character(len=24), dimension(:), pointer :: aname => null() !< variable names
42  integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound
43  real(dp), dimension(:), pointer, contiguous :: hnew => null() !< pointer to model xnew
44  integer(I4B), pointer :: ixt3d => null() !< xt3d flag (0 is off, 1 is lhs, 2 is rhs)
45  integer(I4B), pointer :: ixt3drhs => null() !< xt3d rhs flag, xt3d rhs is set active if 1
46  integer(I4B), pointer :: iperched => null() !< vertical flow corrections if 1
47  integer(I4B), pointer :: ivarcv => null() !< CV is function of water table
48  integer(I4B), pointer :: idewatcv => null() !< CV may be a discontinuous function of water table
49  integer(I4B), pointer :: ithickstrt => null() !< thickstrt option flag
50  integer(I4B), pointer :: igwfnewtonur => null() !< newton head dampening using node bottom option flag
51  integer(I4B), pointer :: icalcspdis => null() !< Calculate specific discharge at cell centers
52  integer(I4B), pointer :: isavspdis => null() !< Save specific discharge at cell centers
53  integer(I4B), pointer :: isavsat => null() !< Save sat to budget file
54  real(dp), pointer :: hnoflo => null() !< default is 1.e30
55  real(dp), pointer :: satomega => null() !< newton-raphson saturation omega
56  integer(I4B), pointer :: irewet => null() !< rewetting (0:off, 1:on)
57  integer(I4B), pointer :: iwetit => null() !< wetting interval (default is 1)
58  integer(I4B), pointer :: ihdwet => null() !< (0 or not 0)
59  integer(I4B), pointer :: icellavg => null() !< harmonic(0), logarithmic(1), or arithmetic thick-log K (2)
60  real(dp), pointer :: wetfct => null() !< wetting factor
61  real(dp), pointer :: hdry => null() !< default is -1.d30
62  integer(I4B), dimension(:), pointer, contiguous :: icelltype => null() !< confined (0) or convertible (1)
63  integer(I4B), dimension(:), pointer, contiguous :: ithickstartflag => null() !< array of flags for handling the thickstrt option
64  !
65  ! K properties
66  real(dp), dimension(:), pointer, contiguous :: k11 => null() !< hydraulic conductivity; if anisotropic, then this is Kx prior to rotation
67  real(dp), dimension(:), pointer, contiguous :: k22 => null() !< hydraulic conductivity; if specified then this is Ky prior to rotation
68  real(dp), dimension(:), pointer, contiguous :: k33 => null() !< hydraulic conductivity; if specified then this is Kz prior to rotation
69  real(dp), dimension(:), pointer, contiguous :: k11input => null() !< hydraulic conductivity originally specified by user prior to TVK or VSC modification
70  real(dp), dimension(:), pointer, contiguous :: k22input => null() !< hydraulic conductivity originally specified by user prior to TVK or VSC modification
71  real(dp), dimension(:), pointer, contiguous :: k33input => null() !< hydraulic conductivity originally specified by user prior to TVK or VSC modification
72  integer(I4B), pointer :: iavgkeff => null() !< effective conductivity averaging (0: harmonic, 1: arithmetic)
73  integer(I4B), pointer :: ik22 => null() !< flag that k22 is specified
74  integer(I4B), pointer :: ik33 => null() !< flag that k33 is specified
75  integer(I4B), pointer :: ik22overk => null() !< flag that k22 is specified as anisotropy ratio
76  integer(I4B), pointer :: ik33overk => null() !< flag that k33 is specified as anisotropy ratio
77  integer(I4B), pointer :: iangle1 => null() !< flag to indicate angle1 was read
78  integer(I4B), pointer :: iangle2 => null() !< flag to indicate angle2 was read
79  integer(I4B), pointer :: iangle3 => null() !< flag to indicate angle3 was read
80  real(dp), dimension(:), pointer, contiguous :: angle1 => null() !< k ellipse rotation in xy plane around z axis (yaw)
81  real(dp), dimension(:), pointer, contiguous :: angle2 => null() !< k ellipse rotation up from xy plane around y axis (pitch)
82  real(dp), dimension(:), pointer, contiguous :: angle3 => null() !< k tensor rotation around x axis (roll)
83  !
84  integer(I4B), pointer :: iwetdry => null() !< flag to indicate angle1 was read
85  real(dp), dimension(:), pointer, contiguous :: wetdry => null() !< wetdry array
86  real(dp), dimension(:), pointer, contiguous :: sat => null() !< saturation (0. to 1.) for each cell
87  real(dp), dimension(:), pointer, contiguous :: condsat => null() !< saturated conductance (symmetric array)
88  integer(I4B), dimension(:), pointer, contiguous :: ibotnode => null() !< bottom node used if igwfnewtonur /= 0
89  ! spdis machinery:
90  real(dp), dimension(:, :), pointer, contiguous :: spdis => null() !< specific discharge : qx, qy, qz (nodes, 3)
91  integer(I4B), pointer :: nedges => null() !< number of cell edges
92  integer(I4B), pointer :: lastedge => null() !< last edge number
93  integer(I4B), dimension(:), pointer, contiguous :: nodedge => null() !< array of node numbers that have edges
94  integer(I4B), dimension(:), pointer, contiguous :: ihcedge => null() !< edge type (horizontal or vertical)
95  real(dp), dimension(:, :), pointer, contiguous :: propsedge => null() !< edge properties (Q, area, nx, ny, distance)
96  integer(I4B), dimension(:), pointer, contiguous :: iedge_ptr => null() !< csr pointer into edge index array
97  integer(I4B), dimension(:), pointer, contiguous :: edge_idxs => null() !< sorted edge indexes for faster lookup
98  type(spdisworkarraytype), pointer :: spdis_wa => null() !< work arrays for spdis calculation
99  !
100  integer(I4B), pointer :: intvk => null() ! TVK (time-varying K) unit number (0 if unused)
101  integer(I4B), pointer :: invsc => null() ! VSC (viscosity) unit number (0 if unused); viscosity leads to time-varying K's
102  type(tvktype), pointer :: tvk => null() ! TVK object
103  integer(I4B), pointer :: kchangeper => null() ! last stress period in which any node K (or K22, or K33) values were changed (0 if unchanged from start of simulation)
104  integer(I4B), pointer :: kchangestp => null() ! last time step in which any node K (or K22, or K33) values were changed (0 if unchanged from start of simulation)
105  integer(I4B), dimension(:), pointer, contiguous :: nodekchange => null() ! grid array of flags indicating for each node whether its K (or K22, or K33) value changed (1) at (kchangeper, kchangestp) or not (0)
106 
107  contains
108 
109  procedure :: npf_df
110  procedure :: npf_ac
111  procedure :: npf_mc
112  procedure :: npf_ar
113  procedure :: npf_rp
114  procedure :: npf_ad
115  procedure :: npf_cf
116  procedure :: npf_fc
117  procedure :: npf_fn
118  procedure :: npf_cq
119  procedure :: npf_save_model_flows
120  procedure :: npf_nur
122  procedure :: npf_da
123  procedure, private :: thksat => sgwf_npf_thksat
124  procedure, private :: qcalc => sgwf_npf_qcalc
125  procedure, private :: wd => sgwf_npf_wetdry
126  procedure, private :: wdmsg => sgwf_npf_wdmsg
127  procedure :: allocate_scalars
128  procedure, private :: store_original_k_arrays
129  procedure, private :: allocate_arrays
130  procedure, private :: source_options
131  procedure, private :: source_griddata
132  procedure, private :: log_options
133  procedure, private :: log_griddata
134  procedure, private :: set_options
135  procedure, private :: check_options
136  procedure, private :: prepcheck
137  procedure, private :: preprocess_input
138  procedure, private :: calc_condsat
139  procedure, private :: calc_initial_sat
140  procedure, public :: rewet_check
141  procedure, public :: hy_eff
142  procedure, public :: calc_spdis
143  procedure, public :: sav_spdis
144  procedure, public :: sav_sat
145  procedure, public :: increase_edge_count
146  procedure, public :: set_edge_properties
147  procedure, public :: calcsatthickness
148  procedure, private :: calc_max_conns
149  procedure, private :: prepare_edge_lookup
150  end type
151 
152 contains
153 
154  !> @brief Create a new NPF object. Pass a inunit value of 0 if npf data will
155  !! initialized from memory
156  !<
157  subroutine npf_cr(npfobj, name_model, input_mempath, inunit, iout)
158  ! -- modules
159  use kindmodule, only: lgp
161  ! -- dummy
162  type(gwfnpftype), pointer :: npfobj
163  character(len=*), intent(in) :: name_model
164  character(len=*), intent(in) :: input_mempath
165  integer(I4B), intent(in) :: inunit
166  integer(I4B), intent(in) :: iout
167  ! -- formats
168  character(len=*), parameter :: fmtheader = &
169  "(1x, /1x, 'NPF -- NODE PROPERTY FLOW PACKAGE, VERSION 1, 3/30/2015', &
170  &' INPUT READ FROM MEMPATH: ', A, /)"
171  !
172  ! -- Create the object
173  allocate (npfobj)
174  !
175  ! -- create name and memory path
176  call npfobj%set_names(1, name_model, 'NPF', 'NPF', input_mempath)
177  !
178  ! -- Allocate scalars
179  call npfobj%allocate_scalars()
180  !
181  ! -- Set variables
182  npfobj%inunit = inunit
183  npfobj%iout = iout
184  !
185  ! -- check if npf is enabled
186  if (inunit > 0) then
187  !
188  ! -- Print a message identifying the node property flow package.
189  write (iout, fmtheader) input_mempath
190  end if
191 
192  ! allocate spdis structure
193  allocate (npfobj%spdis_wa)
194 
195  end subroutine npf_cr
196 
197  !> @brief Define the NPF package instance
198  !!
199  !! This is a hybrid routine: it either reads the options for this package
200  !! from the input file, or the optional argument @param npf_options
201  !! should be passed. A consistency check is performed, and finally
202  !! xt3d_df is called, when enabled.
203  !<
204  subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options)
205  ! -- modules
206  use simmodule, only: store_error
207  use xt3dmodule, only: xt3d_cr
208  ! -- dummy
209  class(gwfnpftype) :: this !< instance of the NPF package
210  class(disbasetype), pointer, intent(inout) :: dis !< the pointer to the discretization
211  type(xt3dtype), pointer :: xt3d !< the pointer to the XT3D 'package'
212  integer(I4B), intent(in) :: ingnc !< ghostnodes enabled? (>0 means yes)
213  integer(I4B), intent(in) :: invsc !< viscosity enabled? (>0 means yes)
214  type(gwfnpfoptionstype), optional, intent(in) :: npf_options !< the optional options, for when not constructing from file
215  !
216  ! -- Set a pointer to dis
217  this%dis => dis
218  !
219  ! -- Set flag signifying whether vsc is active
220  if (invsc > 0) this%invsc = invsc
221  !
222  if (.not. present(npf_options)) then
223  !
224  ! -- source options
225  call this%source_options()
226  !
227  ! -- allocate arrays
228  call this%allocate_arrays(this%dis%nodes, this%dis%njas)
229  !
230  ! -- source griddata, set, and convert/check the input
231  call this%source_griddata()
232  call this%prepcheck()
233  else
234  call this%set_options(npf_options)
235  !
236  ! -- allocate arrays
237  call this%allocate_arrays(this%dis%nodes, this%dis%njas)
238  end if
239  !
240  call this%check_options()
241  !
242  ! -- Save pointer to xt3d object
243  this%xt3d => xt3d
244  if (this%ixt3d /= 0) xt3d%ixt3d = this%ixt3d
245  call this%xt3d%xt3d_df(dis)
246  !
247  ! -- Ensure GNC and XT3D are not both on at the same time
248  if (this%ixt3d /= 0 .and. ingnc > 0) then
249  call store_error('Error in model '//trim(this%name_model)// &
250  '. The XT3D option cannot be used with the GNC &
251  &Package.', terminate=.true.)
252  end if
253  end subroutine npf_df
254 
255  !> @brief Add connections for extended neighbors to the sparse matrix
256  !<
257  subroutine npf_ac(this, moffset, sparse)
258  ! -- modules
259  use sparsemodule, only: sparsematrix
260  ! -- dummy
261  class(gwfnpftype) :: this
262  integer(I4B), intent(in) :: moffset
263  type(sparsematrix), intent(inout) :: sparse
264  !
265  ! -- Add extended neighbors (neighbors of neighbors)
266  if (this%ixt3d /= 0) call this%xt3d%xt3d_ac(moffset, sparse)
267  end subroutine npf_ac
268 
269  !> @brief Map connections and construct iax, jax, and idxglox
270  !<
271  subroutine npf_mc(this, moffset, matrix_sln)
272  ! -- dummy
273  class(gwfnpftype) :: this
274  integer(I4B), intent(in) :: moffset
275  class(matrixbasetype), pointer :: matrix_sln
276  !
277  if (this%ixt3d /= 0) call this%xt3d%xt3d_mc(moffset, matrix_sln)
278  end subroutine npf_mc
279 
280  !> @brief Allocate and read this NPF instance
281  !!
282  !! Allocate remaining package arrays, preprocess the input data and
283  !! call *_ar on xt3d, when active.
284  !<
285  subroutine npf_ar(this, ic, vsc, ibound, hnew)
286  ! -- modules
288  ! -- dummy
289  class(gwfnpftype) :: this !< instance of the NPF package
290  type(gwfictype), pointer, intent(in) :: ic !< initial conditions
291  type(gwfvsctype), pointer, intent(in) :: vsc !< viscosity package
292  integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: ibound !< model ibound array
293  real(DP), dimension(:), pointer, contiguous, intent(inout) :: hnew !< pointer to model head array
294  ! -- local
295  integer(I4B) :: n
296  !
297  ! -- Store pointers to arguments that were passed in
298  this%ic => ic
299  this%ibound => ibound
300  this%hnew => hnew
301  !
302  if (this%icalcspdis == 1) then
303  call mem_reallocate(this%spdis, 3, this%dis%nodes, 'SPDIS', this%memoryPath)
304  call mem_reallocate(this%nodedge, this%nedges, 'NODEDGE', this%memoryPath)
305  call mem_reallocate(this%ihcedge, this%nedges, 'IHCEDGE', this%memoryPath)
306  call mem_reallocate(this%propsedge, 5, this%nedges, 'PROPSEDGE', &
307  this%memoryPath)
308  call mem_reallocate(this%iedge_ptr, this%dis%nodes + 1, &
309  'NREDGESNODE', this%memoryPath)
310  call mem_reallocate(this%edge_idxs, this%nedges, &
311  'EDGEIDXS', this%memoryPath)
312 
313  do n = 1, this%nedges
314  this%edge_idxs(n) = 0
315  end do
316  do n = 1, this%dis%nodes
317  this%iedge_ptr(n) = 0
318  this%spdis(:, n) = dzero
319  end do
320  end if
321  !
322  ! -- Store pointer to VSC if active
323  if (this%invsc /= 0) then
324  this%vsc => vsc
325  end if
326  !
327  ! -- allocate arrays to store original user input in case TVK/VSC modify them
328  if (this%invsc > 0) then
329  !
330  ! -- Reallocate arrays that store user-input values.
331  call mem_reallocate(this%k11input, this%dis%nodes, 'K11INPUT', &
332  this%memoryPath)
333  call mem_reallocate(this%k22input, this%dis%nodes, 'K22INPUT', &
334  this%memoryPath)
335  call mem_reallocate(this%k33input, this%dis%nodes, 'K33INPUT', &
336  this%memoryPath)
337  ! Allocate arrays that will store the original K values. When VSC active,
338  ! the current Kxx arrays carry the viscosity-adjusted K values.
339  ! This approach leverages existing functionality that makes use of K.
340  call this%store_original_k_arrays(this%dis%nodes, this%dis%njas)
341  end if
342  !
343  ! -- preprocess data
344  call this%preprocess_input()
345  !
346  ! -- xt3d
347  if (this%ixt3d /= 0) then
348  call this%xt3d%xt3d_ar(ibound, this%k11, this%ik33, this%k33, &
349  this%sat, this%ik22, this%k22, &
350  this%iangle1, this%iangle2, this%iangle3, &
351  this%angle1, this%angle2, this%angle3, &
352  this%inewton, this%icelltype)
353  end if
354  !
355  ! -- TVK
356  if (this%intvk /= 0) then
357  call this%tvk%ar(this%dis)
358  end if
359  end subroutine npf_ar
360 
361  !> @brief Read and prepare method for package
362  !!
363  !! Read and prepare NPF stress period data.
364  !<
365  subroutine npf_rp(this)
366  implicit none
367  ! -- dummy
368  class(gwfnpftype) :: this
369  !
370  ! -- TVK
371  if (this%intvk /= 0) then
372  call this%tvk%rp()
373  end if
374  end subroutine npf_rp
375 
376  !> @brief Advance
377  !!
378  !! Sets hold (head old) to bot whenever a wettable cell is dry
379  !<
380  subroutine npf_ad(this, nodes, hold, hnew, irestore)
381  ! -- modules
382  use tdismodule, only: kper, kstp
383  !
384  implicit none
385  ! -- dummy
386  class(gwfnpftype) :: this
387  integer(I4B), intent(in) :: nodes
388  real(DP), dimension(nodes), intent(inout) :: hold
389  real(DP), dimension(nodes), intent(inout) :: hnew
390  integer(I4B), intent(in) :: irestore
391  ! -- local
392  integer(I4B) :: n
393  !
394  ! -- loop through all cells and set hold=bot if wettable cell is dry
395  if (this%irewet > 0) then
396  do n = 1, this%dis%nodes
397  if (this%wetdry(n) == dzero) cycle
398  if (this%ibound(n) /= 0) cycle
399  hold(n) = this%dis%bot(n)
400  end do
401  !
402  ! -- if restore state, then set hnew to DRY if it is a dry wettable cell
403  do n = 1, this%dis%nodes
404  if (this%wetdry(n) == dzero) cycle
405  if (this%ibound(n) /= 0) cycle
406  hnew(n) = dhdry
407  end do
408  end if
409  !
410  ! -- TVK
411  if (this%intvk /= 0) then
412  call this%tvk%ad()
413  end if
414  !
415  ! -- VSC
416  ! -- Hit the TVK-updated K's with VSC correction before calling/updating condsat
417  if (this%invsc /= 0) then
418  call this%vsc%update_k_with_vsc()
419  end if
420  !
421  ! -- If any K values have changed, we need to update CONDSAT or XT3D arrays
422  if (this%kchangeper == kper .and. this%kchangestp == kstp) then
423  if (this%ixt3d == 0) then
424  !
425  ! -- Update the saturated conductance for all connections
426  ! -- of the affected nodes
427  do n = 1, this%dis%nodes
428  if (this%nodekchange(n) == 1) then
429  call this%calc_condsat(n, .false.)
430  end if
431  end do
432  else
433  !
434  ! -- Recompute XT3D coefficients for permanently confined connections
435  if (this%xt3d%lamatsaved .and. .not. this%xt3d%ldispersion) then
436  call this%xt3d%xt3d_fcpc(this%dis%nodes, .true.)
437  end if
438  end if
439  end if
440  end subroutine npf_ad
441 
442  !> @brief Routines associated fill coefficients
443  !<
444  subroutine npf_cf(this, kiter, nodes, hnew)
445  ! -- dummy
446  class(gwfnpftype) :: this
447  integer(I4B) :: kiter
448  integer(I4B), intent(in) :: nodes
449  real(DP), intent(inout), dimension(nodes) :: hnew
450  ! -- local
451  integer(I4B) :: n
452  real(DP) :: satn
453  !
454  ! -- Perform wetting and drying
455  if (this%inewton /= 1) then
456  call this%wd(kiter, hnew)
457  end if
458  !
459  ! -- Calculate saturation for convertible cells
460  do n = 1, this%dis%nodes
461  if (this%icelltype(n) /= 0) then
462  if (this%ibound(n) == 0) then
463  satn = dzero
464  else
465  call this%thksat(n, hnew(n), satn)
466  end if
467  this%sat(n) = satn
468  end if
469  end do
470  end subroutine npf_cf
471 
472  !> @brief Formulate coefficients
473  !<
474  subroutine npf_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
475  ! -- modules
476  use constantsmodule, only: done
477  ! -- dummy
478  class(gwfnpftype) :: this
479  integer(I4B) :: kiter
480  class(matrixbasetype), pointer :: matrix_sln
481  integer(I4B), intent(in), dimension(:) :: idxglo
482  real(DP), intent(inout), dimension(:) :: rhs
483  real(DP), intent(inout), dimension(:) :: hnew
484  ! -- local
485  integer(I4B) :: n, m, ii, idiag, ihc
486  integer(I4B) :: isymcon, idiagm
487  real(DP) :: hyn, hym
488  real(DP) :: cond
489  !
490  ! -- Calculate conductance and put into amat
491  !
492  if (this%ixt3d /= 0) then
493  call this%xt3d%xt3d_fc(kiter, matrix_sln, idxglo, rhs, hnew)
494  else
495  do n = 1, this%dis%nodes
496  do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
497  if (this%dis%con%mask(ii) == 0) cycle
498 
499  m = this%dis%con%ja(ii)
500  !
501  ! -- Calculate conductance only for upper triangle but insert into
502  ! upper and lower parts of amat.
503  if (m < n) cycle
504  ihc = this%dis%con%ihc(this%dis%con%jas(ii))
505  hyn = this%hy_eff(n, m, ihc, ipos=ii)
506  hym = this%hy_eff(m, n, ihc, ipos=ii)
507  !
508  ! -- Vertical connection
509  if (ihc == c3d_vertical) then
510  !
511  ! -- Calculate vertical conductance
512  cond = vcond(this%ibound(n), this%ibound(m), &
513  this%icelltype(n), this%icelltype(m), this%inewton, &
514  this%ivarcv, this%idewatcv, &
515  this%condsat(this%dis%con%jas(ii)), hnew(n), hnew(m), &
516  hyn, hym, &
517  this%sat(n), this%sat(m), &
518  this%dis%top(n), this%dis%top(m), &
519  this%dis%bot(n), this%dis%bot(m), &
520  this%dis%con%hwva(this%dis%con%jas(ii)))
521  !
522  ! -- Vertical flow for perched conditions
523  if (this%iperched /= 0) then
524  if (this%icelltype(m) /= 0) then
525  if (hnew(m) < this%dis%top(m)) then
526  !
527  ! -- Fill row n
528  idiag = this%dis%con%ia(n)
529  rhs(n) = rhs(n) - cond * this%dis%bot(n)
530  call matrix_sln%add_value_pos(idxglo(idiag), -cond)
531  !
532  ! -- Fill row m
533  isymcon = this%dis%con%isym(ii)
534  call matrix_sln%add_value_pos(idxglo(isymcon), cond)
535  rhs(m) = rhs(m) + cond * this%dis%bot(n)
536  !
537  ! -- cycle the connection loop
538  cycle
539  end if
540  end if
541  end if
542  !
543  else
544  !
545  ! -- Horizontal conductance
546  cond = hcond(this%ibound(n), this%ibound(m), &
547  this%icelltype(n), this%icelltype(m), &
548  this%inewton, &
549  this%dis%con%ihc(this%dis%con%jas(ii)), &
550  this%icellavg, &
551  this%condsat(this%dis%con%jas(ii)), &
552  hnew(n), hnew(m), this%sat(n), this%sat(m), hyn, hym, &
553  this%dis%top(n), this%dis%top(m), &
554  this%dis%bot(n), this%dis%bot(m), &
555  this%dis%con%cl1(this%dis%con%jas(ii)), &
556  this%dis%con%cl2(this%dis%con%jas(ii)), &
557  this%dis%con%hwva(this%dis%con%jas(ii)))
558  end if
559  !
560  ! -- Fill row n
561  idiag = this%dis%con%ia(n)
562  call matrix_sln%add_value_pos(idxglo(ii), cond)
563  call matrix_sln%add_value_pos(idxglo(idiag), -cond)
564  !
565  ! -- Fill row m
566  isymcon = this%dis%con%isym(ii)
567  idiagm = this%dis%con%ia(m)
568  call matrix_sln%add_value_pos(idxglo(isymcon), cond)
569  call matrix_sln%add_value_pos(idxglo(idiagm), -cond)
570  end do
571  end do
572  !
573  end if
574  end subroutine npf_fc
575 
576  !> @brief Fill newton terms
577  !<
578  subroutine npf_fn(this, kiter, matrix_sln, idxglo, rhs, hnew)
579  ! -- dummy
580  class(gwfnpftype) :: this
581  integer(I4B) :: kiter
582  class(matrixbasetype), pointer :: matrix_sln
583  integer(I4B), intent(in), dimension(:) :: idxglo
584  real(DP), intent(inout), dimension(:) :: rhs
585  real(DP), intent(inout), dimension(:) :: hnew
586  ! -- local
587  integer(I4B) :: nodes, nja
588  integer(I4B) :: n, m, ii, idiag
589  integer(I4B) :: isymcon, idiagm
590  integer(I4B) :: iups
591  integer(I4B) :: idn
592  real(DP) :: cond
593  real(DP) :: consterm
594  real(DP) :: filledterm
595  real(DP) :: derv
596  real(DP) :: hds
597  real(DP) :: term
598  real(DP) :: topup
599  real(DP) :: botup
600  !
601  ! -- add newton terms to solution matrix
602  nodes = this%dis%nodes
603  nja = this%dis%con%nja
604  if (this%ixt3d /= 0) then
605  call this%xt3d%xt3d_fn(kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew)
606  else
607  !
608  do n = 1, nodes
609  idiag = this%dis%con%ia(n)
610  do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
611  if (this%dis%con%mask(ii) == 0) cycle
612 
613  m = this%dis%con%ja(ii)
614  isymcon = this%dis%con%isym(ii)
615  ! work on upper triangle
616  if (m < n) cycle
617  if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0 .and. &
618  this%ivarcv == 0) then
619  !call this%vcond(n,m,hnew(n),hnew(m),ii,cond)
620  ! do nothing
621  else
622  ! determine upstream node
623  iups = m
624  if (hnew(m) < hnew(n)) iups = n
625  idn = n
626  if (iups == n) idn = m
627  !
628  ! -- no newton terms if upstream cell is confined
629  if (this%icelltype(iups) == 0) cycle
630  !
631  ! -- Set the upstream top and bot, and then recalculate for a
632  ! vertically staggered horizontal connection
633  topup = this%dis%top(iups)
634  botup = this%dis%bot(iups)
635  if (this%dis%con%ihc(this%dis%con%jas(ii)) == 2) then
636  topup = min(this%dis%top(n), this%dis%top(m))
637  botup = max(this%dis%bot(n), this%dis%bot(m))
638  end if
639  !
640  ! get saturated conductivity for derivative
641  cond = this%condsat(this%dis%con%jas(ii))
642  !
643  ! compute additional term
644  consterm = -cond * (hnew(iups) - hnew(idn)) !needs to use hwadi instead of hnew(idn)
645  !filledterm = cond
646  filledterm = matrix_sln%get_value_pos(idxglo(ii))
647  derv = squadraticsaturationderivative(topup, botup, hnew(iups), &
648  this%satomega)
649  idiagm = this%dis%con%ia(m)
650  ! fill jacobian for n being the upstream node
651  if (iups == n) then
652  hds = hnew(m)
653  !isymcon = this%dis%con%isym(ii)
654  term = consterm * derv
655  rhs(n) = rhs(n) + term * hnew(n) !+ amat(idxglo(isymcon)) * (dwadi * hds - hds) !need to add dwadi
656  rhs(m) = rhs(m) - term * hnew(n) !- amat(idxglo(isymcon)) * (dwadi * hds - hds) !need to add dwadi
657  ! fill in row of n
658  call matrix_sln%add_value_pos(idxglo(idiag), term)
659  ! fill newton term in off diagonal if active cell
660  if (this%ibound(n) > 0) then
661  filledterm = matrix_sln%get_value_pos(idxglo(ii))
662  call matrix_sln%set_value_pos(idxglo(ii), filledterm) !* dwadi !need to add dwadi
663  end if
664  !fill row of m
665  filledterm = matrix_sln%get_value_pos(idxglo(idiagm))
666  call matrix_sln%set_value_pos(idxglo(idiagm), filledterm) !- filledterm * (dwadi - DONE) !need to add dwadi
667  ! fill newton term in off diagonal if active cell
668  if (this%ibound(m) > 0) then
669  call matrix_sln%add_value_pos(idxglo(isymcon), -term)
670  end if
671  ! fill jacobian for m being the upstream node
672  else
673  hds = hnew(n)
674  term = -consterm * derv
675  rhs(n) = rhs(n) + term * hnew(m) !+ amat(idxglo(ii)) * (dwadi * hds - hds) !need to add dwadi
676  rhs(m) = rhs(m) - term * hnew(m) !- amat(idxglo(ii)) * (dwadi * hds - hds) !need to add dwadi
677  ! fill in row of n
678  filledterm = matrix_sln%get_value_pos(idxglo(idiag))
679  call matrix_sln%set_value_pos(idxglo(idiag), filledterm) !- filledterm * (dwadi - DONE) !need to add dwadi
680  ! fill newton term in off diagonal if active cell
681  if (this%ibound(n) > 0) then
682  call matrix_sln%add_value_pos(idxglo(ii), term)
683  end if
684  !fill row of m
685  call matrix_sln%add_value_pos(idxglo(idiagm), -term)
686  ! fill newton term in off diagonal if active cell
687  if (this%ibound(m) > 0) then
688  filledterm = matrix_sln%get_value_pos(idxglo(isymcon))
689  call matrix_sln%set_value_pos(idxglo(isymcon), filledterm) !* dwadi !need to add dwadi
690  end if
691  end if
692  end if
693 
694  end do
695  end do
696  !
697  end if
698  end subroutine npf_fn
699 
700  !> @brief Under-relaxation
701  !!
702  !! Under-relaxation of Groundwater Flow Model Heads for current outer
703  !! iteration using the cell bottoms at the bottom of the model
704  !<
705  subroutine npf_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax)
706  ! -- dummy
707  class(gwfnpftype) :: this
708  integer(I4B), intent(in) :: neqmod
709  real(DP), dimension(neqmod), intent(inout) :: x
710  real(DP), dimension(neqmod), intent(in) :: xtemp
711  real(DP), dimension(neqmod), intent(inout) :: dx
712  integer(I4B), intent(inout) :: inewtonur
713  real(DP), intent(inout) :: dxmax
714  integer(I4B), intent(inout) :: locmax
715  ! -- local
716  integer(I4B) :: n
717  real(DP) :: botm
718  real(DP) :: xx
719  real(DP) :: dxx
720  !
721  ! -- Newton-Raphson under-relaxation
722  do n = 1, this%dis%nodes
723  if (this%ibound(n) < 1) cycle
724  if (this%icelltype(n) > 0) then
725  botm = this%dis%bot(this%ibotnode(n))
726  ! -- only apply Newton-Raphson under-relaxation if
727  ! solution head is below the bottom of the model
728  if (x(n) < botm) then
729  inewtonur = 1
730  xx = xtemp(n) * (done - dp9) + botm * dp9
731  dxx = x(n) - xx
732  if (abs(dxx) > abs(dxmax)) then
733  locmax = n
734  dxmax = dxx
735  end if
736  x(n) = xx
737  dx(n) = dzero
738  end if
739  end if
740  end do
741  end subroutine npf_nur
742 
743  !> @brief Calculate flowja
744  !<
745  subroutine npf_cq(this, hnew, flowja)
746  ! -- dummy
747  class(gwfnpftype) :: this
748  real(DP), intent(inout), dimension(:) :: hnew
749  real(DP), intent(inout), dimension(:) :: flowja
750  ! -- local
751  integer(I4B) :: n, ipos, m
752  real(DP) :: qnm
753  !
754  ! -- Calculate the flow across each cell face and store in flowja
755  !
756  if (this%ixt3d /= 0) then
757  call this%xt3d%xt3d_flowja(hnew, flowja)
758  else
759  !
760  do n = 1, this%dis%nodes
761  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
762  m = this%dis%con%ja(ipos)
763  if (m < n) cycle
764  call this%qcalc(n, m, hnew(n), hnew(m), ipos, qnm)
765  flowja(ipos) = qnm
766  flowja(this%dis%con%isym(ipos)) = -qnm
767  end do
768  end do
769  !
770  end if
771  end subroutine npf_cq
772 
773  !> @brief Fractional cell saturation
774  !<
775  subroutine sgwf_npf_thksat(this, n, hn, thksat)
776  ! -- dummy
777  class(gwfnpftype) :: this
778  integer(I4B), intent(in) :: n
779  real(DP), intent(in) :: hn
780  real(DP), intent(inout) :: thksat
781  !
782  ! -- Standard Formulation
783  if (hn >= this%dis%top(n)) then
784  thksat = done
785  else
786  thksat = (hn - this%dis%bot(n)) / (this%dis%top(n) - this%dis%bot(n))
787  end if
788  !
789  ! -- Newton-Raphson Formulation
790  if (this%inewton /= 0) then
791  thksat = squadraticsaturation(this%dis%top(n), this%dis%bot(n), hn, &
792  this%satomega)
793  end if
794  end subroutine sgwf_npf_thksat
795 
796  !> @brief Flow between two cells
797  !<
798  subroutine sgwf_npf_qcalc(this, n, m, hn, hm, icon, qnm)
799  ! -- dummy
800  class(gwfnpftype) :: this
801  integer(I4B), intent(in) :: n
802  integer(I4B), intent(in) :: m
803  real(DP), intent(in) :: hn
804  real(DP), intent(in) :: hm
805  integer(I4B), intent(in) :: icon
806  real(DP), intent(inout) :: qnm
807  ! -- local
808  real(DP) :: hyn, hym
809  real(DP) :: condnm
810  real(DP) :: hntemp, hmtemp
811  integer(I4B) :: ihc
812  !
813  ! -- Initialize
814  ihc = this%dis%con%ihc(this%dis%con%jas(icon))
815  hyn = this%hy_eff(n, m, ihc, ipos=icon)
816  hym = this%hy_eff(m, n, ihc, ipos=icon)
817  !
818  ! -- Calculate conductance
819  if (ihc == c3d_vertical) then
820  condnm = vcond(this%ibound(n), this%ibound(m), &
821  this%icelltype(n), this%icelltype(m), this%inewton, &
822  this%ivarcv, this%idewatcv, &
823  this%condsat(this%dis%con%jas(icon)), hn, hm, &
824  hyn, hym, &
825  this%sat(n), this%sat(m), &
826  this%dis%top(n), this%dis%top(m), &
827  this%dis%bot(n), this%dis%bot(m), &
828  this%dis%con%hwva(this%dis%con%jas(icon)))
829  else
830  condnm = hcond(this%ibound(n), this%ibound(m), &
831  this%icelltype(n), this%icelltype(m), &
832  this%inewton, &
833  this%dis%con%ihc(this%dis%con%jas(icon)), &
834  this%icellavg, &
835  this%condsat(this%dis%con%jas(icon)), &
836  hn, hm, this%sat(n), this%sat(m), hyn, hym, &
837  this%dis%top(n), this%dis%top(m), &
838  this%dis%bot(n), this%dis%bot(m), &
839  this%dis%con%cl1(this%dis%con%jas(icon)), &
840  this%dis%con%cl2(this%dis%con%jas(icon)), &
841  this%dis%con%hwva(this%dis%con%jas(icon)))
842  end if
843  !
844  ! -- Initialize hntemp and hmtemp
845  hntemp = hn
846  hmtemp = hm
847  !
848  ! -- Check and adjust for dewatered conditions
849  if (this%iperched /= 0) then
850  if (this%dis%con%ihc(this%dis%con%jas(icon)) == 0) then
851  if (n > m) then
852  if (this%icelltype(n) /= 0) then
853  if (hn < this%dis%top(n)) hntemp = this%dis%bot(m)
854  end if
855  else
856  if (this%icelltype(m) /= 0) then
857  if (hm < this%dis%top(m)) hmtemp = this%dis%bot(n)
858  end if
859  end if
860  end if
861  end if
862  !
863  ! -- Calculate flow positive into cell n
864  qnm = condnm * (hmtemp - hntemp)
865  end subroutine sgwf_npf_qcalc
866 
867  !> @brief Record flowja and calculate specific discharge if requested
868  !<
869  subroutine npf_save_model_flows(this, flowja, icbcfl, icbcun)
870  ! -- dummy
871  class(gwfnpftype) :: this
872  real(DP), dimension(:), intent(in) :: flowja
873  integer(I4B), intent(in) :: icbcfl
874  integer(I4B), intent(in) :: icbcun
875  ! -- local
876  integer(I4B) :: ibinun
877  !
878  ! -- Set unit number for binary output
879  if (this%ipakcb < 0) then
880  ibinun = icbcun
881  elseif (this%ipakcb == 0) then
882  ibinun = 0
883  else
884  ibinun = this%ipakcb
885  end if
886  if (icbcfl == 0) ibinun = 0
887  !
888  ! -- Write the face flows if requested
889  if (ibinun /= 0) then
890  call this%dis%record_connection_array(flowja, ibinun, this%iout)
891  end if
892  !
893  ! -- Calculate specific discharge at cell centers and write, if requested
894  if (this%isavspdis /= 0) then
895  if (ibinun /= 0) call this%sav_spdis(ibinun)
896  end if
897  !
898  ! -- Save saturation, if requested
899  if (this%isavsat /= 0) then
900  if (ibinun /= 0) call this%sav_sat(ibinun)
901  end if
902  end subroutine npf_save_model_flows
903 
904  !> @brief Print budget
905  !<
906  subroutine npf_print_model_flows(this, ibudfl, flowja)
907  ! -- modules
908  use tdismodule, only: kper, kstp
909  use constantsmodule, only: lenbigline
910  ! -- dummy
911  class(gwfnpftype) :: this
912  integer(I4B), intent(in) :: ibudfl
913  real(DP), intent(inout), dimension(:) :: flowja
914  ! -- local
915  character(len=LENBIGLINE) :: line
916  character(len=30) :: tempstr
917  integer(I4B) :: n, ipos, m
918  real(DP) :: qnm
919  ! -- formats
920  character(len=*), parameter :: fmtiprflow = &
921  &"(/,4x,'CALCULATED INTERCELL FLOW FOR PERIOD ', i0, ' STEP ', i0)"
922  !
923  ! -- Write flowja to list file if requested
924  if (ibudfl /= 0 .and. this%iprflow > 0) then
925  write (this%iout, fmtiprflow) kper, kstp
926  do n = 1, this%dis%nodes
927  line = ''
928  call this%dis%noder_to_string(n, tempstr)
929  line = trim(tempstr)//':'
930  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
931  m = this%dis%con%ja(ipos)
932  call this%dis%noder_to_string(m, tempstr)
933  line = trim(line)//' '//trim(tempstr)
934  qnm = flowja(ipos)
935  write (tempstr, '(1pg15.6)') qnm
936  line = trim(line)//' '//trim(adjustl(tempstr))
937  end do
938  write (this%iout, '(a)') trim(line)
939  end do
940  end if
941  end subroutine npf_print_model_flows
942 
943  !> @brief Deallocate variables
944  !<
945  subroutine npf_da(this)
946  ! -- modules
949  ! -- dummy
950  class(gwfnpftype) :: this
951 
952  ! free spdis work structure
953  if (this%icalcspdis == 1) call this%spdis_wa%destroy()
954  deallocate (this%spdis_wa)
955  !
956  ! -- Deallocate input memory
957  call memorystore_remove(this%name_model, 'NPF', idm_context)
958  !
959  ! -- TVK
960  if (this%intvk /= 0) then
961  call this%tvk%da()
962  deallocate (this%tvk)
963  end if
964  !
965  ! -- VSC
966  if (this%invsc /= 0) then
967  nullify (this%vsc)
968  end if
969  !
970  ! -- Strings
971  !
972  ! -- Scalars
973  call mem_deallocate(this%iname)
974  call mem_deallocate(this%ixt3d)
975  call mem_deallocate(this%ixt3drhs)
976  call mem_deallocate(this%satomega)
977  call mem_deallocate(this%hnoflo)
978  call mem_deallocate(this%hdry)
979  call mem_deallocate(this%icellavg)
980  call mem_deallocate(this%iavgkeff)
981  call mem_deallocate(this%ik22)
982  call mem_deallocate(this%ik33)
983  call mem_deallocate(this%iperched)
984  call mem_deallocate(this%ivarcv)
985  call mem_deallocate(this%idewatcv)
986  call mem_deallocate(this%ithickstrt)
987  call mem_deallocate(this%isavspdis)
988  call mem_deallocate(this%isavsat)
989  call mem_deallocate(this%icalcspdis)
990  call mem_deallocate(this%irewet)
991  call mem_deallocate(this%wetfct)
992  call mem_deallocate(this%iwetit)
993  call mem_deallocate(this%ihdwet)
994  call mem_deallocate(this%ibotnode)
995  call mem_deallocate(this%iwetdry)
996  call mem_deallocate(this%iangle1)
997  call mem_deallocate(this%iangle2)
998  call mem_deallocate(this%iangle3)
999  call mem_deallocate(this%nedges)
1000  call mem_deallocate(this%lastedge)
1001  call mem_deallocate(this%ik22overk)
1002  call mem_deallocate(this%ik33overk)
1003  call mem_deallocate(this%intvk)
1004  call mem_deallocate(this%invsc)
1005  call mem_deallocate(this%kchangeper)
1006  call mem_deallocate(this%kchangestp)
1007  !
1008  ! -- Deallocate arrays
1009  deallocate (this%aname)
1010  call mem_deallocate(this%ithickstartflag)
1011  call mem_deallocate(this%icelltype)
1012  call mem_deallocate(this%k11)
1013  call mem_deallocate(this%k22)
1014  call mem_deallocate(this%k33)
1015  call mem_deallocate(this%k11input)
1016  call mem_deallocate(this%k22input)
1017  call mem_deallocate(this%k33input)
1018  call mem_deallocate(this%sat, 'SAT', this%memoryPath)
1019  call mem_deallocate(this%condsat)
1020  call mem_deallocate(this%wetdry)
1021  call mem_deallocate(this%angle1)
1022  call mem_deallocate(this%angle2)
1023  call mem_deallocate(this%angle3)
1024  call mem_deallocate(this%nodedge)
1025  call mem_deallocate(this%ihcedge)
1026  call mem_deallocate(this%propsedge)
1027  call mem_deallocate(this%iedge_ptr)
1028  call mem_deallocate(this%edge_idxs)
1029  call mem_deallocate(this%spdis, 'SPDIS', this%memoryPath)
1030  call mem_deallocate(this%nodekchange)
1031  !
1032  ! -- deallocate parent
1033  call this%NumericalPackageType%da()
1034 
1035  ! pointers
1036  this%hnew => null()
1037 
1038  end subroutine npf_da
1039 
1040  !> @ brief Allocate scalars
1041  !!
1042  !! Allocate and initialize scalars for the VSC package. The base model
1043  !! allocate scalars method is also called.
1044  !<
1045  subroutine allocate_scalars(this)
1046  ! -- modules
1048  ! -- dummy
1049  class(gwfnpftype) :: this
1050  !
1051  ! -- allocate scalars in NumericalPackageType
1052  call this%NumericalPackageType%allocate_scalars()
1053  !
1054  ! -- Allocate scalars
1055  call mem_allocate(this%iname, 'INAME', this%memoryPath)
1056  call mem_allocate(this%ixt3d, 'IXT3D', this%memoryPath)
1057  call mem_allocate(this%ixt3drhs, 'IXT3DRHS', this%memoryPath)
1058  call mem_allocate(this%satomega, 'SATOMEGA', this%memoryPath)
1059  call mem_allocate(this%hnoflo, 'HNOFLO', this%memoryPath)
1060  call mem_allocate(this%hdry, 'HDRY', this%memoryPath)
1061  call mem_allocate(this%icellavg, 'ICELLAVG', this%memoryPath)
1062  call mem_allocate(this%iavgkeff, 'IAVGKEFF', this%memoryPath)
1063  call mem_allocate(this%ik22, 'IK22', this%memoryPath)
1064  call mem_allocate(this%ik33, 'IK33', this%memoryPath)
1065  call mem_allocate(this%ik22overk, 'IK22OVERK', this%memoryPath)
1066  call mem_allocate(this%ik33overk, 'IK33OVERK', this%memoryPath)
1067  call mem_allocate(this%iperched, 'IPERCHED', this%memoryPath)
1068  call mem_allocate(this%ivarcv, 'IVARCV', this%memoryPath)
1069  call mem_allocate(this%idewatcv, 'IDEWATCV', this%memoryPath)
1070  call mem_allocate(this%ithickstrt, 'ITHICKSTRT', this%memoryPath)
1071  call mem_allocate(this%icalcspdis, 'ICALCSPDIS', this%memoryPath)
1072  call mem_allocate(this%isavspdis, 'ISAVSPDIS', this%memoryPath)
1073  call mem_allocate(this%isavsat, 'ISAVSAT', this%memoryPath)
1074  call mem_allocate(this%irewet, 'IREWET', this%memoryPath)
1075  call mem_allocate(this%wetfct, 'WETFCT', this%memoryPath)
1076  call mem_allocate(this%iwetit, 'IWETIT', this%memoryPath)
1077  call mem_allocate(this%ihdwet, 'IHDWET', this%memoryPath)
1078  call mem_allocate(this%iangle1, 'IANGLE1', this%memoryPath)
1079  call mem_allocate(this%iangle2, 'IANGLE2', this%memoryPath)
1080  call mem_allocate(this%iangle3, 'IANGLE3', this%memoryPath)
1081  call mem_allocate(this%iwetdry, 'IWETDRY', this%memoryPath)
1082  call mem_allocate(this%nedges, 'NEDGES', this%memoryPath)
1083  call mem_allocate(this%lastedge, 'LASTEDGE', this%memoryPath)
1084  call mem_allocate(this%intvk, 'INTVK', this%memoryPath)
1085  call mem_allocate(this%invsc, 'INVSC', this%memoryPath)
1086  call mem_allocate(this%kchangeper, 'KCHANGEPER', this%memoryPath)
1087  call mem_allocate(this%kchangestp, 'KCHANGESTP', this%memoryPath)
1088  !
1089  ! -- set pointer to inewtonur
1090  call mem_setptr(this%igwfnewtonur, 'INEWTONUR', &
1091  create_mem_path(this%name_model))
1092  !
1093  ! -- Initialize value
1094  this%iname = 8
1095  this%ixt3d = 0
1096  this%ixt3drhs = 0
1097  this%satomega = dzero
1098  this%hnoflo = dhnoflo !1.d30
1099  this%hdry = dhdry !-1.d30
1100  this%icellavg = ccond_hmean
1101  this%iavgkeff = 0
1102  this%ik22 = 0
1103  this%ik33 = 0
1104  this%ik22overk = 0
1105  this%ik33overk = 0
1106  this%iperched = 0
1107  this%ivarcv = 0
1108  this%idewatcv = 0
1109  this%ithickstrt = 0
1110  this%icalcspdis = 0
1111  this%isavspdis = 0
1112  this%isavsat = 0
1113  this%irewet = 0
1114  this%wetfct = done
1115  this%iwetit = 1
1116  this%ihdwet = 0
1117  this%iangle1 = 0
1118  this%iangle2 = 0
1119  this%iangle3 = 0
1120  this%iwetdry = 0
1121  this%nedges = 0
1122  this%lastedge = 0
1123  this%intvk = 0
1124  this%invsc = 0
1125  this%kchangeper = 0
1126  this%kchangestp = 0
1127  !
1128  ! -- If newton is on, then NPF creates asymmetric matrix
1129  this%iasym = this%inewton
1130  end subroutine allocate_scalars
1131 
1132  !> @ brief Store backup copy of hydraulic conductivity when the VSC
1133  !! package is activate
1134  !!
1135  !! The K arrays (K11, etc.) get multiplied by the viscosity ratio so that
1136  !! subsequent uses of K already take into account the effect of viscosity.
1137  !! Thus the original user-specified K array values are lost unless they are
1138  !! backed up in k11input, for example. In a new stress period/time step,
1139  !! the values in k11input are multiplied by the viscosity ratio, not k11
1140  !! since it contains viscosity-adjusted hydraulic conductivity values.
1141  !<
1142  subroutine store_original_k_arrays(this, ncells, njas)
1143  ! -- modules
1145  ! -- dummy
1146  class(gwfnpftype) :: this
1147  integer(I4B), intent(in) :: ncells
1148  integer(I4B), intent(in) :: njas
1149  ! -- local
1150  integer(I4B) :: n
1151  !
1152  ! -- Retain copy of user-specified K arrays
1153  do n = 1, ncells
1154  this%k11input(n) = this%k11(n)
1155  this%k22input(n) = this%k22(n)
1156  this%k33input(n) = this%k33(n)
1157  end do
1158  end subroutine store_original_k_arrays
1159 
1160  !> @brief Allocate npf arrays
1161  !<
1162  subroutine allocate_arrays(this, ncells, njas)
1163  ! -- dummy
1164  class(gwfnpftype) :: this
1165  integer(I4B), intent(in) :: ncells
1166  integer(I4B), intent(in) :: njas
1167  ! -- local
1168  integer(I4B) :: n
1169  !
1170  call mem_allocate(this%ithickstartflag, ncells, 'ITHICKSTARTFLAG', &
1171  this%memoryPath)
1172  call mem_allocate(this%icelltype, ncells, 'ICELLTYPE', this%memoryPath)
1173  call mem_allocate(this%k11, ncells, 'K11', this%memoryPath)
1174  call mem_allocate(this%sat, ncells, 'SAT', this%memoryPath)
1175  call mem_allocate(this%condsat, njas, 'CONDSAT', this%memoryPath)
1176  !
1177  ! -- Optional arrays dimensioned to full size initially
1178  call mem_allocate(this%k22, ncells, 'K22', this%memoryPath)
1179  call mem_allocate(this%k33, ncells, 'K33', this%memoryPath)
1180  call mem_allocate(this%wetdry, ncells, 'WETDRY', this%memoryPath)
1181  call mem_allocate(this%angle1, ncells, 'ANGLE1', this%memoryPath)
1182  call mem_allocate(this%angle2, ncells, 'ANGLE2', this%memoryPath)
1183  call mem_allocate(this%angle3, ncells, 'ANGLE3', this%memoryPath)
1184  !
1185  ! -- Optional arrays
1186  call mem_allocate(this%ibotnode, 0, 'IBOTNODE', this%memoryPath)
1187  call mem_allocate(this%nodedge, 0, 'NODEDGE', this%memoryPath)
1188  call mem_allocate(this%ihcedge, 0, 'IHCEDGE', this%memoryPath)
1189  call mem_allocate(this%propsedge, 0, 0, 'PROPSEDGE', this%memoryPath)
1190  call mem_allocate(this%iedge_ptr, 0, 'NREDGESNODE', this%memoryPath)
1191  call mem_allocate(this%edge_idxs, 0, 'EDGEIDXS', this%memoryPath)
1192  !
1193  ! -- Optional arrays only needed when vsc package is active
1194  call mem_allocate(this%k11input, 0, 'K11INPUT', this%memoryPath)
1195  call mem_allocate(this%k22input, 0, 'K22INPUT', this%memoryPath)
1196  call mem_allocate(this%k33input, 0, 'K33INPUT', this%memoryPath)
1197  !
1198  ! -- Specific discharge is (re-)allocated when nedges is known
1199  call mem_allocate(this%spdis, 3, 0, 'SPDIS', this%memoryPath)
1200  !
1201  ! -- Time-varying property flag arrays
1202  call mem_allocate(this%nodekchange, ncells, 'NODEKCHANGE', this%memoryPath)
1203  !
1204  ! -- initialize iangle1, iangle2, iangle3, and wetdry
1205  do n = 1, ncells
1206  this%angle1(n) = dzero
1207  this%angle2(n) = dzero
1208  this%angle3(n) = dzero
1209  this%wetdry(n) = dzero
1210  this%nodekchange(n) = dzero
1211  end do
1212  !
1213  ! -- allocate variable names
1214  allocate (this%aname(this%iname))
1215  this%aname = [' ICELLTYPE', ' K', &
1216  ' K33', ' K22', &
1217  ' WETDRY', ' ANGLE1', &
1218  ' ANGLE2', ' ANGLE3']
1219  end subroutine allocate_arrays
1220 
1221  !> @brief Log npf options sourced from the input mempath
1222  !<
1223  subroutine log_options(this, found)
1224  ! -- modules
1225  use kindmodule, only: lgp
1227  ! -- dummy
1228  class(gwfnpftype) :: this
1229  ! -- locals
1230  type(gwfnpfparamfoundtype), intent(in) :: found
1231  !
1232  write (this%iout, '(1x,a)') 'Setting NPF Options'
1233  if (found%iprflow) &
1234  write (this%iout, '(4x,a)') 'Cell-by-cell flow information will be printed &
1235  &to listing file whenever ICBCFL is not zero.'
1236  if (found%ipakcb) &
1237  write (this%iout, '(4x,a)') 'Cell-by-cell flow information will be saved &
1238  &to binary file whenever ICBCFL is not zero.'
1239  if (found%cellavg) &
1240  write (this%iout, '(4x,a,i0)') 'Alternative cell averaging [1=logarithmic, &
1241  &2=AMT-LMK, 3=AMT-HMK] set to: ', &
1242  this%icellavg
1243  if (found%ithickstrt) &
1244  write (this%iout, '(4x,a)') 'THICKSTRT option has been activated.'
1245  if (found%iperched) &
1246  write (this%iout, '(4x,a)') 'Vertical flow will be adjusted for perched &
1247  &conditions.'
1248  if (found%ivarcv) &
1249  write (this%iout, '(4x,a)') 'Vertical conductance varies with water table.'
1250  if (found%idewatcv) &
1251  write (this%iout, '(4x,a)') 'Vertical conductance accounts for dewatered &
1252  &portion of an underlying cell.'
1253  if (found%ixt3d) write (this%iout, '(4x,a)') 'XT3D formulation is selected.'
1254  if (found%ixt3drhs) &
1255  write (this%iout, '(4x,a)') 'XT3D RHS formulation is selected.'
1256  if (found%isavspdis) &
1257  write (this%iout, '(4x,a)') 'Specific discharge will be calculated at cell &
1258  &centers and written to DATA-SPDIS in budget &
1259  &file when requested.'
1260  if (found%isavsat) &
1261  write (this%iout, '(4x,a)') 'Saturation will be written to DATA-SAT in &
1262  &budget file when requested.'
1263  if (found%ik22overk) &
1264  write (this%iout, '(4x,a)') 'Values specified for K22 are anisotropy &
1265  &ratios and will be multiplied by K before &
1266  &being used in calculations.'
1267  if (found%ik33overk) &
1268  write (this%iout, '(4x,a)') 'Values specified for K33 are anisotropy &
1269  &ratios and will be multiplied by K before &
1270  &being used in calculations.'
1271  if (found%inewton) &
1272  write (this%iout, '(4x,a)') 'NEWTON-RAPHSON method disabled for unconfined &
1273  &cells'
1274  if (found%satomega) &
1275  write (this%iout, '(4x,a,1pg15.6)') 'Saturation omega: ', this%satomega
1276  if (found%irewet) &
1277  write (this%iout, '(4x,a)') 'Rewetting is active.'
1278  if (found%wetfct) &
1279  write (this%iout, '(4x,a,1pg15.6)') &
1280  'Wetting factor (WETFCT) has been set to: ', this%wetfct
1281  if (found%iwetit) &
1282  write (this%iout, '(4x,a,i5)') &
1283  'Wetting iteration interval (IWETIT) has been set to: ', this%iwetit
1284  if (found%ihdwet) &
1285  write (this%iout, '(4x,a,i5)') &
1286  'Head rewet equation (IHDWET) has been set to: ', this%ihdwet
1287  write (this%iout, '(1x,a,/)') 'End Setting NPF Options'
1288  end subroutine log_options
1289 
1290  !> @brief Update simulation options from input mempath
1291  !<
1292  subroutine source_options(this)
1293  ! -- modules
1299  use sourcecommonmodule, only: filein_fname
1300  ! -- dummy
1301  class(gwfnpftype) :: this
1302  ! -- locals
1303  character(len=LENVARNAME), dimension(3) :: cellavg_method = &
1304  &[character(len=LENVARNAME) :: 'LOGARITHMIC', 'AMT-LMK', 'AMT-HMK']
1305  type(gwfnpfparamfoundtype) :: found
1306  character(len=LINELENGTH) :: tvk6_filename
1307  !
1308  ! -- update defaults with idm sourced values
1309  call mem_set_value(this%iprflow, 'IPRFLOW', this%input_mempath, found%iprflow)
1310  call mem_set_value(this%ipakcb, 'IPAKCB', this%input_mempath, found%ipakcb)
1311  call mem_set_value(this%icellavg, 'CELLAVG', this%input_mempath, &
1312  cellavg_method, found%cellavg)
1313  call mem_set_value(this%ithickstrt, 'ITHICKSTRT', this%input_mempath, &
1314  found%ithickstrt)
1315  call mem_set_value(this%iperched, 'IPERCHED', this%input_mempath, &
1316  found%iperched)
1317  call mem_set_value(this%ivarcv, 'IVARCV', this%input_mempath, found%ivarcv)
1318  call mem_set_value(this%idewatcv, 'IDEWATCV', this%input_mempath, &
1319  found%idewatcv)
1320  call mem_set_value(this%ixt3d, 'IXT3D', this%input_mempath, found%ixt3d)
1321  call mem_set_value(this%ixt3drhs, 'IXT3DRHS', this%input_mempath, &
1322  found%ixt3drhs)
1323  call mem_set_value(this%isavspdis, 'ISAVSPDIS', this%input_mempath, &
1324  found%isavspdis)
1325  call mem_set_value(this%isavsat, 'ISAVSAT', this%input_mempath, found%isavsat)
1326  call mem_set_value(this%ik22overk, 'IK22OVERK', this%input_mempath, &
1327  found%ik22overk)
1328  call mem_set_value(this%ik33overk, 'IK33OVERK', this%input_mempath, &
1329  found%ik33overk)
1330  call mem_set_value(this%inewton, 'INEWTON', this%input_mempath, found%inewton)
1331  call mem_set_value(this%satomega, 'SATOMEGA', this%input_mempath, &
1332  found%satomega)
1333  call mem_set_value(this%irewet, 'IREWET', this%input_mempath, found%irewet)
1334  call mem_set_value(this%wetfct, 'WETFCT', this%input_mempath, found%wetfct)
1335  call mem_set_value(this%iwetit, 'IWETIT', this%input_mempath, found%iwetit)
1336  call mem_set_value(this%ihdwet, 'IHDWET', this%input_mempath, found%ihdwet)
1337  !
1338  ! -- save flows option active
1339  if (found%ipakcb) this%ipakcb = -1
1340  !
1341  ! -- xt3d active with rhs
1342  if (found%ixt3d .and. found%ixt3drhs) this%ixt3d = 2
1343  !
1344  ! -- save specific discharge active
1345  if (found%isavspdis) this%icalcspdis = this%isavspdis
1346  !
1347  ! -- no newton specified
1348  if (found%inewton) then
1349  this%inewton = 0
1350  this%iasym = 0
1351  end if
1352  !
1353  ! -- enforce 0 or 1 TVK6_FILENAME entries in option block
1354  if (filein_fname(tvk6_filename, 'TVK6_FILENAME', this%input_mempath, &
1355  this%input_fname)) then
1356  call openfile(this%intvk, this%iout, tvk6_filename, 'TVK')
1357  call tvk_cr(this%tvk, this%name_model, this%intvk, this%iout)
1358  end if
1359  !
1360  ! -- log options
1361  if (this%iout > 0) then
1362  call this%log_options(found)
1363  end if
1364  end subroutine source_options
1365 
1366  !> @brief Set options in the NPF object
1367  !<
1368  subroutine set_options(this, options)
1369  ! -- dummy
1370  class(gwfnpftype) :: this
1371  type(gwfnpfoptionstype), intent(in) :: options
1372  !
1373  this%icellavg = options%icellavg
1374  this%ithickstrt = options%ithickstrt
1375  this%iperched = options%iperched
1376  this%ivarcv = options%ivarcv
1377  this%idewatcv = options%idewatcv
1378  this%irewet = options%irewet
1379  this%wetfct = options%wetfct
1380  this%iwetit = options%iwetit
1381  this%ihdwet = options%ihdwet
1382  end subroutine set_options
1383 
1384  !> @brief Check for conflicting NPF options
1385  !<
1386  subroutine check_options(this)
1387  ! -- modules
1389  use constantsmodule, only: linelength
1390  ! -- dummy
1391  class(gwfnpftype) :: this
1392  !
1393  ! -- set omega value used for saturation calculations
1394  if (this%inewton > 0) then
1395  this%satomega = dem6
1396  end if
1397  !
1398  if (this%inewton > 0) then
1399  if (this%iperched > 0) then
1400  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1401  'BE USED WITH PERCHED OPTION.'
1402  call store_error(errmsg)
1403  end if
1404  if (this%ivarcv > 0) then
1405  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1406  'BE USED WITH VARIABLECV OPTION.'
1407  call store_error(errmsg)
1408  end if
1409  if (this%irewet > 0) then
1410  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1411  'BE USED WITH REWET OPTION.'
1412  call store_error(errmsg)
1413  end if
1414  end if
1415  !
1416  if (this%ixt3d /= 0) then
1417  if (this%icellavg > 0) then
1418  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. '// &
1419  'ALTERNATIVE_CELL_AVERAGING OPTION '// &
1420  'CANNOT BE USED WITH XT3D OPTION.'
1421  call store_error(errmsg)
1422  end if
1423  if (this%ithickstrt > 0) then
1424  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. THICKSTRT OPTION '// &
1425  'CANNOT BE USED WITH XT3D OPTION.'
1426  call store_error(errmsg)
1427  end if
1428  if (this%iperched > 0) then
1429  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. PERCHED OPTION '// &
1430  'CANNOT BE USED WITH XT3D OPTION.'
1431  call store_error(errmsg)
1432  end if
1433  if (this%ivarcv > 0) then
1434  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. VARIABLECV OPTION '// &
1435  'CANNOT BE USED WITH XT3D OPTION.'
1436  call store_error(errmsg)
1437  end if
1438  end if
1439  !
1440  ! -- Terminate if errors
1441  if (count_errors() > 0) then
1442  call store_error_filename(this%input_fname)
1443  end if
1444  end subroutine check_options
1445 
1446  !> @brief Write dimensions to list file
1447  !<
1448  subroutine log_griddata(this, found)
1449  ! -- modules
1451  ! -- dummy
1452  class(gwfnpftype) :: this
1453  type(gwfnpfparamfoundtype), intent(in) :: found
1454  !
1455  write (this%iout, '(1x,a)') 'Setting NPF Griddata'
1456  !
1457  if (found%icelltype) then
1458  write (this%iout, '(4x,a)') 'ICELLTYPE set from input file'
1459  end if
1460  !
1461  if (found%k) then
1462  write (this%iout, '(4x,a)') 'K set from input file'
1463  end if
1464  !
1465  if (found%k33) then
1466  write (this%iout, '(4x,a)') 'K33 set from input file'
1467  else
1468  write (this%iout, '(4x,a)') 'K33 not provided. Setting K33 = K.'
1469  end if
1470  !
1471  if (found%k22) then
1472  write (this%iout, '(4x,a)') 'K22 set from input file'
1473  else
1474  write (this%iout, '(4x,a)') 'K22 not provided. Setting K22 = K.'
1475  end if
1476  !
1477  if (found%wetdry) then
1478  write (this%iout, '(4x,a)') 'WETDRY set from input file'
1479  end if
1480  !
1481  if (found%angle1) then
1482  write (this%iout, '(4x,a)') 'ANGLE1 set from input file'
1483  end if
1484  !
1485  if (found%angle2) then
1486  write (this%iout, '(4x,a)') 'ANGLE2 set from input file'
1487  end if
1488  !
1489  if (found%angle3) then
1490  write (this%iout, '(4x,a)') 'ANGLE3 set from input file'
1491  end if
1492  !
1493  write (this%iout, '(1x,a,/)') 'End Setting NPF Griddata'
1494  end subroutine log_griddata
1495 
1496  !> @brief Update simulation griddata from input mempath
1497  !<
1498  subroutine source_griddata(this)
1499  ! -- modules
1500  use simmodule, only: count_errors, store_error
1504  ! -- dummy
1505  class(gwfnpftype) :: this
1506  ! -- locals
1507  character(len=LINELENGTH) :: errmsg
1508  type(gwfnpfparamfoundtype) :: found
1509  logical, dimension(2) :: afound
1510  integer(I4B), dimension(:), pointer, contiguous :: map
1511  !
1512  ! -- set map to convert user input data into reduced data
1513  map => null()
1514  if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
1515  !
1516  ! -- update defaults with idm sourced values
1517  call mem_set_value(this%icelltype, 'ICELLTYPE', this%input_mempath, map, &
1518  found%icelltype)
1519  call mem_set_value(this%k11, 'K', this%input_mempath, map, found%k)
1520  call mem_set_value(this%k33, 'K33', this%input_mempath, map, found%k33)
1521  call mem_set_value(this%k22, 'K22', this%input_mempath, map, found%k22)
1522  call mem_set_value(this%wetdry, 'WETDRY', this%input_mempath, map, &
1523  found%wetdry)
1524  call mem_set_value(this%angle1, 'ANGLE1', this%input_mempath, map, &
1525  found%angle1)
1526  call mem_set_value(this%angle2, 'ANGLE2', this%input_mempath, map, &
1527  found%angle2)
1528  call mem_set_value(this%angle3, 'ANGLE3', this%input_mempath, map, &
1529  found%angle3)
1530  !
1531  ! -- ensure ICELLTYPE was found
1532  if (.not. found%icelltype) then
1533  write (errmsg, '(a)') 'Error in GRIDDATA block: ICELLTYPE not found.'
1534  call store_error(errmsg)
1535  end if
1536  !
1537  ! -- ensure K was found
1538  if (.not. found%k) then
1539  write (errmsg, '(a)') 'Error in GRIDDATA block: K not found.'
1540  call store_error(errmsg)
1541  end if
1542  !
1543  ! -- set error if ik33overk set with no k33
1544  if (.not. found%k33 .and. this%ik33overk /= 0) then
1545  write (errmsg, '(a)') 'K33OVERK option specified but K33 not specified.'
1546  call store_error(errmsg)
1547  end if
1548  !
1549  ! -- set error if ik22overk set with no k22
1550  if (.not. found%k22 .and. this%ik22overk /= 0) then
1551  write (errmsg, '(a)') 'K22OVERK option specified but K22 not specified.'
1552  call store_error(errmsg)
1553  end if
1554  !
1555  ! -- handle found side effects
1556  if (found%k33) this%ik33 = 1
1557  if (found%k22) this%ik22 = 1
1558  if (found%wetdry) this%iwetdry = 1
1559  if (found%angle1) this%iangle1 = 1
1560  if (found%angle2) this%iangle2 = 1
1561  if (found%angle3) this%iangle3 = 1
1562  !
1563  ! -- handle not found side effects
1564  if (.not. found%k33) then
1565  call mem_set_value(this%k33, 'K', this%input_mempath, map, afound(1))
1566  end if
1567  if (.not. found%k22) then
1568  call mem_set_value(this%k22, 'K', this%input_mempath, map, afound(2))
1569  end if
1570  if (.not. found%wetdry) call mem_reallocate(this%wetdry, 1, 'WETDRY', &
1571  trim(this%memoryPath))
1572  if (.not. found%angle1 .and. this%ixt3d == 0) &
1573  call mem_reallocate(this%angle1, 0, 'ANGLE1', trim(this%memoryPath))
1574  if (.not. found%angle2 .and. this%ixt3d == 0) &
1575  call mem_reallocate(this%angle2, 0, 'ANGLE2', trim(this%memoryPath))
1576  if (.not. found%angle3 .and. this%ixt3d == 0) &
1577  call mem_reallocate(this%angle3, 0, 'ANGLE3', trim(this%memoryPath))
1578  !
1579  ! -- log griddata
1580  if (this%iout > 0) then
1581  call this%log_griddata(found)
1582  end if
1583  end subroutine source_griddata
1584 
1585  !> @brief Initialize and check NPF data
1586  !<
1587  subroutine prepcheck(this)
1588  ! -- modules
1589  use constantsmodule, only: linelength, dpio180
1591  ! -- dummy
1592  class(gwfnpftype) :: this
1593  ! -- local
1594  character(len=24), dimension(:), pointer :: aname
1595  character(len=LINELENGTH) :: cellstr, errmsg
1596  integer(I4B) :: nerr, n
1597  ! -- format
1598  character(len=*), parameter :: fmtkerr = &
1599  &"(1x, 'Hydraulic property ',a,' is <= 0 for cell ',a, ' ', 1pg15.6)"
1600  character(len=*), parameter :: fmtkerr2 = &
1601  &"(1x, '... ', i0,' additional errors not shown for ',a)"
1602  !
1603  ! -- initialize
1604  aname => this%aname
1605  !
1606  ! -- check k11
1607  nerr = 0
1608  do n = 1, size(this%k11)
1609  if (this%k11(n) <= dzero) then
1610  nerr = nerr + 1
1611  if (nerr <= 20) then
1612  call this%dis%noder_to_string(n, cellstr)
1613  write (errmsg, fmtkerr) trim(adjustl(aname(2))), trim(cellstr), &
1614  this%k11(n)
1615  call store_error(errmsg)
1616  end if
1617  end if
1618  end do
1619  if (nerr > 20) then
1620  write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(2)))
1621  call store_error(errmsg)
1622  end if
1623  !
1624  ! -- check k33 because it was read
1625  if (this%ik33 /= 0) then
1626  !
1627  ! -- Check to make sure values are greater than or equal to zero
1628  nerr = 0
1629  do n = 1, size(this%k33)
1630  if (this%ik33overk /= 0) this%k33(n) = this%k33(n) * this%k11(n)
1631  if (this%k33(n) <= dzero) then
1632  nerr = nerr + 1
1633  if (nerr <= 20) then
1634  call this%dis%noder_to_string(n, cellstr)
1635  write (errmsg, fmtkerr) trim(adjustl(aname(3))), trim(cellstr), &
1636  this%k33(n)
1637  call store_error(errmsg)
1638  end if
1639  end if
1640  end do
1641  if (nerr > 20) then
1642  write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(3)))
1643  call store_error(errmsg)
1644  end if
1645  end if
1646  !
1647  ! -- check k22 because it was read
1648  if (this%ik22 /= 0) then
1649  !
1650  ! -- Check to make sure that angles are available
1651  if (this%dis%con%ianglex == 0) then
1652  write (errmsg, '(a)') 'Error. ANGLDEGX not provided in '// &
1653  'discretization file, but K22 was specified. '
1654  call store_error(errmsg)
1655  end if
1656  !
1657  ! -- Check to make sure values are greater than or equal to zero
1658  nerr = 0
1659  do n = 1, size(this%k22)
1660  if (this%ik22overk /= 0) this%k22(n) = this%k22(n) * this%k11(n)
1661  if (this%k22(n) <= dzero) then
1662  nerr = nerr + 1
1663  if (nerr <= 20) then
1664  call this%dis%noder_to_string(n, cellstr)
1665  write (errmsg, fmtkerr) trim(adjustl(aname(4))), trim(cellstr), &
1666  this%k22(n)
1667  call store_error(errmsg)
1668  end if
1669  end if
1670  end do
1671  if (nerr > 20) then
1672  write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(4)))
1673  call store_error(errmsg)
1674  end if
1675  end if
1676  !
1677  ! -- check for wetdry conflicts
1678  if (this%irewet == 1) then
1679  if (this%iwetdry == 0) then
1680  write (errmsg, '(a, a, a)') 'Error in GRIDDATA block: ', &
1681  trim(adjustl(aname(5))), ' not found.'
1682  call store_error(errmsg)
1683  end if
1684  end if
1685  !
1686  ! -- Check for angle conflicts
1687  if (this%iangle1 /= 0) then
1688  do n = 1, size(this%angle1)
1689  this%angle1(n) = this%angle1(n) * dpio180
1690  end do
1691  else
1692  if (this%ixt3d /= 0) then
1693  this%iangle1 = 1
1694  write (this%iout, '(a)') 'XT3D IN USE, BUT ANGLE1 NOT SPECIFIED. '// &
1695  'SETTING ANGLE1 TO ZERO.'
1696  do n = 1, size(this%angle1)
1697  this%angle1(n) = dzero
1698  end do
1699  end if
1700  end if
1701  if (this%iangle2 /= 0) then
1702  if (this%iangle1 == 0) then
1703  write (errmsg, '(a)') 'ANGLE2 SPECIFIED BUT NOT ANGLE1. '// &
1704  'ANGLE2 REQUIRES ANGLE1. '
1705  call store_error(errmsg)
1706  end if
1707  if (this%iangle3 == 0) then
1708  write (errmsg, '(a)') 'ANGLE2 SPECIFIED BUT NOT ANGLE3. '// &
1709  'SPECIFY BOTH OR NEITHER ONE. '
1710  call store_error(errmsg)
1711  end if
1712  do n = 1, size(this%angle2)
1713  this%angle2(n) = this%angle2(n) * dpio180
1714  end do
1715  end if
1716  if (this%iangle3 /= 0) then
1717  if (this%iangle1 == 0) then
1718  write (errmsg, '(a)') 'ANGLE3 SPECIFIED BUT NOT ANGLE1. '// &
1719  'ANGLE3 REQUIRES ANGLE1. '
1720  call store_error(errmsg)
1721  end if
1722  if (this%iangle2 == 0) then
1723  write (errmsg, '(a)') 'ANGLE3 SPECIFIED BUT NOT ANGLE2. '// &
1724  'SPECIFY BOTH OR NEITHER ONE. '
1725  call store_error(errmsg)
1726  end if
1727  do n = 1, size(this%angle3)
1728  this%angle3(n) = this%angle3(n) * dpio180
1729  end do
1730  end if
1731  !
1732  ! -- terminate if data errors
1733  if (count_errors() > 0) then
1734  call store_error_filename(this%input_fname)
1735  end if
1736  end subroutine prepcheck
1737 
1738  !> @brief preprocess the NPF input data
1739  !!
1740  !! This routine consists of the following steps:
1741  !!
1742  !! 1. convert cells to noflow when all transmissive parameters equal zero
1743  !! 2. perform initial wetting and drying
1744  !! 3. initialize cell saturation
1745  !! 4. calculate saturated conductance (when not xt3d)
1746  !! 5. If NEWTON under-relaxation, determine lower most node
1747  !<
1748  subroutine preprocess_input(this)
1749  ! -- modules
1750  use constantsmodule, only: linelength
1752  ! -- dummy
1753  class(gwfnpftype) :: this !< the instance of the NPF package
1754  ! -- local
1755  integer(I4B) :: n, m, ii, nn
1756  real(DP) :: hyn, hym
1757  real(DP) :: satn, topn, botn
1758  integer(I4B) :: nextn
1759  real(DP) :: minbot, botm
1760  logical :: finished
1761  character(len=LINELENGTH) :: cellstr, errmsg
1762  ! -- format
1763  character(len=*), parameter :: fmtcnv = &
1764  "(1X,'CELL ', A, &
1765  &' ELIMINATED BECAUSE ALL HYDRAULIC CONDUCTIVITIES TO NODE ARE 0.')"
1766  character(len=*), parameter :: fmtnct = &
1767  &"(1X,'Negative cell thickness at cell ', A)"
1768  character(len=*), parameter :: fmtihbe = &
1769  &"(1X,'Initial head, bottom elevation:',1P,2G13.5)"
1770  character(len=*), parameter :: fmttebe = &
1771  &"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
1772  !
1773  do n = 1, this%dis%nodes
1774  this%ithickstartflag(n) = 0
1775  end do
1776  !
1777  ! -- Insure that each cell has at least one non-zero transmissive parameter
1778  ! Note that a cell can be deactivated even if it has a valid connection
1779  ! to another model.
1780  nodeloop: do n = 1, this%dis%nodes
1781  !
1782  ! -- Skip if already inactive
1783  if (this%ibound(n) == 0) then
1784  if (this%irewet /= 0) then
1785  if (this%wetdry(n) == dzero) cycle nodeloop
1786  else
1787  cycle nodeloop
1788  end if
1789  end if
1790  !
1791  ! -- Cycle if k11 is not zero
1792  if (this%k11(n) /= dzero) cycle nodeloop
1793  !
1794  ! -- Cycle if at least one vertical connection has non-zero k33
1795  ! for n and m
1796  do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1797  m = this%dis%con%ja(ii)
1798  if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0) then
1799  hyn = this%k11(n)
1800  if (this%ik33 /= 0) hyn = this%k33(n)
1801  if (hyn /= dzero) then
1802  hym = this%k11(m)
1803  if (this%ik33 /= 0) hym = this%k33(m)
1804  if (hym /= dzero) cycle
1805  end if
1806  end if
1807  end do
1808  !
1809  ! -- If this part of the loop is reached, then all connections have
1810  ! zero transmissivity, so convert to noflow.
1811  this%ibound(n) = 0
1812  this%hnew(n) = this%hnoflo
1813  if (this%irewet /= 0) this%wetdry(n) = dzero
1814  call this%dis%noder_to_string(n, cellstr)
1815  write (this%iout, fmtcnv) trim(adjustl(cellstr))
1816  !
1817  end do nodeloop
1818  !
1819  ! -- Preprocess cell status and heads based on initial conditions
1820  if (this%inewton == 0) then
1821  !
1822  ! -- For standard formulation (non-Newton) call wetdry routine
1823  call this%wd(0, this%hnew)
1824  else
1825  !
1826  ! -- Newton formulation, so adjust heads to be above bottom
1827  ! (Not used in present formulation because variable cv
1828  ! cannot be used with Newton)
1829  if (this%ivarcv == 1) then
1830  do n = 1, this%dis%nodes
1831  if (this%hnew(n) < this%dis%bot(n)) then
1832  this%hnew(n) = this%dis%bot(n) + dem6
1833  end if
1834  end do
1835  end if
1836  end if
1837  !
1838  ! -- If THCKSTRT is not active, then loop through icelltype and replace
1839  ! any negative values with 1.
1840  if (this%ithickstrt == 0) then
1841  do n = 1, this%dis%nodes
1842  if (this%icelltype(n) < 0) then
1843  this%icelltype(n) = 1
1844  end if
1845  end do
1846  end if
1847  !
1848  ! -- Initialize sat to zero for ibound=0 cells, unless the cell can
1849  ! rewet. Initialize sat to the saturated fraction based on strt
1850  ! if icelltype is negative and the THCKSTRT option is in effect.
1851  ! Initialize sat to 1.0 for all other cells in order to calculate
1852  ! condsat in next section.
1853  do n = 1, this%dis%nodes
1854  if (this%ibound(n) == 0) then
1855  this%sat(n) = done
1856  if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0) then
1857  this%ithickstartflag(n) = 1
1858  this%icelltype(n) = 0
1859  end if
1860  else
1861  topn = this%dis%top(n)
1862  botn = this%dis%bot(n)
1863  if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0) then
1864  call this%thksat(n, this%ic%strt(n), satn)
1865  if (botn > this%ic%strt(n)) then
1866  call this%dis%noder_to_string(n, cellstr)
1867  write (errmsg, fmtnct) trim(adjustl(cellstr))
1868  call store_error(errmsg)
1869  write (errmsg, fmtihbe) this%ic%strt(n), botn
1870  call store_error(errmsg)
1871  end if
1872  this%ithickstartflag(n) = 1
1873  this%icelltype(n) = 0
1874  else
1875  satn = done
1876  if (botn > topn) then
1877  call this%dis%noder_to_string(n, cellstr)
1878  write (errmsg, fmtnct) trim(adjustl(cellstr))
1879  call store_error(errmsg)
1880  write (errmsg, fmttebe) topn, botn
1881  call store_error(errmsg)
1882  end if
1883  end if
1884  this%sat(n) = satn
1885  end if
1886  end do
1887  if (count_errors() > 0) then
1888  call store_error_filename(this%input_fname)
1889  end if
1890  !
1891  ! -- Calculate condsat, but only if xt3d is not active. If xt3d is
1892  ! active, then condsat is allocated to size of zero.
1893  if (this%ixt3d == 0) then
1894  !
1895  ! -- Calculate the saturated conductance for all connections assuming
1896  ! that saturation is 1 (except for case where icelltype was entered
1897  ! as a negative value and THCKSTRT option in effect)
1898  do n = 1, this%dis%nodes
1899  call this%calc_condsat(n, .true.)
1900  end do
1901  !
1902  end if
1903  !
1904  ! -- Determine the lower most node
1905  if (this%igwfnewtonur /= 0) then
1906  call mem_reallocate(this%ibotnode, this%dis%nodes, 'IBOTNODE', &
1907  trim(this%memoryPath))
1908  do n = 1, this%dis%nodes
1909  !
1910  minbot = this%dis%bot(n)
1911  nn = n
1912  finished = .false.
1913  do while (.not. finished)
1914  nextn = 0
1915  !
1916  ! -- Go through the connecting cells
1917  do ii = this%dis%con%ia(nn) + 1, this%dis%con%ia(nn + 1) - 1
1918  !
1919  ! -- Set the m cell number
1920  m = this%dis%con%ja(ii)
1921  botm = this%dis%bot(m)
1922  !
1923  ! -- select vertical connections: ihc == 0
1924  if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0) then
1925  if (m > nn .and. botm < minbot) then
1926  nextn = m
1927  minbot = botm
1928  end if
1929  end if
1930  end do
1931  if (nextn > 0) then
1932  nn = nextn
1933  else
1934  finished = .true.
1935  end if
1936  end do
1937  this%ibotnode(n) = nn
1938  end do
1939  end if
1940  !
1941  ! -- nullify unneeded gwf pointers
1942  this%igwfnewtonur => null()
1943  end subroutine preprocess_input
1944 
1945  !> @brief Calculate CONDSAT array entries for the given node
1946  !!
1947  !! Calculate saturated conductances for all connections of the given node,
1948  !! or optionally for the upper portion of the matrix only.
1949  !<
1950  subroutine calc_condsat(this, node, upperOnly)
1951  ! -- dummy variables
1952  class(gwfnpftype) :: this
1953  integer(I4B), intent(in) :: node
1954  logical, intent(in) :: upperOnly
1955  ! -- local variables
1956  integer(I4B) :: ii, m, n, ihc, jj
1957  real(DP) :: topm, topn, topnode, botm, botn, botnode, satm, satn, satnode
1958  real(DP) :: hyn, hym, hn, hm, fawidth, csat
1959  !
1960  satnode = this%calc_initial_sat(node)
1961  !
1962  topnode = this%dis%top(node)
1963  botnode = this%dis%bot(node)
1964  !
1965  ! -- Go through the connecting cells
1966  do ii = this%dis%con%ia(node) + 1, this%dis%con%ia(node + 1) - 1
1967  !
1968  ! -- Set the m cell number and cycle if lower triangle connection and
1969  ! -- we're not updating both upper and lower matrix parts for this node
1970  m = this%dis%con%ja(ii)
1971  jj = this%dis%con%jas(ii)
1972  if (m < node) then
1973  if (upperonly) cycle
1974  ! m => node, n => neighbour
1975  n = m
1976  m = node
1977  topm = topnode
1978  botm = botnode
1979  satm = satnode
1980  topn = this%dis%top(n)
1981  botn = this%dis%bot(n)
1982  satn = this%calc_initial_sat(n)
1983  else
1984  ! n => node, m => neighbour
1985  n = node
1986  topn = topnode
1987  botn = botnode
1988  satn = satnode
1989  topm = this%dis%top(m)
1990  botm = this%dis%bot(m)
1991  satm = this%calc_initial_sat(m)
1992  end if
1993  !
1994  ihc = this%dis%con%ihc(jj)
1995  hyn = this%hy_eff(n, m, ihc, ipos=ii)
1996  hym = this%hy_eff(m, n, ihc, ipos=ii)
1997  if (this%ithickstartflag(n) == 0) then
1998  hn = topn
1999  else
2000  hn = this%ic%strt(n)
2001  end if
2002  if (this%ithickstartflag(m) == 0) then
2003  hm = topm
2004  else
2005  hm = this%ic%strt(m)
2006  end if
2007  !
2008  ! -- Calculate conductance depending on whether connection is
2009  ! vertical (0), horizontal (1), or staggered horizontal (2)
2010  if (ihc == c3d_vertical) then
2011  !
2012  ! -- Vertical conductance for fully saturated conditions
2013  csat = vcond(1, 1, 1, 1, 0, 1, 1, done, &
2014  botn, botm, &
2015  hyn, hym, &
2016  satn, satm, &
2017  topn, topm, &
2018  botn, botm, &
2019  this%dis%con%hwva(jj))
2020  else
2021  !
2022  ! -- Horizontal conductance for fully saturated conditions
2023  fawidth = this%dis%con%hwva(jj)
2024  csat = hcond(1, 1, 1, 1, 0, &
2025  ihc, &
2026  this%icellavg, &
2027  done, &
2028  hn, hm, satn, satm, hyn, hym, &
2029  topn, topm, &
2030  botn, botm, &
2031  this%dis%con%cl1(jj), &
2032  this%dis%con%cl2(jj), &
2033  fawidth)
2034  end if
2035  this%condsat(jj) = csat
2036  end do
2037  end subroutine calc_condsat
2038 
2039  !> @brief Calculate initial saturation for the given node
2040  !!
2041  !! Calculate saturation as a fraction of thickness for the given node, used
2042  !! for saturated conductance calculations: full thickness by default (1.0) or
2043  !! saturation based on initial conditions if the THICKSTRT option is used.
2044  !!
2045  !<
2046  function calc_initial_sat(this, n) result(satn)
2047  ! -- dummy variables
2048  class(gwfnpftype) :: this
2049  integer(I4B), intent(in) :: n
2050  ! -- Return
2051  real(dp) :: satn
2052  !
2053  satn = done
2054  if (this%ibound(n) /= 0 .and. this%ithickstartflag(n) /= 0) then
2055  call this%thksat(n, this%ic%strt(n), satn)
2056  end if
2057  end function calc_initial_sat
2058 
2059  !> @brief Perform wetting and drying
2060  !<
2061  subroutine sgwf_npf_wetdry(this, kiter, hnew)
2062  ! -- modules
2063  use tdismodule, only: kstp, kper
2065  use constantsmodule, only: linelength
2066  ! -- dummy
2067  class(gwfnpftype) :: this
2068  integer(I4B), intent(in) :: kiter
2069  real(DP), intent(inout), dimension(:) :: hnew
2070  ! -- local
2071  integer(I4B) :: n, m, ii, ihc
2072  real(DP) :: ttop, bbot, thick
2073  integer(I4B) :: ncnvrt, ihdcnv
2074  character(len=30), dimension(5) :: nodcnvrt
2075  character(len=30) :: nodestr
2076  character(len=3), dimension(5) :: acnvrt
2077  character(len=LINELENGTH) :: errmsg
2078  integer(I4B) :: irewet
2079  ! -- formats
2080  character(len=*), parameter :: fmtnct = &
2081  "(1X,/1X,'Negative cell thickness at (layer,row,col)', &
2082  &I4,',',I5,',',I5)"
2083  character(len=*), parameter :: fmttopbot = &
2084  &"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
2085  character(len=*), parameter :: fmttopbotthk = &
2086  &"(1X,'Top elevation, bottom elevation, thickness:',1P,3G13.5)"
2087  character(len=*), parameter :: fmtdrychd = &
2088  &"(1X,/1X,'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED')"
2089  character(len=*), parameter :: fmtni = &
2090  &"(1X,'CELLID=',a,' ITERATION=',I0,' TIME STEP=',I0,' STRESS PERIOD=',I0)"
2091  !
2092  ! -- Initialize
2093  ncnvrt = 0
2094  ihdcnv = 0
2095  !
2096  ! -- Convert dry cells to wet
2097  do n = 1, this%dis%nodes
2098  do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2099  m = this%dis%con%ja(ii)
2100  ihc = this%dis%con%ihc(this%dis%con%jas(ii))
2101  call this%rewet_check(kiter, n, hnew(m), this%ibound(m), ihc, hnew, &
2102  irewet)
2103  if (irewet == 1) then
2104  call this%wdmsg(2, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2105  end if
2106  end do
2107  end do
2108  !
2109  ! -- Perform drying
2110  do n = 1, this%dis%nodes
2111  !
2112  ! -- cycle if inactive or confined
2113  if (this%ibound(n) == 0) cycle
2114  if (this%icelltype(n) == 0) cycle
2115  !
2116  ! -- check for negative cell thickness
2117  bbot = this%dis%bot(n)
2118  ttop = this%dis%top(n)
2119  if (bbot > ttop) then
2120  write (errmsg, fmtnct) n
2121  call store_error(errmsg)
2122  write (errmsg, fmttopbot) ttop, bbot
2123  call store_error(errmsg)
2124  call store_error_filename(this%input_fname)
2125  end if
2126  !
2127  ! -- Calculate saturated thickness
2128  if (this%icelltype(n) /= 0) then
2129  if (hnew(n) < ttop) ttop = hnew(n)
2130  end if
2131  thick = ttop - bbot
2132  !
2133  ! -- If thick<0 print message, set hnew, and ibound
2134  if (thick <= dzero) then
2135  call this%wdmsg(1, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2136  hnew(n) = this%hdry
2137  if (this%ibound(n) < 0) then
2138  errmsg = 'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED'
2139  call store_error(errmsg)
2140  write (errmsg, fmttopbotthk) ttop, bbot, thick
2141  call store_error(errmsg)
2142  call this%dis%noder_to_string(n, nodestr)
2143  write (errmsg, fmtni) trim(adjustl(nodestr)), kiter, kstp, kper
2144  call store_error(errmsg)
2145  call store_error_filename(this%input_fname)
2146  end if
2147  this%ibound(n) = 0
2148  end if
2149  end do
2150  !
2151  ! -- Print remaining cell conversions
2152  call this%wdmsg(0, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2153  !
2154  ! -- Change ibound from 30000 to 1
2155  do n = 1, this%dis%nodes
2156  if (this%ibound(n) == 30000) this%ibound(n) = 1
2157  end do
2158  end subroutine sgwf_npf_wetdry
2159 
2160  !> @brief Determine if a cell should rewet
2161  !!
2162  !! This method can be called from any external object that has a head that
2163  !! can be used to rewet the GWF cell node. The ihc value is used to
2164  !! determine if it is a vertical or horizontal connection, which can operate
2165  !! differently depending on user settings.
2166  !<
2167  subroutine rewet_check(this, kiter, node, hm, ibdm, ihc, hnew, irewet)
2168  ! -- dummy
2169  class(gwfnpftype) :: this
2170  integer(I4B), intent(in) :: kiter
2171  integer(I4B), intent(in) :: node
2172  real(DP), intent(in) :: hm
2173  integer(I4B), intent(in) :: ibdm
2174  integer(I4B), intent(in) :: ihc
2175  real(DP), intent(inout), dimension(:) :: hnew
2176  integer(I4B), intent(out) :: irewet
2177  ! -- local
2178  integer(I4B) :: itflg
2179  real(DP) :: wd, awd, turnon, bbot
2180  !
2181  irewet = 0
2182  !
2183  ! -- Convert a dry cell to wet if it meets the criteria
2184  if (this%irewet > 0) then
2185  itflg = mod(kiter, this%iwetit)
2186  if (itflg == 0) then
2187  if (this%ibound(node) == 0 .and. this%wetdry(node) /= dzero) then
2188  !
2189  ! -- Calculate wetting elevation
2190  bbot = this%dis%bot(node)
2191  wd = this%wetdry(node)
2192  awd = wd
2193  if (wd < 0) awd = -wd
2194  turnon = bbot + awd
2195  !
2196  ! -- Check head in adjacent cells to see if wetting elevation has
2197  ! been reached
2198  if (ihc == c3d_vertical) then
2199  !
2200  ! -- check cell below
2201  if (ibdm > 0 .and. hm >= turnon) irewet = 1
2202  else
2203  if (wd > dzero) then
2204  !
2205  ! -- check horizontally adjacent cells
2206  if (ibdm > 0 .and. hm >= turnon) irewet = 1
2207  end if
2208  end if
2209  !
2210  if (irewet == 1) then
2211  ! -- rewet cell; use equation 3a if ihdwet=0; use equation 3b if
2212  ! ihdwet is not 0.
2213  if (this%ihdwet == 0) then
2214  hnew(node) = bbot + this%wetfct * (hm - bbot)
2215  else
2216  hnew(node) = bbot + this%wetfct * awd !(hm - bbot)
2217  end if
2218  this%ibound(node) = 30000
2219  end if
2220  end if
2221  end if
2222  end if
2223  end subroutine rewet_check
2224 
2225  !> @brief Print wet/dry message
2226  !<
2227  subroutine sgwf_npf_wdmsg(this, icode, ncnvrt, nodcnvrt, acnvrt, ihdcnv, &
2228  kiter, n)
2229  ! -- modules
2230  use tdismodule, only: kstp, kper
2231  ! -- dummy
2232  class(gwfnpftype) :: this
2233  integer(I4B), intent(in) :: icode
2234  integer(I4B), intent(inout) :: ncnvrt
2235  character(len=30), dimension(5), intent(inout) :: nodcnvrt
2236  character(len=3), dimension(5), intent(inout) :: acnvrt
2237  integer(I4B), intent(inout) :: ihdcnv
2238  integer(I4B), intent(in) :: kiter
2239  integer(I4B), intent(in) :: n
2240  ! -- local
2241  integer(I4B) :: l
2242  ! -- formats
2243  character(len=*), parameter :: fmtcnvtn = &
2244  "(1X,/1X,'CELL CONVERSIONS FOR ITER.=',I0, &
2245  &' STEP=',I0,' PERIOD=',I0,' (NODE or LRC)')"
2246  character(len=*), parameter :: fmtnode = "(1X,3X,5(A4, A20))"
2247  !
2248  ! -- Keep track of cell conversions
2249  if (icode > 0) then
2250  ncnvrt = ncnvrt + 1
2251  call this%dis%noder_to_string(n, nodcnvrt(ncnvrt))
2252  if (icode == 1) then
2253  acnvrt(ncnvrt) = 'DRY'
2254  else
2255  acnvrt(ncnvrt) = 'WET'
2256  end if
2257  end if
2258  !
2259  ! -- Print a line if 5 conversions have occurred or if icode indicates that a
2260  ! partial line should be printed
2261  if (ncnvrt == 5 .or. (icode == 0 .and. ncnvrt > 0)) then
2262  if (ihdcnv == 0) write (this%iout, fmtcnvtn) kiter, kstp, kper
2263  ihdcnv = 1
2264  write (this%iout, fmtnode) &
2265  (acnvrt(l), trim(adjustl(nodcnvrt(l))), l=1, ncnvrt)
2266  ncnvrt = 0
2267  end if
2268  end subroutine sgwf_npf_wdmsg
2269 
2270  !> @brief Calculate the effective hydraulic conductivity for the n-m connection
2271  !!
2272  !! n is primary node node number
2273  !! m is connected node (not used if vg is provided)
2274  !! ihc is horizontal indicator (0 vertical, 1 horizontal, 2 vertically
2275  !! staggered)
2276  !! ipos_opt is position of connection in ja array
2277  !! vg is the global unit vector that expresses the direction from which to
2278  !! calculate an effective hydraulic conductivity.
2279  !<
2280  function hy_eff(this, n, m, ihc, ipos, vg) result(hy)
2281  ! -- return
2282  real(dp) :: hy
2283  ! -- dummy
2284  class(gwfnpftype) :: this
2285  integer(I4B), intent(in) :: n
2286  integer(I4B), intent(in) :: m
2287  integer(I4B), intent(in) :: ihc
2288  integer(I4B), intent(in), optional :: ipos
2289  real(dp), dimension(3), intent(in), optional :: vg
2290  ! -- local
2291  integer(I4B) :: iipos
2292  real(dp) :: hy11, hy22, hy33
2293  real(dp) :: ang1, ang2, ang3
2294  real(dp) :: vg1, vg2, vg3
2295  !
2296  ! -- Initialize
2297  iipos = 0
2298  if (present(ipos)) iipos = ipos
2299  hy11 = this%k11(n)
2300  hy22 = this%k11(n)
2301  hy33 = this%k11(n)
2302  hy22 = this%k22(n)
2303  hy33 = this%k33(n)
2304  !
2305  ! -- Calculate effective K based on whether connection is vertical
2306  ! or horizontal
2307  if (ihc == c3d_vertical) then
2308  !
2309  ! -- Handle rotated anisotropy case that would affect the effective
2310  ! vertical hydraulic conductivity
2311  hy = hy33
2312  if (this%iangle2 > 0) then
2313  if (present(vg)) then
2314  vg1 = vg(1)
2315  vg2 = vg(2)
2316  vg3 = vg(3)
2317  else
2318  call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2319  end if
2320  ang1 = this%angle1(n)
2321  ang2 = this%angle2(n)
2322  ang3 = dzero
2323  if (this%iangle3 > 0) ang3 = this%angle3(n)
2324  hy = hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2325  this%iavgkeff)
2326  end if
2327  !
2328  else
2329  !
2330  ! -- Handle horizontal case
2331  hy = hy11
2332  if (this%ik22 > 0) then
2333  if (present(vg)) then
2334  vg1 = vg(1)
2335  vg2 = vg(2)
2336  vg3 = vg(3)
2337  else
2338  call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2339  end if
2340  ang1 = dzero
2341  ang2 = dzero
2342  ang3 = dzero
2343  if (this%iangle1 > 0) then
2344  ang1 = this%angle1(n)
2345  if (this%iangle2 > 0) then
2346  ang2 = this%angle2(n)
2347  if (this%iangle3 > 0) ang3 = this%angle3(n)
2348  end if
2349  end if
2350  hy = hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2351  this%iavgkeff)
2352  end if
2353  !
2354  end if
2355  end function hy_eff
2356 
2357  !> @brief Calculate the 3 components of specific discharge at the cell center
2358  !<
2359  subroutine calc_spdis(this, flowja)
2360  ! -- modules
2361  use simmodule, only: store_error
2362  ! -- dummy
2363  class(gwfnpftype) :: this
2364  real(DP), intent(in), dimension(:) :: flowja
2365  ! -- local
2366  integer(I4B) :: n
2367  integer(I4B) :: m
2368  integer(I4B) :: ipos
2369  integer(I4B) :: iedge
2370  integer(I4B) :: isympos
2371  integer(I4B) :: ihc
2372  integer(I4B) :: ic
2373  integer(I4B) :: iz
2374  integer(I4B) :: nc
2375  integer(I4B) :: ncz
2376  real(DP) :: qz
2377  real(DP) :: vx
2378  real(DP) :: vy
2379  real(DP) :: vz
2380  real(DP) :: xn
2381  real(DP) :: yn
2382  real(DP) :: zn
2383  real(DP) :: xc
2384  real(DP) :: yc
2385  real(DP) :: zc
2386  real(DP) :: cl1
2387  real(DP) :: cl2
2388  real(DP) :: dltot
2389  real(DP) :: ooclsum
2390  real(DP) :: dsumx
2391  real(DP) :: dsumy
2392  real(DP) :: dsumz
2393  real(DP) :: denom
2394  real(DP) :: area
2395  real(DP) :: dz
2396  real(DP) :: axy
2397  real(DP) :: ayx
2398  logical :: nozee = .true.
2399  type(spdisworkarraytype), pointer :: swa => null() !< pointer to spdis work arrays structure
2400  !
2401  ! -- Ensure dis has necessary information
2402  if (this%icalcspdis /= 0 .and. this%dis%con%ianglex == 0) then
2403  call store_error('Error. ANGLDEGX not provided in '// &
2404  'discretization file. ANGLDEGX required for '// &
2405  'calculation of specific discharge.', terminate=.true.)
2406  end if
2407 
2408  swa => this%spdis_wa
2409  if (.not. swa%is_created()) then
2410  ! prepare work arrays
2411  call this%spdis_wa%create(this%calc_max_conns())
2412 
2413  ! prepare lookup table
2414  if (this%nedges > 0) call this%prepare_edge_lookup()
2415  end if
2416  !
2417  ! -- Go through each cell and calculate specific discharge
2418  do n = 1, this%dis%nodes
2419  !
2420  ! -- first calculate geometric properties for x and y directions and
2421  ! the specific discharge at a face (vi)
2422  ic = 0
2423  iz = 0
2424 
2425  ! reset work arrays
2426  call swa%reset()
2427 
2428  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2429  m = this%dis%con%ja(ipos)
2430  isympos = this%dis%con%jas(ipos)
2431  ihc = this%dis%con%ihc(isympos)
2432  area = this%dis%con%hwva(isympos)
2433  if (ihc == c3d_vertical) then
2434  !
2435  ! -- vertical connection
2436  iz = iz + 1
2437  !call this%dis%connection_normal(n, m, ihc, xn, yn, zn, ipos)
2438  call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2439  ihc, xc, yc, zc, dltot)
2440  cl1 = this%dis%con%cl1(isympos)
2441  cl2 = this%dis%con%cl2(isympos)
2442  if (m < n) then
2443  cl1 = this%dis%con%cl2(isympos)
2444  cl2 = this%dis%con%cl1(isympos)
2445  end if
2446  ooclsum = done / (cl1 + cl2)
2447  swa%diz(iz) = dltot * cl1 * ooclsum
2448  qz = flowja(ipos)
2449  if (n > m) qz = -qz
2450  swa%viz(iz) = qz / area
2451  else
2452  !
2453  ! -- horizontal connection
2454  ic = ic + 1
2455  dz = thksatnm(this%ibound(n), this%ibound(m), &
2456  this%icelltype(n), this%icelltype(m), &
2457  this%inewton, ihc, &
2458  this%hnew(n), this%hnew(m), this%sat(n), this%sat(m), &
2459  this%dis%top(n), this%dis%top(m), this%dis%bot(n), &
2460  this%dis%bot(m))
2461  area = area * dz
2462  call this%dis%connection_normal(n, m, ihc, xn, yn, zn, ipos)
2463  call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2464  ihc, xc, yc, zc, dltot)
2465  cl1 = this%dis%con%cl1(isympos)
2466  cl2 = this%dis%con%cl2(isympos)
2467  if (m < n) then
2468  cl1 = this%dis%con%cl2(isympos)
2469  cl2 = this%dis%con%cl1(isympos)
2470  end if
2471  ooclsum = done / (cl1 + cl2)
2472  swa%nix(ic) = -xn
2473  swa%niy(ic) = -yn
2474  swa%di(ic) = dltot * cl1 * ooclsum
2475  if (area > dzero) then
2476  swa%vi(ic) = flowja(ipos) / area
2477  else
2478  swa%vi(ic) = dzero
2479  end if
2480  end if
2481  end do
2482 
2483  ! add contribution from edge flows (i.e. from exchanges)
2484  if (this%nedges > 0) then
2485  do ipos = this%iedge_ptr(n), this%iedge_ptr(n + 1) - 1
2486  iedge = this%edge_idxs(ipos)
2487 
2488  ! propsedge: (Q, area, nx, ny, distance)
2489  ihc = this%ihcedge(iedge)
2490  area = this%propsedge(2, iedge)
2491  if (ihc == c3d_vertical) then
2492  iz = iz + 1
2493  swa%viz(iz) = this%propsedge(1, iedge) / area
2494  swa%diz(iz) = this%propsedge(5, iedge)
2495  else
2496  ic = ic + 1
2497  swa%nix(ic) = -this%propsedge(3, iedge)
2498  swa%niy(ic) = -this%propsedge(4, iedge)
2499  swa%di(ic) = this%propsedge(5, iedge)
2500  if (area > dzero) then
2501  swa%vi(ic) = this%propsedge(1, iedge) / area
2502  else
2503  swa%vi(ic) = dzero
2504  end if
2505  end if
2506  end do
2507  end if
2508  !
2509  ! -- Assign number of vertical and horizontal connections
2510  ncz = iz
2511  nc = ic
2512  !
2513  ! -- calculate z weight (wiz) and z velocity
2514  if (ncz == 1) then
2515  swa%wiz(1) = done
2516  else
2517  dsumz = dzero
2518  do iz = 1, ncz
2519  dsumz = dsumz + swa%diz(iz)
2520  end do
2521  denom = (ncz - done)
2522  if (denom < dzero) denom = dzero
2523  dsumz = dsumz + dem10 * dsumz
2524  do iz = 1, ncz
2525  if (dsumz > dzero) swa%wiz(iz) = done - swa%diz(iz) / dsumz
2526  if (denom > 0) then
2527  swa%wiz(iz) = swa%wiz(iz) / denom
2528  else
2529  swa%wiz(iz) = dzero
2530  end if
2531  end do
2532  end if
2533  vz = dzero
2534  do iz = 1, ncz
2535  vz = vz + swa%wiz(iz) * swa%viz(iz)
2536  end do
2537  !
2538  ! -- distance-based weighting
2539  nc = ic
2540  dsumx = dzero
2541  dsumy = dzero
2542  dsumz = dzero
2543  do ic = 1, nc
2544  swa%wix(ic) = swa%di(ic) * abs(swa%nix(ic))
2545  swa%wiy(ic) = swa%di(ic) * abs(swa%niy(ic))
2546  dsumx = dsumx + swa%wix(ic)
2547  dsumy = dsumy + swa%wiy(ic)
2548  end do
2549  !
2550  ! -- Finish computing omega weights. Add a tiny bit
2551  ! to dsum so that the normalized omega weight later
2552  ! evaluates to (essentially) 1 in the case of a single
2553  ! relevant connection, avoiding 0/0.
2554  dsumx = dsumx + dem10 * dsumx
2555  dsumy = dsumy + dem10 * dsumy
2556  do ic = 1, nc
2557  swa%wix(ic) = (dsumx - swa%wix(ic)) * abs(swa%nix(ic))
2558  swa%wiy(ic) = (dsumy - swa%wiy(ic)) * abs(swa%niy(ic))
2559  end do
2560  !
2561  ! -- compute B weights
2562  dsumx = dzero
2563  dsumy = dzero
2564  do ic = 1, nc
2565  swa%bix(ic) = swa%wix(ic) * sign(done, swa%nix(ic))
2566  swa%biy(ic) = swa%wiy(ic) * sign(done, swa%niy(ic))
2567  dsumx = dsumx + swa%wix(ic) * abs(swa%nix(ic))
2568  dsumy = dsumy + swa%wiy(ic) * abs(swa%niy(ic))
2569  end do
2570  if (dsumx > dzero) dsumx = done / dsumx
2571  if (dsumy > dzero) dsumy = done / dsumy
2572  axy = dzero
2573  ayx = dzero
2574  do ic = 1, nc
2575  swa%bix(ic) = swa%bix(ic) * dsumx
2576  swa%biy(ic) = swa%biy(ic) * dsumy
2577  axy = axy + swa%bix(ic) * swa%niy(ic)
2578  ayx = ayx + swa%biy(ic) * swa%nix(ic)
2579  end do
2580  !
2581  ! -- Calculate specific discharge. The divide by zero checking below
2582  ! is problematic for cells with only one flow, such as can happen
2583  ! with triangular cells in corners. In this case, the resulting
2584  ! cell velocity will be calculated as zero. The method should be
2585  ! improved so that edge flows of zero are included in these
2586  ! calculations. But this needs to be done with consideration for LGR
2587  ! cases in which flows are submitted from an exchange.
2588  vx = dzero
2589  vy = dzero
2590  do ic = 1, nc
2591  vx = vx + (swa%bix(ic) - axy * swa%biy(ic)) * swa%vi(ic)
2592  vy = vy + (swa%biy(ic) - ayx * swa%bix(ic)) * swa%vi(ic)
2593  end do
2594  denom = done - axy * ayx
2595  if (denom /= dzero) then
2596  vx = vx / denom
2597  vy = vy / denom
2598  end if
2599  !
2600  this%spdis(1, n) = vx
2601  this%spdis(2, n) = vy
2602  this%spdis(3, n) = vz
2603  !
2604  end do
2605 
2606  end subroutine calc_spdis
2607 
2608  !> @brief Save specific discharge in binary format to ibinun
2609  !<
2610  subroutine sav_spdis(this, ibinun)
2611  ! -- dummy
2612  class(gwfnpftype) :: this
2613  integer(I4B), intent(in) :: ibinun
2614  ! -- local
2615  character(len=16) :: text
2616  character(len=16), dimension(3) :: auxtxt
2617  integer(I4B) :: n
2618  integer(I4B) :: naux
2619  !
2620  ! -- Write the header
2621  text = ' DATA-SPDIS'
2622  naux = 3
2623  auxtxt(:) = [' qx', ' qy', ' qz']
2624  call this%dis%record_srcdst_list_header(text, this%name_model, &
2625  this%packName, this%name_model, &
2626  this%packName, naux, auxtxt, ibinun, &
2627  this%dis%nodes, this%iout)
2628  !
2629  ! -- Write a zero for Q, and then write qx, qy, qz as aux variables
2630  do n = 1, this%dis%nodes
2631  call this%dis%record_mf6_list_entry(ibinun, n, n, dzero, naux, &
2632  this%spdis(:, n))
2633  end do
2634  end subroutine sav_spdis
2635 
2636  !> @brief Save saturation in binary format to ibinun
2637  !<
2638  subroutine sav_sat(this, ibinun)
2639  ! -- dummy
2640  class(gwfnpftype) :: this
2641  integer(I4B), intent(in) :: ibinun
2642  ! -- local
2643  character(len=16) :: text
2644  character(len=16), dimension(1) :: auxtxt
2645  real(DP), dimension(1) :: a
2646  integer(I4B) :: n
2647  integer(I4B) :: naux
2648  !
2649  ! -- Write the header
2650  text = ' DATA-SAT'
2651  naux = 1
2652  auxtxt(:) = [' sat']
2653  call this%dis%record_srcdst_list_header(text, this%name_model, &
2654  this%packName, this%name_model, &
2655  this%packName, naux, auxtxt, ibinun, &
2656  this%dis%nodes, this%iout)
2657  !
2658  ! -- Write a zero for Q, and then write saturation as an aux variables
2659  do n = 1, this%dis%nodes
2660  a(1) = this%sat(n)
2661  call this%dis%record_mf6_list_entry(ibinun, n, n, dzero, naux, a)
2662  end do
2663  end subroutine sav_sat
2664 
2665  !> @brief Reserve space for nedges cells that have an edge on them.
2666  !!
2667  !! This must be called before the npf%allocate_arrays routine, which is
2668  !! called from npf%ar.
2669  !<
2670  subroutine increase_edge_count(this, nedges)
2671  ! -- dummy
2672  class(gwfnpftype) :: this
2673  integer(I4B), intent(in) :: nedges
2674  !
2675  this%nedges = this%nedges + nedges
2676  end subroutine increase_edge_count
2677 
2678  !> @brief Calculate the maximum number of connections for any cell
2679  !<
2680  function calc_max_conns(this) result(max_conns)
2681  class(gwfnpftype) :: this
2682  integer(I4B) :: max_conns
2683  ! local
2684  integer(I4B) :: n, m, ic
2685 
2686  max_conns = 0
2687  do n = 1, this%dis%nodes
2688 
2689  ! Count internal model connections
2690  ic = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
2691 
2692  ! Add edge connections
2693  do m = 1, this%nedges
2694  if (this%nodedge(m) == n) then
2695  ic = ic + 1
2696  end if
2697  end do
2698 
2699  ! Set max number of connections for any cell
2700  if (ic > max_conns) max_conns = ic
2701  end do
2702 
2703  end function calc_max_conns
2704 
2705  !> @brief Provide the npf package with edge properties
2706  !<
2707  subroutine set_edge_properties(this, nodedge, ihcedge, q, area, nx, ny, &
2708  distance)
2709  ! -- dummy
2710  class(gwfnpftype) :: this
2711  integer(I4B), intent(in) :: nodedge
2712  integer(I4B), intent(in) :: ihcedge
2713  real(DP), intent(in) :: q
2714  real(DP), intent(in) :: area
2715  real(DP), intent(in) :: nx
2716  real(DP), intent(in) :: ny
2717  real(DP), intent(in) :: distance
2718  ! -- local
2719  integer(I4B) :: lastedge
2720  !
2721  this%lastedge = this%lastedge + 1
2722  lastedge = this%lastedge
2723  this%nodedge(lastedge) = nodedge
2724  this%ihcedge(lastedge) = ihcedge
2725  this%propsedge(1, lastedge) = q
2726  this%propsedge(2, lastedge) = area
2727  this%propsedge(3, lastedge) = nx
2728  this%propsedge(4, lastedge) = ny
2729  this%propsedge(5, lastedge) = distance
2730  !
2731  ! -- If this is the last edge, then the next call must be starting a new
2732  ! edge properties assignment loop, so need to reset lastedge to 0
2733  if (this%lastedge == this%nedges) this%lastedge = 0
2734  end subroutine set_edge_properties
2735 
2736  subroutine prepare_edge_lookup(this)
2737  class(gwfnpftype) :: this
2738  ! local
2739  integer(I4B) :: i, inode, iedge
2740  integer(I4B) :: n, start, end
2741  integer(I4B) :: prev_cnt, strt_idx, ipos
2742 
2743  do i = 1, size(this%iedge_ptr)
2744  this%iedge_ptr(i) = 0
2745  end do
2746  do i = 1, size(this%edge_idxs)
2747  this%edge_idxs(i) = 0
2748  end do
2749 
2750  ! count
2751  do iedge = 1, this%nedges
2752  n = this%nodedge(iedge)
2753  this%iedge_ptr(n) = this%iedge_ptr(n) + 1
2754  end do
2755 
2756  ! determine start indexes
2757  prev_cnt = this%iedge_ptr(1)
2758  this%iedge_ptr(1) = 1
2759  do inode = 2, this%dis%nodes + 1
2760  strt_idx = this%iedge_ptr(inode - 1) + prev_cnt
2761  prev_cnt = this%iedge_ptr(inode)
2762  this%iedge_ptr(inode) = strt_idx
2763  end do
2764 
2765  ! loop over edges to fill lookup table
2766  do iedge = 1, this%nedges
2767  n = this%nodedge(iedge)
2768  start = this%iedge_ptr(n)
2769  end = this%iedge_ptr(n + 1) - 1
2770  do ipos = start, end
2771  if (this%edge_idxs(ipos) > 0) cycle ! go to next
2772  this%edge_idxs(ipos) = iedge
2773  exit
2774  end do
2775  end do
2776 
2777  end subroutine prepare_edge_lookup
2778 
2779  !> Calculate saturated thickness between cell n and m
2780  !<
2781  function calcsatthickness(this, n, m, ihc) result(satThickness)
2782  ! -- dummy
2783  class(gwfnpftype) :: this !< this NPF instance
2784  integer(I4B) :: n !< node n
2785  integer(I4B) :: m !< node m
2786  integer(I4B) :: ihc !< 1 = horizontal connection, 0 for vertical
2787  ! -- return
2788  real(dp) :: satthickness !< saturated thickness
2789  !
2790  satthickness = thksatnm(this%ibound(n), &
2791  this%ibound(m), &
2792  this%icelltype(n), &
2793  this%icelltype(m), &
2794  this%inewton, &
2795  ihc, &
2796  this%hnew(n), &
2797  this%hnew(m), &
2798  this%sat(n), &
2799  this%sat(m), &
2800  this%dis%top(n), &
2801  this%dis%top(m), &
2802  this%dis%bot(n), &
2803  this%dis%bot(m))
2804  end function calcsatthickness
2805 
2806 end module gwfnpfmodule
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
@ c3d_vertical
vertical connection
Definition: Constants.f90:222
real(dp), parameter dhdry
real dry cell constant
Definition: Constants.f90:94
real(dp), parameter dp9
real constant 9/10
Definition: Constants.f90:72
real(dp), parameter dem10
real constant 1e-10
Definition: Constants.f90:113
real(dp), parameter dem7
real constant 1e-7
Definition: Constants.f90:110
real(dp), parameter dem8
real constant 1e-8
Definition: Constants.f90:111
integer(i4b), parameter lenbigline
maximum length of a big line
Definition: Constants.f90:15
real(dp), parameter dhnoflo
real no flow constant
Definition: Constants.f90:93
integer(i4b), parameter lenvarname
maximum length of a variable name
Definition: Constants.f90:17
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
real(dp), parameter dpio180
real constant
Definition: Constants.f90:130
real(dp), parameter dem6
real constant 1e-6
Definition: Constants.f90:109
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter dem9
real constant 1e-9
Definition: Constants.f90:112
real(dp), parameter dem2
real constant 1e-2
Definition: Constants.f90:105
real(dp), parameter dtwo
real constant 2
Definition: Constants.f90:79
integer(i4b), parameter lenmempath
maximum length of the memory path
Definition: Constants.f90:27
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
This module contains stateless conductance functions.
real(dp) function, public thksatnm(ibdn, ibdm, ictn, ictm, iupstream, ihc, hn, hm, satn, satm, topn, topm, botn, botm)
Calculate wetted cell thickness at interface between two cells.
real(dp) function, public condmean(k1, k2, thick1, thick2, cl1, cl2, width, iavgmeth)
Calculate the conductance between two cells.
real(dp) function, public hcond(ibdn, ibdm, ictn, ictm, iupstream, ihc, icellavg, condsat, hn, hm, satn, satm, hkn, hkm, topn, topm, botn, botm, cln, clm, fawidth)
Horizontal conductance between two cells.
@, public ccond_hmean
Harmonic mean.
real(dp) function, public vcond(ibdn, ibdm, ictn, ictm, inewton, ivarcv, idewatcv, condsat, hn, hm, vkn, vkm, satn, satm, topn, topm, botn, botm, flowarea)
Vertical conductance between two cells.
subroutine set_options(this, options)
Set options in the NPF object.
Definition: gwf-npf.f90:1369
subroutine npf_save_model_flows(this, flowja, icbcfl, icbcun)
Record flowja and calculate specific discharge if requested.
Definition: gwf-npf.f90:870
subroutine source_options(this)
Update simulation options from input mempath.
Definition: gwf-npf.f90:1293
real(dp) function calc_initial_sat(this, n)
Calculate initial saturation for the given node.
Definition: gwf-npf.f90:2047
subroutine calc_condsat(this, node, upperOnly)
Calculate CONDSAT array entries for the given node.
Definition: gwf-npf.f90:1951
subroutine rewet_check(this, kiter, node, hm, ibdm, ihc, hnew, irewet)
Determine if a cell should rewet.
Definition: gwf-npf.f90:2168
integer(i4b) function calc_max_conns(this)
Calculate the maximum number of connections for any cell.
Definition: gwf-npf.f90:2681
subroutine npf_mc(this, moffset, matrix_sln)
Map connections and construct iax, jax, and idxglox.
Definition: gwf-npf.f90:272
subroutine npf_fc(this, kiter, matrix_sln, idxglo, rhs, hnew)
Formulate coefficients.
Definition: gwf-npf.f90:475
subroutine npf_ac(this, moffset, sparse)
Add connections for extended neighbors to the sparse matrix.
Definition: gwf-npf.f90:258
subroutine sgwf_npf_wetdry(this, kiter, hnew)
Perform wetting and drying.
Definition: gwf-npf.f90:2062
subroutine source_griddata(this)
Update simulation griddata from input mempath.
Definition: gwf-npf.f90:1499
subroutine npf_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax)
Under-relaxation.
Definition: gwf-npf.f90:706
subroutine sav_spdis(this, ibinun)
Save specific discharge in binary format to ibinun.
Definition: gwf-npf.f90:2611
subroutine sgwf_npf_thksat(this, n, hn, thksat)
Fractional cell saturation.
Definition: gwf-npf.f90:776
subroutine preprocess_input(this)
preprocess the NPF input data
Definition: gwf-npf.f90:1749
subroutine set_edge_properties(this, nodedge, ihcedge, q, area, nx, ny, distance)
Provide the npf package with edge properties.
Definition: gwf-npf.f90:2709
subroutine prepare_edge_lookup(this)
Definition: gwf-npf.f90:2737
subroutine npf_da(this)
Deallocate variables.
Definition: gwf-npf.f90:946
subroutine log_griddata(this, found)
Write dimensions to list file.
Definition: gwf-npf.f90:1449
real(dp) function hy_eff(this, n, m, ihc, ipos, vg)
Calculate the effective hydraulic conductivity for the n-m connection.
Definition: gwf-npf.f90:2281
subroutine calc_spdis(this, flowja)
Calculate the 3 components of specific discharge at the cell center.
Definition: gwf-npf.f90:2360
subroutine, public npf_cr(npfobj, name_model, input_mempath, inunit, iout)
Create a new NPF object. Pass a inunit value of 0 if npf data will initialized from memory.
Definition: gwf-npf.f90:158
subroutine increase_edge_count(this, nedges)
Reserve space for nedges cells that have an edge on them.
Definition: gwf-npf.f90:2671
subroutine npf_fn(this, kiter, matrix_sln, idxglo, rhs, hnew)
Fill newton terms.
Definition: gwf-npf.f90:579
subroutine log_options(this, found)
Log npf options sourced from the input mempath.
Definition: gwf-npf.f90:1224
subroutine npf_print_model_flows(this, ibudfl, flowja)
Print budget.
Definition: gwf-npf.f90:907
subroutine allocate_arrays(this, ncells, njas)
Allocate npf arrays.
Definition: gwf-npf.f90:1163
subroutine sav_sat(this, ibinun)
Save saturation in binary format to ibinun.
Definition: gwf-npf.f90:2639
subroutine prepcheck(this)
Initialize and check NPF data.
Definition: gwf-npf.f90:1588
subroutine npf_df(this, dis, xt3d, ingnc, invsc, npf_options)
Define the NPF package instance.
Definition: gwf-npf.f90:205
subroutine npf_ad(this, nodes, hold, hnew, irestore)
Advance.
Definition: gwf-npf.f90:381
subroutine sgwf_npf_wdmsg(this, icode, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
Print wet/dry message.
Definition: gwf-npf.f90:2229
subroutine sgwf_npf_qcalc(this, n, m, hn, hm, icon, qnm)
Flow between two cells.
Definition: gwf-npf.f90:799
subroutine npf_rp(this)
Read and prepare method for package.
Definition: gwf-npf.f90:366
real(dp) function calcsatthickness(this, n, m, ihc)
Calculate saturated thickness between cell n and m.
Definition: gwf-npf.f90:2782
subroutine store_original_k_arrays(this, ncells, njas)
@ brief Store backup copy of hydraulic conductivity when the VSC package is activate
Definition: gwf-npf.f90:1143
subroutine allocate_scalars(this)
@ brief Allocate scalars
Definition: gwf-npf.f90:1046
subroutine npf_cq(this, hnew, flowja)
Calculate flowja.
Definition: gwf-npf.f90:746
subroutine npf_ar(this, ic, vsc, ibound, hnew)
Allocate and read this NPF instance.
Definition: gwf-npf.f90:286
subroutine check_options(this)
Check for conflicting NPF options.
Definition: gwf-npf.f90:1387
subroutine npf_cf(this, kiter, nodes, hnew)
Routines associated fill coefficients.
Definition: gwf-npf.f90:445
General-purpose hydrogeologic functions.
Definition: HGeoUtil.f90:2
real(dp) function, public hyeff(k11, k22, k33, ang1, ang2, ang3, vg1, vg2, vg3, iavgmeth)
Calculate the effective horizontal hydraulic conductivity from an ellipse using a specified direction...
Definition: HGeoUtil.f90:31
integer(i4b) function, public getunit()
Get a free unit number.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
This module defines variable data types.
Definition: kind.f90:8
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
subroutine, public memorystore_remove(component, subcomponent, context)
subroutine, public get_isize(name, mem_path, isize)
@ brief Get the number of elements for this variable
This module contains the base numerical package type.
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
Definition: Sim.f90:203
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
character(len=linelength) idm_context
real(dp) function squadraticsaturation(top, bot, x, eps)
@ brief sQuadraticSaturation
real(dp) function squadraticsaturationderivative(top, bot, x, eps)
@ brief Derivative of the quadratic saturation function
This module contains the SourceCommonModule.
Definition: SourceCommon.f90:7
logical(lgp) function, public filein_fname(filename, tagname, input_mempath, input_fname)
enforce and set a single input filename provided via FILEIN keyword
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
This module contains time-varying conductivity package methods.
Definition: gwf-tvk.f90:8
subroutine, public tvk_cr(tvk, name_model, inunit, iout)
Create a new TvkType object.
Definition: gwf-tvk.f90:53
subroutine, public xt3d_cr(xt3dobj, name_model, inunit, iout, ldispopt)
Create a new xt3d object.
This class is used to store a single deferred-length character string. It was designed to work in an ...
Definition: CharString.f90:23
Data structure and helper methods for passing NPF options into npf_df, as an alternative to reading t...
Helper class with work arrays for the SPDIS calculation in NPF.