MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
swf-dfw.f90
Go to the documentation of this file.
1 !> @brief Stream Network Flow (SWF) Diffusive Wave (DFW) Module
2 !!
3 !! This module solves one-dimensional flow routing using a diffusive
4 !< wave approach.
5 module swfdfwmodule
6 
7  use kindmodule, only: dp, i4b, lgp
9  dzero, dhalf, done, dtwo, &
11  dprec, dem10
15  use simvariablesmodule, only: errmsg
19  use basedismodule, only: disbasetype
20  use swfcxsmodule, only: swfcxstype
21  use obsmodule, only: obstype, obs_cr
22  use observemodule, only: observetype
24 
25  implicit none
26  private
27  public :: swfdfwtype, dfw_cr
28 
29  type, extends(numericalpackagetype) :: swfdfwtype
30 
31  ! user-provided input
32  integer(I4B), pointer :: is2d => null() !< flag to indicate this model is 2D overland flow and not 1d channel flow
33  integer(I4B), pointer :: icentral => null() !< flag to use central in space weighting (default is upstream weighting)
34  integer(I4B), pointer :: iswrcond => null() !< flag to activate the dev SWR conductance formulation
35  real(DP), pointer :: unitconv !< conversion factor used in mannings equation; calculated from timeconv and lengthconv
36  real(DP), pointer :: timeconv !< conversion factor from model length units to meters (1.0 if model uses meters for length)
37  real(DP), pointer :: lengthconv !< conversion factor from model time units to seconds (1.0 if model uses seconds for time)
38  real(DP), dimension(:), pointer, contiguous :: hnew => null() !< pointer to model xnew
39  real(DP), dimension(:), pointer, contiguous :: manningsn => null() !< mannings roughness for each reach
40  integer(I4B), dimension(:), pointer, contiguous :: idcxs => null() !< cross section id for each reach
41  integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound
42  integer(I4B), dimension(:), pointer, contiguous :: icelltype => null() !< set to 1 and is accessed by chd for checking
43 
44  ! velocity
45  integer(I4B), pointer :: icalcvelocity => null() !< flag to indicate velocity will be calculated (always on)
46  integer(I4B), pointer :: isavvelocity => null() !< flag to indicate velocity will be saved
47  real(DP), dimension(:, :), pointer, contiguous :: vcomp => null() !< velocity components: vx, vy, vz (nodes, 3)
48  real(DP), dimension(:), pointer, contiguous :: vmag => null() !< velocity magnitude (of size nodes)
49  integer(I4B), pointer :: nedges => null() !< number of cell edges
50  integer(I4B), pointer :: lastedge => null() !< last edge number
51  integer(I4B), dimension(:), pointer, contiguous :: nodedge => null() !< array of node numbers that have edges
52  integer(I4B), dimension(:), pointer, contiguous :: ihcedge => null() !< edge type (horizontal or vertical)
53  real(DP), dimension(:, :), pointer, contiguous :: propsedge => null() !< edge properties (Q, area, nx, ny, distance)
54  real(DP), dimension(:), pointer, contiguous :: grad_dhds_mag => null() !< magnitude of the gradient (of size nodes)
55  real(DP), dimension(:), pointer, contiguous :: dhdsja => null() !< gradient for each connection (of size njas)
56 
57  ! observation data
58  integer(I4B), pointer :: inobspkg => null() !< unit number for obs package
59  type(ObsType), pointer :: obs => null() !< observation package
60 
61  ! pointer to cross section data
62  type(SwfCxsType), pointer :: cxs
63 
64  contains
65 
66  procedure :: dfw_df
67  procedure :: allocate_scalars
68  procedure :: allocate_arrays
69  procedure :: dfw_load
70  procedure :: source_options
71  procedure :: log_options
72  procedure :: source_griddata
73  procedure :: log_griddata
74  procedure :: dfw_ar
75  procedure :: dfw_rp
76  procedure :: dfw_ad
77  procedure :: dfw_fc
78  procedure :: dfw_qnm_fc_nr
79  !procedure :: dfw_qnm_fc
80  procedure :: dfw_fn
81  procedure :: dfw_nur
82  procedure :: dfw_cq
83  procedure :: dfw_bd
84  procedure :: dfw_save_model_flows
85  procedure :: dfw_print_model_flows
86  procedure :: dfw_da
87  procedure :: dfw_df_obs
88  procedure :: dfw_rp_obs
89  procedure :: dfw_bd_obs
90  procedure :: qcalc
91  procedure :: get_cond
92  procedure :: get_cond_swr
93  procedure :: get_cond_n
94  procedure :: get_flow_area_nm
95  procedure :: calc_velocity
96  procedure :: sav_velocity
97  procedure, public :: increase_edge_count
98  procedure, public :: set_edge_properties
99  procedure :: calc_dhds
100  procedure :: write_cxs_tables
101 
102  end type swfdfwtype
103 
104 contains
105 
106  !> @brief create package
107  !<
108  subroutine dfw_cr(dfwobj, name_model, input_mempath, inunit, iout, &
109  cxs)
110  ! modules
112  ! dummy
113  type(SwfDfwType), pointer :: dfwobj !< object to create
114  character(len=*), intent(in) :: name_model !< name of the SWF model
115  character(len=*), intent(in) :: input_mempath !< memory path
116  integer(I4B), intent(in) :: inunit !< flag to indicate if package is active
117  integer(I4B), intent(in) :: iout !< unit number for output
118  type(SwfCxsType), pointer, intent(in) :: cxs !< the pointer to the cxs package
119  ! locals
120  logical(LGP) :: found_fname
121  ! formats
122  character(len=*), parameter :: fmtheader = &
123  "(1x, /1x, 'DFW -- DIFFUSIVE WAVE (DFW) PACKAGE, VERSION 1, 9/25/2023', &
124  &' INPUT READ FROM MEMPATH: ', A, /)"
125  !
126  ! Create the object
127  allocate (dfwobj)
128 
129  ! create name and memory path
130  call dfwobj%set_names(1, name_model, 'DFW', 'DFW')
131 
132  ! Allocate scalars
133  call dfwobj%allocate_scalars()
134 
135  ! Set variables
136  dfwobj%input_mempath = input_mempath
137  dfwobj%inunit = inunit
138  dfwobj%iout = iout
139 
140  ! set name of input file
141  call mem_set_value(dfwobj%input_fname, 'INPUT_FNAME', dfwobj%input_mempath, &
142  found_fname)
143 
144  ! Set a pointers to passed in objects
145  dfwobj%cxs => cxs
146 
147  ! create obs package
148  call obs_cr(dfwobj%obs, dfwobj%inobspkg)
149 
150  ! check if dfw is enabled
151  if (inunit > 0) then
152 
153  ! Print a message identifying the package.
154  write (iout, fmtheader) input_mempath
155 
156  end if
157 
158  end subroutine dfw_cr
159 
160  !> @brief load data from IDM to package
161  !<
162  subroutine dfw_df(this, dis)
163  ! dummy
164  class(SwfDfwType) :: this !< this instance
165  class(DisBaseType), pointer, intent(inout) :: dis !< the pointer to the discretization
166 
167  ! Set a pointers to passed in objects
168  this%dis => dis
169 
170  ! Set the distype (either DISV1D or DIS2D)
171  if (this%dis%is_2d()) then
172  this%is2d = 1
173  end if
174 
175  ! check if dfw is enabled
176  ! this will need to become if (.not. present(dfw_options)) then
177  !if (inunit > 0) then
178 
179  ! allocate arrays
180  call this%allocate_arrays()
181 
182  ! load dfw
183  call this%dfw_load()
184 
185  !end if
186 
187  end subroutine dfw_df
188 
189  !> @ brief Allocate scalars
190  !!
191  !! Allocate and initialize scalars for the package. The base model
192  !! allocate scalars method is also called.
193  !!
194  !<
195  subroutine allocate_scalars(this)
196  ! modules
197  ! dummy
198  class(SwfDfwtype) :: this !< this instance
199  !
200  ! allocate scalars in NumericalPackageType
201  call this%NumericalPackageType%allocate_scalars()
202  !
203  ! Allocate scalars
204  call mem_allocate(this%is2d, 'IS2D', this%memoryPath)
205  call mem_allocate(this%icentral, 'ICENTRAL', this%memoryPath)
206  call mem_allocate(this%iswrcond, 'ISWRCOND', this%memoryPath)
207  call mem_allocate(this%unitconv, 'UNITCONV', this%memoryPath)
208  call mem_allocate(this%lengthconv, 'LENGTHCONV', this%memoryPath)
209  call mem_allocate(this%timeconv, 'TIMECONV', this%memoryPath)
210  call mem_allocate(this%inobspkg, 'INOBSPKG', this%memoryPath)
211  call mem_allocate(this%icalcvelocity, 'ICALCVELOCITY', this%memoryPath)
212  call mem_allocate(this%isavvelocity, 'ISAVVELOCITY', this%memoryPath)
213  call mem_allocate(this%nedges, 'NEDGES', this%memoryPath)
214  call mem_allocate(this%lastedge, 'LASTEDGE', this%memoryPath)
215 
216  this%is2d = 0
217  this%icentral = 0
218  this%iswrcond = 0
219  this%unitconv = done
220  this%lengthconv = done
221  this%timeconv = done
222  this%inobspkg = 0
223  this%icalcvelocity = 0
224  this%isavvelocity = 0
225  this%nedges = 0
226  this%lastedge = 0
227 
228  end subroutine allocate_scalars
229 
230  !> @brief allocate memory for arrays
231  !<
232  subroutine allocate_arrays(this)
233  ! dummy
234  class(SwfDfwType) :: this !< this instance
235  ! locals
236  integer(I4B) :: n
237  !
238  ! user-provided input
239  call mem_allocate(this%manningsn, this%dis%nodes, &
240  'MANNINGSN', this%memoryPath)
241  call mem_allocate(this%idcxs, this%dis%nodes, &
242  'IDCXS', this%memoryPath)
243  call mem_allocate(this%icelltype, this%dis%nodes, &
244  'ICELLTYPE', this%memoryPath)
245 
246  ! optional arrays
247  call mem_allocate(this%nodedge, 0, 'NODEDGE', this%memoryPath)
248  call mem_allocate(this%ihcedge, 0, 'IHCEDGE', this%memoryPath)
249  call mem_allocate(this%propsedge, 0, 0, 'PROPSEDGE', this%memoryPath)
250 
251  ! Specific discharge is (re-)allocated when nedges is known
252  call mem_allocate(this%vcomp, 3, 0, 'VCOMP', this%memoryPath)
253  call mem_allocate(this%vmag, 0, 'VMAG', this%memoryPath)
254 
255  do n = 1, this%dis%nodes
256  this%manningsn(n) = dzero
257  this%idcxs(n) = 0
258  this%icelltype(n) = 1
259  end do
260 
261  ! for 2d models, need to calculate and store dhds magnitude
262  if (this%is2d == 1) then
263  call mem_allocate(this%grad_dhds_mag, this%dis%nodes, &
264  'GRAD_DHDS_MAG', this%memoryPath)
265  call mem_allocate(this%dhdsja, this%dis%njas, &
266  'DHDSJA', this%memoryPath)
267  do n = 1, this%dis%nodes
268  this%grad_dhds_mag(n) = dzero
269  end do
270  do n = 1, this%dis%njas
271  this%dhdsja(n) = dzero
272  end do
273  end if
274 
275  end subroutine allocate_arrays
276 
277  !> @brief load data from IDM to package
278  !<
279  subroutine dfw_load(this)
280  ! dummy
281  class(SwfDfwType) :: this !< this instance
282 
283  ! source input data
284  call this%source_options()
285  call this%source_griddata()
286 
287  end subroutine dfw_load
288 
289  !> @brief Copy options from IDM into package
290  !<
291  subroutine source_options(this)
292  ! modules
293  use kindmodule, only: lgp
298  ! dummy
299  class(SwfDfwType) :: this !< this instance
300  ! locals
301  integer(I4B) :: isize
302  type(SwfDfwParamFoundType) :: found
303  type(CharacterStringType), dimension(:), pointer, &
304  contiguous :: obs6_fnames
305 
306  ! update defaults with idm sourced values
307  call mem_set_value(this%icentral, 'ICENTRAL', &
308  this%input_mempath, found%icentral)
309  call mem_set_value(this%iswrcond, 'ISWRCOND', &
310  this%input_mempath, found%iswrcond)
311  call mem_set_value(this%lengthconv, 'LENGTHCONV', &
312  this%input_mempath, found%lengthconv)
313  call mem_set_value(this%timeconv, 'TIMECONV', &
314  this%input_mempath, found%timeconv)
315  call mem_set_value(this%iprflow, 'IPRFLOW', &
316  this%input_mempath, found%iprflow)
317  call mem_set_value(this%ipakcb, 'IPAKCB', &
318  this%input_mempath, found%ipakcb)
319  call mem_set_value(this%isavvelocity, 'ISAVVELOCITY', &
320  this%input_mempath, found%isavvelocity)
321 
322  ! save flows option active
323  if (found%icentral) this%icentral = 1
324  if (found%ipakcb) this%ipakcb = -1
325 
326  ! calculate unit conversion
327  this%unitconv = this%lengthconv**donethird
328  this%unitconv = this%unitconv * this%timeconv
329 
330  ! save velocity active
331  if (found%isavvelocity) this%icalcvelocity = this%isavvelocity
332 
333  ! check for obs6_filename
334  call get_isize('OBS6_FILENAME', this%input_mempath, isize)
335  if (isize > 0) then
336  !
337  if (isize /= 1) then
338  errmsg = 'Multiple OBS6 keywords detected in OPTIONS block.'// &
339  ' Only one OBS6 entry allowed.'
340  call store_error(errmsg)
341  call store_error_filename(this%input_fname)
342  end if
343 
344  call mem_setptr(obs6_fnames, 'OBS6_FILENAME', this%input_mempath)
345 
346  found%obs6_filename = .true.
347  this%obs%inputFilename = obs6_fnames(1)
348  this%obs%active = .true.
349  this%inobspkg = getunit()
350  this%obs%inUnitObs = this%inobspkg
351  call openfile(this%inobspkg, this%iout, this%obs%inputFilename, 'OBS')
352  call this%obs%obs_df(this%iout, this%packName, this%filtyp, this%dis)
353  call this%dfw_df_obs()
354  end if
355 
356  ! log values to list file
357  if (this%iout > 0) then
358  call this%log_options(found)
359  end if
360 
361  end subroutine source_options
362 
363  !> @brief Write user options to list file
364  !<
365  subroutine log_options(this, found)
367  class(SwfDfwType) :: this !< this instance
368  type(SwfDfwParamFoundType), intent(in) :: found
369 
370  write (this%iout, '(1x,a)') 'Setting DFW Options'
371 
372  if (found%lengthconv) then
373  write (this%iout, '(4x,a, G0)') 'Mannings length conversion value &
374  &specified as ', this%lengthconv
375  end if
376 
377  if (found%timeconv) then
378  write (this%iout, '(4x,a, G0)') 'Mannings time conversion value &
379  &specified as ', this%timeconv
380  end if
381 
382  if (found%lengthconv .or. found%timeconv) then
383  write (this%iout, '(4x,a, G0)') 'Mannings conversion value calculated &
384  &from user-provided length_conversion and &
385  &time_conversion is ', this%unitconv
386  end if
387 
388  if (found%iprflow) then
389  write (this%iout, '(4x,a)') 'Cell-by-cell flow information will be printed &
390  &to listing file whenever ICBCFL is not zero.'
391  end if
392 
393  if (found%ipakcb) then
394  write (this%iout, '(4x,a)') 'Cell-by-cell flow information will be printed &
395  &to listing file whenever ICBCFL is not zero.'
396  end if
397 
398  if (found%obs6_filename) then
399  write (this%iout, '(4x,a)') 'Observation package is active.'
400  end if
401 
402  if (found%isavvelocity) &
403  write (this%iout, '(4x,a)') 'Velocity will be calculated at cell &
404  &centers and written to DATA-VCOMP in budget &
405  &file when requested.'
406 
407  if (found%iswrcond) then
408  write (this%iout, '(4x,a, G0)') 'Conductance will be calculated using &
409  &the SWR development option.'
410  end if
411 
412  write (this%iout, '(1x,a,/)') 'End Setting DFW Options'
413 
414  end subroutine log_options
415 
416  !> @brief copy griddata from IDM to package
417  !<
418  subroutine source_griddata(this)
419  ! modules
425  ! dummy
426  class(SwfDfwType) :: this !< this instance
427  ! locals
428  character(len=LENMEMPATH) :: idmMemoryPath
429  type(SwfDfwParamFoundType) :: found
430  integer(I4B), dimension(:), pointer, contiguous :: map
431 
432  ! set memory path
433  idmmemorypath = create_mem_path(this%name_model, 'DFW', idm_context)
434 
435  ! set map to convert user input data into reduced data
436  map => null()
437  if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
438 
439  ! update defaults with idm sourced values
440  call mem_set_value(this%manningsn, 'MANNINGSN', &
441  idmmemorypath, map, found%manningsn)
442  call mem_set_value(this%idcxs, 'IDCXS', idmmemorypath, map, found%idcxs)
443 
444  ! ensure MANNINGSN was found
445  if (.not. found%manningsn) then
446  write (errmsg, '(a)') 'Error in GRIDDATA block: MANNINGSN not found.'
447  call store_error(errmsg)
448  end if
449 
450  if (count_errors() > 0) then
451  call store_error_filename(this%input_fname)
452  end if
453 
454  ! log griddata
455  if (this%iout > 0) then
456  call this%log_griddata(found)
457  end if
458 
459  end subroutine source_griddata
460 
461  !> @brief log griddata to list file
462  !<
463  subroutine log_griddata(this, found)
465  class(SwfDfwType) :: this !< this instance
466  type(SwfDfwParamFoundType), intent(in) :: found
467 
468  write (this%iout, '(1x,a)') 'Setting DFW Griddata'
469 
470  if (found%manningsn) then
471  write (this%iout, '(4x,a)') 'MANNINGSN set from input file'
472  end if
473 
474  if (found%idcxs) then
475  write (this%iout, '(4x,a)') 'IDCXS set from input file'
476  end if
477 
478  call this%write_cxs_tables()
479 
480  write (this%iout, '(1x,a,/)') 'End Setting DFW Griddata'
481 
482  end subroutine log_griddata
483 
484  subroutine write_cxs_tables(this)
485  ! modules
486  ! dummy
487  class(SwfDfwType) :: this !< this instance
488  ! local
489  ! integer(I4B) :: idcxs
490  ! integer(I4B) :: n
491 
492  !-- TODO: write cross section tables
493  ! do n = 1, this%dis%nodes
494  ! idcxs = this%idcxs(n)
495  ! if (idcxs > 0) then
496  ! call this%cxs%write_cxs_table(idcxs, this%width(n), this%slope(n), &
497  ! this%manningsn(n), this%unitconv)
498  ! end if
499  ! end do
500  end subroutine write_cxs_tables
501 
502  !> @brief allocate memory
503  !<
504  subroutine dfw_ar(this, ibound, hnew)
505  ! modules
506  ! dummy
507  class(SwfDfwType) :: this !< this instance
508  integer(I4B), dimension(:), pointer, contiguous :: ibound !< model ibound array
509  real(DP), dimension(:), pointer, contiguous, intent(inout) :: hnew !< pointer to model head array
510  ! local
511  integer(I4B) :: n
512 
513  ! store pointer to ibound
514  this%ibound => ibound
515  this%hnew => hnew
516 
517  if (this%icalcvelocity == 1) then
518  call mem_reallocate(this%vcomp, 3, this%dis%nodes, 'VCOMP', this%memoryPath)
519  call mem_reallocate(this%vmag, this%dis%nodes, 'VMAG', this%memoryPath)
520  call mem_reallocate(this%nodedge, this%nedges, 'NODEDGE', this%memoryPath)
521  call mem_reallocate(this%ihcedge, this%nedges, 'IHCEDGE', this%memoryPath)
522  call mem_reallocate(this%propsedge, 5, this%nedges, 'PROPSEDGE', &
523  this%memoryPath)
524  do n = 1, this%dis%nodes
525  this%vcomp(:, n) = dzero
526  this%vmag(n) = dzero
527  end do
528  end if
529 
530  ! observation data
531  call this%obs%obs_ar()
532 
533  end subroutine dfw_ar
534 
535  !> @brief allocate memory
536  !<
537  subroutine dfw_rp(this)
538  ! modules
539  ! dummy
540  class(SwfDfwType) :: this !< this instance
541 
542  ! read observations
543  call this%dfw_rp_obs()
544 
545  end subroutine dfw_rp
546 
547  !> @brief advance
548  !<
549  subroutine dfw_ad(this, irestore)
550  class(SwfDfwType) :: this !< this instance
551  integer(I4B), intent(in) :: irestore !< ATS flag for retrying time step (1) or advancing (0)
552 
553  ! Push simulated values to preceding time/subtime step
554  call this%obs%obs_ad()
555 
556  end subroutine dfw_ad
557 
558  !> @brief fill coefficients
559  !!
560  !! The DFW Package is entirely Newton based. All matrix and rhs terms
561  !! are added from thish routine.
562  !!
563  !<
564  subroutine dfw_fc(this, kiter, matrix_sln, idxglo, rhs, stage, stage_old)
565  ! modules
566  ! dummy
567  class(SwfDfwType) :: this !< this instance
568  integer(I4B) :: kiter
569  class(MatrixBaseType), pointer :: matrix_sln
570  integer(I4B), intent(in), dimension(:) :: idxglo
571  real(DP), intent(inout), dimension(:) :: rhs
572  real(DP), intent(inout), dimension(:) :: stage
573  real(DP), intent(inout), dimension(:) :: stage_old
574  ! local
575 
576  ! calculate dhds at cell center for 2d case
577  if (this%is2d == 1) then
578  call this%calc_dhds()
579  end if
580 
581  ! add qnm contributions to matrix equations
582  call this%dfw_qnm_fc_nr(kiter, matrix_sln, idxglo, rhs, stage, stage_old)
583 
584  end subroutine dfw_fc
585 
586  !> @brief fill coefficients
587  !!
588  !< Add qnm contributions to matrix equations
589  subroutine dfw_qnm_fc_nr(this, kiter, matrix_sln, idxglo, rhs, stage, stage_old)
590  ! modules
592  ! dummy
593  class(SwfDfwType) :: this !< this instance
594  integer(I4B) :: kiter
595  class(MatrixBaseType), pointer :: matrix_sln
596  integer(I4B), intent(in), dimension(:) :: idxglo
597  real(DP), intent(inout), dimension(:) :: rhs
598  real(DP), intent(inout), dimension(:) :: stage
599  real(DP), intent(inout), dimension(:) :: stage_old
600  ! local
601  integer(I4B) :: n, m, ii, idiag
602  real(DP) :: qnm
603  real(DP) :: qeps
604  real(DP) :: eps
605  real(DP) :: derv
606 
607  ! Calculate conductance and put into amat
608  do n = 1, this%dis%nodes
609 
610  ! Find diagonal position for row n
611  idiag = this%dis%con%ia(n)
612 
613  ! Loop through connections adding matrix terms
614  do ii = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
615 
616  ! skip for masked cells
617  if (this%dis%con%mask(ii) == 0) cycle
618 
619  ! connection variables
620  m = this%dis%con%ja(ii)
621 
622  ! Fill the qnm term on the right-hand side
623  qnm = this%qcalc(n, m, stage(n), stage(m), ii)
624  rhs(n) = rhs(n) - qnm
625 
626  ! Derivative calculation and fill of n terms
627  eps = get_perturbation(stage(n))
628  qeps = this%qcalc(n, m, stage(n) + eps, stage(m), ii)
629  derv = (qeps - qnm) / eps
630  call matrix_sln%add_value_pos(idxglo(idiag), derv)
631  rhs(n) = rhs(n) + derv * stage(n)
632 
633  ! Derivative calculation and fill of m terms
634  eps = get_perturbation(stage(m))
635  qeps = this%qcalc(n, m, stage(n), stage(m) + eps, ii)
636  derv = (qeps - qnm) / eps
637  call matrix_sln%add_value_pos(idxglo(ii), derv)
638  rhs(n) = rhs(n) + derv * stage(m)
639 
640  end do
641  end do
642 
643  end subroutine dfw_qnm_fc_nr
644 
645  !> @brief fill newton
646  !<
647  subroutine dfw_fn(this, kiter, matrix_sln, idxglo, rhs, stage)
648  ! dummy
649  class(SwfDfwType) :: this !< this instance
650  integer(I4B) :: kiter
651  class(MatrixBaseType), pointer :: matrix_sln
652  integer(I4B), intent(in), dimension(:) :: idxglo
653  real(DP), intent(inout), dimension(:) :: rhs
654  real(DP), intent(inout), dimension(:) :: stage
655  ! local
656 
657  ! add newton terms to solution matrix
658  ! todo: add newton terms here instead?
659  ! this routine is probably not necessary as method is fully newton
660 
661  end subroutine dfw_fn
662 
663  !> @brief calculate flow between cells n and m
664  !<
665  function qcalc(this, n, m, stage_n, stage_m, ipos) result(qnm)
666  ! dummy
667  class(SwfDfwType) :: this !< this instance
668  integer(I4B), intent(in) :: n !< number for cell n
669  integer(I4B), intent(in) :: m !< number for cell m
670  real(DP), intent(in) :: stage_n !< stage in reach n
671  real(DP), intent(in) :: stage_m !< stage in reach m
672  integer(I4B), intent(in) :: ipos !< connection number
673  ! local
674  integer(I4B) :: isympos
675  real(DP) :: qnm
676  real(DP) :: cond
677  real(DP) :: cl1
678  real(DP) :: cl2
679 
680  ! Set connection lengths
681  isympos = this%dis%con%jas(ipos)
682  if (n < m) then
683  cl1 = this%dis%con%cl1(isympos)
684  cl2 = this%dis%con%cl2(isympos)
685  else
686  cl1 = this%dis%con%cl2(isympos)
687  cl2 = this%dis%con%cl1(isympos)
688  end if
689 
690  ! Calculate conductance
691  if (this%iswrcond == 0) then
692  cond = this%get_cond(n, m, ipos, stage_n, stage_m, cl1, cl2)
693  else if (this%iswrcond == 1) then
694  cond = this%get_cond_swr(n, m, ipos, stage_n, stage_m, cl1, cl2)
695  end if
696 
697  ! calculate flow between n and m
698  qnm = cond * (stage_m - stage_n)
699 
700  end function qcalc
701 
702  !> @brief calculate effective conductance between cells n and m
703  !!
704  !! Calculate half-cell conductances for cell n and cell m and then use
705  !! harmonic averaging to calculate the effective conductance between the
706  !< two cells.
707  function get_cond(this, n, m, ipos, stage_n, stage_m, cln, clm) result(cond)
708  ! modules
709  use smoothingmodule, only: squadratic
710  ! dummy
711  class(SwfDfwType) :: this !< this instance
712  integer(I4B), intent(in) :: n !< number for cell n
713  integer(I4B), intent(in) :: m !< number for cell m
714  integer(I4B), intent(in) :: ipos !< connection number
715  real(DP), intent(in) :: stage_n !< stage in reach n
716  real(DP), intent(in) :: stage_m !< stage in reach m
717  real(DP), intent(in) :: cln !< distance from cell n to shared face with m
718  real(DP), intent(in) :: clm !< distance from cell m to shared face with n
719  ! local
720  real(DP) :: depth_n
721  real(DP) :: depth_m
722  real(DP) :: dhds_n
723  real(DP) :: dhds_m
724  real(DP) :: width_n
725  real(DP) :: width_m
726  real(DP) :: range = 1.d-6
727  real(DP) :: dydx
728  real(DP) :: smooth_factor
729  real(DP) :: length_nm
730  real(DP) :: cond
731  real(DP) :: cn
732  real(DP) :: cm
733 
734  ! we are using a harmonic conductance approach here; however
735  ! the SWR Process for MODFLOW-2005/NWT uses length-weighted
736  ! average areas and hydraulic radius instead.
737  length_nm = cln + clm
738  cond = dzero
739  if (length_nm > dprec) then
740 
741  ! Calculate depth in each reach
742  depth_n = stage_n - this%dis%bot(n)
743  depth_m = stage_m - this%dis%bot(m)
744 
745  ! assign gradients
746  if (this%is2d == 0) then
747  dhds_n = abs(stage_m - stage_n) / (cln + clm)
748  dhds_m = dhds_n
749  else
750  dhds_n = this%grad_dhds_mag(n)
751  dhds_m = this%grad_dhds_mag(m)
752  end if
753 
754  ! Assign upstream depth, if not central
755  if (this%icentral == 0) then
756  ! use upstream weighting
757  if (stage_n > stage_m) then
758  depth_m = depth_n
759  else
760  depth_n = depth_m
761  end if
762  end if
763 
764  ! Calculate a smoothed depth that goes to zero over
765  ! the specified range
766  call squadratic(depth_n, range, dydx, smooth_factor)
767  depth_n = depth_n * smooth_factor
768  call squadratic(depth_m, range, dydx, smooth_factor)
769  depth_m = depth_m * smooth_factor
770 
771  ! Get the flow widths for n and m from dis package
772  call this%dis%get_flow_width(n, m, ipos, width_n, width_m)
773 
774  ! Calculate half-cell conductance for reach
775  ! n and m
776  cn = this%get_cond_n(n, depth_n, cln, width_n, dhds_n)
777  cm = this%get_cond_n(m, depth_m, clm, width_m, dhds_m)
778 
779  ! Use harmonic mean to calculate weighted
780  ! conductance between the centers of reaches
781  ! n and m
782  if ((cn + cm) > dprec) then
783  cond = cn * cm / (cn + cm)
784  else
785  cond = dzero
786  end if
787 
788  end if
789 
790  end function get_cond
791 
792  !> @brief Calculate half cell conductance
793  !!
794  !! Calculate half-cell conductance for cell n
795  !< using conveyance and Manning's equation
796  function get_cond_n(this, n, depth, dx, width, dhds) result(c)
797  ! modules
798  ! dummy
799  class(SwfDfwType) :: this !< this instance
800  integer(I4B), intent(in) :: n !< reach number
801  real(DP), intent(in) :: depth !< simulated depth (stage - elevation) in reach n for this iteration
802  real(DP), intent(in) :: dx !< half-cell distance
803  real(DP), intent(in) :: width !< width of the reach perpendicular to flow
804  real(DP), intent(in) :: dhds !< gradient
805  ! return
806  real(DP) :: c
807  ! local
808  real(DP) :: rough
809  real(DP) :: dhds_sqr
810  real(DP) :: conveyance
811 
812  ! Calculate conveyance, which is a * r**DTWOTHIRDS / roughc
813  rough = this%manningsn(n)
814  conveyance = this%cxs%get_conveyance(this%idcxs(n), width, depth, rough)
815  dhds_sqr = dhds**dhalf
816  if (dhds_sqr < dem10) then
817  dhds_sqr = dem10
818  end if
819 
820  ! Multiply by unitconv and divide conveyance by sqrt of friction slope and dx
821  c = this%unitconv * conveyance / dx / dhds_sqr
822 
823  end function get_cond_n
824 
825  !> @brief Calculate effective conductance for cells n and m using SWR method
826  !!
827  !! The SWR Process for MODFLOW uses average cell parameters from cell n and
828  !! m to calculate an effective conductance. This is different from the
829  !! default approach used in SWF, which uses harmonic averaging on two half-
830  !< cell conductances.
831  function get_cond_swr(this, n, m, ipos, stage_n, stage_m, cln, clm) result(cond)
832  ! modules
833  use smoothingmodule, only: squadratic
834  ! dummy
835  class(SwfDfwType) :: this !< this instance
836  integer(I4B), intent(in) :: n !< number for cell n
837  integer(I4B), intent(in) :: m !< number for cell m
838  integer(I4B), intent(in) :: ipos !< connection number
839  real(DP), intent(in) :: stage_n !< stage in reach n
840  real(DP), intent(in) :: stage_m !< stage in reach m
841  real(DP), intent(in) :: cln !< distance from cell n to shared face with m
842  real(DP), intent(in) :: clm !< distance from cell m to shared face with n
843  ! local
844  real(DP) :: depth_n
845  real(DP) :: depth_m
846  real(DP) :: dhds_n
847  real(DP) :: dhds_m
848  real(DP) :: dhds_nm
849  real(DP) :: dhds_sqr
850  real(DP) :: width_n
851  real(DP) :: width_m
852  real(DP) :: range = 1.d-6
853  real(DP) :: dydx
854  real(DP) :: smooth_factor
855  real(DP) :: length_nm
856  real(DP) :: cond
857  real(DP) :: ravg
858  real(DP) :: rinv_avg
859  real(DP) :: area_n, area_m, area_avg
860  real(DP) :: rhn, rhm, rhavg
861  real(DP) :: weight_n
862  real(DP) :: weight_m
863  real(DP) :: rough_n
864  real(DP) :: rough_m
865 
866  ! Use harmonic weighting for 1/manningsn, but using length-weighted
867  ! averaging for other terms
868  length_nm = cln + clm
869  cond = dzero
870  if (length_nm > dprec) then
871 
872  ! Calculate depth in each reach
873  depth_n = stage_n - this%dis%bot(n)
874  depth_m = stage_m - this%dis%bot(m)
875 
876  ! Assign upstream depth, if not central
877  if (this%icentral == 0) then
878  ! use upstream weighting
879  if (stage_n > stage_m) then
880  depth_m = depth_n
881  else
882  depth_n = depth_m
883  end if
884  end if
885 
886  ! Calculate a smoothed depth that goes to zero over
887  ! the specified range
888  call squadratic(depth_n, range, dydx, smooth_factor)
889  depth_n = depth_n * smooth_factor
890  call squadratic(depth_m, range, dydx, smooth_factor)
891  depth_m = depth_m * smooth_factor
892 
893  ! Get the flow widths for n and m from dis package
894  call this%dis%get_flow_width(n, m, ipos, width_n, width_m)
895 
896  ! linear weight toward node closer to shared face
897  weight_n = clm / length_nm
898  weight_m = done - weight_n
899 
900  ! average cross sectional flow area
901  area_n = this%cxs%get_area(this%idcxs(n), width_n, depth_n)
902  area_m = this%cxs%get_area(this%idcxs(m), width_m, depth_m)
903  area_avg = weight_n * area_n + weight_m * area_m
904 
905  ! average hydraulic radius
906  if (this%is2d == 0) then
907  rhn = this%cxs%get_hydraulic_radius(this%idcxs(n), width_n, &
908  depth_n, area_n)
909  rhm = this%cxs%get_hydraulic_radius(this%idcxs(m), width_m, &
910  depth_m, area_m)
911  rhavg = weight_n * rhn + weight_m * rhm
912  else
913  rhavg = area_avg / width_n
914  end if
915  rhavg = rhavg**dtwothirds
916 
917  ! average gradient
918  if (this%is2d == 0) then
919  dhds_nm = abs(stage_m - stage_n) / (length_nm)
920  else
921  dhds_n = this%grad_dhds_mag(n)
922  dhds_m = this%grad_dhds_mag(m)
923  dhds_nm = weight_n * dhds_n + weight_m * dhds_m
924  end if
925  dhds_sqr = dhds_nm**dhalf
926  if (dhds_sqr < dem10) then
927  dhds_sqr = dem10
928  end if
929 
930  ! weighted harmonic mean for inverse mannings value
931  weight_n = cln / length_nm
932  weight_m = done - weight_n
933  rough_n = this%cxs%get_roughness(this%idcxs(n), width_n, depth_n, &
934  this%manningsn(n), dhds_nm)
935  rough_m = this%cxs%get_roughness(this%idcxs(m), width_m, depth_m, &
936  this%manningsn(m), dhds_nm)
937  ravg = (weight_n + weight_m) / &
938  (weight_n / rough_n + weight_m / rough_m)
939  rinv_avg = done / ravg
940 
941  ! calculate conductance using averaged values
942  cond = this%unitconv * rinv_avg * area_avg * rhavg / dhds_sqr / length_nm
943 
944  end if
945 
946  end function get_cond_swr
947 
948  !> @brief Calculate flow area between cell n and m
949  !!
950  !! Calculate an average flow area between cell n and m.
951  !! First calculate a flow area for cell n and then for
952  !! cell m and linearly weight the areas using the connection
953  !< distances.
954  function get_flow_area_nm(this, n, m, stage_n, stage_m, cln, clm, &
955  ipos) result(area_avg)
956  ! module
957  use smoothingmodule, only: squadratic
958  ! dummy
959  class(SwfDfwType) :: this !< this instance
960  integer(I4B), intent(in) :: n
961  integer(I4B), intent(in) :: m
962  real(DP), intent(in) :: stage_n
963  real(DP), intent(in) :: stage_m
964  real(DP), intent(in) :: cln
965  real(DP), intent(in) :: clm
966  integer(I4B), intent(in) :: ipos
967  ! local
968  real(DP) :: depth_n
969  real(DP) :: depth_m
970  real(DP) :: width_n
971  real(DP) :: width_m
972  real(DP) :: area_n
973  real(DP) :: area_m
974  real(DP) :: weight_n
975  real(DP) :: weight_m
976  real(DP) :: length_nm
977  real(DP) :: range = 1.d-6
978  real(DP) :: dydx
979  real(DP) :: smooth_factor
980  ! return
981  real(DP) :: area_avg
982 
983  ! depths
984  depth_n = stage_n - this%dis%bot(n)
985  depth_m = stage_m - this%dis%bot(m)
986 
987  ! Assign upstream depth, if not central
988  if (this%icentral == 0) then
989  ! use upstream weighting
990  if (stage_n > stage_m) then
991  depth_m = depth_n
992  else
993  depth_n = depth_m
994  end if
995  end if
996 
997  ! Calculate a smoothed depth that goes to zero over
998  ! the specified range
999  call squadratic(depth_n, range, dydx, smooth_factor)
1000  depth_n = depth_n * smooth_factor
1001  call squadratic(depth_m, range, dydx, smooth_factor)
1002  depth_m = depth_m * smooth_factor
1003 
1004  ! Get the flow widths for n and m from dis package
1005  call this%dis%get_flow_width(n, m, ipos, width_n, width_m)
1006 
1007  ! linear weight toward node closer to shared face
1008  length_nm = cln + clm
1009  weight_n = clm / length_nm
1010  weight_m = done - weight_n
1011 
1012  ! average cross sectional flow area
1013  area_n = this%cxs%get_area(this%idcxs(n), width_n, depth_n)
1014  area_m = this%cxs%get_area(this%idcxs(m), width_m, depth_m)
1015  area_avg = weight_n * area_n + weight_m * area_m
1016 
1017  end function get_flow_area_nm
1018 
1019  !> @brief Calculate average hydraulic gradient magnitude for each cell
1020  !!
1021  !! Go through each cell and calculate the average hydraulic gradient using
1022  !! an xt3d-style gradient interpolation. This is used for 2D grids in the
1023  !! calculation of an effective conductance. For 1D grids, gradients are
1024  !< calculated between cell centers.
1025  subroutine calc_dhds(this)
1026  ! modules
1028  ! dummy
1029  class(SwfDfwType) :: this !< this instance
1030  ! local
1031  integer(I4B) :: n
1032  integer(I4B) :: m
1033  integer(I4B) :: ipos
1034  integer(I4B) :: isympos
1035  real(DP) :: cl1
1036  real(DP) :: cl2
1037 
1038  do n = 1, this%dis%nodes
1039  this%grad_dhds_mag(n) = dzero
1040  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1041  m = this%dis%con%ja(ipos)
1042  isympos = this%dis%con%jas(ipos)
1043 
1044  ! determine cl1 and cl2
1045  if (n < m) then
1046  cl1 = this%dis%con%cl1(isympos)
1047  cl2 = this%dis%con%cl2(isympos)
1048  else
1049  cl1 = this%dis%con%cl2(isympos)
1050  cl2 = this%dis%con%cl1(isympos)
1051  end if
1052 
1053  ! store for n < m in upper right triangular part of symmetric dhdsja array
1054  if (n < m) then
1055  if (cl1 + cl2 > dprec) then
1056  this%dhdsja(isympos) = (this%hnew(m) - this%hnew(n)) / (cl1 + cl2)
1057  else
1058  this%dhdsja(isympos) = dzero
1059  end if
1060  end if
1061  end do
1062  end do
1063 
1064  ! pass dhdsja into the vector interpolation to get the components
1065  ! of the gradient at the cell center
1066  call vector_interpolation_2d(this%dis, this%dhdsja, vmag=this%grad_dhds_mag)
1067 
1068  end subroutine calc_dhds
1069 
1070  !> @ brief Newton under relaxation
1071  !!
1072  !< If stage is below the bottom, then pull it up a bit
1073  subroutine dfw_nur(this, neqmod, x, xtemp, dx, inewtonur, dxmax, locmax)
1074  ! dummy
1075  class(SwfDfwType) :: this !< this instance
1076  integer(I4B), intent(in) :: neqmod !< number of equations
1077  real(DP), dimension(neqmod), intent(inout) :: x !< dependent variable
1078  real(DP), dimension(neqmod), intent(in) :: xtemp !< temporary dependent variable
1079  real(DP), dimension(neqmod), intent(inout) :: dx !< change in dependent variable
1080  integer(I4B), intent(inout) :: inewtonur !< flag to indication relaxation was applied
1081  real(DP), intent(inout) :: dxmax !< max change in x
1082  integer(I4B), intent(inout) :: locmax !< location of max change
1083  ! local
1084  integer(I4B) :: n
1085  real(DP) :: botm
1086  real(DP) :: xx
1087  real(DP) :: dxx
1088 
1089  ! Newton-Raphson under-relaxation
1090  do n = 1, this%dis%nodes
1091  if (this%ibound(n) < 1) cycle
1092  if (this%icelltype(n) > 0) then
1093  botm = this%dis%bot(n)
1094  ! only apply Newton-Raphson under-relaxation if
1095  ! solution head is below the bottom of the model
1096  if (x(n) < botm) then
1097  inewtonur = 1
1098  xx = xtemp(n) * (done - dp9) + botm * dp9
1099  dxx = x(n) - xx
1100  if (abs(dxx) > abs(dxmax)) then
1101  locmax = n
1102  dxmax = dxx
1103  end if
1104  x(n) = xx
1105  dx(n) = dzero
1106  end if
1107  end if
1108  end do
1109 
1110  end subroutine dfw_nur
1111 
1112  !> @ brief Calculate flow for each connection and store in flowja
1113  !<
1114  subroutine dfw_cq(this, stage, stage_old, flowja)
1115  ! dummy
1116  class(SwfDfwType) :: this !< this instance
1117  real(DP), intent(inout), dimension(:) :: stage !< calculated head
1118  real(DP), intent(inout), dimension(:) :: stage_old !< calculated head from previous time step
1119  real(DP), intent(inout), dimension(:) :: flowja !< vector of flows in CSR format
1120  ! local
1121  integer(I4B) :: n, ipos, m
1122  real(DP) :: qnm
1123 
1124  do n = 1, this%dis%nodes
1125  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1126  m = this%dis%con%ja(ipos)
1127  if (m < n) cycle
1128  qnm = this%qcalc(n, m, stage(n), stage(m), ipos)
1129  flowja(ipos) = qnm
1130  flowja(this%dis%con%isym(ipos)) = -qnm
1131  end do
1132  end do
1133 
1134  end subroutine dfw_cq
1135 
1136  !> @ brief Model budget calculation for package
1137  !<
1138  subroutine dfw_bd(this, isuppress_output, model_budget)
1139  ! modules
1140  use budgetmodule, only: budgettype
1141  ! dummy variables
1142  class(SwfDfwType) :: this !< this instance
1143  integer(I4B), intent(in) :: isuppress_output !< flag to suppress model output
1144  type(BudgetType), intent(inout) :: model_budget !< model budget object
1145 
1146  ! Add any DFW budget terms
1147  ! none
1148 
1149  end subroutine dfw_bd
1150 
1151  !> @ brief save flows for package
1152  !<
1153  subroutine dfw_save_model_flows(this, flowja, icbcfl, icbcun)
1154  ! dummy
1155  class(SwfDfwType) :: this !< this instance
1156  real(DP), dimension(:), intent(in) :: flowja !< vector of flows in CSR format
1157  integer(I4B), intent(in) :: icbcfl !< flag to indicate if flows should be saved
1158  integer(I4B), intent(in) :: icbcun !< unit number for flow output
1159  ! local
1160  integer(I4B) :: ibinun
1161 
1162  ! Set unit number for binary output
1163  if (this%ipakcb < 0) then
1164  ibinun = icbcun
1165  elseif (this%ipakcb == 0) then
1166  ibinun = 0
1167  else
1168  ibinun = this%ipakcb
1169  end if
1170  if (icbcfl == 0) ibinun = 0
1171 
1172  ! Write the face flows if requested
1173  if (ibinun /= 0) then
1174  ! flowja
1175  call this%dis%record_connection_array(flowja, ibinun, this%iout)
1176  end if
1177 
1178  ! Calculate velocities at cell centers and write, if requested
1179  if (this%isavvelocity /= 0) then
1180  if (ibinun /= 0) call this%sav_velocity(ibinun)
1181  end if
1182 
1183  end subroutine dfw_save_model_flows
1184 
1185  !> @ brief print flows for package
1186  !<
1187  subroutine dfw_print_model_flows(this, ibudfl, flowja)
1188  ! modules
1189  use tdismodule, only: kper, kstp
1190  use constantsmodule, only: lenbigline
1191  ! dummy
1192  class(SwfDfwType) :: this !< this instance
1193  integer(I4B), intent(in) :: ibudfl !< print flag
1194  real(DP), intent(inout), dimension(:) :: flowja !< vector of flows in CSR format
1195  ! local
1196  character(len=LENBIGLINE) :: line
1197  character(len=30) :: tempstr
1198  integer(I4B) :: n, ipos, m
1199  real(DP) :: qnm
1200  ! formats
1201  character(len=*), parameter :: fmtiprflow = &
1202  &"(/,4x,'CALCULATED INTERCELL FLOW FOR PERIOD ', i0, ' STEP ', i0)"
1203 
1204  ! Write flowja to list file if requested
1205  if (ibudfl /= 0 .and. this%iprflow > 0) then
1206  write (this%iout, fmtiprflow) kper, kstp
1207  do n = 1, this%dis%nodes
1208  line = ''
1209  call this%dis%noder_to_string(n, tempstr)
1210  line = trim(tempstr)//':'
1211  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1212  m = this%dis%con%ja(ipos)
1213  call this%dis%noder_to_string(m, tempstr)
1214  line = trim(line)//' '//trim(tempstr)
1215  qnm = flowja(ipos)
1216  write (tempstr, '(1pg15.6)') qnm
1217  line = trim(line)//' '//trim(adjustl(tempstr))
1218  end do
1219  write (this%iout, '(a)') trim(line)
1220  end do
1221  end if
1222 
1223  end subroutine dfw_print_model_flows
1224 
1225  !> @brief deallocate memory
1226  !<
1227  subroutine dfw_da(this)
1228  ! modules
1231  use simvariablesmodule, only: idm_context
1232  ! dummy
1233  class(SwfDfwType) :: this !< this instance
1234 
1235  ! Deallocate input memory
1236  call memorystore_remove(this%name_model, 'DFW', idm_context)
1237 
1238  ! Deallocate arrays
1239  call mem_deallocate(this%manningsn)
1240  call mem_deallocate(this%idcxs)
1241  call mem_deallocate(this%icelltype)
1242  call mem_deallocate(this%nodedge)
1243  call mem_deallocate(this%ihcedge)
1244  call mem_deallocate(this%propsedge)
1245  call mem_deallocate(this%vcomp)
1246  call mem_deallocate(this%vmag)
1247  if (this%is2d == 1) then
1248  call mem_deallocate(this%grad_dhds_mag)
1249  call mem_deallocate(this%dhdsja)
1250  end if
1251 
1252  ! Scalars
1253  call mem_deallocate(this%is2d)
1254  call mem_deallocate(this%icentral)
1255  call mem_deallocate(this%iswrcond)
1256  call mem_deallocate(this%unitconv)
1257  call mem_deallocate(this%lengthconv)
1258  call mem_deallocate(this%timeconv)
1259  call mem_deallocate(this%isavvelocity)
1260  call mem_deallocate(this%icalcvelocity)
1261  call mem_deallocate(this%nedges)
1262  call mem_deallocate(this%lastedge)
1263 
1264  ! obs package
1265  call mem_deallocate(this%inobspkg)
1266  call this%obs%obs_da()
1267  deallocate (this%obs)
1268  nullify (this%obs)
1269  nullify (this%cxs)
1270 
1271  ! deallocate parent
1272  call this%NumericalPackageType%da()
1273 
1274  ! pointers
1275  this%hnew => null()
1276 
1277  end subroutine dfw_da
1278 
1279  !> @brief Calculate the 3 components of velocity at the cell center
1280  !!
1281  !! todo: duplicated from NPF; should consolidate
1282  !<
1283  subroutine calc_velocity(this, flowja)
1284  ! modules
1285  ! dummy
1286  class(SwfDfwType) :: this !< this instance
1287  real(DP), intent(in), dimension(:) :: flowja !< vector of flows in CSR format
1288  ! local
1289  integer(I4B) :: n
1290  integer(I4B) :: m
1291  integer(I4B) :: ipos
1292  integer(I4B) :: isympos
1293  integer(I4B) :: ihc
1294  integer(I4B) :: ic
1295  integer(I4B) :: iz
1296  integer(I4B) :: nc
1297  integer(I4B) :: ncz
1298  real(DP) :: vx
1299  real(DP) :: vy
1300  real(DP) :: vz
1301  real(DP) :: xn
1302  real(DP) :: yn
1303  real(DP) :: zn
1304  real(DP) :: xc
1305  real(DP) :: yc
1306  real(DP) :: zc
1307  real(DP) :: cl1
1308  real(DP) :: cl2
1309  real(DP) :: dltot
1310  real(DP) :: ooclsum
1311  real(DP) :: dsumx
1312  real(DP) :: dsumy
1313  real(DP) :: dsumz
1314  real(DP) :: denom
1315  real(DP) :: area
1316  real(DP) :: axy
1317  real(DP) :: ayx
1318  real(DP), allocatable, dimension(:) :: vi
1319  real(DP), allocatable, dimension(:) :: di
1320  real(DP), allocatable, dimension(:) :: viz
1321  real(DP), allocatable, dimension(:) :: diz
1322  real(DP), allocatable, dimension(:) :: nix
1323  real(DP), allocatable, dimension(:) :: niy
1324  real(DP), allocatable, dimension(:) :: wix
1325  real(DP), allocatable, dimension(:) :: wiy
1326  real(DP), allocatable, dimension(:) :: wiz
1327  real(DP), allocatable, dimension(:) :: bix
1328  real(DP), allocatable, dimension(:) :: biy
1329  logical :: nozee = .true.
1330 
1331  ! Ensure dis has necessary information
1332  ! todo: do we need this for SWF?
1333  if (this%icalcvelocity /= 0 .and. this%dis%con%ianglex == 0) then
1334  call store_error('Error. ANGLDEGX not provided in '// &
1335  'discretization file. ANGLDEGX required for '// &
1336  'calculation of velocity.', terminate=.true.)
1337  end if
1338 
1339  ! Find max number of connections and allocate weight arrays
1340  nc = 0
1341  do n = 1, this%dis%nodes
1342 
1343  ! Count internal model connections
1344  ic = this%dis%con%ia(n + 1) - this%dis%con%ia(n) - 1
1345 
1346  ! Count edge connections
1347  do m = 1, this%nedges
1348  if (this%nodedge(m) == n) then
1349  ic = ic + 1
1350  end if
1351  end do
1352 
1353  ! Set max number of connections for any cell
1354  if (ic > nc) nc = ic
1355  end do
1356 
1357  ! Allocate storage arrays needed for cell-centered calculation
1358  allocate (vi(nc))
1359  allocate (di(nc))
1360  allocate (viz(nc))
1361  allocate (diz(nc))
1362  allocate (nix(nc))
1363  allocate (niy(nc))
1364  allocate (wix(nc))
1365  allocate (wiy(nc))
1366  allocate (wiz(nc))
1367  allocate (bix(nc))
1368  allocate (biy(nc))
1369 
1370  ! Go through each cell and calculate specific discharge
1371  do n = 1, this%dis%nodes
1372 
1373  ! first calculate geometric properties for x and y directions and
1374  ! the specific discharge at a face (vi)
1375  ic = 0
1376  iz = 0
1377  vi(:) = dzero
1378  di(:) = dzero
1379  viz(:) = dzero
1380  diz(:) = dzero
1381  nix(:) = dzero
1382  niy(:) = dzero
1383  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
1384  m = this%dis%con%ja(ipos)
1385  isympos = this%dis%con%jas(ipos)
1386  ihc = this%dis%con%ihc(isympos)
1387  ic = ic + 1
1388  call this%dis%connection_normal(n, m, ihc, xn, yn, zn, ipos)
1389  call this%dis%connection_vector(n, m, nozee, done, done, &
1390  ihc, xc, yc, zc, dltot)
1391  cl1 = this%dis%con%cl1(isympos)
1392  cl2 = this%dis%con%cl2(isympos)
1393  if (m < n) then
1394  cl1 = this%dis%con%cl2(isympos)
1395  cl2 = this%dis%con%cl1(isympos)
1396  end if
1397  ooclsum = done / (cl1 + cl2)
1398  nix(ic) = -xn
1399  niy(ic) = -yn
1400  di(ic) = dltot * cl1 * ooclsum
1401  area = this%get_flow_area_nm(n, m, this%hnew(n), this%hnew(m), &
1402  cl1, cl2, ipos)
1403  if (area > dzero) then
1404  vi(ic) = flowja(ipos) / area
1405  else
1406  vi(ic) = dzero
1407  end if
1408 
1409  end do
1410 
1411  ! Look through edge flows that may have been provided by an exchange
1412  ! and incorporate them into the averaging arrays
1413  do m = 1, this%nedges
1414  if (this%nodedge(m) == n) then
1415 
1416  ! propsedge: (Q, area, nx, ny, distance)
1417  ihc = this%ihcedge(m)
1418  area = this%propsedge(2, m)
1419 
1420  ic = ic + 1
1421  nix(ic) = -this%propsedge(3, m)
1422  niy(ic) = -this%propsedge(4, m)
1423  di(ic) = this%propsedge(5, m)
1424  if (area > dzero) then
1425  vi(ic) = this%propsedge(1, m) / area
1426  else
1427  vi(ic) = dzero
1428  end if
1429 
1430  end if
1431  end do
1432 
1433  ! Assign number of vertical and horizontal connections
1434  ncz = iz
1435  nc = ic
1436 
1437  ! calculate z weight (wiz) and z velocity
1438  if (ncz == 1) then
1439  wiz(1) = done
1440  else
1441  dsumz = dzero
1442  do iz = 1, ncz
1443  dsumz = dsumz + diz(iz)
1444  end do
1445  denom = (ncz - done)
1446  if (denom < dzero) denom = dzero
1447  dsumz = dsumz + dem10 * dsumz
1448  do iz = 1, ncz
1449  if (dsumz > dzero) wiz(iz) = done - diz(iz) / dsumz
1450  if (denom > 0) then
1451  wiz(iz) = wiz(iz) / denom
1452  else
1453  wiz(iz) = dzero
1454  end if
1455  end do
1456  end if
1457  vz = dzero
1458  do iz = 1, ncz
1459  vz = vz + wiz(iz) * viz(iz)
1460  end do
1461 
1462  ! distance-based weighting
1463  nc = ic
1464  dsumx = dzero
1465  dsumy = dzero
1466  dsumz = dzero
1467  do ic = 1, nc
1468  wix(ic) = di(ic) * abs(nix(ic))
1469  wiy(ic) = di(ic) * abs(niy(ic))
1470  dsumx = dsumx + wix(ic)
1471  dsumy = dsumy + wiy(ic)
1472  end do
1473 
1474  ! Finish computing omega weights. Add a tiny bit
1475  ! to dsum so that the normalized omega weight later
1476  ! evaluates to (essentially) 1 in the case of a single
1477  ! relevant connection, avoiding 0/0.
1478  dsumx = dsumx + dem10 * dsumx
1479  dsumy = dsumy + dem10 * dsumy
1480  do ic = 1, nc
1481  wix(ic) = (dsumx - wix(ic)) * abs(nix(ic))
1482  wiy(ic) = (dsumy - wiy(ic)) * abs(niy(ic))
1483  end do
1484 
1485  ! compute B weights
1486  dsumx = dzero
1487  dsumy = dzero
1488  do ic = 1, nc
1489  bix(ic) = wix(ic) * sign(done, nix(ic))
1490  biy(ic) = wiy(ic) * sign(done, niy(ic))
1491  dsumx = dsumx + wix(ic) * abs(nix(ic))
1492  dsumy = dsumy + wiy(ic) * abs(niy(ic))
1493  end do
1494  if (dsumx > dzero) dsumx = done / dsumx
1495  if (dsumy > dzero) dsumy = done / dsumy
1496  axy = dzero
1497  ayx = dzero
1498  do ic = 1, nc
1499  bix(ic) = bix(ic) * dsumx
1500  biy(ic) = biy(ic) * dsumy
1501  axy = axy + bix(ic) * niy(ic)
1502  ayx = ayx + biy(ic) * nix(ic)
1503  end do
1504 
1505  ! Calculate specific discharge. The divide by zero checking below
1506  ! is problematic for cells with only one flow, such as can happen
1507  ! with triangular cells in corners. In this case, the resulting
1508  ! cell velocity will be calculated as zero. The method should be
1509  ! improved so that edge flows of zero are included in these
1510  ! calculations. But this needs to be done with consideration for LGR
1511  ! cases in which flows are submitted from an exchange.
1512  vx = dzero
1513  vy = dzero
1514  do ic = 1, nc
1515  vx = vx + (bix(ic) - axy * biy(ic)) * vi(ic)
1516  vy = vy + (biy(ic) - ayx * bix(ic)) * vi(ic)
1517  end do
1518  denom = done - axy * ayx
1519  if (denom /= dzero) then
1520  vx = vx / denom
1521  vy = vy / denom
1522  end if
1523 
1524  this%vcomp(1, n) = vx
1525  this%vcomp(2, n) = vy
1526  this%vcomp(3, n) = vz
1527  this%vmag(n) = sqrt(vx**2 + vy**2 + vz**2)
1528 
1529  end do
1530 
1531  ! cleanup
1532  deallocate (vi)
1533  deallocate (di)
1534  deallocate (nix)
1535  deallocate (niy)
1536  deallocate (wix)
1537  deallocate (wiy)
1538  deallocate (wiz)
1539  deallocate (bix)
1540  deallocate (biy)
1541 
1542  end subroutine calc_velocity
1543 
1544  !> @brief Reserve space for nedges cells that have an edge on them.
1545  !!
1546  !! todo: duplicated from NPF; should consolidate
1547  !! This must be called before the swf%allocate_arrays routine, which is
1548  !< called from swf%ar.
1549  subroutine increase_edge_count(this, nedges)
1550  ! dummy
1551  class(SwfDfwType) :: this !< this instance
1552  integer(I4B), intent(in) :: nedges
1553 
1554  this%nedges = this%nedges + nedges
1555 
1556  end subroutine increase_edge_count
1557 
1558  !> @brief Provide the swf package with edge properties
1559  !!
1560  !! todo: duplicated from NPF; should consolidate
1561  !<
1562  subroutine set_edge_properties(this, nodedge, ihcedge, q, area, nx, ny, &
1563  distance)
1564  ! dummy
1565  class(SwfDfwType) :: this !< this instance
1566  integer(I4B), intent(in) :: nodedge
1567  integer(I4B), intent(in) :: ihcedge
1568  real(DP), intent(in) :: q
1569  real(DP), intent(in) :: area
1570  real(DP), intent(in) :: nx
1571  real(DP), intent(in) :: ny
1572  real(DP), intent(in) :: distance
1573  ! local
1574  integer(I4B) :: lastedge
1575 
1576  this%lastedge = this%lastedge + 1
1577  lastedge = this%lastedge
1578  this%nodedge(lastedge) = nodedge
1579  this%ihcedge(lastedge) = ihcedge
1580  this%propsedge(1, lastedge) = q
1581  this%propsedge(2, lastedge) = area
1582  this%propsedge(3, lastedge) = nx
1583  this%propsedge(4, lastedge) = ny
1584  this%propsedge(5, lastedge) = distance
1585 
1586  ! If this is the last edge, then the next call must be starting a new
1587  ! edge properties assignment loop, so need to reset lastedge to 0
1588  if (this%lastedge == this%nedges) this%lastedge = 0
1589 
1590  end subroutine set_edge_properties
1591 
1592  !> @brief Save specific discharge in binary format to ibinun
1593  !!
1594  !! todo: should write 2d velocity; what about for 1D channel?
1595  !<
1596  subroutine sav_velocity(this, ibinun)
1597  ! dummy
1598  class(SwfDfwType) :: this !< this instance
1599  integer(I4B), intent(in) :: ibinun
1600  ! local
1601  character(len=16) :: text
1602  character(len=16), dimension(3) :: auxtxt
1603  integer(I4B) :: n
1604  integer(I4B) :: naux
1605 
1606  ! Write the header
1607  text = ' DATA-VCOMP'
1608  naux = 3
1609  auxtxt(:) = [' vx', ' vy', ' vz']
1610  call this%dis%record_srcdst_list_header(text, this%name_model, &
1611  this%packName, this%name_model, &
1612  this%packName, naux, auxtxt, ibinun, &
1613  this%dis%nodes, this%iout)
1614 
1615  ! Write a zero for Q, and then write qx, qy, qz as aux variables
1616  do n = 1, this%dis%nodes
1617  call this%dis%record_mf6_list_entry(ibinun, n, n, dzero, naux, &
1618  this%vcomp(:, n))
1619  end do
1620 
1621  end subroutine sav_velocity
1622 
1623  !> @brief Define the observation types available in the package
1624  !!
1625  !< Method to define the observation types available in the package.
1626  subroutine dfw_df_obs(this)
1627  ! dummy variables
1628  class(SwfDfwType) :: this !< this instance
1629  ! local variables
1630  integer(I4B) :: indx
1631 
1632  ! Store obs type and assign procedure pointer
1633  ! for ext-outflow observation type.
1634  call this%obs%StoreObsType('ext-outflow', .true., indx)
1635  this%obs%obsData(indx)%ProcessIdPtr => dfwobsidprocessor
1636 
1637  end subroutine dfw_df_obs
1638 
1639  subroutine dfwobsidprocessor(obsrv, dis, inunitobs, iout)
1640  ! dummy
1641  type(ObserveType), intent(inout) :: obsrv
1642  class(DisBaseType), intent(in) :: dis
1643  integer(I4B), intent(in) :: inunitobs
1644  integer(I4B), intent(in) :: iout
1645  ! local
1646  integer(I4B) :: n
1647  character(len=LINELENGTH) :: string
1648 
1649  ! Initialize variables
1650  string = obsrv%IDstring
1651  read (string, *) n
1652 
1653  if (n > 0) then
1654  obsrv%NodeNumber = n
1655  else
1656  errmsg = 'Error reading data from ID string'
1657  call store_error(errmsg)
1658  call store_error_unit(inunitobs)
1659  end if
1660 
1661  end subroutine dfwobsidprocessor
1662 
1663  !> @brief Save observations for the package
1664  !!
1665  !< Method to save simulated values for the package.
1666  subroutine dfw_bd_obs(this)
1667  ! dummy variables
1668  class(SwfDfwType) :: this !< this instance
1669  ! local variables
1670  integer(I4B) :: i
1671  integer(I4B) :: j
1672  integer(I4B) :: n
1673  real(DP) :: v
1674  character(len=100) :: msg
1675  type(ObserveType), pointer :: obsrv => null()
1676 
1677  ! Write simulated values for all observations
1678  if (this%obs%npakobs > 0) then
1679  call this%obs%obs_bd_clear()
1680  do i = 1, this%obs%npakobs
1681  obsrv => this%obs%pakobs(i)%obsrv
1682  do j = 1, obsrv%indxbnds_count
1683  n = obsrv%indxbnds(j)
1684  v = dzero
1685  select case (obsrv%ObsTypeId)
1686  case default
1687  msg = 'Unrecognized observation type: '//trim(obsrv%ObsTypeId)
1688  call store_error(msg)
1689  end select
1690  call this%obs%SaveOneSimval(obsrv, v)
1691  end do
1692  end do
1693 
1694  ! write summary of package error messages
1695  if (count_errors() > 0) then
1696  call store_error_filename(this%input_fname)
1697  end if
1698  end if
1699 
1700  end subroutine dfw_bd_obs
1701 
1702  !> @brief Read and prepare observations for a package
1703  !!
1704  !< Method to read and prepare observations for a package.
1705  subroutine dfw_rp_obs(this)
1706  ! modules
1707  use tdismodule, only: kper
1708  ! dummy
1709  class(SwfDfwType), intent(inout) :: this !< this instance
1710  ! local
1711  integer(I4B) :: i
1712  integer(I4B) :: j
1713  integer(I4B) :: nn1
1714  class(ObserveType), pointer :: obsrv => null()
1715 
1716  ! process each package observation
1717  ! only done the first stress period since boundaries are fixed
1718  ! for the simulation
1719  if (kper == 1) then
1720  do i = 1, this%obs%npakobs
1721  obsrv => this%obs%pakobs(i)%obsrv
1722 
1723  ! get node number 1
1724  nn1 = obsrv%NodeNumber
1725  if (nn1 < 1 .or. nn1 > this%dis%nodes) then
1726  write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
1727  trim(adjustl(obsrv%ObsTypeId)), &
1728  'reach must be greater than 0 and less than or equal to', &
1729  this%dis%nodes, '(specified value is ', nn1, ')'
1730  call store_error(errmsg)
1731  else
1732  if (obsrv%indxbnds_count == 0) then
1733  call obsrv%AddObsIndex(nn1)
1734  else
1735  errmsg = 'Programming error in dfw_rp_obs'
1736  call store_error(errmsg)
1737  end if
1738  end if
1739 
1740  ! check that node number 1 is valid; call store_error if not
1741  do j = 1, obsrv%indxbnds_count
1742  nn1 = obsrv%indxbnds(j)
1743  if (nn1 < 1 .or. nn1 > this%dis%nodes) then
1744  write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') &
1745  trim(adjustl(obsrv%ObsTypeId)), &
1746  'reach must be greater than 0 and less than or equal to', &
1747  this%dis%nodes, '(specified value is ', nn1, ')'
1748  call store_error(errmsg)
1749  end if
1750  end do
1751  end do
1752 
1753  ! evaluate if there are any observation errors
1754  if (count_errors() > 0) then
1755  call store_error_filename(this%input_fname)
1756  end if
1757 
1758  end if
1759 
1760  end subroutine dfw_rp_obs
1761 
1762 end module swfdfwmodule
This module contains the BudgetModule.
Definition: Budget.f90:20
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
real(dp), parameter dtwothirds
real constant 2/3
Definition: Constants.f90:70
real(dp), parameter dp9
real constant 9/10
Definition: Constants.f90:72
real(dp), parameter dem10
real constant 1e-10
Definition: Constants.f90:113
real(dp), parameter donethird
real constant 1/3
Definition: Constants.f90:67
integer(i4b), parameter lenbigline
maximum length of a big line
Definition: Constants.f90:15
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter dprec
real constant machine precision
Definition: Constants.f90:120
real(dp), parameter dtwo
real constant 2
Definition: Constants.f90:79
integer(i4b), parameter lenmempath
maximum length of the memory path
Definition: Constants.f90:27
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
integer(i4b) function, public getunit()
Get a free unit number.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
This module defines variable data types.
Definition: kind.f90:8
real(dp) function, public get_perturbation(x)
Calculate a numerical perturbation given the value of x.
Definition: MathUtil.f90:372
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
subroutine, public memorystore_remove(component, subcomponent, context)
subroutine, public get_isize(name, mem_path, isize)
@ brief Get the number of elements for this variable
This module contains the base numerical package type.
This module contains the derived types ObserveType and ObsDataType.
Definition: Observe.f90:15
This module contains the derived type ObsType.
Definition: Obs.f90:127
subroutine, public obs_cr(obs, inobs)
@ brief Create a new ObsType object
Definition: Obs.f90:225
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
Definition: Sim.f90:203
subroutine, public store_error_unit(iunit, terminate)
Store the file unit number.
Definition: Sim.f90:168
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
character(len=linelength) idm_context
subroutine squadratic(x, range, dydx, y)
@ brief sQuadratic
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
subroutine, public vector_interpolation_2d(dis, flowja, nedges, nodedge, propsedge, vcomp, vmag, flowareaja)
Interpolate 2D vector components at cell center.
Derived type for the Budget object.
Definition: Budget.f90:39
This class is used to store a single deferred-length character string. It was designed to work in an ...
Definition: CharString.f90:23