MODFLOW 6  version 6.7.0.dev2
USGS Modular Hydrologic Model
tspadvmodule Module Reference

Data Types

type  tspadvtype
 

Functions/Subroutines

subroutine, public adv_cr (advobj, name_model, input_mempath, inunit, iout, fmi, eqnsclfac)
 @ brief Create a new ADV object More...
 
subroutine adv_df (this, adv_options)
 Define ADV object. More...
 
subroutine adv_ar (this, dis, ibound)
 Allocate and read method for package. More...
 
subroutine adv_dt (this, dtmax, msg, thetam)
 Calculate maximum time step length. More...
 
subroutine adv_fc (this, nodes, matrix_sln, idxglo, cnew, rhs)
 Fill coefficient method for ADV package. More...
 
subroutine adv_cq (this, cnew, flowja)
 Calculate advection contribution to flowja. More...
 
subroutine adv_da (this)
 Deallocate memory. More...
 
subroutine allocate_scalars (this)
 Allocate scalars specific to the streamflow energy transport (SFE) package. More...
 
subroutine source_options (this)
 Source input options. More...
 

Function/Subroutine Documentation

◆ adv_ar()

subroutine tspadvmodule::adv_ar ( class(tspadvtype this,
class(disbasetype), intent(in), pointer  dis,
integer(i4b), dimension(:), intent(in), pointer, contiguous  ibound 
)
private

Method to allocate and read static data for the ADV package.

Definition at line 115 of file tsp-adv.f90.

116  ! -- modules
117  use simmodule, only: store_error
118  ! -- dummy
119  class(TspAdvType) :: this
120  class(DisBaseType), pointer, intent(in) :: dis
121  integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ibound
122  ! -- local
123  integer(I4B) :: iadvwt_value
124  !
125  ! -- adv pointers to arguments that were passed in
126  this%dis => dis
127  this%ibound => ibound
128  !
129  ! -- Create interpolation scheme
130  iadvwt_value = this%iadvwt ! Dereference iadvwt to work with case statement
131  select case (iadvwt_value)
132  case (adv_scheme_upstream)
133  this%face_interpolation = &
134  upstreamschemetype(this%dis, this%fmi)
135  case (adv_scheme_central)
136  this%face_interpolation = &
137  centraldifferenceschemetype(this%dis, this%fmi)
138  case (adv_scheme_tvd)
139  this%face_interpolation = &
140  tvdschemetype(this%dis, this%fmi, this%ibound)
141  case (adv_scheme_utvd)
142  this%gradient = leastsquaresgradienttype(this%dis)
143  this%face_interpolation = &
144  utvdschemetype(this%dis, this%fmi, this%gradient)
145  case default
146  call store_error("Unknown advection scheme", terminate=.true.)
147  end select
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
Here is the call graph for this function:

◆ adv_cq()

subroutine tspadvmodule::adv_cq ( class(tspadvtype this,
real(dp), dimension(:), intent(in)  cnew,
real(dp), dimension(:), intent(inout)  flowja 
)
private

Definition at line 254 of file tsp-adv.f90.

255  ! -- modules
256  ! -- dummy
257  class(TspAdvType) :: this
258  real(DP), intent(in), dimension(:) :: cnew
259  real(DP), intent(inout), dimension(:) :: flowja
260  ! -- local
261  integer(I4B) :: nodes
262  integer(I4B) :: n, m, ipos
263  real(DP) :: qnm
264  type(CoefficientsType) :: coefficients
265  !
266  ! -- Calculate advection and add to flowja. qnm is the volumetric flow
267  ! rate and has dimensions of L^/T.
268  nodes = this%dis%nodes
269 
270  do n = 1, nodes
271  if (this%ibound(n) == 0) cycle
272  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
273  m = this%dis%con%ja(ipos)
274  if (this%ibound(m) == 0) cycle
275  qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
276 
277  coefficients = this%face_interpolation%compute(n, m, ipos, cnew)
278  flowja(ipos) = flowja(ipos) &
279  + qnm * coefficients%c_n * cnew(n) &
280  + qnm * coefficients%c_m * cnew(m) &
281  - qnm * coefficients%rhs
282  end do
283  end do
284 

◆ adv_cr()

subroutine, public tspadvmodule::adv_cr ( type(tspadvtype), pointer  advobj,
character(len=*), intent(in)  name_model,
character(len=*), intent(in)  input_mempath,
integer(i4b), intent(in)  inunit,
integer(i4b), intent(in)  iout,
type(tspfmitype), intent(in), target  fmi,
real(dp), intent(in), pointer  eqnsclfac 
)

Create a new ADV package

Parameters
[in]eqnsclfacgoverning equation scale factor

Definition at line 56 of file tsp-adv.f90.

58  ! -- dummy
59  type(TspAdvType), pointer :: advobj
60  character(len=*), intent(in) :: name_model
61  character(len=*), intent(in) :: input_mempath
62  integer(I4B), intent(in) :: inunit
63  integer(I4B), intent(in) :: iout
64  type(TspFmiType), intent(in), target :: fmi
65  real(DP), intent(in), pointer :: eqnsclfac !< governing equation scale factor
66  !
67  ! -- Create the object
68  allocate (advobj)
69  !
70  ! -- create name and memory path
71  call advobj%set_names(1, name_model, 'ADV', 'ADV', input_mempath)
72  !
73  ! -- Allocate scalars
74  call advobj%allocate_scalars()
75  !
76  ! -- Set variables
77  advobj%inunit = inunit
78  advobj%iout = iout
79  advobj%fmi => fmi
80  advobj%eqnsclfac => eqnsclfac
Here is the caller graph for this function:

◆ adv_da()

subroutine tspadvmodule::adv_da ( class(tspadvtype this)
private

Definition at line 289 of file tsp-adv.f90.

290  ! -- modules
292  ! -- dummy
293  class(TspAdvType) :: this
294  !
295  ! -- Deallocate arrays if package was active
296  if (this%inunit > 0) then
297  end if
298  !
299  ! -- nullify pointers
300  this%ibound => null()
301  !
302  ! -- Scalars
303  call mem_deallocate(this%iadvwt)
304  call mem_deallocate(this%ats_percel)
305  !
306  ! -- deallocate parent
307  call this%NumericalPackageType%da()

◆ adv_df()

subroutine tspadvmodule::adv_df ( class(tspadvtype this,
type(tspadvoptionstype), intent(in), optional  adv_options 
)
private

Define the ADV package

Parameters
[in]adv_optionsthe optional options, for when not constructing from file

Definition at line 87 of file tsp-adv.f90.

88  ! -- dummy
89  class(TspAdvType) :: this
90  type(TspAdvOptionsType), optional, intent(in) :: adv_options !< the optional options, for when not constructing from file
91  ! -- local
92  character(len=*), parameter :: fmtadv = &
93  "(1x,/1x,'ADV -- ADVECTION PACKAGE, VERSION 1, 8/25/2017', &
94  &' INPUT READ FROM MEMPATH: ', A, //)"
95  !
96  ! -- Read or set advection options
97  if (.not. present(adv_options)) then
98  !
99  ! --print a message identifying the advection package.
100  write (this%iout, fmtadv) this%input_mempath
101  !
102  ! --read options from file
103  call this%source_options()
104  else
105  !
106  ! --set options from input arg
107  this%iadvwt = adv_options%iAdvScheme
108  end if

◆ adv_dt()

subroutine tspadvmodule::adv_dt ( class(tspadvtype this,
real(dp), intent(out)  dtmax,
character(len=*), intent(inout)  msg,
real(dp), dimension(:), intent(in)  thetam 
)

Return the largest time step that meets stability constraints

Parameters
thisthis instance
[out]dtmaxmaximum allowable dt subject to stability constraint
[in,out]msgpackage/cell dt constraint message
[in]thetamporosity

Definition at line 154 of file tsp-adv.f90.

155  ! dummy
156  class(TspAdvType) :: this !< this instance
157  real(DP), intent(out) :: dtmax !< maximum allowable dt subject to stability constraint
158  character(len=*), intent(inout) :: msg !< package/cell dt constraint message
159  real(DP), dimension(:), intent(in) :: thetam !< porosity
160  ! local
161  integer(I4B) :: n
162  integer(I4B) :: m
163  integer(I4B) :: ipos
164  integer(I4B) :: nrmax
165  character(len=LINELENGTH) :: cellstr
166  real(DP) :: dt
167  real(DP) :: flowmax
168  real(DP) :: flowsumpos
169  real(DP) :: flowsumneg
170  real(DP) :: flownm
171  real(DP) :: cell_volume
172  dtmax = dnodata
173  nrmax = 0
174  msg = ''
175 
176  ! If ats_percel not specified by user, then return without making
177  ! the courant time step calculation
178  if (this%ats_percel == dnodata) then
179  return
180  end if
181 
182  ! Calculate time step lengths based on stability constraint for each cell
183  ! and store the smallest one
184  do n = 1, this%dis%nodes
185  if (this%ibound(n) == 0) cycle
186  flowsumneg = dzero
187  flowsumpos = dzero
188  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
189  if (this%dis%con%mask(ipos) == 0) cycle
190  m = this%dis%con%ja(ipos)
191  if (this%ibound(m) == 0) cycle
192  flownm = this%fmi%gwfflowja(ipos)
193  if (flownm < dzero) then
194  flowsumneg = flowsumneg - flownm
195  else
196  flowsumpos = flowsumpos + flownm
197  end if
198  end do
199  flowmax = max(flowsumneg, flowsumpos)
200  if (flowmax < dprec) cycle
201  cell_volume = this%dis%get_cell_volume(n, this%dis%top(n))
202  dt = cell_volume * this%fmi%gwfsat(n) * thetam(n) / flowmax
203  dt = dt * this%ats_percel
204  if (dt < dtmax) then
205  dtmax = dt
206  nrmax = n
207  end if
208  end do
209  if (nrmax > 0) then
210  call this%dis%noder_to_string(nrmax, cellstr)
211  write (msg, *) adjustl(trim(this%memoryPath))//'-'//trim(cellstr)
212  end if

◆ adv_fc()

subroutine tspadvmodule::adv_fc ( class(tspadvtype this,
integer(i4b), intent(in)  nodes,
class(matrixbasetype), pointer  matrix_sln,
integer(i4b), dimension(:), intent(in)  idxglo,
real(dp), dimension(:), intent(in)  cnew,
real(dp), dimension(:), intent(inout)  rhs 
)
private

Method to calculate coefficients and fill amat and rhs.

Parameters
thisthis instance
[in]nodesnumber of nodes
matrix_slnpointer to solution matrix
[in]idxgloglobal indices for matrix
[in]cnewnew concentration/temperature values
[in,out]rhsright-hand side vector

Definition at line 219 of file tsp-adv.f90.

220  ! -- modules
221  ! -- dummy
222  class(TspAdvType) :: this !< this instance
223  integer(I4B), intent(in) :: nodes !< number of nodes
224  class(MatrixBaseType), pointer :: matrix_sln !< pointer to solution matrix
225  integer(I4B), intent(in), dimension(:) :: idxglo !< global indices for matrix
226  real(DP), intent(in), dimension(:) :: cnew !< new concentration/temperature values
227  real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector
228  ! -- local
229  integer(I4B) :: n, m, idiag, ipos
230  real(DP) :: qnm !< volumetric flow rate
231  type(CoefficientsType) :: coefficients
232 
233  ! Calculate internal domain fluxes and add to matrix_sln and rhs.
234  do n = 1, nodes
235  if (this%ibound(n) == 0) cycle ! skip inactive nodes
236  idiag = this%dis%con%ia(n)
237  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
238  m = this%dis%con%ja(ipos)
239  if (this%ibound(m) == 0) cycle ! skip inactive nodes
240  if (this%dis%con%mask(ipos) == 0) cycle ! skip masked connections
241 
242  qnm = this%fmi%gwfflowja(ipos) * this%eqnsclfac
243  coefficients = this%face_interpolation%compute(n, m, ipos, cnew)
244 
245  call matrix_sln%add_value_pos(idxglo(idiag), qnm * coefficients%c_n)
246  call matrix_sln%add_value_pos(idxglo(ipos), qnm * coefficients%c_m)
247  rhs(n) = rhs(n) + qnm * coefficients%rhs
248  end do
249  end do

◆ allocate_scalars()

subroutine tspadvmodule::allocate_scalars ( class(tspadvtype this)

Definition at line 313 of file tsp-adv.f90.

314  ! -- modules
316  ! -- dummy
317  class(TspAdvType) :: this
318  ! -- local
319  !
320  ! -- allocate scalars in NumericalPackageType
321  call this%NumericalPackageType%allocate_scalars()
322  !
323  ! -- Allocate
324  call mem_allocate(this%iadvwt, 'IADVWT', this%memoryPath)
325  call mem_allocate(this%ats_percel, 'ATS_PERCEL', this%memoryPath)
326  !
327  ! -- Initialize
328  this%iadvwt = adv_scheme_upstream
329  this%ats_percel = dnodata
330  !
331  ! -- Advection creates an asymmetric coefficient matrix
332  this%iasym = 1

◆ source_options()

subroutine tspadvmodule::source_options ( class(tspadvtype this)

Definition at line 337 of file tsp-adv.f90.

338  ! -- modules
339  use kindmodule, only: lgp
340  use constantsmodule, only: lenvarname
341  use simvariablesmodule, only: errmsg
344  ! -- dummy
345  class(TspAdvType) :: this
346  ! -- locals
347  character(len=LENVARNAME), dimension(4) :: scheme = &
348  &[character(len=LENVARNAME) :: 'UPSTREAM', 'CENTRAL', 'TVD', 'UTVD']
349  logical(lgp) :: found_scheme, found_atspercel
350  ! -- formats
351  character(len=*), parameter :: fmtiadvwt = &
352  &"(4x,'ADVECTION WEIGHTING SCHEME HAS BEEN SET TO: ', a)"
353 
354  ! update defaults with input sourced values
355  call mem_set_value(this%iadvwt, 'SCHEME', this%input_mempath, &
356  scheme, found_scheme)
357  call mem_set_value(this%ats_percel, 'ATS_PERCEL', this%input_mempath, &
358  found_atspercel)
359 
360  if (found_scheme) then
361  ! should currently be set to index of scheme names
362  if (this%iadvwt == 0) then
363  write (errmsg, '(a, a)') &
364  'Unknown scheme, must be "UPSTREAM", "CENTRAL", "TVD" or "UTVD"'
365  call store_error(errmsg)
366  call store_error_filename(this%input_fname)
367  else
368  ! scheme parameters are 0 based
369  this%iadvwt = this%iadvwt - 1
370  end if
371  end if
372  if (found_atspercel) then
373  if (this%ats_percel == dzero) this%ats_percel = dnodata
374  end if
375 
376  ! log options
377  write (this%iout, '(1x,a)') 'PROCESSING ADVECTION OPTIONS'
378  if (found_scheme) then
379  write (this%iout, fmtiadvwt) trim(scheme(this%iadvwt + 1))
380  end if
381  if (found_atspercel) then
382  write (this%iout, '(4x,a,1pg15.6)') &
383  'User-specified fractional cell distance for adaptive time &
384  &steps: ', this%ats_percel
385  end if
386  write (this%iout, '(1x,a)') 'END OF ADVECTION OPTIONS'
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter lenvarname
maximum length of a variable name
Definition: Constants.f90:17
This module defines variable data types.
Definition: kind.f90:8
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
Definition: Sim.f90:203
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
Here is the call graph for this function: