MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
Data Types | Functions/Subroutines
gwfnpfmodule Module Reference

Data Types

type  gwfnpftype
 

Functions/Subroutines

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. More...
 
subroutine npf_df (this, dis, xt3d, ingnc, invsc, npf_options)
 Define the NPF package instance. More...
 
subroutine npf_ac (this, moffset, sparse)
 Add connections for extended neighbors to the sparse matrix. More...
 
subroutine npf_mc (this, moffset, matrix_sln)
 Map connections and construct iax, jax, and idxglox. More...
 
subroutine npf_ar (this, ic, vsc, ibound, hnew)
 Allocate and read this NPF instance. More...
 
subroutine npf_rp (this)
 Read and prepare method for package. More...
 
subroutine npf_ad (this, nodes, hold, hnew, irestore)
 Advance. More...
 
subroutine npf_cf (this, kiter, nodes, hnew)
 Routines associated fill coefficients. More...
 
subroutine npf_fc (this, kiter, matrix_sln, idxglo, rhs, hnew)
 Formulate coefficients. More...
 
subroutine npf_fn (this, kiter, matrix_sln, idxglo, rhs, hnew)
 Fill newton terms. More...
 
subroutine npf_nur (this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax)
 Under-relaxation. More...
 
subroutine npf_cq (this, hnew, flowja)
 Calculate flowja. More...
 
subroutine sgwf_npf_thksat (this, n, hn, thksat)
 Fractional cell saturation. More...
 
subroutine sgwf_npf_qcalc (this, n, m, hn, hm, icon, qnm)
 Flow between two cells. More...
 
subroutine npf_save_model_flows (this, flowja, icbcfl, icbcun)
 Record flowja and calculate specific discharge if requested. More...
 
subroutine npf_print_model_flows (this, ibudfl, flowja)
 Print budget. More...
 
subroutine npf_da (this)
 Deallocate variables. More...
 
subroutine allocate_scalars (this)
 @ brief Allocate scalars More...
 
subroutine store_original_k_arrays (this, ncells, njas)
 @ brief Store backup copy of hydraulic conductivity when the VSC package is activate More...
 
subroutine allocate_arrays (this, ncells, njas)
 Allocate npf arrays. More...
 
subroutine log_options (this, found)
 Log npf options sourced from the input mempath. More...
 
subroutine source_options (this)
 Update simulation options from input mempath. More...
 
subroutine set_options (this, options)
 Set options in the NPF object. More...
 
subroutine check_options (this)
 Check for conflicting NPF options. More...
 
subroutine log_griddata (this, found)
 Write dimensions to list file. More...
 
subroutine source_griddata (this)
 Update simulation griddata from input mempath. More...
 
subroutine prepcheck (this)
 Initialize and check NPF data. More...
 
subroutine preprocess_input (this)
 preprocess the NPF input data More...
 
subroutine calc_condsat (this, node, upperOnly)
 Calculate CONDSAT array entries for the given node. More...
 
real(dp) function calc_initial_sat (this, n)
 Calculate initial saturation for the given node. More...
 
subroutine sgwf_npf_wetdry (this, kiter, hnew)
 Perform wetting and drying. More...
 
subroutine rewet_check (this, kiter, node, hm, ibdm, ihc, hnew, irewet)
 Determine if a cell should rewet. More...
 
subroutine sgwf_npf_wdmsg (this, icode, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
 Print wet/dry message. More...
 
real(dp) function hy_eff (this, n, m, ihc, ipos, vg)
 Calculate the effective hydraulic conductivity for the n-m connection. More...
 
subroutine calc_spdis (this, flowja)
 Calculate the 3 components of specific discharge at the cell center. More...
 
subroutine sav_spdis (this, ibinun)
 Save specific discharge in binary format to ibinun. More...
 
subroutine sav_sat (this, ibinun)
 Save saturation in binary format to ibinun. More...
 
subroutine increase_edge_count (this, nedges)
 Reserve space for nedges cells that have an edge on them. More...
 
subroutine set_edge_properties (this, nodedge, ihcedge, q, area, nx, ny, distance)
 Provide the npf package with edge properties. More...
 
real(dp) function calcsatthickness (this, n, m, ihc)
 Calculate saturated thickness between cell n and m. More...
 

Function/Subroutine Documentation

◆ allocate_arrays()

subroutine gwfnpfmodule::allocate_arrays ( class(gwfnpftype)  this,
integer(i4b), intent(in)  ncells,
integer(i4b), intent(in)  njas 
)

Definition at line 1138 of file gwf-npf.f90.

1139  ! -- dummy
1140  class(GwfNpftype) :: this
1141  integer(I4B), intent(in) :: ncells
1142  integer(I4B), intent(in) :: njas
1143  ! -- local
1144  integer(I4B) :: n
1145  !
1146  call mem_allocate(this%ithickstartflag, ncells, 'ITHICKSTARTFLAG', &
1147  this%memoryPath)
1148  call mem_allocate(this%icelltype, ncells, 'ICELLTYPE', this%memoryPath)
1149  call mem_allocate(this%k11, ncells, 'K11', this%memoryPath)
1150  call mem_allocate(this%sat, ncells, 'SAT', this%memoryPath)
1151  call mem_allocate(this%condsat, njas, 'CONDSAT', this%memoryPath)
1152  !
1153  ! -- Optional arrays dimensioned to full size initially
1154  call mem_allocate(this%k22, ncells, 'K22', this%memoryPath)
1155  call mem_allocate(this%k33, ncells, 'K33', this%memoryPath)
1156  call mem_allocate(this%wetdry, ncells, 'WETDRY', this%memoryPath)
1157  call mem_allocate(this%angle1, ncells, 'ANGLE1', this%memoryPath)
1158  call mem_allocate(this%angle2, ncells, 'ANGLE2', this%memoryPath)
1159  call mem_allocate(this%angle3, ncells, 'ANGLE3', this%memoryPath)
1160  !
1161  ! -- Optional arrays
1162  call mem_allocate(this%ibotnode, 0, 'IBOTNODE', this%memoryPath)
1163  call mem_allocate(this%nodedge, 0, 'NODEDGE', this%memoryPath)
1164  call mem_allocate(this%ihcedge, 0, 'IHCEDGE', this%memoryPath)
1165  call mem_allocate(this%propsedge, 0, 0, 'PROPSEDGE', this%memoryPath)
1166  !
1167  ! -- Optional arrays only needed when vsc package is active
1168  call mem_allocate(this%k11input, 0, 'K11INPUT', this%memoryPath)
1169  call mem_allocate(this%k22input, 0, 'K22INPUT', this%memoryPath)
1170  call mem_allocate(this%k33input, 0, 'K33INPUT', this%memoryPath)
1171  !
1172  ! -- Specific discharge is (re-)allocated when nedges is known
1173  call mem_allocate(this%spdis, 3, 0, 'SPDIS', this%memoryPath)
1174  !
1175  ! -- Time-varying property flag arrays
1176  call mem_allocate(this%nodekchange, ncells, 'NODEKCHANGE', this%memoryPath)
1177  !
1178  ! -- initialize iangle1, iangle2, iangle3, and wetdry
1179  do n = 1, ncells
1180  this%angle1(n) = dzero
1181  this%angle2(n) = dzero
1182  this%angle3(n) = dzero
1183  this%wetdry(n) = dzero
1184  this%nodekchange(n) = dzero
1185  end do
1186  !
1187  ! -- allocate variable names
1188  allocate (this%aname(this%iname))
1189  this%aname = [' ICELLTYPE', ' K', &
1190  ' K33', ' K22', &
1191  ' WETDRY', ' ANGLE1', &
1192  ' ANGLE2', ' ANGLE3']

◆ allocate_scalars()

subroutine gwfnpfmodule::allocate_scalars ( class(gwfnpftype)  this)

Allocate and initialize scalars for the VSC package. The base model allocate scalars method is also called.

Definition at line 1021 of file gwf-npf.f90.

1022  ! -- modules
1023  use memoryhelpermodule, only: create_mem_path
1024  ! -- dummy
1025  class(GwfNpftype) :: this
1026  !
1027  ! -- allocate scalars in NumericalPackageType
1028  call this%NumericalPackageType%allocate_scalars()
1029  !
1030  ! -- Allocate scalars
1031  call mem_allocate(this%iname, 'INAME', this%memoryPath)
1032  call mem_allocate(this%ixt3d, 'IXT3D', this%memoryPath)
1033  call mem_allocate(this%ixt3drhs, 'IXT3DRHS', this%memoryPath)
1034  call mem_allocate(this%satomega, 'SATOMEGA', this%memoryPath)
1035  call mem_allocate(this%hnoflo, 'HNOFLO', this%memoryPath)
1036  call mem_allocate(this%hdry, 'HDRY', this%memoryPath)
1037  call mem_allocate(this%icellavg, 'ICELLAVG', this%memoryPath)
1038  call mem_allocate(this%iavgkeff, 'IAVGKEFF', this%memoryPath)
1039  call mem_allocate(this%ik22, 'IK22', this%memoryPath)
1040  call mem_allocate(this%ik33, 'IK33', this%memoryPath)
1041  call mem_allocate(this%ik22overk, 'IK22OVERK', this%memoryPath)
1042  call mem_allocate(this%ik33overk, 'IK33OVERK', this%memoryPath)
1043  call mem_allocate(this%iperched, 'IPERCHED', this%memoryPath)
1044  call mem_allocate(this%ivarcv, 'IVARCV', this%memoryPath)
1045  call mem_allocate(this%idewatcv, 'IDEWATCV', this%memoryPath)
1046  call mem_allocate(this%ithickstrt, 'ITHICKSTRT', this%memoryPath)
1047  call mem_allocate(this%icalcspdis, 'ICALCSPDIS', this%memoryPath)
1048  call mem_allocate(this%isavspdis, 'ISAVSPDIS', this%memoryPath)
1049  call mem_allocate(this%isavsat, 'ISAVSAT', this%memoryPath)
1050  call mem_allocate(this%irewet, 'IREWET', this%memoryPath)
1051  call mem_allocate(this%wetfct, 'WETFCT', this%memoryPath)
1052  call mem_allocate(this%iwetit, 'IWETIT', this%memoryPath)
1053  call mem_allocate(this%ihdwet, 'IHDWET', this%memoryPath)
1054  call mem_allocate(this%iangle1, 'IANGLE1', this%memoryPath)
1055  call mem_allocate(this%iangle2, 'IANGLE2', this%memoryPath)
1056  call mem_allocate(this%iangle3, 'IANGLE3', this%memoryPath)
1057  call mem_allocate(this%iwetdry, 'IWETDRY', this%memoryPath)
1058  call mem_allocate(this%nedges, 'NEDGES', this%memoryPath)
1059  call mem_allocate(this%lastedge, 'LASTEDGE', this%memoryPath)
1060  call mem_allocate(this%intvk, 'INTVK', this%memoryPath)
1061  call mem_allocate(this%invsc, 'INVSC', this%memoryPath)
1062  call mem_allocate(this%kchangeper, 'KCHANGEPER', this%memoryPath)
1063  call mem_allocate(this%kchangestp, 'KCHANGESTP', this%memoryPath)
1064  !
1065  ! -- set pointer to inewtonur
1066  call mem_setptr(this%igwfnewtonur, 'INEWTONUR', &
1067  create_mem_path(this%name_model))
1068  !
1069  ! -- Initialize value
1070  this%iname = 8
1071  this%ixt3d = 0
1072  this%ixt3drhs = 0
1073  this%satomega = dzero
1074  this%hnoflo = dhnoflo !1.d30
1075  this%hdry = dhdry !-1.d30
1076  this%icellavg = ccond_hmean
1077  this%iavgkeff = 0
1078  this%ik22 = 0
1079  this%ik33 = 0
1080  this%ik22overk = 0
1081  this%ik33overk = 0
1082  this%iperched = 0
1083  this%ivarcv = 0
1084  this%idewatcv = 0
1085  this%ithickstrt = 0
1086  this%icalcspdis = 0
1087  this%isavspdis = 0
1088  this%isavsat = 0
1089  this%irewet = 0
1090  this%wetfct = done
1091  this%iwetit = 1
1092  this%ihdwet = 0
1093  this%iangle1 = 0
1094  this%iangle2 = 0
1095  this%iangle3 = 0
1096  this%iwetdry = 0
1097  this%nedges = 0
1098  this%lastedge = 0
1099  this%intvk = 0
1100  this%invsc = 0
1101  this%kchangeper = 0
1102  this%kchangestp = 0
1103  !
1104  ! -- If newton is on, then NPF creates asymmetric matrix
1105  this%iasym = this%inewton
memoryhelpermodule::create_mem_path
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
Definition: MemoryHelper.f90:22
Here is the call graph for this function:

◆ calc_condsat()

subroutine gwfnpfmodule::calc_condsat ( class(gwfnpftype)  this,
integer(i4b), intent(in)  node,
logical, intent(in)  upperOnly 
)

Calculate saturated conductances for all connections of the given node, or optionally for the upper portion of the matrix only.

Definition at line 1924 of file gwf-npf.f90.

1925  ! -- dummy variables
1926  class(GwfNpfType) :: this
1927  integer(I4B), intent(in) :: node
1928  logical, intent(in) :: upperOnly
1929  ! -- local variables
1930  integer(I4B) :: ii, m, n, ihc, jj
1931  real(DP) :: topm, topn, topnode, botm, botn, botnode, satm, satn, satnode
1932  real(DP) :: hyn, hym, hn, hm, fawidth, csat
1933  !
1934  satnode = this%calc_initial_sat(node)
1935  !
1936  topnode = this%dis%top(node)
1937  botnode = this%dis%bot(node)
1938  !
1939  ! -- Go through the connecting cells
1940  do ii = this%dis%con%ia(node) + 1, this%dis%con%ia(node + 1) - 1
1941  !
1942  ! -- Set the m cell number and cycle if lower triangle connection and
1943  ! -- we're not updating both upper and lower matrix parts for this node
1944  m = this%dis%con%ja(ii)
1945  jj = this%dis%con%jas(ii)
1946  if (m < node) then
1947  if (upperonly) cycle
1948  ! m => node, n => neighbour
1949  n = m
1950  m = node
1951  topm = topnode
1952  botm = botnode
1953  satm = satnode
1954  topn = this%dis%top(n)
1955  botn = this%dis%bot(n)
1956  satn = this%calc_initial_sat(n)
1957  else
1958  ! n => node, m => neighbour
1959  n = node
1960  topn = topnode
1961  botn = botnode
1962  satn = satnode
1963  topm = this%dis%top(m)
1964  botm = this%dis%bot(m)
1965  satm = this%calc_initial_sat(m)
1966  end if
1967  !
1968  ihc = this%dis%con%ihc(jj)
1969  hyn = this%hy_eff(n, m, ihc, ipos=ii)
1970  hym = this%hy_eff(m, n, ihc, ipos=ii)
1971  if (this%ithickstartflag(n) == 0) then
1972  hn = topn
1973  else
1974  hn = this%ic%strt(n)
1975  end if
1976  if (this%ithickstartflag(m) == 0) then
1977  hm = topm
1978  else
1979  hm = this%ic%strt(m)
1980  end if
1981  !
1982  ! -- Calculate conductance depending on whether connection is
1983  ! vertical (0), horizontal (1), or staggered horizontal (2)
1984  if (ihc == c3d_vertical) then
1985  !
1986  ! -- Vertical conductance for fully saturated conditions
1987  csat = vcond(1, 1, 1, 1, 0, 1, 1, done, &
1988  botn, botm, &
1989  hyn, hym, &
1990  satn, satm, &
1991  topn, topm, &
1992  botn, botm, &
1993  this%dis%con%hwva(jj))
1994  else
1995  !
1996  ! -- Horizontal conductance for fully saturated conditions
1997  fawidth = this%dis%con%hwva(jj)
1998  csat = hcond(1, 1, 1, 1, 0, &
1999  ihc, &
2000  this%icellavg, &
2001  done, &
2002  hn, hm, satn, satm, hyn, hym, &
2003  topn, topm, &
2004  botn, botm, &
2005  this%dis%con%cl1(jj), &
2006  this%dis%con%cl2(jj), &
2007  fawidth)
2008  end if
2009  this%condsat(jj) = csat
2010  end do
Here is the call graph for this function:

◆ calc_initial_sat()

real(dp) function gwfnpfmodule::calc_initial_sat ( class(gwfnpftype)  this,
integer(i4b), intent(in)  n 
)
private

Calculate saturation as a fraction of thickness for the given node, used for saturated conductance calculations: full thickness by default (1.0) or saturation based on initial conditions if the THICKSTRT option is used.

Definition at line 2020 of file gwf-npf.f90.

2021  ! -- dummy variables
2022  class(GwfNpfType) :: this
2023  integer(I4B), intent(in) :: n
2024  ! -- Return
2025  real(DP) :: satn
2026  !
2027  satn = done
2028  if (this%ibound(n) /= 0 .and. this%ithickstartflag(n) /= 0) then
2029  call this%thksat(n, this%ic%strt(n), satn)
2030  end if

◆ calc_spdis()

subroutine gwfnpfmodule::calc_spdis ( class(gwfnpftype)  this,
real(dp), dimension(:), intent(in)  flowja 
)
private

Definition at line 2333 of file gwf-npf.f90.

2334  ! -- modules
2335  use simmodule, only: store_error
2336  ! -- dummy
2337  class(GwfNpfType) :: this
2338  real(DP), intent(in), dimension(:) :: flowja
2339  ! -- local
2340  integer(I4B) :: n
2341  integer(I4B) :: m
2342  integer(I4B) :: ipos
2343  integer(I4B) :: isympos
2344  integer(I4B) :: ihc
2345  integer(I4B) :: ic
2346  integer(I4B) :: iz
2347  integer(I4B) :: nc
2348  integer(I4B) :: ncz
2349  real(DP) :: qz
2350  real(DP) :: vx
2351  real(DP) :: vy
2352  real(DP) :: vz
2353  real(DP) :: xn
2354  real(DP) :: yn
2355  real(DP) :: zn
2356  real(DP) :: xc
2357  real(DP) :: yc
2358  real(DP) :: zc
2359  real(DP) :: cl1
2360  real(DP) :: cl2
2361  real(DP) :: dltot
2362  real(DP) :: ooclsum
2363  real(DP) :: dsumx
2364  real(DP) :: dsumy
2365  real(DP) :: dsumz
2366  real(DP) :: denom
2367  real(DP) :: area
2368  real(DP) :: dz
2369  real(DP) :: axy
2370  real(DP) :: ayx
2371  real(DP), allocatable, dimension(:) :: vi
2372  real(DP), allocatable, dimension(:) :: di
2373  real(DP), allocatable, dimension(:) :: viz
2374  real(DP), allocatable, dimension(:) :: diz
2375  real(DP), allocatable, dimension(:) :: nix
2376  real(DP), allocatable, dimension(:) :: niy
2377  real(DP), allocatable, dimension(:) :: wix
2378  real(DP), allocatable, dimension(:) :: wiy
2379  real(DP), allocatable, dimension(:) :: wiz
2380  real(DP), allocatable, dimension(:) :: bix
2381  real(DP), allocatable, dimension(:) :: biy
2382  logical :: nozee = .true.
2383  !
2384  ! -- Ensure dis has necessary information
2385  if (this%icalcspdis /= 0 .and. this%dis%con%ianglex == 0) then
2386  call store_error('Error. ANGLDEGX not provided in '// &
2387  'discretization file. ANGLDEGX required for '// &
2388  'calculation of specific discharge.', terminate=.true.)
2389  end if
2390  !
2391  ! -- Find max number of connections and allocate weight arrays
2392  nc = 0
2393  do n = 1, this%dis%nodes
2394  !
2395  ! -- Count internal model connections
2396  ic = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
2397  !
2398  ! -- Count edge connections
2399  do m = 1, this%nedges
2400  if (this%nodedge(m) == n) then
2401  ic = ic + 1
2402  end if
2403  end do
2404  !
2405  ! -- Set max number of connections for any cell
2406  if (ic > nc) nc = ic
2407  end do
2408  !
2409  ! -- Allocate storage arrays needed for cell-centered spdis calculation
2410  allocate (vi(nc))
2411  allocate (di(nc))
2412  allocate (viz(nc))
2413  allocate (diz(nc))
2414  allocate (nix(nc))
2415  allocate (niy(nc))
2416  allocate (wix(nc))
2417  allocate (wiy(nc))
2418  allocate (wiz(nc))
2419  allocate (bix(nc))
2420  allocate (biy(nc))
2421  !
2422  ! -- Go through each cell and calculate specific discharge
2423  do n = 1, this%dis%nodes
2424  !
2425  ! -- first calculate geometric properties for x and y directions and
2426  ! the specific discharge at a face (vi)
2427  ic = 0
2428  iz = 0
2429  vi(:) = dzero
2430  di(:) = dzero
2431  viz(:) = dzero
2432  diz(:) = dzero
2433  nix(:) = dzero
2434  niy(:) = dzero
2435  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2436  m = this%dis%con%ja(ipos)
2437  isympos = this%dis%con%jas(ipos)
2438  ihc = this%dis%con%ihc(isympos)
2439  area = this%dis%con%hwva(isympos)
2440  if (ihc == c3d_vertical) then
2441  !
2442  ! -- vertical connection
2443  iz = iz + 1
2444  !call this%dis%connection_normal(n, m, ihc, xn, yn, zn, ipos)
2445  call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2446  ihc, xc, yc, zc, dltot)
2447  cl1 = this%dis%con%cl1(isympos)
2448  cl2 = this%dis%con%cl2(isympos)
2449  if (m < n) then
2450  cl1 = this%dis%con%cl2(isympos)
2451  cl2 = this%dis%con%cl1(isympos)
2452  end if
2453  ooclsum = done / (cl1 + cl2)
2454  diz(iz) = dltot * cl1 * ooclsum
2455  qz = flowja(ipos)
2456  if (n > m) qz = -qz
2457  viz(iz) = qz / area
2458  else
2459  !
2460  ! -- horizontal connection
2461  ic = ic + 1
2462  dz = thksatnm(this%ibound(n), this%ibound(m), &
2463  this%icelltype(n), this%icelltype(m), &
2464  this%inewton, ihc, &
2465  this%hnew(n), this%hnew(m), this%sat(n), this%sat(m), &
2466  this%dis%top(n), this%dis%top(m), this%dis%bot(n), &
2467  this%dis%bot(m))
2468  area = area * dz
2469  call this%dis%connection_normal(n, m, ihc, xn, yn, zn, ipos)
2470  call this%dis%connection_vector(n, m, nozee, this%sat(n), this%sat(m), &
2471  ihc, xc, yc, zc, dltot)
2472  cl1 = this%dis%con%cl1(isympos)
2473  cl2 = this%dis%con%cl2(isympos)
2474  if (m < n) then
2475  cl1 = this%dis%con%cl2(isympos)
2476  cl2 = this%dis%con%cl1(isympos)
2477  end if
2478  ooclsum = done / (cl1 + cl2)
2479  nix(ic) = -xn
2480  niy(ic) = -yn
2481  di(ic) = dltot * cl1 * ooclsum
2482  if (area > dzero) then
2483  vi(ic) = flowja(ipos) / area
2484  else
2485  vi(ic) = dzero
2486  end if
2487  end if
2488  end do
2489  !
2490  ! -- Look through edge flows that may have been provided by an exchange
2491  ! and incorporate them into the averaging arrays
2492  do m = 1, this%nedges
2493  if (this%nodedge(m) == n) then
2494  !
2495  ! -- propsedge: (Q, area, nx, ny, distance)
2496  ihc = this%ihcedge(m)
2497  area = this%propsedge(2, m)
2498  if (ihc == c3d_vertical) then
2499  iz = iz + 1
2500  viz(iz) = this%propsedge(1, m) / area
2501  diz(iz) = this%propsedge(5, m)
2502  else
2503  ic = ic + 1
2504  nix(ic) = -this%propsedge(3, m)
2505  niy(ic) = -this%propsedge(4, m)
2506  di(ic) = this%propsedge(5, m)
2507  if (area > dzero) then
2508  vi(ic) = this%propsedge(1, m) / area
2509  else
2510  vi(ic) = dzero
2511  end if
2512  end if
2513  end if
2514  end do
2515  !
2516  ! -- Assign number of vertical and horizontal connections
2517  ncz = iz
2518  nc = ic
2519  !
2520  ! -- calculate z weight (wiz) and z velocity
2521  if (ncz == 1) then
2522  wiz(1) = done
2523  else
2524  dsumz = dzero
2525  do iz = 1, ncz
2526  dsumz = dsumz + diz(iz)
2527  end do
2528  denom = (ncz - done)
2529  if (denom < dzero) denom = dzero
2530  dsumz = dsumz + dem10 * dsumz
2531  do iz = 1, ncz
2532  if (dsumz > dzero) wiz(iz) = done - diz(iz) / dsumz
2533  if (denom > 0) then
2534  wiz(iz) = wiz(iz) / denom
2535  else
2536  wiz(iz) = dzero
2537  end if
2538  end do
2539  end if
2540  vz = dzero
2541  do iz = 1, ncz
2542  vz = vz + wiz(iz) * viz(iz)
2543  end do
2544  !
2545  ! -- distance-based weighting
2546  nc = ic
2547  dsumx = dzero
2548  dsumy = dzero
2549  dsumz = dzero
2550  do ic = 1, nc
2551  wix(ic) = di(ic) * abs(nix(ic))
2552  wiy(ic) = di(ic) * abs(niy(ic))
2553  dsumx = dsumx + wix(ic)
2554  dsumy = dsumy + wiy(ic)
2555  end do
2556  !
2557  ! -- Finish computing omega weights. Add a tiny bit
2558  ! to dsum so that the normalized omega weight later
2559  ! evaluates to (essentially) 1 in the case of a single
2560  ! relevant connection, avoiding 0/0.
2561  dsumx = dsumx + dem10 * dsumx
2562  dsumy = dsumy + dem10 * dsumy
2563  do ic = 1, nc
2564  wix(ic) = (dsumx - wix(ic)) * abs(nix(ic))
2565  wiy(ic) = (dsumy - wiy(ic)) * abs(niy(ic))
2566  end do
2567  !
2568  ! -- compute B weights
2569  dsumx = dzero
2570  dsumy = dzero
2571  do ic = 1, nc
2572  bix(ic) = wix(ic) * sign(done, nix(ic))
2573  biy(ic) = wiy(ic) * sign(done, niy(ic))
2574  dsumx = dsumx + wix(ic) * abs(nix(ic))
2575  dsumy = dsumy + wiy(ic) * abs(niy(ic))
2576  end do
2577  if (dsumx > dzero) dsumx = done / dsumx
2578  if (dsumy > dzero) dsumy = done / dsumy
2579  axy = dzero
2580  ayx = dzero
2581  do ic = 1, nc
2582  bix(ic) = bix(ic) * dsumx
2583  biy(ic) = biy(ic) * dsumy
2584  axy = axy + bix(ic) * niy(ic)
2585  ayx = ayx + biy(ic) * nix(ic)
2586  end do
2587  !
2588  ! -- Calculate specific discharge. The divide by zero checking below
2589  ! is problematic for cells with only one flow, such as can happen
2590  ! with triangular cells in corners. In this case, the resulting
2591  ! cell velocity will be calculated as zero. The method should be
2592  ! improved so that edge flows of zero are included in these
2593  ! calculations. But this needs to be done with consideration for LGR
2594  ! cases in which flows are submitted from an exchange.
2595  vx = dzero
2596  vy = dzero
2597  do ic = 1, nc
2598  vx = vx + (bix(ic) - axy * biy(ic)) * vi(ic)
2599  vy = vy + (biy(ic) - ayx * bix(ic)) * vi(ic)
2600  end do
2601  denom = done - axy * ayx
2602  if (denom /= dzero) then
2603  vx = vx / denom
2604  vy = vy / denom
2605  end if
2606  !
2607  this%spdis(1, n) = vx
2608  this%spdis(2, n) = vy
2609  this%spdis(3, n) = vz
2610  !
2611  end do
2612  !
2613  ! -- cleanup
2614  deallocate (vi)
2615  deallocate (di)
2616  deallocate (nix)
2617  deallocate (niy)
2618  deallocate (wix)
2619  deallocate (wiy)
2620  deallocate (wiz)
2621  deallocate (bix)
2622  deallocate (biy)
simmodule
This module contains simulation methods.
Definition: Sim.f90:10
simmodule::store_error
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
Here is the call graph for this function:

◆ calcsatthickness()

real(dp) function gwfnpfmodule::calcsatthickness ( class(gwfnpftype)  this,
integer(i4b)  n,
integer(i4b)  m,
integer(i4b)  ihc 
)
private
Parameters
thisthis NPF instance
nnode n
mnode m
ihc1 = horizontal connection, 0 for vertical
Returns
saturated thickness

Definition at line 2728 of file gwf-npf.f90.

2729  ! -- dummy
2730  class(GwfNpfType) :: this !< this NPF instance
2731  integer(I4B) :: n !< node n
2732  integer(I4B) :: m !< node m
2733  integer(I4B) :: ihc !< 1 = horizontal connection, 0 for vertical
2734  ! -- return
2735  real(DP) :: satThickness !< saturated thickness
2736  !
2737  satthickness = thksatnm(this%ibound(n), &
2738  this%ibound(m), &
2739  this%icelltype(n), &
2740  this%icelltype(m), &
2741  this%inewton, &
2742  ihc, &
2743  this%hnew(n), &
2744  this%hnew(m), &
2745  this%sat(n), &
2746  this%sat(m), &
2747  this%dis%top(n), &
2748  this%dis%top(m), &
2749  this%dis%bot(n), &
2750  this%dis%bot(m))
Here is the call graph for this function:

◆ check_options()

subroutine gwfnpfmodule::check_options ( class(gwfnpftype)  this)
private

Definition at line 1360 of file gwf-npf.f90.

1361  ! -- modules
1362  use simmodule, only: store_error, count_errors, store_error_filename
1363  use constantsmodule, only: linelength
1364  ! -- dummy
1365  class(GwfNpftype) :: this
1366  !
1367  ! -- set omega value used for saturation calculations
1368  if (this%inewton > 0) then
1369  this%satomega = dem6
1370  end if
1371  !
1372  if (this%inewton > 0) then
1373  if (this%iperched > 0) then
1374  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1375  'BE USED WITH PERCHED OPTION.'
1376  call store_error(errmsg)
1377  end if
1378  if (this%ivarcv > 0) then
1379  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1380  'BE USED WITH VARIABLECV OPTION.'
1381  call store_error(errmsg)
1382  end if
1383  if (this%irewet > 0) then
1384  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. NEWTON OPTION CANNOT '// &
1385  'BE USED WITH REWET OPTION.'
1386  call store_error(errmsg)
1387  end if
1388  end if
1389  !
1390  if (this%ixt3d /= 0) then
1391  if (this%icellavg > 0) then
1392  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. '// &
1393  'ALTERNATIVE_CELL_AVERAGING OPTION '// &
1394  'CANNOT BE USED WITH XT3D OPTION.'
1395  call store_error(errmsg)
1396  end if
1397  if (this%ithickstrt > 0) then
1398  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. THICKSTRT OPTION '// &
1399  'CANNOT BE USED WITH XT3D OPTION.'
1400  call store_error(errmsg)
1401  end if
1402  if (this%iperched > 0) then
1403  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. PERCHED OPTION '// &
1404  'CANNOT BE USED WITH XT3D OPTION.'
1405  call store_error(errmsg)
1406  end if
1407  if (this%ivarcv > 0) then
1408  write (errmsg, '(a)') 'ERROR IN NPF OPTIONS. VARIABLECV OPTION '// &
1409  'CANNOT BE USED WITH XT3D OPTION.'
1410  call store_error(errmsg)
1411  end if
1412  end if
1413  !
1414  ! -- Terminate if errors
1415  if (count_errors() > 0) then
1416  call store_error_filename(this%input_fname)
1417  end if
constantsmodule
This module contains simulation constants.
Definition: Constants.f90:9
constantsmodule::linelength
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
simmodule::count_errors
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
simmodule::store_error_filename
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
Definition: Sim.f90:203
Here is the call graph for this function:

◆ hy_eff()

real(dp) function gwfnpfmodule::hy_eff ( class(gwfnpftype)  this,
integer(i4b), intent(in)  n,
integer(i4b), intent(in)  m,
integer(i4b), intent(in)  ihc,
integer(i4b), intent(in), optional  ipos,
real(dp), dimension(3), intent(in), optional  vg 
)

n is primary node node number m is connected node (not used if vg is provided) ihc is horizontal indicator (0 vertical, 1 horizontal, 2 vertically staggered) ipos_opt is position of connection in ja array vg is the global unit vector that expresses the direction from which to calculate an effective hydraulic conductivity.

Definition at line 2254 of file gwf-npf.f90.

2255  ! -- return
2256  real(DP) :: hy
2257  ! -- dummy
2258  class(GwfNpfType) :: this
2259  integer(I4B), intent(in) :: n
2260  integer(I4B), intent(in) :: m
2261  integer(I4B), intent(in) :: ihc
2262  integer(I4B), intent(in), optional :: ipos
2263  real(DP), dimension(3), intent(in), optional :: vg
2264  ! -- local
2265  integer(I4B) :: iipos
2266  real(DP) :: hy11, hy22, hy33
2267  real(DP) :: ang1, ang2, ang3
2268  real(DP) :: vg1, vg2, vg3
2269  !
2270  ! -- Initialize
2271  iipos = 0
2272  if (present(ipos)) iipos = ipos
2273  hy11 = this%k11(n)
2274  hy22 = this%k11(n)
2275  hy33 = this%k11(n)
2276  hy22 = this%k22(n)
2277  hy33 = this%k33(n)
2278  !
2279  ! -- Calculate effective K based on whether connection is vertical
2280  ! or horizontal
2281  if (ihc == c3d_vertical) then
2282  !
2283  ! -- Handle rotated anisotropy case that would affect the effective
2284  ! vertical hydraulic conductivity
2285  hy = hy33
2286  if (this%iangle2 > 0) then
2287  if (present(vg)) then
2288  vg1 = vg(1)
2289  vg2 = vg(2)
2290  vg3 = vg(3)
2291  else
2292  call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2293  end if
2294  ang1 = this%angle1(n)
2295  ang2 = this%angle2(n)
2296  ang3 = dzero
2297  if (this%iangle3 > 0) ang3 = this%angle3(n)
2298  hy = hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2299  this%iavgkeff)
2300  end if
2301  !
2302  else
2303  !
2304  ! -- Handle horizontal case
2305  hy = hy11
2306  if (this%ik22 > 0) then
2307  if (present(vg)) then
2308  vg1 = vg(1)
2309  vg2 = vg(2)
2310  vg3 = vg(3)
2311  else
2312  call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, iipos)
2313  end if
2314  ang1 = dzero
2315  ang2 = dzero
2316  ang3 = dzero
2317  if (this%iangle1 > 0) then
2318  ang1 = this%angle1(n)
2319  if (this%iangle2 > 0) then
2320  ang2 = this%angle2(n)
2321  if (this%iangle3 > 0) ang3 = this%angle3(n)
2322  end if
2323  end if
2324  hy = hyeff(hy11, hy22, hy33, ang1, ang2, ang3, vg1, vg2, vg3, &
2325  this%iavgkeff)
2326  end if
2327  !
2328  end if
Here is the call graph for this function:

◆ increase_edge_count()

subroutine gwfnpfmodule::increase_edge_count ( class(gwfnpftype)  this,
integer(i4b), intent(in)  nedges 
)
private

This must be called before the npfallocate_arrays routine, which is called from npfar.

Definition at line 2687 of file gwf-npf.f90.

2688  ! -- dummy
2689  class(GwfNpfType) :: this
2690  integer(I4B), intent(in) :: nedges
2691  !
2692  this%nedges = this%nedges + nedges

◆ log_griddata()

subroutine gwfnpfmodule::log_griddata ( class(gwfnpftype)  this,
type(gwfnpfparamfoundtype), intent(in)  found 
)

Definition at line 1422 of file gwf-npf.f90.

1423  ! -- modules
1424  use gwfnpfinputmodule, only: gwfnpfparamfoundtype
1425  ! -- dummy
1426  class(GwfNpfType) :: this
1427  type(GwfNpfParamFoundType), intent(in) :: found
1428  !
1429  write (this%iout, '(1x,a)') 'Setting NPF Griddata'
1430  !
1431  if (found%icelltype) then
1432  write (this%iout, '(4x,a)') 'ICELLTYPE set from input file'
1433  end if
1434  !
1435  if (found%k) then
1436  write (this%iout, '(4x,a)') 'K set from input file'
1437  end if
1438  !
1439  if (found%k33) then
1440  write (this%iout, '(4x,a)') 'K33 set from input file'
1441  else
1442  write (this%iout, '(4x,a)') 'K33 not provided. Setting K33 = K.'
1443  end if
1444  !
1445  if (found%k22) then
1446  write (this%iout, '(4x,a)') 'K22 set from input file'
1447  else
1448  write (this%iout, '(4x,a)') 'K22 not provided. Setting K22 = K.'
1449  end if
1450  !
1451  if (found%wetdry) then
1452  write (this%iout, '(4x,a)') 'WETDRY set from input file'
1453  end if
1454  !
1455  if (found%angle1) then
1456  write (this%iout, '(4x,a)') 'ANGLE1 set from input file'
1457  end if
1458  !
1459  if (found%angle2) then
1460  write (this%iout, '(4x,a)') 'ANGLE2 set from input file'
1461  end if
1462  !
1463  if (found%angle3) then
1464  write (this%iout, '(4x,a)') 'ANGLE3 set from input file'
1465  end if
1466  !
1467  write (this%iout, '(1x,a,/)') 'End Setting NPF Griddata'

◆ log_options()

subroutine gwfnpfmodule::log_options ( class(gwfnpftype)  this,
type(gwfnpfparamfoundtype), intent(in)  found 
)
private

Definition at line 1197 of file gwf-npf.f90.

1198  ! -- modules
1199  use kindmodule, only: lgp
1200  use gwfnpfinputmodule, only: gwfnpfparamfoundtype
1201  ! -- dummy
1202  class(GwfNpftype) :: this
1203  ! -- locals
1204  type(GwfNpfParamFoundType), intent(in) :: found
1205  !
1206  write (this%iout, '(1x,a)') 'Setting NPF Options'
1207  if (found%iprflow) &
1208  write (this%iout, '(4x,a)') 'Cell-by-cell flow information will be printed &
1209  &to listing file whenever ICBCFL is not zero.'
1210  if (found%ipakcb) &
1211  write (this%iout, '(4x,a)') 'Cell-by-cell flow information will be saved &
1212  &to binary file whenever ICBCFL is not zero.'
1213  if (found%cellavg) &
1214  write (this%iout, '(4x,a,i0)') 'Alternative cell averaging [1=logarithmic, &
1215  &2=AMT-LMK, 3=AMT-HMK] set to: ', &
1216  this%icellavg
1217  if (found%ithickstrt) &
1218  write (this%iout, '(4x,a)') 'THICKSTRT option has been activated.'
1219  if (found%iperched) &
1220  write (this%iout, '(4x,a)') 'Vertical flow will be adjusted for perched &
1221  &conditions.'
1222  if (found%ivarcv) &
1223  write (this%iout, '(4x,a)') 'Vertical conductance varies with water table.'
1224  if (found%idewatcv) &
1225  write (this%iout, '(4x,a)') 'Vertical conductance accounts for dewatered &
1226  &portion of an underlying cell.'
1227  if (found%ixt3d) write (this%iout, '(4x,a)') 'XT3D formulation is selected.'
1228  if (found%ixt3drhs) &
1229  write (this%iout, '(4x,a)') 'XT3D RHS formulation is selected.'
1230  if (found%isavspdis) &
1231  write (this%iout, '(4x,a)') 'Specific discharge will be calculated at cell &
1232  &centers and written to DATA-SPDIS in budget &
1233  &file when requested.'
1234  if (found%isavsat) &
1235  write (this%iout, '(4x,a)') 'Saturation will be written to DATA-SAT in &
1236  &budget file when requested.'
1237  if (found%ik22overk) &
1238  write (this%iout, '(4x,a)') 'Values specified for K22 are anisotropy &
1239  &ratios and will be multiplied by K before &
1240  &being used in calculations.'
1241  if (found%ik33overk) &
1242  write (this%iout, '(4x,a)') 'Values specified for K33 are anisotropy &
1243  &ratios and will be multiplied by K before &
1244  &being used in calculations.'
1245  if (found%inewton) &
1246  write (this%iout, '(4x,a)') 'NEWTON-RAPHSON method disabled for unconfined &
1247  &cells'
1248  if (found%satomega) &
1249  write (this%iout, '(4x,a,1pg15.6)') 'Saturation omega: ', this%satomega
1250  if (found%irewet) &
1251  write (this%iout, '(4x,a)') 'Rewetting is active.'
1252  if (found%wetfct) &
1253  write (this%iout, '(4x,a,1pg15.6)') &
1254  'Wetting factor (WETFCT) has been set to: ', this%wetfct
1255  if (found%iwetit) &
1256  write (this%iout, '(4x,a,i5)') &
1257  'Wetting iteration interval (IWETIT) has been set to: ', this%iwetit
1258  if (found%ihdwet) &
1259  write (this%iout, '(4x,a,i5)') &
1260  'Head rewet equation (IHDWET) has been set to: ', this%ihdwet
1261  write (this%iout, '(1x,a,/)') 'End Setting NPF Options'
kindmodule
This module defines variable data types.
Definition: kind.f90:8

◆ npf_ac()

subroutine gwfnpfmodule::npf_ac ( class(gwfnpftype)  this,
integer(i4b), intent(in)  moffset,
type(sparsematrix), intent(inout)  sparse 
)

Definition at line 248 of file gwf-npf.f90.

249  ! -- modules
250  use sparsemodule, only: sparsematrix
251  ! -- dummy
252  class(GwfNpftype) :: this
253  integer(I4B), intent(in) :: moffset
254  type(sparsematrix), intent(inout) :: sparse
255  !
256  ! -- Add extended neighbors (neighbors of neighbors)
257  if (this%ixt3d /= 0) call this%xt3d%xt3d_ac(moffset, sparse)

◆ npf_ad()

subroutine gwfnpfmodule::npf_ad ( class(gwfnpftype)  this,
integer(i4b), intent(in)  nodes,
real(dp), dimension(nodes), intent(inout)  hold,
real(dp), dimension(nodes), intent(inout)  hnew,
integer(i4b), intent(in)  irestore 
)
private

Sets hold (head old) to bot whenever a wettable cell is dry

Definition at line 362 of file gwf-npf.f90.

363  ! -- modules
364  use tdismodule, only: kper, kstp
365  !
366  implicit none
367  ! -- dummy
368  class(GwfNpfType) :: this
369  integer(I4B), intent(in) :: nodes
370  real(DP), dimension(nodes), intent(inout) :: hold
371  real(DP), dimension(nodes), intent(inout) :: hnew
372  integer(I4B), intent(in) :: irestore
373  ! -- local
374  integer(I4B) :: n
375  !
376  ! -- loop through all cells and set hold=bot if wettable cell is dry
377  if (this%irewet > 0) then
378  do n = 1, this%dis%nodes
379  if (this%wetdry(n) == dzero) cycle
380  if (this%ibound(n) /= 0) cycle
381  hold(n) = this%dis%bot(n)
382  end do
383  !
384  ! -- if restore state, then set hnew to DRY if it is a dry wettable cell
385  do n = 1, this%dis%nodes
386  if (this%wetdry(n) == dzero) cycle
387  if (this%ibound(n) /= 0) cycle
388  hnew(n) = dhdry
389  end do
390  end if
391  !
392  ! -- TVK
393  if (this%intvk /= 0) then
394  call this%tvk%ad()
395  end if
396  !
397  ! -- VSC
398  ! -- Hit the TVK-updated K's with VSC correction before calling/updating condsat
399  if (this%invsc /= 0) then
400  call this%vsc%update_k_with_vsc()
401  end if
402  !
403  ! -- If any K values have changed, we need to update CONDSAT or XT3D arrays
404  if (this%kchangeper == kper .and. this%kchangestp == kstp) then
405  if (this%ixt3d == 0) then
406  !
407  ! -- Update the saturated conductance for all connections
408  ! -- of the affected nodes
409  do n = 1, this%dis%nodes
410  if (this%nodekchange(n) == 1) then
411  call this%calc_condsat(n, .false.)
412  end if
413  end do
414  else
415  !
416  ! -- Recompute XT3D coefficients for permanently confined connections
417  if (this%xt3d%lamatsaved .and. .not. this%xt3d%ldispersion) then
418  call this%xt3d%xt3d_fcpc(this%dis%nodes, .true.)
419  end if
420  end if
421  end if
tdismodule::kstp
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
tdismodule::kper
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23

◆ npf_ar()

subroutine gwfnpfmodule::npf_ar ( class(gwfnpftype)  this,
type(gwfictype), intent(in), pointer  ic,
type(gwfvsctype), intent(in), pointer  vsc,
integer(i4b), dimension(:), intent(inout), pointer, contiguous  ibound,
real(dp), dimension(:), intent(inout), pointer, contiguous  hnew 
)
private

Allocate remaining package arrays, preprocess the input data and call *_ar on xt3d, when active.

Parameters
thisinstance of the NPF package
[in]icinitial conditions
[in]vscviscosity package
[in,out]iboundmodel ibound array
[in,out]hnewpointer to model head array

Definition at line 276 of file gwf-npf.f90.

277  ! -- modules
278  use memorymanagermodule, only: mem_reallocate
279  ! -- dummy
280  class(GwfNpftype) :: this !< instance of the NPF package
281  type(GwfIcType), pointer, intent(in) :: ic !< initial conditions
282  type(GwfVscType), pointer, intent(in) :: vsc !< viscosity package
283  integer(I4B), dimension(:), pointer, contiguous, intent(inout) :: ibound !< model ibound array
284  real(DP), dimension(:), pointer, contiguous, intent(inout) :: hnew !< pointer to model head array
285  ! -- local
286  integer(I4B) :: n
287  !
288  ! -- Store pointers to arguments that were passed in
289  this%ic => ic
290  this%ibound => ibound
291  this%hnew => hnew
292  !
293  if (this%icalcspdis == 1) then
294  call mem_reallocate(this%spdis, 3, this%dis%nodes, 'SPDIS', this%memoryPath)
295  call mem_reallocate(this%nodedge, this%nedges, 'NODEDGE', this%memoryPath)
296  call mem_reallocate(this%ihcedge, this%nedges, 'IHCEDGE', this%memoryPath)
297  call mem_reallocate(this%propsedge, 5, this%nedges, 'PROPSEDGE', &
298  this%memoryPath)
299  do n = 1, this%dis%nodes
300  this%spdis(:, n) = dzero
301  end do
302  end if
303  !
304  ! -- Store pointer to VSC if active
305  if (this%invsc /= 0) then
306  this%vsc => vsc
307  end if
308  !
309  ! -- allocate arrays to store original user input in case TVK/VSC modify them
310  if (this%invsc > 0) then
311  !
312  ! -- Reallocate arrays that store user-input values.
313  call mem_reallocate(this%k11input, this%dis%nodes, 'K11INPUT', &
314  this%memoryPath)
315  call mem_reallocate(this%k22input, this%dis%nodes, 'K22INPUT', &
316  this%memoryPath)
317  call mem_reallocate(this%k33input, this%dis%nodes, 'K33INPUT', &
318  this%memoryPath)
319  ! Allocate arrays that will store the original K values. When VSC active,
320  ! the current Kxx arrays carry the viscosity-adjusted K values.
321  ! This approach leverages existing functionality that makes use of K.
322  call this%store_original_k_arrays(this%dis%nodes, this%dis%njas)
323  end if
324  !
325  ! -- preprocess data
326  call this%preprocess_input()
327  !
328  ! -- xt3d
329  if (this%ixt3d /= 0) then
330  call this%xt3d%xt3d_ar(ibound, this%k11, this%ik33, this%k33, &
331  this%sat, this%ik22, this%k22, &
332  this%iangle1, this%iangle2, this%iangle3, &
333  this%angle1, this%angle2, this%angle3, &
334  this%inewton, this%icelltype)
335  end if
336  !
337  ! -- TVK
338  if (this%intvk /= 0) then
339  call this%tvk%ar(this%dis)
340  end if
memorymanagermodule::mem_reallocate
Definition: MemoryManager.f90:79

◆ npf_cf()

subroutine gwfnpfmodule::npf_cf ( class(gwfnpftype)  this,
integer(i4b)  kiter,
integer(i4b), intent(in)  nodes,
real(dp), dimension(nodes), intent(inout)  hnew 
)

Definition at line 426 of file gwf-npf.f90.

427  ! -- dummy
428  class(GwfNpfType) :: this
429  integer(I4B) :: kiter
430  integer(I4B), intent(in) :: nodes
431  real(DP), intent(inout), dimension(nodes) :: hnew
432  ! -- local
433  integer(I4B) :: n
434  real(DP) :: satn
435  !
436  ! -- Perform wetting and drying
437  if (this%inewton /= 1) then
438  call this%wd(kiter, hnew)
439  end if
440  !
441  ! -- Calculate saturation for convertible cells
442  do n = 1, this%dis%nodes
443  if (this%icelltype(n) /= 0) then
444  if (this%ibound(n) == 0) then
445  satn = dzero
446  else
447  call this%thksat(n, hnew(n), satn)
448  end if
449  this%sat(n) = satn
450  end if
451  end do

◆ npf_cq()

subroutine gwfnpfmodule::npf_cq ( class(gwfnpftype)  this,
real(dp), dimension(:), intent(inout)  hnew,
real(dp), dimension(:), intent(inout)  flowja 
)
private

Definition at line 727 of file gwf-npf.f90.

728  ! -- dummy
729  class(GwfNpfType) :: this
730  real(DP), intent(inout), dimension(:) :: hnew
731  real(DP), intent(inout), dimension(:) :: flowja
732  ! -- local
733  integer(I4B) :: n, ipos, m
734  real(DP) :: qnm
735  !
736  ! -- Calculate the flow across each cell face and store in flowja
737  !
738  if (this%ixt3d /= 0) then
739  call this%xt3d%xt3d_flowja(hnew, flowja)
740  else
741  !
742  do n = 1, this%dis%nodes
743  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
744  m = this%dis%con%ja(ipos)
745  if (m < n) cycle
746  call this%qcalc(n, m, hnew(n), hnew(m), ipos, qnm)
747  flowja(ipos) = qnm
748  flowja(this%dis%con%isym(ipos)) = -qnm
749  end do
750  end do
751  !
752  end if

◆ npf_cr()

subroutine, public gwfnpfmodule::npf_cr ( type(gwfnpftype), pointer  npfobj,
character(len=*), intent(in)  name_model,
character(len=*), intent(in)  input_mempath,
integer(i4b), intent(in)  inunit,
integer(i4b), intent(in)  iout 
)

Definition at line 152 of file gwf-npf.f90.

153  ! -- modules
154  use kindmodule, only: lgp
155  use memorymanagerextmodule, only: mem_set_value
156  ! -- dummy
157  type(GwfNpfType), pointer :: npfobj
158  character(len=*), intent(in) :: name_model
159  character(len=*), intent(in) :: input_mempath
160  integer(I4B), intent(in) :: inunit
161  integer(I4B), intent(in) :: iout
162  ! -- formats
163  character(len=*), parameter :: fmtheader = &
164  "(1x, /1x, 'NPF -- NODE PROPERTY FLOW PACKAGE, VERSION 1, 3/30/2015', &
165  &' INPUT READ FROM MEMPATH: ', A, /)"
166  !
167  ! -- Create the object
168  allocate (npfobj)
169  !
170  ! -- create name and memory path
171  call npfobj%set_names(1, name_model, 'NPF', 'NPF', input_mempath)
172  !
173  ! -- Allocate scalars
174  call npfobj%allocate_scalars()
175  !
176  ! -- Set variables
177  npfobj%inunit = inunit
178  npfobj%iout = iout
179  !
180  ! -- check if npf is enabled
181  if (inunit > 0) then
182  !
183  ! -- Print a message identifying the node property flow package.
184  write (iout, fmtheader) input_mempath
185  end if
memorymanagerextmodule::mem_set_value
Definition: MemoryManagerExt.f90:14
Here is the caller graph for this function:

◆ npf_da()

subroutine gwfnpfmodule::npf_da ( class(gwfnpftype)  this)

Definition at line 927 of file gwf-npf.f90.

928  ! -- modules
929  use memorymanagerextmodule, only: memorystore_remove
930  use simvariablesmodule, only: idm_context
931  ! -- dummy
932  class(GwfNpftype) :: this
933  !
934  ! -- Deallocate input memory
935  call memorystore_remove(this%name_model, 'NPF', idm_context)
936  !
937  ! -- TVK
938  if (this%intvk /= 0) then
939  call this%tvk%da()
940  deallocate (this%tvk)
941  end if
942  !
943  ! -- VSC
944  if (this%invsc /= 0) then
945  nullify (this%vsc)
946  end if
947  !
948  ! -- Strings
949  !
950  ! -- Scalars
951  call mem_deallocate(this%iname)
952  call mem_deallocate(this%ixt3d)
953  call mem_deallocate(this%ixt3drhs)
954  call mem_deallocate(this%satomega)
955  call mem_deallocate(this%hnoflo)
956  call mem_deallocate(this%hdry)
957  call mem_deallocate(this%icellavg)
958  call mem_deallocate(this%iavgkeff)
959  call mem_deallocate(this%ik22)
960  call mem_deallocate(this%ik33)
961  call mem_deallocate(this%iperched)
962  call mem_deallocate(this%ivarcv)
963  call mem_deallocate(this%idewatcv)
964  call mem_deallocate(this%ithickstrt)
965  call mem_deallocate(this%isavspdis)
966  call mem_deallocate(this%isavsat)
967  call mem_deallocate(this%icalcspdis)
968  call mem_deallocate(this%irewet)
969  call mem_deallocate(this%wetfct)
970  call mem_deallocate(this%iwetit)
971  call mem_deallocate(this%ihdwet)
972  call mem_deallocate(this%ibotnode)
973  call mem_deallocate(this%iwetdry)
974  call mem_deallocate(this%iangle1)
975  call mem_deallocate(this%iangle2)
976  call mem_deallocate(this%iangle3)
977  call mem_deallocate(this%nedges)
978  call mem_deallocate(this%lastedge)
979  call mem_deallocate(this%ik22overk)
980  call mem_deallocate(this%ik33overk)
981  call mem_deallocate(this%intvk)
982  call mem_deallocate(this%invsc)
983  call mem_deallocate(this%kchangeper)
984  call mem_deallocate(this%kchangestp)
985  !
986  ! -- Deallocate arrays
987  deallocate (this%aname)
988  call mem_deallocate(this%ithickstartflag)
989  call mem_deallocate(this%icelltype)
990  call mem_deallocate(this%k11)
991  call mem_deallocate(this%k22)
992  call mem_deallocate(this%k33)
993  call mem_deallocate(this%k11input)
994  call mem_deallocate(this%k22input)
995  call mem_deallocate(this%k33input)
996  call mem_deallocate(this%sat, 'SAT', this%memoryPath)
997  call mem_deallocate(this%condsat)
998  call mem_deallocate(this%wetdry)
999  call mem_deallocate(this%angle1)
1000  call mem_deallocate(this%angle2)
1001  call mem_deallocate(this%angle3)
1002  call mem_deallocate(this%nodedge)
1003  call mem_deallocate(this%ihcedge)
1004  call mem_deallocate(this%propsedge)
1005  call mem_deallocate(this%spdis, 'SPDIS', this%memoryPath)
1006  call mem_deallocate(this%nodekchange)
1007  !
1008  ! -- deallocate parent
1009  call this%NumericalPackageType%da()
1010 
1011  ! pointers
1012  this%hnew => null()
1013 
memorymanagerextmodule::memorystore_remove
subroutine, public memorystore_remove(component, subcomponent, context)
Definition: MemoryManagerExt.f90:27
simvariablesmodule
This module contains simulation variables.
Definition: SimVariables.f90:9
simvariablesmodule::idm_context
character(len=linelength) idm_context
Definition: SimVariables.f90:18
Here is the call graph for this function:

◆ npf_df()

subroutine gwfnpfmodule::npf_df ( class(gwfnpftype)  this,
class(disbasetype), intent(inout), pointer  dis,
type(xt3dtype), pointer  xt3d,
integer(i4b), intent(in)  ingnc,
integer(i4b), intent(in)  invsc,
type(gwfnpfoptionstype), intent(in), optional  npf_options 
)

This is a hybrid routine: it either reads the options for this package from the input file, or the optional argument

Parameters
npf_optionsshould be passed. A consistency check is performed, and finally xt3d_df is called, when enabled.
thisinstance of the NPF package
[in,out]disthe pointer to the discretization
xt3dthe pointer to the XT3D 'package'
[in]ingncghostnodes enabled? (>0 means yes)
[in]invscviscosity enabled? (>0 means yes)
[in]npf_optionsthe optional options, for when not constructing from file

Definition at line 195 of file gwf-npf.f90.

196  ! -- modules
197  use simmodule, only: store_error
198  use xt3dmodule, only: xt3d_cr
199  ! -- dummy
200  class(GwfNpftype) :: this !< instance of the NPF package
201  class(DisBaseType), pointer, intent(inout) :: dis !< the pointer to the discretization
202  type(Xt3dType), pointer :: xt3d !< the pointer to the XT3D 'package'
203  integer(I4B), intent(in) :: ingnc !< ghostnodes enabled? (>0 means yes)
204  integer(I4B), intent(in) :: invsc !< viscosity enabled? (>0 means yes)
205  type(GwfNpfOptionsType), optional, intent(in) :: npf_options !< the optional options, for when not constructing from file
206  !
207  ! -- Set a pointer to dis
208  this%dis => dis
209  !
210  ! -- Set flag signifying whether vsc is active
211  if (invsc > 0) this%invsc = invsc
212  !
213  if (.not. present(npf_options)) then
214  !
215  ! -- source options
216  call this%source_options()
217  !
218  ! -- allocate arrays
219  call this%allocate_arrays(this%dis%nodes, this%dis%njas)
220  !
221  ! -- source griddata, set, and convert/check the input
222  call this%source_griddata()
223  call this%prepcheck()
224  else
225  call this%set_options(npf_options)
226  !
227  ! -- allocate arrays
228  call this%allocate_arrays(this%dis%nodes, this%dis%njas)
229  end if
230  !
231  call this%check_options()
232  !
233  ! -- Save pointer to xt3d object
234  this%xt3d => xt3d
235  if (this%ixt3d /= 0) xt3d%ixt3d = this%ixt3d
236  call this%xt3d%xt3d_df(dis)
237  !
238  ! -- Ensure GNC and XT3D are not both on at the same time
239  if (this%ixt3d /= 0 .and. ingnc > 0) then
240  call store_error('Error in model '//trim(this%name_model)// &
241  '. The XT3D option cannot be used with the GNC &
242  &Package.', terminate=.true.)
243  end if
xt3dmodule::xt3d_cr
subroutine, public xt3d_cr(xt3dobj, name_model, inunit, iout, ldispopt)
Create a new xt3d object.
Definition: Xt3dInterface.f90:90
Here is the call graph for this function:

◆ npf_fc()

subroutine gwfnpfmodule::npf_fc ( class(gwfnpftype)  this,
integer(i4b)  kiter,
class(matrixbasetype), pointer  matrix_sln,
integer(i4b), dimension(:), intent(in)  idxglo,
real(dp), dimension(:), intent(inout)  rhs,
real(dp), dimension(:), intent(inout)  hnew 
)
private

Definition at line 456 of file gwf-npf.f90.

457  ! -- modules
458  use constantsmodule, only: done
459  ! -- dummy
460  class(GwfNpfType) :: this
461  integer(I4B) :: kiter
462  class(MatrixBaseType), pointer :: matrix_sln
463  integer(I4B), intent(in), dimension(:) :: idxglo
464  real(DP), intent(inout), dimension(:) :: rhs
465  real(DP), intent(inout), dimension(:) :: hnew
466  ! -- local
467  integer(I4B) :: n, m, ii, idiag, ihc
468  integer(I4B) :: isymcon, idiagm
469  real(DP) :: hyn, hym
470  real(DP) :: cond
471  !
472  ! -- Calculate conductance and put into amat
473  !
474  if (this%ixt3d /= 0) then
475  call this%xt3d%xt3d_fc(kiter, matrix_sln, idxglo, rhs, hnew)
476  else
477  do n = 1, this%dis%nodes
478  do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
479  if (this%dis%con%mask(ii) == 0) cycle
480 
481  m = this%dis%con%ja(ii)
482  !
483  ! -- Calculate conductance only for upper triangle but insert into
484  ! upper and lower parts of amat.
485  if (m < n) cycle
486  ihc = this%dis%con%ihc(this%dis%con%jas(ii))
487  hyn = this%hy_eff(n, m, ihc, ipos=ii)
488  hym = this%hy_eff(m, n, ihc, ipos=ii)
489  !
490  ! -- Vertical connection
491  if (ihc == c3d_vertical) then
492  !
493  ! -- Calculate vertical conductance
494  cond = vcond(this%ibound(n), this%ibound(m), &
495  this%icelltype(n), this%icelltype(m), this%inewton, &
496  this%ivarcv, this%idewatcv, &
497  this%condsat(this%dis%con%jas(ii)), hnew(n), hnew(m), &
498  hyn, hym, &
499  this%sat(n), this%sat(m), &
500  this%dis%top(n), this%dis%top(m), &
501  this%dis%bot(n), this%dis%bot(m), &
502  this%dis%con%hwva(this%dis%con%jas(ii)))
503  !
504  ! -- Vertical flow for perched conditions
505  if (this%iperched /= 0) then
506  if (this%icelltype(m) /= 0) then
507  if (hnew(m) < this%dis%top(m)) then
508  !
509  ! -- Fill row n
510  idiag = this%dis%con%ia(n)
511  rhs(n) = rhs(n) - cond * this%dis%bot(n)
512  call matrix_sln%add_value_pos(idxglo(idiag), -cond)
513  !
514  ! -- Fill row m
515  isymcon = this%dis%con%isym(ii)
516  call matrix_sln%add_value_pos(idxglo(isymcon), cond)
517  rhs(m) = rhs(m) + cond * this%dis%bot(n)
518  !
519  ! -- cycle the connection loop
520  cycle
521  end if
522  end if
523  end if
524  !
525  else
526  !
527  ! -- Horizontal conductance
528  cond = hcond(this%ibound(n), this%ibound(m), &
529  this%icelltype(n), this%icelltype(m), &
530  this%inewton, &
531  this%dis%con%ihc(this%dis%con%jas(ii)), &
532  this%icellavg, &
533  this%condsat(this%dis%con%jas(ii)), &
534  hnew(n), hnew(m), this%sat(n), this%sat(m), hyn, hym, &
535  this%dis%top(n), this%dis%top(m), &
536  this%dis%bot(n), this%dis%bot(m), &
537  this%dis%con%cl1(this%dis%con%jas(ii)), &
538  this%dis%con%cl2(this%dis%con%jas(ii)), &
539  this%dis%con%hwva(this%dis%con%jas(ii)))
540  end if
541  !
542  ! -- Fill row n
543  idiag = this%dis%con%ia(n)
544  call matrix_sln%add_value_pos(idxglo(ii), cond)
545  call matrix_sln%add_value_pos(idxglo(idiag), -cond)
546  !
547  ! -- Fill row m
548  isymcon = this%dis%con%isym(ii)
549  idiagm = this%dis%con%ia(m)
550  call matrix_sln%add_value_pos(idxglo(isymcon), cond)
551  call matrix_sln%add_value_pos(idxglo(idiagm), -cond)
552  end do
553  end do
554  !
555  end if
constantsmodule::done
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
Here is the call graph for this function:

◆ npf_fn()

subroutine gwfnpfmodule::npf_fn ( class(gwfnpftype)  this,
integer(i4b)  kiter,
class(matrixbasetype), pointer  matrix_sln,
integer(i4b), dimension(:), intent(in)  idxglo,
real(dp), dimension(:), intent(inout)  rhs,
real(dp), dimension(:), intent(inout)  hnew 
)

Definition at line 560 of file gwf-npf.f90.

561  ! -- dummy
562  class(GwfNpfType) :: this
563  integer(I4B) :: kiter
564  class(MatrixBaseType), pointer :: matrix_sln
565  integer(I4B), intent(in), dimension(:) :: idxglo
566  real(DP), intent(inout), dimension(:) :: rhs
567  real(DP), intent(inout), dimension(:) :: hnew
568  ! -- local
569  integer(I4B) :: nodes, nja
570  integer(I4B) :: n, m, ii, idiag
571  integer(I4B) :: isymcon, idiagm
572  integer(I4B) :: iups
573  integer(I4B) :: idn
574  real(DP) :: cond
575  real(DP) :: consterm
576  real(DP) :: filledterm
577  real(DP) :: derv
578  real(DP) :: hds
579  real(DP) :: term
580  real(DP) :: topup
581  real(DP) :: botup
582  !
583  ! -- add newton terms to solution matrix
584  nodes = this%dis%nodes
585  nja = this%dis%con%nja
586  if (this%ixt3d /= 0) then
587  call this%xt3d%xt3d_fn(kiter, nodes, nja, matrix_sln, idxglo, rhs, hnew)
588  else
589  !
590  do n = 1, nodes
591  idiag = this%dis%con%ia(n)
592  do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
593  if (this%dis%con%mask(ii) == 0) cycle
594 
595  m = this%dis%con%ja(ii)
596  isymcon = this%dis%con%isym(ii)
597  ! work on upper triangle
598  if (m < n) cycle
599  if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0 .and. &
600  this%ivarcv == 0) then
601  !call this%vcond(n,m,hnew(n),hnew(m),ii,cond)
602  ! do nothing
603  else
604  ! determine upstream node
605  iups = m
606  if (hnew(m) < hnew(n)) iups = n
607  idn = n
608  if (iups == n) idn = m
609  !
610  ! -- no newton terms if upstream cell is confined
611  if (this%icelltype(iups) == 0) cycle
612  !
613  ! -- Set the upstream top and bot, and then recalculate for a
614  ! vertically staggered horizontal connection
615  topup = this%dis%top(iups)
616  botup = this%dis%bot(iups)
617  if (this%dis%con%ihc(this%dis%con%jas(ii)) == 2) then
618  topup = min(this%dis%top(n), this%dis%top(m))
619  botup = max(this%dis%bot(n), this%dis%bot(m))
620  end if
621  !
622  ! get saturated conductivity for derivative
623  cond = this%condsat(this%dis%con%jas(ii))
624  !
625  ! compute additional term
626  consterm = -cond * (hnew(iups) - hnew(idn)) !needs to use hwadi instead of hnew(idn)
627  !filledterm = cond
628  filledterm = matrix_sln%get_value_pos(idxglo(ii))
629  derv = squadraticsaturationderivative(topup, botup, hnew(iups), &
630  this%satomega)
631  idiagm = this%dis%con%ia(m)
632  ! fill jacobian for n being the upstream node
633  if (iups == n) then
634  hds = hnew(m)
635  !isymcon = this%dis%con%isym(ii)
636  term = consterm * derv
637  rhs(n) = rhs(n) + term * hnew(n) !+ amat(idxglo(isymcon)) * (dwadi * hds - hds) !need to add dwadi
638  rhs(m) = rhs(m) - term * hnew(n) !- amat(idxglo(isymcon)) * (dwadi * hds - hds) !need to add dwadi
639  ! fill in row of n
640  call matrix_sln%add_value_pos(idxglo(idiag), term)
641  ! fill newton term in off diagonal if active cell
642  if (this%ibound(n) > 0) then
643  filledterm = matrix_sln%get_value_pos(idxglo(ii))
644  call matrix_sln%set_value_pos(idxglo(ii), filledterm) !* dwadi !need to add dwadi
645  end if
646  !fill row of m
647  filledterm = matrix_sln%get_value_pos(idxglo(idiagm))
648  call matrix_sln%set_value_pos(idxglo(idiagm), filledterm) !- filledterm * (dwadi - DONE) !need to add dwadi
649  ! fill newton term in off diagonal if active cell
650  if (this%ibound(m) > 0) then
651  call matrix_sln%add_value_pos(idxglo(isymcon), -term)
652  end if
653  ! fill jacobian for m being the upstream node
654  else
655  hds = hnew(n)
656  term = -consterm * derv
657  rhs(n) = rhs(n) + term * hnew(m) !+ amat(idxglo(ii)) * (dwadi * hds - hds) !need to add dwadi
658  rhs(m) = rhs(m) - term * hnew(m) !- amat(idxglo(ii)) * (dwadi * hds - hds) !need to add dwadi
659  ! fill in row of n
660  filledterm = matrix_sln%get_value_pos(idxglo(idiag))
661  call matrix_sln%set_value_pos(idxglo(idiag), filledterm) !- filledterm * (dwadi - DONE) !need to add dwadi
662  ! fill newton term in off diagonal if active cell
663  if (this%ibound(n) > 0) then
664  call matrix_sln%add_value_pos(idxglo(ii), term)
665  end if
666  !fill row of m
667  call matrix_sln%add_value_pos(idxglo(idiagm), -term)
668  ! fill newton term in off diagonal if active cell
669  if (this%ibound(m) > 0) then
670  filledterm = matrix_sln%get_value_pos(idxglo(isymcon))
671  call matrix_sln%set_value_pos(idxglo(isymcon), filledterm) !* dwadi !need to add dwadi
672  end if
673  end if
674  end if
675 
676  end do
677  end do
678  !
679  end if
Here is the call graph for this function:

◆ npf_mc()

subroutine gwfnpfmodule::npf_mc ( class(gwfnpftype)  this,
integer(i4b), intent(in)  moffset,
class(matrixbasetype), pointer  matrix_sln 
)

Definition at line 262 of file gwf-npf.f90.

263  ! -- dummy
264  class(GwfNpftype) :: this
265  integer(I4B), intent(in) :: moffset
266  class(MatrixBaseType), pointer :: matrix_sln
267  !
268  if (this%ixt3d /= 0) call this%xt3d%xt3d_mc(moffset, matrix_sln)

◆ npf_nur()

subroutine gwfnpfmodule::npf_nur ( class(gwfnpftype)  this,
integer(i4b), intent(in)  neqmod,
real(dp), dimension(neqmod), intent(inout)  x,
real(dp), dimension(neqmod), intent(in)  xtemp,
real(dp), dimension(neqmod), intent(inout)  dx,
integer(i4b), intent(inout)  inewtonur,
real(dp), intent(inout)  dxmax,
integer(i4b), intent(inout)  locmax 
)
private

Under-relaxation of Groundwater Flow Model Heads for current outer iteration using the cell bottoms at the bottom of the model

Definition at line 687 of file gwf-npf.f90.

688  ! -- dummy
689  class(GwfNpfType) :: this
690  integer(I4B), intent(in) :: neqmod
691  real(DP), dimension(neqmod), intent(inout) :: x
692  real(DP), dimension(neqmod), intent(in) :: xtemp
693  real(DP), dimension(neqmod), intent(inout) :: dx
694  integer(I4B), intent(inout) :: inewtonur
695  real(DP), intent(inout) :: dxmax
696  integer(I4B), intent(inout) :: locmax
697  ! -- local
698  integer(I4B) :: n
699  real(DP) :: botm
700  real(DP) :: xx
701  real(DP) :: dxx
702  !
703  ! -- Newton-Raphson under-relaxation
704  do n = 1, this%dis%nodes
705  if (this%ibound(n) < 1) cycle
706  if (this%icelltype(n) > 0) then
707  botm = this%dis%bot(this%ibotnode(n))
708  ! -- only apply Newton-Raphson under-relaxation if
709  ! solution head is below the bottom of the model
710  if (x(n) < botm) then
711  inewtonur = 1
712  xx = xtemp(n) * (done - dp9) + botm * dp9
713  dxx = x(n) - xx
714  if (abs(dxx) > abs(dxmax)) then
715  locmax = n
716  dxmax = dxx
717  end if
718  x(n) = xx
719  dx(n) = dzero
720  end if
721  end if
722  end do

◆ npf_print_model_flows()

subroutine gwfnpfmodule::npf_print_model_flows ( class(gwfnpftype)  this,
integer(i4b), intent(in)  ibudfl,
real(dp), dimension(:), intent(inout)  flowja 
)
private

Definition at line 888 of file gwf-npf.f90.

889  ! -- modules
890  use tdismodule, only: kper, kstp
891  use constantsmodule, only: lenbigline
892  ! -- dummy
893  class(GwfNpfType) :: this
894  integer(I4B), intent(in) :: ibudfl
895  real(DP), intent(inout), dimension(:) :: flowja
896  ! -- local
897  character(len=LENBIGLINE) :: line
898  character(len=30) :: tempstr
899  integer(I4B) :: n, ipos, m
900  real(DP) :: qnm
901  ! -- formats
902  character(len=*), parameter :: fmtiprflow = &
903  &"(/,4x,'CALCULATED INTERCELL FLOW FOR PERIOD ', i0, ' STEP ', i0)"
904  !
905  ! -- Write flowja to list file if requested
906  if (ibudfl /= 0 .and. this%iprflow > 0) then
907  write (this%iout, fmtiprflow) kper, kstp
908  do n = 1, this%dis%nodes
909  line = ''
910  call this%dis%noder_to_string(n, tempstr)
911  line = trim(tempstr)//':'
912  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
913  m = this%dis%con%ja(ipos)
914  call this%dis%noder_to_string(m, tempstr)
915  line = trim(line)//' '//trim(tempstr)
916  qnm = flowja(ipos)
917  write (tempstr, '(1pg15.6)') qnm
918  line = trim(line)//' '//trim(adjustl(tempstr))
919  end do
920  write (this%iout, '(a)') trim(line)
921  end do
922  end if
constantsmodule::lenbigline
integer(i4b), parameter lenbigline
maximum length of a big line
Definition: Constants.f90:15

◆ npf_rp()

subroutine gwfnpfmodule::npf_rp ( class(gwfnpftype)  this)

Read and prepare NPF stress period data.

Definition at line 347 of file gwf-npf.f90.

348  implicit none
349  ! -- dummy
350  class(GwfNpfType) :: this
351  !
352  ! -- TVK
353  if (this%intvk /= 0) then
354  call this%tvk%rp()
355  end if

◆ npf_save_model_flows()

subroutine gwfnpfmodule::npf_save_model_flows ( class(gwfnpftype)  this,
real(dp), dimension(:), intent(in)  flowja,
integer(i4b), intent(in)  icbcfl,
integer(i4b), intent(in)  icbcun 
)
private

Definition at line 851 of file gwf-npf.f90.

852  ! -- dummy
853  class(GwfNpfType) :: this
854  real(DP), dimension(:), intent(in) :: flowja
855  integer(I4B), intent(in) :: icbcfl
856  integer(I4B), intent(in) :: icbcun
857  ! -- local
858  integer(I4B) :: ibinun
859  !
860  ! -- Set unit number for binary output
861  if (this%ipakcb < 0) then
862  ibinun = icbcun
863  elseif (this%ipakcb == 0) then
864  ibinun = 0
865  else
866  ibinun = this%ipakcb
867  end if
868  if (icbcfl == 0) ibinun = 0
869  !
870  ! -- Write the face flows if requested
871  if (ibinun /= 0) then
872  call this%dis%record_connection_array(flowja, ibinun, this%iout)
873  end if
874  !
875  ! -- Calculate specific discharge at cell centers and write, if requested
876  if (this%isavspdis /= 0) then
877  if (ibinun /= 0) call this%sav_spdis(ibinun)
878  end if
879  !
880  ! -- Save saturation, if requested
881  if (this%isavsat /= 0) then
882  if (ibinun /= 0) call this%sav_sat(ibinun)
883  end if

◆ prepcheck()

subroutine gwfnpfmodule::prepcheck ( class(gwfnpftype)  this)

Definition at line 1561 of file gwf-npf.f90.

1562  ! -- modules
1563  use constantsmodule, only: linelength, dpio180
1564  use simmodule, only: store_error, count_errors, store_error_filename
1565  ! -- dummy
1566  class(GwfNpfType) :: this
1567  ! -- local
1568  character(len=24), dimension(:), pointer :: aname
1569  character(len=LINELENGTH) :: cellstr, errmsg
1570  integer(I4B) :: nerr, n
1571  ! -- format
1572  character(len=*), parameter :: fmtkerr = &
1573  &"(1x, 'Hydraulic property ',a,' is <= 0 for cell ',a, ' ', 1pg15.6)"
1574  character(len=*), parameter :: fmtkerr2 = &
1575  &"(1x, '... ', i0,' additional errors not shown for ',a)"
1576  !
1577  ! -- initialize
1578  aname => this%aname
1579  !
1580  ! -- check k11
1581  nerr = 0
1582  do n = 1, size(this%k11)
1583  if (this%k11(n) <= dzero) then
1584  nerr = nerr + 1
1585  if (nerr <= 20) then
1586  call this%dis%noder_to_string(n, cellstr)
1587  write (errmsg, fmtkerr) trim(adjustl(aname(2))), trim(cellstr), &
1588  this%k11(n)
1589  call store_error(errmsg)
1590  end if
1591  end if
1592  end do
1593  if (nerr > 20) then
1594  write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(2)))
1595  call store_error(errmsg)
1596  end if
1597  !
1598  ! -- check k33 because it was read
1599  if (this%ik33 /= 0) then
1600  !
1601  ! -- Check to make sure values are greater than or equal to zero
1602  nerr = 0
1603  do n = 1, size(this%k33)
1604  if (this%ik33overk /= 0) this%k33(n) = this%k33(n) * this%k11(n)
1605  if (this%k33(n) <= dzero) then
1606  nerr = nerr + 1
1607  if (nerr <= 20) then
1608  call this%dis%noder_to_string(n, cellstr)
1609  write (errmsg, fmtkerr) trim(adjustl(aname(3))), trim(cellstr), &
1610  this%k33(n)
1611  call store_error(errmsg)
1612  end if
1613  end if
1614  end do
1615  if (nerr > 20) then
1616  write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(3)))
1617  call store_error(errmsg)
1618  end if
1619  end if
1620  !
1621  ! -- check k22 because it was read
1622  if (this%ik22 /= 0) then
1623  !
1624  ! -- Check to make sure that angles are available
1625  if (this%dis%con%ianglex == 0) then
1626  write (errmsg, '(a)') 'Error. ANGLDEGX not provided in '// &
1627  'discretization file, but K22 was specified. '
1628  call store_error(errmsg)
1629  end if
1630  !
1631  ! -- Check to make sure values are greater than or equal to zero
1632  nerr = 0
1633  do n = 1, size(this%k22)
1634  if (this%ik22overk /= 0) this%k22(n) = this%k22(n) * this%k11(n)
1635  if (this%k22(n) <= dzero) then
1636  nerr = nerr + 1
1637  if (nerr <= 20) then
1638  call this%dis%noder_to_string(n, cellstr)
1639  write (errmsg, fmtkerr) trim(adjustl(aname(4))), trim(cellstr), &
1640  this%k22(n)
1641  call store_error(errmsg)
1642  end if
1643  end if
1644  end do
1645  if (nerr > 20) then
1646  write (errmsg, fmtkerr2) nerr, trim(adjustl(aname(4)))
1647  call store_error(errmsg)
1648  end if
1649  end if
1650  !
1651  ! -- check for wetdry conflicts
1652  if (this%irewet == 1) then
1653  if (this%iwetdry == 0) then
1654  write (errmsg, '(a, a, a)') 'Error in GRIDDATA block: ', &
1655  trim(adjustl(aname(5))), ' not found.'
1656  call store_error(errmsg)
1657  end if
1658  end if
1659  !
1660  ! -- Check for angle conflicts
1661  if (this%iangle1 /= 0) then
1662  do n = 1, size(this%angle1)
1663  this%angle1(n) = this%angle1(n) * dpio180
1664  end do
1665  else
1666  if (this%ixt3d /= 0) then
1667  this%iangle1 = 1
1668  write (this%iout, '(a)') 'XT3D IN USE, BUT ANGLE1 NOT SPECIFIED. '// &
1669  'SETTING ANGLE1 TO ZERO.'
1670  do n = 1, size(this%angle1)
1671  this%angle1(n) = dzero
1672  end do
1673  end if
1674  end if
1675  if (this%iangle2 /= 0) then
1676  if (this%iangle1 == 0) then
1677  write (errmsg, '(a)') 'ANGLE2 SPECIFIED BUT NOT ANGLE1. '// &
1678  'ANGLE2 REQUIRES ANGLE1. '
1679  call store_error(errmsg)
1680  end if
1681  if (this%iangle3 == 0) then
1682  write (errmsg, '(a)') 'ANGLE2 SPECIFIED BUT NOT ANGLE3. '// &
1683  'SPECIFY BOTH OR NEITHER ONE. '
1684  call store_error(errmsg)
1685  end if
1686  do n = 1, size(this%angle2)
1687  this%angle2(n) = this%angle2(n) * dpio180
1688  end do
1689  end if
1690  if (this%iangle3 /= 0) then
1691  if (this%iangle1 == 0) then
1692  write (errmsg, '(a)') 'ANGLE3 SPECIFIED BUT NOT ANGLE1. '// &
1693  'ANGLE3 REQUIRES ANGLE1. '
1694  call store_error(errmsg)
1695  end if
1696  if (this%iangle2 == 0) then
1697  write (errmsg, '(a)') 'ANGLE3 SPECIFIED BUT NOT ANGLE2. '// &
1698  'SPECIFY BOTH OR NEITHER ONE. '
1699  call store_error(errmsg)
1700  end if
1701  do n = 1, size(this%angle3)
1702  this%angle3(n) = this%angle3(n) * dpio180
1703  end do
1704  end if
1705  !
1706  ! -- terminate if data errors
1707  if (count_errors() > 0) then
1708  call store_error_filename(this%input_fname)
1709  end if
constantsmodule::dpio180
real(dp), parameter dpio180
real constant
Definition: Constants.f90:130
Here is the call graph for this function:

◆ preprocess_input()

subroutine gwfnpfmodule::preprocess_input ( class(gwfnpftype)  this)

This routine consists of the following steps:

  1. convert cells to noflow when all transmissive parameters equal zero
  2. perform initial wetting and drying
  3. initialize cell saturation
  4. calculate saturated conductance (when not xt3d)
  5. If NEWTON under-relaxation, determine lower most node
    Parameters
    thisthe instance of the NPF package

Definition at line 1722 of file gwf-npf.f90.

1723  ! -- modules
1724  use constantsmodule, only: linelength
1725  use simmodule, only: store_error, count_errors, store_error_filename
1726  ! -- dummy
1727  class(GwfNpfType) :: this !< the instance of the NPF package
1728  ! -- local
1729  integer(I4B) :: n, m, ii, nn
1730  real(DP) :: hyn, hym
1731  real(DP) :: satn, topn, botn
1732  integer(I4B) :: nextn
1733  real(DP) :: minbot, botm
1734  logical :: finished
1735  character(len=LINELENGTH) :: cellstr, errmsg
1736  ! -- format
1737  character(len=*), parameter :: fmtcnv = &
1738  "(1X,'CELL ', A, &
1739  &' ELIMINATED BECAUSE ALL HYDRAULIC CONDUCTIVITIES TO NODE ARE 0.')"
1740  character(len=*), parameter :: fmtnct = &
1741  &"(1X,'Negative cell thickness at cell ', A)"
1742  character(len=*), parameter :: fmtihbe = &
1743  &"(1X,'Initial head, bottom elevation:',1P,2G13.5)"
1744  character(len=*), parameter :: fmttebe = &
1745  &"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
1746  !
1747  do n = 1, this%dis%nodes
1748  this%ithickstartflag(n) = 0
1749  end do
1750  !
1751  ! -- Insure that each cell has at least one non-zero transmissive parameter
1752  ! Note that a cell can be deactivated even if it has a valid connection
1753  ! to another model.
1754  nodeloop: do n = 1, this%dis%nodes
1755  !
1756  ! -- Skip if already inactive
1757  if (this%ibound(n) == 0) then
1758  if (this%irewet /= 0) then
1759  if (this%wetdry(n) == dzero) cycle nodeloop
1760  else
1761  cycle nodeloop
1762  end if
1763  end if
1764  !
1765  ! -- Cycle if k11 is not zero
1766  if (this%k11(n) /= dzero) cycle nodeloop
1767  !
1768  ! -- Cycle if at least one vertical connection has non-zero k33
1769  ! for n and m
1770  do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1771  m = this%dis%con%ja(ii)
1772  if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0) then
1773  hyn = this%k11(n)
1774  if (this%ik33 /= 0) hyn = this%k33(n)
1775  if (hyn /= dzero) then
1776  hym = this%k11(m)
1777  if (this%ik33 /= 0) hym = this%k33(m)
1778  if (hym /= dzero) cycle
1779  end if
1780  end if
1781  end do
1782  !
1783  ! -- If this part of the loop is reached, then all connections have
1784  ! zero transmissivity, so convert to noflow.
1785  this%ibound(n) = 0
1786  this%hnew(n) = this%hnoflo
1787  if (this%irewet /= 0) this%wetdry(n) = dzero
1788  call this%dis%noder_to_string(n, cellstr)
1789  write (this%iout, fmtcnv) trim(adjustl(cellstr))
1790  !
1791  end do nodeloop
1792  !
1793  ! -- Preprocess cell status and heads based on initial conditions
1794  if (this%inewton == 0) then
1795  !
1796  ! -- For standard formulation (non-Newton) call wetdry routine
1797  call this%wd(0, this%hnew)
1798  else
1799  !
1800  ! -- Newton formulation, so adjust heads to be above bottom
1801  ! (Not used in present formulation because variable cv
1802  ! cannot be used with Newton)
1803  if (this%ivarcv == 1) then
1804  do n = 1, this%dis%nodes
1805  if (this%hnew(n) < this%dis%bot(n)) then
1806  this%hnew(n) = this%dis%bot(n) + dem6
1807  end if
1808  end do
1809  end if
1810  end if
1811  !
1812  ! -- If THCKSTRT is not active, then loop through icelltype and replace
1813  ! any negative values with 1.
1814  if (this%ithickstrt == 0) then
1815  do n = 1, this%dis%nodes
1816  if (this%icelltype(n) < 0) then
1817  this%icelltype(n) = 1
1818  end if
1819  end do
1820  end if
1821  !
1822  ! -- Initialize sat to zero for ibound=0 cells, unless the cell can
1823  ! rewet. Initialize sat to the saturated fraction based on strt
1824  ! if icelltype is negative and the THCKSTRT option is in effect.
1825  ! Initialize sat to 1.0 for all other cells in order to calculate
1826  ! condsat in next section.
1827  do n = 1, this%dis%nodes
1828  if (this%ibound(n) == 0) then
1829  this%sat(n) = done
1830  if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0) then
1831  this%ithickstartflag(n) = 1
1832  this%icelltype(n) = 0
1833  end if
1834  else
1835  topn = this%dis%top(n)
1836  botn = this%dis%bot(n)
1837  if (this%icelltype(n) < 0 .and. this%ithickstrt /= 0) then
1838  call this%thksat(n, this%ic%strt(n), satn)
1839  if (botn > this%ic%strt(n)) then
1840  call this%dis%noder_to_string(n, cellstr)
1841  write (errmsg, fmtnct) trim(adjustl(cellstr))
1842  call store_error(errmsg)
1843  write (errmsg, fmtihbe) this%ic%strt(n), botn
1844  call store_error(errmsg)
1845  end if
1846  this%ithickstartflag(n) = 1
1847  this%icelltype(n) = 0
1848  else
1849  satn = done
1850  if (botn > topn) then
1851  call this%dis%noder_to_string(n, cellstr)
1852  write (errmsg, fmtnct) trim(adjustl(cellstr))
1853  call store_error(errmsg)
1854  write (errmsg, fmttebe) topn, botn
1855  call store_error(errmsg)
1856  end if
1857  end if
1858  this%sat(n) = satn
1859  end if
1860  end do
1861  if (count_errors() > 0) then
1862  call store_error_filename(this%input_fname)
1863  end if
1864  !
1865  ! -- Calculate condsat, but only if xt3d is not active. If xt3d is
1866  ! active, then condsat is allocated to size of zero.
1867  if (this%ixt3d == 0) then
1868  !
1869  ! -- Calculate the saturated conductance for all connections assuming
1870  ! that saturation is 1 (except for case where icelltype was entered
1871  ! as a negative value and THCKSTRT option in effect)
1872  do n = 1, this%dis%nodes
1873  call this%calc_condsat(n, .true.)
1874  end do
1875  !
1876  end if
1877  !
1878  ! -- Determine the lower most node
1879  if (this%igwfnewtonur /= 0) then
1880  call mem_reallocate(this%ibotnode, this%dis%nodes, 'IBOTNODE', &
1881  trim(this%memoryPath))
1882  do n = 1, this%dis%nodes
1883  !
1884  minbot = this%dis%bot(n)
1885  nn = n
1886  finished = .false.
1887  do while (.not. finished)
1888  nextn = 0
1889  !
1890  ! -- Go through the connecting cells
1891  do ii = this%dis%con%ia(nn) + 1, this%dis%con%ia(nn + 1) - 1
1892  !
1893  ! -- Set the m cell number
1894  m = this%dis%con%ja(ii)
1895  botm = this%dis%bot(m)
1896  !
1897  ! -- select vertical connections: ihc == 0
1898  if (this%dis%con%ihc(this%dis%con%jas(ii)) == 0) then
1899  if (m > nn .and. botm < minbot) then
1900  nextn = m
1901  minbot = botm
1902  end if
1903  end if
1904  end do
1905  if (nextn > 0) then
1906  nn = nextn
1907  else
1908  finished = .true.
1909  end if
1910  end do
1911  this%ibotnode(n) = nn
1912  end do
1913  end if
1914  !
1915  ! -- nullify unneeded gwf pointers
1916  this%igwfnewtonur => null()
Here is the call graph for this function:

◆ rewet_check()

subroutine gwfnpfmodule::rewet_check ( class(gwfnpftype)  this,
integer(i4b), intent(in)  kiter,
integer(i4b), intent(in)  node,
real(dp), intent(in)  hm,
integer(i4b), intent(in)  ibdm,
integer(i4b), intent(in)  ihc,
real(dp), dimension(:), intent(inout)  hnew,
integer(i4b), intent(out)  irewet 
)

This method can be called from any external object that has a head that can be used to rewet the GWF cell node. The ihc value is used to determine if it is a vertical or horizontal connection, which can operate differently depending on user settings.

Definition at line 2141 of file gwf-npf.f90.

2142  ! -- dummy
2143  class(GwfNpfType) :: this
2144  integer(I4B), intent(in) :: kiter
2145  integer(I4B), intent(in) :: node
2146  real(DP), intent(in) :: hm
2147  integer(I4B), intent(in) :: ibdm
2148  integer(I4B), intent(in) :: ihc
2149  real(DP), intent(inout), dimension(:) :: hnew
2150  integer(I4B), intent(out) :: irewet
2151  ! -- local
2152  integer(I4B) :: itflg
2153  real(DP) :: wd, awd, turnon, bbot
2154  !
2155  irewet = 0
2156  !
2157  ! -- Convert a dry cell to wet if it meets the criteria
2158  if (this%irewet > 0) then
2159  itflg = mod(kiter, this%iwetit)
2160  if (itflg == 0) then
2161  if (this%ibound(node) == 0 .and. this%wetdry(node) /= dzero) then
2162  !
2163  ! -- Calculate wetting elevation
2164  bbot = this%dis%bot(node)
2165  wd = this%wetdry(node)
2166  awd = wd
2167  if (wd < 0) awd = -wd
2168  turnon = bbot + awd
2169  !
2170  ! -- Check head in adjacent cells to see if wetting elevation has
2171  ! been reached
2172  if (ihc == c3d_vertical) then
2173  !
2174  ! -- check cell below
2175  if (ibdm > 0 .and. hm >= turnon) irewet = 1
2176  else
2177  if (wd > dzero) then
2178  !
2179  ! -- check horizontally adjacent cells
2180  if (ibdm > 0 .and. hm >= turnon) irewet = 1
2181  end if
2182  end if
2183  !
2184  if (irewet == 1) then
2185  ! -- rewet cell; use equation 3a if ihdwet=0; use equation 3b if
2186  ! ihdwet is not 0.
2187  if (this%ihdwet == 0) then
2188  hnew(node) = bbot + this%wetfct * (hm - bbot)
2189  else
2190  hnew(node) = bbot + this%wetfct * awd !(hm - bbot)
2191  end if
2192  this%ibound(node) = 30000
2193  end if
2194  end if
2195  end if
2196  end if

◆ sav_sat()

subroutine gwfnpfmodule::sav_sat ( class(gwfnpftype)  this,
integer(i4b), intent(in)  ibinun 
)
private

Definition at line 2655 of file gwf-npf.f90.

2656  ! -- dummy
2657  class(GwfNpfType) :: this
2658  integer(I4B), intent(in) :: ibinun
2659  ! -- local
2660  character(len=16) :: text
2661  character(len=16), dimension(1) :: auxtxt
2662  real(DP), dimension(1) :: a
2663  integer(I4B) :: n
2664  integer(I4B) :: naux
2665  !
2666  ! -- Write the header
2667  text = ' DATA-SAT'
2668  naux = 1
2669  auxtxt(:) = [' sat']
2670  call this%dis%record_srcdst_list_header(text, this%name_model, &
2671  this%packName, this%name_model, &
2672  this%packName, naux, auxtxt, ibinun, &
2673  this%dis%nodes, this%iout)
2674  !
2675  ! -- Write a zero for Q, and then write saturation as an aux variables
2676  do n = 1, this%dis%nodes
2677  a(1) = this%sat(n)
2678  call this%dis%record_mf6_list_entry(ibinun, n, n, dzero, naux, a)
2679  end do

◆ sav_spdis()

subroutine gwfnpfmodule::sav_spdis ( class(gwfnpftype)  this,
integer(i4b), intent(in)  ibinun 
)

Definition at line 2627 of file gwf-npf.f90.

2628  ! -- dummy
2629  class(GwfNpfType) :: this
2630  integer(I4B), intent(in) :: ibinun
2631  ! -- local
2632  character(len=16) :: text
2633  character(len=16), dimension(3) :: auxtxt
2634  integer(I4B) :: n
2635  integer(I4B) :: naux
2636  !
2637  ! -- Write the header
2638  text = ' DATA-SPDIS'
2639  naux = 3
2640  auxtxt(:) = [' qx', ' qy', ' qz']
2641  call this%dis%record_srcdst_list_header(text, this%name_model, &
2642  this%packName, this%name_model, &
2643  this%packName, naux, auxtxt, ibinun, &
2644  this%dis%nodes, this%iout)
2645  !
2646  ! -- Write a zero for Q, and then write qx, qy, qz as aux variables
2647  do n = 1, this%dis%nodes
2648  call this%dis%record_mf6_list_entry(ibinun, n, n, dzero, naux, &
2649  this%spdis(:, n))
2650  end do

◆ set_edge_properties()

subroutine gwfnpfmodule::set_edge_properties ( class(gwfnpftype)  this,
integer(i4b), intent(in)  nodedge,
integer(i4b), intent(in)  ihcedge,
real(dp), intent(in)  q,
real(dp), intent(in)  area,
real(dp), intent(in)  nx,
real(dp), intent(in)  ny,
real(dp), intent(in)  distance 
)
private

Definition at line 2697 of file gwf-npf.f90.

2699  ! -- dummy
2700  class(GwfNpfType) :: this
2701  integer(I4B), intent(in) :: nodedge
2702  integer(I4B), intent(in) :: ihcedge
2703  real(DP), intent(in) :: q
2704  real(DP), intent(in) :: area
2705  real(DP), intent(in) :: nx
2706  real(DP), intent(in) :: ny
2707  real(DP), intent(in) :: distance
2708  ! -- local
2709  integer(I4B) :: lastedge
2710  !
2711  this%lastedge = this%lastedge + 1
2712  lastedge = this%lastedge
2713  this%nodedge(lastedge) = nodedge
2714  this%ihcedge(lastedge) = ihcedge
2715  this%propsedge(1, lastedge) = q
2716  this%propsedge(2, lastedge) = area
2717  this%propsedge(3, lastedge) = nx
2718  this%propsedge(4, lastedge) = ny
2719  this%propsedge(5, lastedge) = distance
2720  !
2721  ! -- If this is the last edge, then the next call must be starting a new
2722  ! edge properties assignment loop, so need to reset lastedge to 0
2723  if (this%lastedge == this%nedges) this%lastedge = 0

◆ set_options()

subroutine gwfnpfmodule::set_options ( class(gwfnpftype)  this,
type(gwfnpfoptionstype), intent(in)  options 
)

Definition at line 1342 of file gwf-npf.f90.

1343  ! -- dummy
1344  class(GwfNpftype) :: this
1345  type(GwfNpfOptionsType), intent(in) :: options
1346  !
1347  this%icellavg = options%icellavg
1348  this%ithickstrt = options%ithickstrt
1349  this%iperched = options%iperched
1350  this%ivarcv = options%ivarcv
1351  this%idewatcv = options%idewatcv
1352  this%irewet = options%irewet
1353  this%wetfct = options%wetfct
1354  this%iwetit = options%iwetit
1355  this%ihdwet = options%ihdwet

◆ sgwf_npf_qcalc()

subroutine gwfnpfmodule::sgwf_npf_qcalc ( class(gwfnpftype)  this,
integer(i4b), intent(in)  n,
integer(i4b), intent(in)  m,
real(dp), intent(in)  hn,
real(dp), intent(in)  hm,
integer(i4b), intent(in)  icon,
real(dp), intent(inout)  qnm 
)
private

Definition at line 780 of file gwf-npf.f90.

781  ! -- dummy
782  class(GwfNpfType) :: this
783  integer(I4B), intent(in) :: n
784  integer(I4B), intent(in) :: m
785  real(DP), intent(in) :: hn
786  real(DP), intent(in) :: hm
787  integer(I4B), intent(in) :: icon
788  real(DP), intent(inout) :: qnm
789  ! -- local
790  real(DP) :: hyn, hym
791  real(DP) :: condnm
792  real(DP) :: hntemp, hmtemp
793  integer(I4B) :: ihc
794  !
795  ! -- Initialize
796  ihc = this%dis%con%ihc(this%dis%con%jas(icon))
797  hyn = this%hy_eff(n, m, ihc, ipos=icon)
798  hym = this%hy_eff(m, n, ihc, ipos=icon)
799  !
800  ! -- Calculate conductance
801  if (ihc == c3d_vertical) then
802  condnm = vcond(this%ibound(n), this%ibound(m), &
803  this%icelltype(n), this%icelltype(m), this%inewton, &
804  this%ivarcv, this%idewatcv, &
805  this%condsat(this%dis%con%jas(icon)), hn, hm, &
806  hyn, hym, &
807  this%sat(n), this%sat(m), &
808  this%dis%top(n), this%dis%top(m), &
809  this%dis%bot(n), this%dis%bot(m), &
810  this%dis%con%hwva(this%dis%con%jas(icon)))
811  else
812  condnm = hcond(this%ibound(n), this%ibound(m), &
813  this%icelltype(n), this%icelltype(m), &
814  this%inewton, &
815  this%dis%con%ihc(this%dis%con%jas(icon)), &
816  this%icellavg, &
817  this%condsat(this%dis%con%jas(icon)), &
818  hn, hm, this%sat(n), this%sat(m), hyn, hym, &
819  this%dis%top(n), this%dis%top(m), &
820  this%dis%bot(n), this%dis%bot(m), &
821  this%dis%con%cl1(this%dis%con%jas(icon)), &
822  this%dis%con%cl2(this%dis%con%jas(icon)), &
823  this%dis%con%hwva(this%dis%con%jas(icon)))
824  end if
825  !
826  ! -- Initialize hntemp and hmtemp
827  hntemp = hn
828  hmtemp = hm
829  !
830  ! -- Check and adjust for dewatered conditions
831  if (this%iperched /= 0) then
832  if (this%dis%con%ihc(this%dis%con%jas(icon)) == 0) then
833  if (n > m) then
834  if (this%icelltype(n) /= 0) then
835  if (hn < this%dis%top(n)) hntemp = this%dis%bot(m)
836  end if
837  else
838  if (this%icelltype(m) /= 0) then
839  if (hm < this%dis%top(m)) hmtemp = this%dis%bot(n)
840  end if
841  end if
842  end if
843  end if
844  !
845  ! -- Calculate flow positive into cell n
846  qnm = condnm * (hmtemp - hntemp)
Here is the call graph for this function:

◆ sgwf_npf_thksat()

subroutine gwfnpfmodule::sgwf_npf_thksat ( class(gwfnpftype)  this,
integer(i4b), intent(in)  n,
real(dp), intent(in)  hn,
real(dp), intent(inout)  thksat 
)
private

Definition at line 757 of file gwf-npf.f90.

758  ! -- dummy
759  class(GwfNpfType) :: this
760  integer(I4B), intent(in) :: n
761  real(DP), intent(in) :: hn
762  real(DP), intent(inout) :: thksat
763  !
764  ! -- Standard Formulation
765  if (hn >= this%dis%top(n)) then
766  thksat = done
767  else
768  thksat = (hn - this%dis%bot(n)) / (this%dis%top(n) - this%dis%bot(n))
769  end if
770  !
771  ! -- Newton-Raphson Formulation
772  if (this%inewton /= 0) then
773  thksat = squadraticsaturation(this%dis%top(n), this%dis%bot(n), hn, &
774  this%satomega)
775  end if
Here is the call graph for this function:

◆ sgwf_npf_wdmsg()

subroutine gwfnpfmodule::sgwf_npf_wdmsg ( class(gwfnpftype)  this,
integer(i4b), intent(in)  icode,
integer(i4b), intent(inout)  ncnvrt,
character(len=30), dimension(5), intent(inout)  nodcnvrt,
character(len=3), dimension(5), intent(inout)  acnvrt,
integer(i4b), intent(inout)  ihdcnv,
integer(i4b), intent(in)  kiter,
integer(i4b), intent(in)  n 
)
private

Definition at line 2201 of file gwf-npf.f90.

2203  ! -- modules
2204  use tdismodule, only: kstp, kper
2205  ! -- dummy
2206  class(GwfNpfType) :: this
2207  integer(I4B), intent(in) :: icode
2208  integer(I4B), intent(inout) :: ncnvrt
2209  character(len=30), dimension(5), intent(inout) :: nodcnvrt
2210  character(len=3), dimension(5), intent(inout) :: acnvrt
2211  integer(I4B), intent(inout) :: ihdcnv
2212  integer(I4B), intent(in) :: kiter
2213  integer(I4B), intent(in) :: n
2214  ! -- local
2215  integer(I4B) :: l
2216  ! -- formats
2217  character(len=*), parameter :: fmtcnvtn = &
2218  "(1X,/1X,'CELL CONVERSIONS FOR ITER.=',I0, &
2219  &' STEP=',I0,' PERIOD=',I0,' (NODE or LRC)')"
2220  character(len=*), parameter :: fmtnode = "(1X,3X,5(A4, A20))"
2221  !
2222  ! -- Keep track of cell conversions
2223  if (icode > 0) then
2224  ncnvrt = ncnvrt + 1
2225  call this%dis%noder_to_string(n, nodcnvrt(ncnvrt))
2226  if (icode == 1) then
2227  acnvrt(ncnvrt) = 'DRY'
2228  else
2229  acnvrt(ncnvrt) = 'WET'
2230  end if
2231  end if
2232  !
2233  ! -- Print a line if 5 conversions have occurred or if icode indicates that a
2234  ! partial line should be printed
2235  if (ncnvrt == 5 .or. (icode == 0 .and. ncnvrt > 0)) then
2236  if (ihdcnv == 0) write (this%iout, fmtcnvtn) kiter, kstp, kper
2237  ihdcnv = 1
2238  write (this%iout, fmtnode) &
2239  (acnvrt(l), trim(adjustl(nodcnvrt(l))), l=1, ncnvrt)
2240  ncnvrt = 0
2241  end if

◆ sgwf_npf_wetdry()

subroutine gwfnpfmodule::sgwf_npf_wetdry ( class(gwfnpftype)  this,
integer(i4b), intent(in)  kiter,
real(dp), dimension(:), intent(inout)  hnew 
)
private

Definition at line 2035 of file gwf-npf.f90.

2036  ! -- modules
2037  use tdismodule, only: kstp, kper
2038  use simmodule, only: store_error, store_error_filename
2039  use constantsmodule, only: linelength
2040  ! -- dummy
2041  class(GwfNpfType) :: this
2042  integer(I4B), intent(in) :: kiter
2043  real(DP), intent(inout), dimension(:) :: hnew
2044  ! -- local
2045  integer(I4B) :: n, m, ii, ihc
2046  real(DP) :: ttop, bbot, thick
2047  integer(I4B) :: ncnvrt, ihdcnv
2048  character(len=30), dimension(5) :: nodcnvrt
2049  character(len=30) :: nodestr
2050  character(len=3), dimension(5) :: acnvrt
2051  character(len=LINELENGTH) :: errmsg
2052  integer(I4B) :: irewet
2053  ! -- formats
2054  character(len=*), parameter :: fmtnct = &
2055  "(1X,/1X,'Negative cell thickness at (layer,row,col)', &
2056  &I4,',',I5,',',I5)"
2057  character(len=*), parameter :: fmttopbot = &
2058  &"(1X,'Top elevation, bottom elevation:',1P,2G13.5)"
2059  character(len=*), parameter :: fmttopbotthk = &
2060  &"(1X,'Top elevation, bottom elevation, thickness:',1P,3G13.5)"
2061  character(len=*), parameter :: fmtdrychd = &
2062  &"(1X,/1X,'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED')"
2063  character(len=*), parameter :: fmtni = &
2064  &"(1X,'CELLID=',a,' ITERATION=',I0,' TIME STEP=',I0,' STRESS PERIOD=',I0)"
2065  !
2066  ! -- Initialize
2067  ncnvrt = 0
2068  ihdcnv = 0
2069  !
2070  ! -- Convert dry cells to wet
2071  do n = 1, this%dis%nodes
2072  do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
2073  m = this%dis%con%ja(ii)
2074  ihc = this%dis%con%ihc(this%dis%con%jas(ii))
2075  call this%rewet_check(kiter, n, hnew(m), this%ibound(m), ihc, hnew, &
2076  irewet)
2077  if (irewet == 1) then
2078  call this%wdmsg(2, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2079  end if
2080  end do
2081  end do
2082  !
2083  ! -- Perform drying
2084  do n = 1, this%dis%nodes
2085  !
2086  ! -- cycle if inactive or confined
2087  if (this%ibound(n) == 0) cycle
2088  if (this%icelltype(n) == 0) cycle
2089  !
2090  ! -- check for negative cell thickness
2091  bbot = this%dis%bot(n)
2092  ttop = this%dis%top(n)
2093  if (bbot > ttop) then
2094  write (errmsg, fmtnct) n
2095  call store_error(errmsg)
2096  write (errmsg, fmttopbot) ttop, bbot
2097  call store_error(errmsg)
2098  call store_error_filename(this%input_fname)
2099  end if
2100  !
2101  ! -- Calculate saturated thickness
2102  if (this%icelltype(n) /= 0) then
2103  if (hnew(n) < ttop) ttop = hnew(n)
2104  end if
2105  thick = ttop - bbot
2106  !
2107  ! -- If thick<0 print message, set hnew, and ibound
2108  if (thick <= dzero) then
2109  call this%wdmsg(1, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2110  hnew(n) = this%hdry
2111  if (this%ibound(n) < 0) then
2112  errmsg = 'CONSTANT-HEAD CELL WENT DRY -- SIMULATION ABORTED'
2113  call store_error(errmsg)
2114  write (errmsg, fmttopbotthk) ttop, bbot, thick
2115  call store_error(errmsg)
2116  call this%dis%noder_to_string(n, nodestr)
2117  write (errmsg, fmtni) trim(adjustl(nodestr)), kiter, kstp, kper
2118  call store_error(errmsg)
2119  call store_error_filename(this%input_fname)
2120  end if
2121  this%ibound(n) = 0
2122  end if
2123  end do
2124  !
2125  ! -- Print remaining cell conversions
2126  call this%wdmsg(0, ncnvrt, nodcnvrt, acnvrt, ihdcnv, kiter, n)
2127  !
2128  ! -- Change ibound from 30000 to 1
2129  do n = 1, this%dis%nodes
2130  if (this%ibound(n) == 30000) this%ibound(n) = 1
2131  end do
Here is the call graph for this function:

◆ source_griddata()

subroutine gwfnpfmodule::source_griddata ( class(gwfnpftype)  this)

Definition at line 1472 of file gwf-npf.f90.

1473  ! -- modules
1474  use simmodule, only: count_errors, store_error
1475  use memorymanagermodule, only: mem_reallocate
1476  use memorymanagerextmodule, only: mem_set_value
1477  use gwfnpfinputmodule, only: gwfnpfparamfoundtype
1478  ! -- dummy
1479  class(GwfNpftype) :: this
1480  ! -- locals
1481  character(len=LINELENGTH) :: errmsg
1482  type(GwfNpfParamFoundType) :: found
1483  logical, dimension(2) :: afound
1484  integer(I4B), dimension(:), pointer, contiguous :: map
1485  !
1486  ! -- set map to convert user input data into reduced data
1487  map => null()
1488  if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
1489  !
1490  ! -- update defaults with idm sourced values
1491  call mem_set_value(this%icelltype, 'ICELLTYPE', this%input_mempath, map, &
1492  found%icelltype)
1493  call mem_set_value(this%k11, 'K', this%input_mempath, map, found%k)
1494  call mem_set_value(this%k33, 'K33', this%input_mempath, map, found%k33)
1495  call mem_set_value(this%k22, 'K22', this%input_mempath, map, found%k22)
1496  call mem_set_value(this%wetdry, 'WETDRY', this%input_mempath, map, &
1497  found%wetdry)
1498  call mem_set_value(this%angle1, 'ANGLE1', this%input_mempath, map, &
1499  found%angle1)
1500  call mem_set_value(this%angle2, 'ANGLE2', this%input_mempath, map, &
1501  found%angle2)
1502  call mem_set_value(this%angle3, 'ANGLE3', this%input_mempath, map, &
1503  found%angle3)
1504  !
1505  ! -- ensure ICELLTYPE was found
1506  if (.not. found%icelltype) then
1507  write (errmsg, '(a)') 'Error in GRIDDATA block: ICELLTYPE not found.'
1508  call store_error(errmsg)
1509  end if
1510  !
1511  ! -- ensure K was found
1512  if (.not. found%k) then
1513  write (errmsg, '(a)') 'Error in GRIDDATA block: K not found.'
1514  call store_error(errmsg)
1515  end if
1516  !
1517  ! -- set error if ik33overk set with no k33
1518  if (.not. found%k33 .and. this%ik33overk /= 0) then
1519  write (errmsg, '(a)') 'K33OVERK option specified but K33 not specified.'
1520  call store_error(errmsg)
1521  end if
1522  !
1523  ! -- set error if ik22overk set with no k22
1524  if (.not. found%k22 .and. this%ik22overk /= 0) then
1525  write (errmsg, '(a)') 'K22OVERK option specified but K22 not specified.'
1526  call store_error(errmsg)
1527  end if
1528  !
1529  ! -- handle found side effects
1530  if (found%k33) this%ik33 = 1
1531  if (found%k22) this%ik22 = 1
1532  if (found%wetdry) this%iwetdry = 1
1533  if (found%angle1) this%iangle1 = 1
1534  if (found%angle2) this%iangle2 = 1
1535  if (found%angle3) this%iangle3 = 1
1536  !
1537  ! -- handle not found side effects
1538  if (.not. found%k33) then
1539  call mem_set_value(this%k33, 'K', this%input_mempath, map, afound(1))
1540  end if
1541  if (.not. found%k22) then
1542  call mem_set_value(this%k22, 'K', this%input_mempath, map, afound(2))
1543  end if
1544  if (.not. found%wetdry) call mem_reallocate(this%wetdry, 1, 'WETDRY', &
1545  trim(this%memoryPath))
1546  if (.not. found%angle1 .and. this%ixt3d == 0) &
1547  call mem_reallocate(this%angle1, 0, 'ANGLE1', trim(this%memoryPath))
1548  if (.not. found%angle2 .and. this%ixt3d == 0) &
1549  call mem_reallocate(this%angle2, 0, 'ANGLE2', trim(this%memoryPath))
1550  if (.not. found%angle3 .and. this%ixt3d == 0) &
1551  call mem_reallocate(this%angle3, 0, 'ANGLE3', trim(this%memoryPath))
1552  !
1553  ! -- log griddata
1554  if (this%iout > 0) then
1555  call this%log_griddata(found)
1556  end if
Here is the call graph for this function:

◆ source_options()

subroutine gwfnpfmodule::source_options ( class(gwfnpftype)  this)

Definition at line 1266 of file gwf-npf.f90.

1267  ! -- modules
1268  use simmodule, only: store_error, store_error_filename
1269  use memorymanagermodule, only: mem_setptr, get_isize
1270  use memorymanagerextmodule, only: mem_set_value
1271  use characterstringmodule, only: characterstringtype
1272  use gwfnpfinputmodule, only: gwfnpfparamfoundtype
1273  use sourcecommonmodule, only: filein_fname
1274  ! -- dummy
1275  class(GwfNpftype) :: this
1276  ! -- locals
1277  character(len=LENVARNAME), dimension(3) :: cellavg_method = &
1278  &[character(len=LENVARNAME) :: 'LOGARITHMIC', 'AMT-LMK', 'AMT-HMK']
1279  type(gwfnpfparamfoundtype) :: found
1280  character(len=LINELENGTH) :: tvk6_filename
1281  !
1282  ! -- update defaults with idm sourced values
1283  call mem_set_value(this%iprflow, 'IPRFLOW', this%input_mempath, found%iprflow)
1284  call mem_set_value(this%ipakcb, 'IPAKCB', this%input_mempath, found%ipakcb)
1285  call mem_set_value(this%icellavg, 'CELLAVG', this%input_mempath, &
1286  cellavg_method, found%cellavg)
1287  call mem_set_value(this%ithickstrt, 'ITHICKSTRT', this%input_mempath, &
1288  found%ithickstrt)
1289  call mem_set_value(this%iperched, 'IPERCHED', this%input_mempath, &
1290  found%iperched)
1291  call mem_set_value(this%ivarcv, 'IVARCV', this%input_mempath, found%ivarcv)
1292  call mem_set_value(this%idewatcv, 'IDEWATCV', this%input_mempath, &
1293  found%idewatcv)
1294  call mem_set_value(this%ixt3d, 'IXT3D', this%input_mempath, found%ixt3d)
1295  call mem_set_value(this%ixt3drhs, 'IXT3DRHS', this%input_mempath, &
1296  found%ixt3drhs)
1297  call mem_set_value(this%isavspdis, 'ISAVSPDIS', this%input_mempath, &
1298  found%isavspdis)
1299  call mem_set_value(this%isavsat, 'ISAVSAT', this%input_mempath, found%isavsat)
1300  call mem_set_value(this%ik22overk, 'IK22OVERK', this%input_mempath, &
1301  found%ik22overk)
1302  call mem_set_value(this%ik33overk, 'IK33OVERK', this%input_mempath, &
1303  found%ik33overk)
1304  call mem_set_value(this%inewton, 'INEWTON', this%input_mempath, found%inewton)
1305  call mem_set_value(this%satomega, 'SATOMEGA', this%input_mempath, &
1306  found%satomega)
1307  call mem_set_value(this%irewet, 'IREWET', this%input_mempath, found%irewet)
1308  call mem_set_value(this%wetfct, 'WETFCT', this%input_mempath, found%wetfct)
1309  call mem_set_value(this%iwetit, 'IWETIT', this%input_mempath, found%iwetit)
1310  call mem_set_value(this%ihdwet, 'IHDWET', this%input_mempath, found%ihdwet)
1311  !
1312  ! -- save flows option active
1313  if (found%ipakcb) this%ipakcb = -1
1314  !
1315  ! -- xt3d active with rhs
1316  if (found%ixt3d .and. found%ixt3drhs) this%ixt3d = 2
1317  !
1318  ! -- save specific discharge active
1319  if (found%isavspdis) this%icalcspdis = this%isavspdis
1320  !
1321  ! -- no newton specified
1322  if (found%inewton) then
1323  this%inewton = 0
1324  this%iasym = 0
1325  end if
1326  !
1327  ! -- enforce 0 or 1 TVK6_FILENAME entries in option block
1328  if (filein_fname(tvk6_filename, 'TVK6_FILENAME', this%input_mempath, &
1329  this%input_fname)) then
1330  call openfile(this%intvk, this%iout, tvk6_filename, 'TVK')
1331  call tvk_cr(this%tvk, this%name_model, this%intvk, this%iout)
1332  end if
1333  !
1334  ! -- log options
1335  if (this%iout > 0) then
1336  call this%log_options(found)
1337  end if
memorymanagermodule::mem_setptr
Definition: MemoryManager.f90:89
memorymanagermodule::get_isize
subroutine, public get_isize(name, mem_path, isize)
@ brief Get the number of elements for this variable
Definition: MemoryManager.f90:265
sourcecommonmodule
This module contains the SourceCommonModule.
Definition: SourceCommon.f90:7
sourcecommonmodule::filein_fname
logical(lgp) function, public filein_fname(filename, tagname, input_mempath, input_fname)
enforce and set a single input filename provided via FILEIN keyword
Definition: SourceCommon.f90:452
characterstringmodule::characterstringtype
This class is used to store a single deferred-length character string. It was designed to work in an ...
Definition: CharString.f90:23
Here is the call graph for this function:

◆ store_original_k_arrays()

subroutine gwfnpfmodule::store_original_k_arrays ( class(gwfnpftype)  this,
integer(i4b), intent(in)  ncells,
integer(i4b), intent(in)  njas 
)

The K arrays (K11, etc.) get multiplied by the viscosity ratio so that subsequent uses of K already take into account the effect of viscosity. Thus the original user-specified K array values are lost unless they are backed up in k11input, for example. In a new stress period/time step, the values in k11input are multiplied by the viscosity ratio, not k11 since it contains viscosity-adjusted hydraulic conductivity values.

Definition at line 1118 of file gwf-npf.f90.

1119  ! -- modules
1120  use memorymanagermodule, only: mem_allocate
1121  ! -- dummy
1122  class(GwfNpftype) :: this
1123  integer(I4B), intent(in) :: ncells
1124  integer(I4B), intent(in) :: njas
1125  ! -- local
1126  integer(I4B) :: n
1127  !
1128  ! -- Retain copy of user-specified K arrays
1129  do n = 1, ncells
1130  this%k11input(n) = this%k11(n)
1131  this%k22input(n) = this%k22(n)
1132  this%k33input(n) = this%k33(n)
1133  end do
memorymanagermodule::mem_allocate
Definition: MemoryManager.f90:54