MODFLOW 6  version 6.8.0.dev0
USGS Modular Hydrologic Model
swf-cxs.f90
Go to the documentation of this file.
1 ! The SwfCxsType package is assigned for a model
2 ! and can be used to calculate wetted area, wetted
3 ! perimeter, hydraulic radius, composite roughness, etc.
4 ! even if the user doesn't specify a CXS Package.
6 
7  use kindmodule, only: dp, i4b, lgp
11  use simvariablesmodule, only: errmsg
14  use basedismodule, only: disbasetype
15 
16  implicit none
17  private
18  public :: swfcxstype, cxs_cr
19 
20  !> @brief flags indicating which CXS options were found in the input
22  logical(LGP) :: iprpak = .false.
23  end type cxsoptionsfoundtype
24 
25  !> @brief flags indicating which CXS dimensions were found in the input
27  logical(LGP) :: nsections = .false.
28  logical(LGP) :: npoints = .false.
29  end type cxsdimensionsfoundtype
30 
31  !> @brief flags indicating which CXS packagedata were found in the input
33  logical(LGP) :: idcxs = .false.
34  logical(LGP) :: nxspoints = .false.
36 
37  !> @brief flags indicating which CXS crosssectiondata were found in the input
39  logical(LGP) :: ifno = .false.
40  logical(LGP) :: xfraction = .false.
41  logical(LGP) :: height = .false.
42  logical(LGP) :: manfraction = .false.
44 
45  type, extends(numericalpackagetype) :: swfcxstype
46 
47  ! provided as input
48  integer(I4B), pointer :: nsections => null() !< number of cross section
49  integer(I4B), pointer :: npoints => null() !< total number of cross-section points
50  integer(I4B), dimension(:), pointer, contiguous :: idcxs => null() !< cross section id number, size nsections
51  integer(I4B), dimension(:), pointer, contiguous :: nxspoints => null() !< number of cross section points for section, size nsections
52  real(dp), dimension(:), pointer, contiguous :: xfraction => null() !< cross-section relative x distance, of size npoints
53  real(dp), dimension(:), pointer, contiguous :: height => null() !< cross-section heights, of size npoints
54  real(dp), dimension(:), pointer, contiguous :: manfraction => null() !< cross-section roughness data, of size npoints
55 
56  ! calculated from input
57  integer(I4B), dimension(:), pointer, contiguous :: iacross => null() !< pointers to cross-section data for each section, of size nsections + 1
58 
59  contains
60 
61  procedure :: allocate_scalars
62  procedure :: allocate_arrays
63  procedure :: source_options
64  procedure :: log_options
65  procedure :: source_dimensions
66  procedure :: log_dimensions
67  procedure :: source_packagedata
68  procedure :: log_packagedata
69  procedure :: check_packagedata
72  procedure :: log_crosssectiondata
73  procedure :: cxs_da
75  procedure :: get_area
76  procedure :: get_wetted_perimeter => cxs_wetted_perimeter
77  procedure :: get_roughness
78  procedure :: get_conveyance => cxs_conveyance
79  procedure :: get_hydraulic_radius
80  procedure :: get_wetted_top_width
81  procedure :: get_maximum_top_width
82  procedure :: write_cxs_table
83 
84  end type swfcxstype
85 
86 contains
87 
88  !> @brief create package
89  !<
90  subroutine cxs_cr(pobj, name_model, input_mempath, inunit, iout, dis)
91  ! -- modules
93  ! -- dummy
94  type(swfcxstype), pointer :: pobj
95  character(len=*), intent(in) :: name_model
96  character(len=*), intent(in) :: input_mempath
97  integer(I4B), intent(in) :: inunit
98  integer(I4B), intent(in) :: iout
99  class(disbasetype), pointer, intent(inout) :: dis !< the pointer to the discretization
100  ! -- locals
101  logical(LGP) :: found_fname
102  ! -- formats
103  character(len=*), parameter :: fmtheader = &
104  "(1x, /1x, 'CXS -- CROSS SECTION PACKAGE, VERSION 1, 5/24/2023', &
105  &' INPUT READ FROM MEMPATH: ', A, /)"
106  !
107  ! -- Create the object
108  allocate (pobj)
109 
110  ! -- create name and memory path
111  call pobj%set_names(1, name_model, 'CXS', 'CXS')
112 
113  ! -- Allocate scalars
114  call pobj%allocate_scalars()
115 
116  ! -- Set variables
117  pobj%input_mempath = input_mempath
118  pobj%inunit = inunit
119  pobj%iout = iout
120  pobj%dis => dis
121 
122  ! -- set name of input file
123  call mem_set_value(pobj%input_fname, 'INPUT_FNAME', pobj%input_mempath, &
124  found_fname)
125 
126  ! -- check if package is enabled
127  if (inunit > 0) then
128 
129  ! -- Print a message identifying the package.
130  write (iout, fmtheader) input_mempath
131 
132  ! -- source options
133  call pobj%source_options()
134 
135  ! -- source dimensions
136  call pobj%source_dimensions()
137 
138  ! -- allocate arrays
139  call pobj%allocate_arrays()
140 
141  ! -- source dimensions
142  call pobj%source_packagedata()
143 
144  ! -- source dimensions
145  call pobj%source_crosssectiondata()
146 
147  end if
148  end subroutine cxs_cr
149 
150  !> @ brief Allocate scalars
151  !!
152  !! Allocate and initialize scalars for the package. The base model
153  !! allocate scalars method is also called.
154  !!
155  !<
156  subroutine allocate_scalars(this)
157  ! -- modules
158  ! -- dummy
159  class(swfcxstype) :: this
160  !
161  ! -- allocate scalars in NumericalPackageType
162  call this%NumericalPackageType%allocate_scalars()
163  !
164  ! -- Allocate scalars
165  call mem_allocate(this%nsections, 'NSECTIONS', this%memoryPath)
166  call mem_allocate(this%npoints, 'NPOINTS', this%memoryPath)
167 
168  ! -- initialize
169  this%nsections = 0
170  this%npoints = 0
171  end subroutine allocate_scalars
172 
173  !> @brief Copy options from IDM into package
174  !<
175  subroutine source_options(this)
176  ! -- modules
177  use kindmodule, only: lgp
180  ! -- dummy
181  class(swfcxstype) :: this
182  ! -- locals
183  character(len=LENMEMPATH) :: idmMemoryPath
184  type(cxsoptionsfoundtype) :: found
185  !
186  ! -- set memory path
187  idmmemorypath = create_mem_path(this%name_model, 'CXS', idm_context)
188  !
189  ! -- update defaults with idm sourced values
190  call mem_set_value(this%iprpak, 'PRINT_INPUT', idmmemorypath, &
191  found%iprpak)
192  !
193  ! -- log values to list file
194  if (this%iout > 0) then
195  call this%log_options(found)
196  end if
197  end subroutine source_options
198 
199  !> @brief Write user options to list file
200  !<
201  subroutine log_options(this, found)
202  class(swfcxstype) :: this
203  type(cxsoptionsfoundtype), intent(in) :: found
204 
205  write (this%iout, '(1x,a)') 'Setting CXS Options'
206 
207  if (found%iprpak) then
208  write (this%iout, '(4x,a)') 'Package information will be printed.'
209  end if
210 
211  write (this%iout, '(1x,a,/)') 'End Setting CXS Options'
212 
213  end subroutine log_options
214 
215  !> @brief Copy options from IDM into package
216  !<
217  subroutine source_dimensions(this)
218  ! -- modules
219  use kindmodule, only: lgp
222  ! -- dummy
223  class(swfcxstype) :: this
224  ! -- locals
225  character(len=LENMEMPATH) :: idmMemoryPath
226  type(cxsdimensionsfoundtype) :: found
227  !
228  ! -- set memory path
229  idmmemorypath = create_mem_path(this%name_model, 'CXS', idm_context)
230  !
231  ! -- update defaults with idm sourced values
232  call mem_set_value(this%nsections, 'NSECTIONS', idmmemorypath, &
233  found%nsections)
234  call mem_set_value(this%npoints, 'NPOINTS', idmmemorypath, &
235  found%npoints)
236  !
237  ! -- ensure nsections was found
238  if (.not. found%nsections) then
239  write (errmsg, '(a)') 'Error in DIMENSIONS block: NSECTIONS not found.'
240  call store_error(errmsg)
241  end if
242  !
243  ! -- ensure npoints was found
244  if (.not. found%npoints) then
245  write (errmsg, '(a)') 'Error in DIMENSIONS block: NPOINTS not found.'
246  call store_error(errmsg)
247  end if
248  !
249  ! -- log values to list file
250  if (this%iout > 0) then
251  call this%log_dimensions(found)
252  end if
253  end subroutine source_dimensions
254 
255  !> @brief Write user options to list file
256  !<
257  subroutine log_dimensions(this, found)
258  class(swfcxstype) :: this
259  type(cxsdimensionsfoundtype), intent(in) :: found
260 
261  write (this%iout, '(1x,a)') 'Setting CXS Dimensions'
262 
263  if (found%nsections) then
264  write (this%iout, '(4x,a)') 'NSECTIONS set from input file.'
265  end if
266 
267  if (found%npoints) then
268  write (this%iout, '(4x,a)') 'NPOINTS set from input file.'
269  end if
270 
271  write (this%iout, '(1x,a,/)') 'End Setting CXS Dimensions'
272 
273  end subroutine log_dimensions
274 
275  !> @brief allocate memory for arrays
276  !<
277  subroutine allocate_arrays(this)
278  ! -- dummy
279  class(swfcxstype) :: this
280  ! -- locals
281  integer(I4B) :: n
282  !
283  ! -- arrays allocation
284  call mem_allocate(this%idcxs, this%nsections, &
285  'IDCXS', this%memoryPath)
286  call mem_allocate(this%nxspoints, this%nsections, &
287  'NXSPOINTS', this%memoryPath)
288  call mem_allocate(this%xfraction, this%npoints, &
289  'XFRACTION', this%memoryPath)
290  call mem_allocate(this%height, this%npoints, &
291  'HEIGHT', this%memoryPath)
292  call mem_allocate(this%manfraction, this%npoints, &
293  'MANFRACTION', this%memoryPath)
294  call mem_allocate(this%iacross, this%nsections + 1, &
295  'IACROSS', this%memoryPath)
296 
297  ! -- initialization
298  do n = 1, this%nsections
299  this%idcxs(n) = 0
300  this%nxspoints(n) = 0
301  end do
302  do n = 1, this%npoints
303  this%xfraction(n) = dzero
304  this%height(n) = dzero
305  this%manfraction(n) = dzero
306  end do
307  do n = 1, this%nsections + 1
308  this%iacross(n) = 0
309  end do
310  end subroutine allocate_arrays
311 
312  !> @brief Copy options from IDM into package
313  !<
314  subroutine source_packagedata(this)
315  ! modules
316  use kindmodule, only: lgp
319  ! dummy
320  class(swfcxstype) :: this
321  ! locals
322  character(len=LENMEMPATH) :: idmMemoryPath
323  type(cxspackagedatafoundtype) :: found
324 
325  ! set memory path
326  idmmemorypath = create_mem_path(this%name_model, 'CXS', idm_context)
327 
328  ! update defaults with idm sourced values
329  call mem_set_value(this%idcxs, 'IFNO_PKGDATA', idmmemorypath, &
330  found%idcxs)
331  call mem_set_value(this%nxspoints, 'NXSPOINTS', idmmemorypath, &
332  found%nxspoints)
333 
334  ! ensure idcxs was found
335  if (.not. found%idcxs) then
336  write (errmsg, '(a)') 'Error in PACKAGEDATA block: IFNO not found.'
337  call store_error(errmsg)
338  end if
339 
340  ! ensure nxspoints was found
341  if (.not. found%nxspoints) then
342  write (errmsg, '(a)') 'Error in PACKAGEDATA block: NXSPOINTS not found.'
343  call store_error(errmsg)
344  end if
345 
346  ! log values to list file
347  if (this%iout > 0) then
348  call this%log_packagedata(found)
349  end if
350 
351  ! Check to make sure package data is valid
352  call this%check_packagedata()
353 
354  ! Calculate the iacross index array using nxspoints
355  call calc_iacross(this%nxspoints, this%iacross)
356  end subroutine source_packagedata
357 
358  !> @brief Calculate index pointer array iacross from nxspoints
359  !<
360  subroutine calc_iacross(nxspoints, iacross)
361  integer(I4B), dimension(:), intent(in) :: nxspoints
362  integer(I4B), dimension(:), intent(inout) :: iacross
363  integer(I4B) :: n
364  iacross(1) = 1
365  do n = 1, size(nxspoints)
366  iacross(n + 1) = iacross(n) + nxspoints(n)
367  end do
368  end subroutine calc_iacross
369 
370  !> @brief Check packagedata
371  !<
372  subroutine check_packagedata(this)
373  ! dummy arguments
374  class(swfcxstype) :: this !< this instance
375  ! local variables
376  integer(I4B) :: i
377 
378  ! Check that all cross section IDs are in range
379  do i = 1, size(this%idcxs)
380  if (this%idcxs(i) <= 0 .or. this%idcxs(i) > this%nsections) then
381  write (errmsg, '(a, i0, a)') &
382  'IFNO values must be greater than 0 and less than NSECTIONS. &
383  &Found ', this%idcxs(i), '.'
384  call store_error(errmsg)
385  end if
386  end do
387 
388  ! Check that nxspoints are greater than one
389  do i = 1, size(this%nxspoints)
390  if (this%nxspoints(i) <= 1) then
391  write (errmsg, '(a, i0, a, i0, a)') &
392  'NXSPOINTS values must be greater than 1 for each cross section. &
393  &Found ', this%nxspoints(i), ' for cross section ', this%idcxs(i), '.'
394  call store_error(errmsg)
395  end if
396  end do
397 
398  ! write summary of package error messages
399  if (count_errors() > 0) then
400  call store_error_filename(this%input_fname)
401  end if
402 
403  end subroutine check_packagedata
404 
405  !> @brief Write user packagedata to list file
406  !<
407  subroutine log_packagedata(this, found)
408  class(swfcxstype) :: this
409  type(cxspackagedatafoundtype), intent(in) :: found
410 
411  write (this%iout, '(1x,a)') 'Setting CXS Package Data'
412 
413  if (found%idcxs) then
414  write (this%iout, '(4x,a)') 'IFNO set from input file.'
415  end if
416 
417  if (found%nxspoints) then
418  write (this%iout, '(4x,a)') 'NXSPOINTS set from input file.'
419  end if
420 
421  write (this%iout, '(1x,a,/)') 'End Setting CXS Package Data'
422 
423  end subroutine log_packagedata
424 
425  !> @brief Copy options from IDM into package
426  !<
427  subroutine source_crosssectiondata(this)
428  ! -- modules
429  use kindmodule, only: lgp
432  ! -- dummy
433  class(swfcxstype) :: this
434  ! -- locals
435  character(len=LENMEMPATH) :: idmMemoryPath
436  type(cxscrosssectiondatafoundtype) :: found
437  integer(I4B), dimension(:), pointer, contiguous :: ifno => null()
438  !
439  ! -- set memory path
440  idmmemorypath = create_mem_path(this%name_model, 'CXS', idm_context)
441  !
442  ! -- update defaults with idm sourced values
443  allocate (ifno(this%npoints))
444  call mem_set_value(ifno, 'IFNO', idmmemorypath, &
445  found%ifno)
446  call mem_set_value(this%xfraction, 'XFRACTION', idmmemorypath, &
447  found%xfraction)
448  call mem_set_value(this%height, 'HEIGHT', idmmemorypath, &
449  found%height)
450  call mem_set_value(this%manfraction, 'MANFRACTION', idmmemorypath, &
451  found%manfraction)
452  !
453  ! -- ensure ifno was found
454  if (.not. found%ifno) then
455  write (errmsg, '(a)') &
456  'Error in CROSSSECTIONDATA block: IFNO not found.'
457  call store_error(errmsg)
458  end if
459  !
460  ! -- ensure xfraction was found
461  if (.not. found%xfraction) then
462  write (errmsg, '(a)') &
463  'Error in CROSSSECTIONDATA block: xfraction not found.'
464  call store_error(errmsg)
465  end if
466  !
467  ! -- ensure height was found
468  if (.not. found%height) then
469  write (errmsg, '(a)') &
470  'Error in CROSSSECTIONDATA block: HEIGHT not found.'
471  call store_error(errmsg)
472  end if
473  !
474  ! -- ensure manfraction was found
475  if (.not. found%manfraction) then
476  write (errmsg, '(a)') &
477  'Error in CROSSSECTIONDATA block: MANFRACTION not found.'
478  call store_error(errmsg)
479  end if
480  !
481  ! -- check that ifno is consistent with the section groupings
482  ! implied by PACKAGEDATA's NXSPOINTS
483  if (found%ifno) then
484  call this%check_crosssectiondata(ifno)
485  end if
486  deallocate (ifno)
487  !
488  ! -- log values to list file
489  if (this%iout > 0) then
490  call this%log_crosssectiondata(found)
491  end if
492  end subroutine source_crosssectiondata
493 
494  !> @brief Check crosssectiondata IFNO against PACKAGEDATA section groupings
495  !<
496  subroutine check_crosssectiondata(this, ifno)
497  ! -- dummy
498  class(swfcxstype) :: this !< this instance
499  integer(I4B), dimension(:), intent(in) :: ifno !< cross section number for each point
500  ! -- local
501  integer(I4B) :: n
502  integer(I4B) :: i
503  integer(I4B) :: i0
504  integer(I4B) :: i1
505 
506  do n = 1, this%nsections
507  i0 = this%iacross(n)
508  i1 = this%iacross(n + 1) - 1
509  do i = i0, i1
510  if (ifno(i) /= this%idcxs(n)) then
511  write (errmsg, '(a, i0, a, i0, a, i0, a)') &
512  'IFNO value in CROSSSECTIONDATA does not match the cross &
513  &section implied by PACKAGEDATA. Found IFNO = ', ifno(i), &
514  ' at point ', i, ' but expected ', this%idcxs(n), '.'
515  call store_error(errmsg)
516  end if
517  end do
518  end do
519 
520  if (count_errors() > 0) then
521  call store_error_filename(this%input_fname)
522  end if
523  end subroutine check_crosssectiondata
524 
525  !> @brief Write user packagedata to list file
526  !<
527  subroutine log_crosssectiondata(this, found)
528  class(swfcxstype) :: this
529  type(cxscrosssectiondatafoundtype), intent(in) :: found
530 
531  write (this%iout, '(1x,a)') 'Setting CXS Cross Section Data'
532 
533  if (found%ifno) then
534  write (this%iout, '(4x,a)') 'IFNO set from input file.'
535  end if
536 
537  if (found%xfraction) then
538  write (this%iout, '(4x,a)') 'XFRACTION set from input file.'
539  end if
540 
541  if (found%height) then
542  write (this%iout, '(4x,a)') 'HEIGHT set from input file.'
543  end if
544 
545  if (found%manfraction) then
546  write (this%iout, '(4x,a)') 'MANFRACTION set from input file.'
547  end if
548 
549  write (this%iout, '(1x,a,/)') 'End Setting CXS Cross Section Data'
550 
551  end subroutine log_crosssectiondata
552 
553  subroutine write_cxs_table(this, idcxs, width, slope, rough, unitconv)
554  ! -- module
555  use sortmodule, only: qsort, unique_values
556  ! -- dummy
557  class(swfcxstype) :: this
558  integer(I4B), intent(in) :: idcxs
559  real(DP), intent(in) :: width
560  real(DP), intent(in) :: slope
561  real(DP), intent(in) :: rough
562  real(DP), intent(in) :: unitconv
563  ! -- local
564  integer(I4B) :: ipt
565  real(DP) :: d
566  real(DP) :: a
567  real(DP) :: rh
568  real(DP) :: wp
569  real(DP) :: r
570  real(DP) :: c
571  real(DP) :: q
572  integer(I4B) :: i0
573  integer(I4B) :: i1
574  integer(I4B) :: npts
575  integer(I4B) :: icalcmeth
576  real(DP), dimension(:), allocatable :: depths
577  real(DP), dimension(:), allocatable :: depths_unique
578  integer(I4B), dimension(:), allocatable :: indx
579 
580  call this%get_cross_section_info(idcxs, i0, i1, npts, icalcmeth)
581 
582  if (npts > 0) then
583 
584  write (this%iout, *) 'Processing information for cross section ', idcxs
585  write (this%iout, *) 'Depth Area WettedP HydRad Rough Conveyance Q'
586 
587  allocate (depths(npts))
588  allocate (indx(size(depths)))
589 
590  depths(:) = this%height(:)
591  call qsort(indx, depths)
592  call unique_values(depths, depths_unique)
593 
594  do ipt = 1, size(depths_unique)
595  d = depths_unique(ipt)
596  a = this%get_area(idcxs, width, d)
597  wp = this%get_wetted_perimeter(idcxs, width, d)
598  rh = this%get_hydraulic_radius(idcxs, width, d, a)
599  r = this%get_roughness(idcxs, width, d, rough, slope)
600  c = this%get_conveyance(idcxs, width, d, rough)
601  if (slope > dzero) then
602  q = unitconv * c * sqrt(slope)
603  else
604  q = dzero
605  end if
606  write (this%iout, *) d, a, wp, rh, r, c, q
607  end do
608 
609  deallocate (depths)
610  deallocate (depths_unique)
611  write (this%iout, *) 'Done processing information for cross section ', idcxs
612 
613  end if
614  end subroutine write_cxs_table
615 
616  !> @brief deallocate memory
617  !<
618  subroutine cxs_da(this)
619  ! -- modules
623  ! -- dummy
624  class(swfcxstype) :: this
625  !
626  ! -- Deallocate input memory
627  call memorystore_remove(this%name_model, 'CXS', idm_context)
628  !
629  ! -- Scalars
630  call mem_deallocate(this%nsections)
631  call mem_deallocate(this%npoints)
632  !
633  ! -- Deallocate arrays if the package was created
634  ! from an input file
635  if (this%inunit > 0) then
636  call mem_deallocate(this%idcxs)
637  call mem_deallocate(this%nxspoints)
638  call mem_deallocate(this%xfraction)
639  call mem_deallocate(this%height)
640  call mem_deallocate(this%manfraction)
641  call mem_deallocate(this%iacross)
642  end if
643  !
644  ! -- deallocate parent
645  call this%NumericalPackageType%da()
646  end subroutine cxs_da
647 
648  subroutine get_cross_section_info(this, idcxs, i0, i1, npts, icalcmeth)
649  ! -- dummy
650  class(swfcxstype) :: this
651  integer(I4B), intent(in) :: idcxs !< cross section id number
652  integer(I4B), intent(inout) :: i0 !< starting cross section point number
653  integer(I4B), intent(inout) :: i1 !< ending cross section point number
654  integer(I4B), intent(inout) :: npts !< number of points in cross section
655  integer(I4B), intent(inout) :: icalcmeth !< calculation method for mannings roughness
656  ! -- local
657  !
658  ! -- Return npts = 0 if this package does not have input file
659  if (this%inunit == 0 .or. idcxs == 0) then
660  npts = 0
661  i0 = 1
662  i1 = 1
663  icalcmeth = 0
664  else
665  !
666  ! -- If the cross section id is 0, then it is a hydraulically wide channel,
667  ! and only width and rough are needed (not xfraction, height, and manfraction)
668  if (idcxs > 0) then
669  i0 = this%iacross(idcxs)
670  i1 = this%iacross(idcxs + 1) - 1
671  else
672  i0 = 1
673  i1 = 1
674  end if
675  ! set icalcmeth based on number of cross section points
676  npts = i1 - i0 + 1
677  icalcmeth = 0 ! linear composite mannings resistance
678  if (npts > 4) then
679  icalcmeth = 0 ! sum q by cross section segments
680  end if
681  end if
682  end subroutine get_cross_section_info
683 
684  function get_area(this, idcxs, width, depth) result(area)
685  ! -- modules
687  ! -- dummy
688  class(swfcxstype) :: this
689  integer(I4B), intent(in) :: idcxs !< cross section id
690  real(dp), intent(in) :: width !< width in reach
691  real(dp), intent(in) :: depth !< stage in reach
692  ! -- local
693  real(dp) :: area
694  integer(I4B) :: i0
695  integer(I4B) :: i1
696  integer(I4B) :: npts
697  integer(I4B) :: icalcmeth
698  call this%get_cross_section_info(idcxs, i0, i1, npts, icalcmeth)
699  if (npts == 0) then
700  area = width * depth
701  else
702  area = get_cross_section_area(npts, &
703  this%xfraction(i0:i1), &
704  this%height(i0:i1), &
705  width, depth)
706  end if
707  end function get_area
708 
709  function cxs_wetted_perimeter(this, idcxs, width, depth) result(wp)
710  ! -- modules
712  ! -- dummy
713  class(swfcxstype) :: this
714  integer(I4B), intent(in) :: idcxs !< cross section id
715  real(dp), intent(in) :: width !< width in reach
716  real(dp), intent(in) :: depth !< stage in reach
717  ! -- local
718  real(dp) :: wp
719  integer(I4B) :: i0
720  integer(I4B) :: i1
721  integer(I4B) :: npts
722  integer(I4B) :: icalcmeth
723  call this%get_cross_section_info(idcxs, i0, i1, npts, icalcmeth)
724  if (npts == 0) then
725  wp = width
726  else
727  wp = get_wetted_perimeter(npts, &
728  this%xfraction(i0:i1), &
729  this%height(i0:i1), &
730  width, depth)
731  end if
732  end function cxs_wetted_perimeter
733 
734  function get_roughness(this, idcxs, width, depth, rough, &
735  slope) result(roughc)
736  ! -- modules
738  ! -- dummy
739  class(swfcxstype) :: this
740  integer(I4B), intent(in) :: idcxs !< cross section id
741  real(dp), intent(in) :: width !< width in reach
742  real(dp), intent(in) :: depth !< stage in reach
743  real(dp), intent(in) :: rough !< mannings value provided for the reach
744  real(dp), intent(in) :: slope !< slope value provided for the reach
745  ! -- local
746  real(dp) :: roughc !< calculated composite roughness
747  integer(I4B) :: i0
748  integer(I4B) :: i1
749  integer(I4B) :: npts
750  integer(I4B) :: icalcmeth
751  call this%get_cross_section_info(idcxs, i0, i1, npts, icalcmeth)
752  if (npts == 0) then
753  roughc = rough
754  else
755  roughc = calc_composite_roughness(npts, &
756  depth, &
757  width, &
758  rough, &
759  slope, &
760  this%xfraction(i0:i1), &
761  this%height(i0:i1), &
762  this%manfraction(i0:i1), &
763  icalcmeth)
764  end if
765  end function get_roughness
766 
767  !> @brief Calculate and return conveyance
768  !!
769  !! Conveyance = area * hydraulic_radius ** (2/3) / mannings_roughness
770  !! If idcxs = 0 (no cross section specified) then reach is
771  !< hydraulically wide and hydraulic radius is equal to depth.
772  function cxs_conveyance(this, idcxs, width, depth, &
773  rough) result(conveyance)
774  ! -- modules
776  ! -- dummy
777  class(swfcxstype) :: this
778  integer(I4B), intent(in) :: idcxs !< cross section id
779  real(dp), intent(in) :: width !< width in reach
780  real(dp), intent(in) :: depth !< stage in reach
781  real(dp), intent(in) :: rough !< mannings value provided for the reach
782  ! -- return
783  real(dp) :: conveyance !< calculated composite roughness
784  ! -- local
785  real(dp) :: a
786  real(dp) :: rh
787  integer(I4B) :: i0
788  integer(I4B) :: i1
789  integer(I4B) :: npts
790  integer(I4B) :: icalcmeth
791  call this%get_cross_section_info(idcxs, i0, i1, npts, icalcmeth)
792  if (npts == 0) then
793  a = depth * width
794  rh = depth
795  conveyance = a * rh**dtwothirds / rough
796  else
797  conveyance = get_conveyance(npts, &
798  this%xfraction(i0:i1), &
799  this%height(i0:i1), &
800  this%manfraction(i0:i1), &
801  width, rough, depth)
802  end if
803  end function cxs_conveyance
804 
805  function get_hydraulic_radius(this, idcxs, width, depth, area) result(r)
806  ! -- modules
808  ! -- dummy
809  class(swfcxstype) :: this
810  integer(I4B), intent(in) :: idcxs !< cross section id
811  real(dp), intent(in) :: width !< width in reach
812  real(dp), intent(in) :: depth !< stage in reach
813  real(dp), intent(in), optional :: area !< area of the reach
814  ! -- local
815  real(dp) :: r !< calculated hydraulic radius
816  real(dp) :: a
817  integer(I4B) :: i0
818  integer(I4B) :: i1
819  integer(I4B) :: npts
820  integer(I4B) :: icalcmeth
821  call this%get_cross_section_info(idcxs, i0, i1, npts, icalcmeth)
822  if (present(area)) then
823  a = area
824  else
825  a = this%get_area(idcxs, width, depth)
826  end if
827  if (npts == 0) then
828  r = a / width
829  else
830  r = get_hydraulic_radius_xf(npts, &
831  this%xfraction(i0:i1), &
832  this%height(i0:i1), &
833  width, depth)
834  end if
835  end function get_hydraulic_radius
836 
837  function get_wetted_top_width(this, idcxs, width, depth) result(r)
838  ! modules
840  ! dummy
841  class(swfcxstype) :: this
842  integer(I4B), intent(in) :: idcxs !< cross section id
843  real(dp), intent(in) :: width !< width in reach
844  real(dp), intent(in) :: depth !< stage in reach
845  ! local
846  real(dp) :: r !< calculated hydraulic radius
847  integer(I4B) :: i0
848  integer(I4B) :: i1
849  integer(I4B) :: npts
850  integer(I4B) :: icalcmeth
851  call this%get_cross_section_info(idcxs, i0, i1, npts, icalcmeth)
852  if (npts == 0) then
853  r = width
854  else
855  r = get_wetted_topwidth(npts, this%xfraction(i0:i1), &
856  this%height(i0:i1), width, depth)
857  end if
858  end function get_wetted_top_width
859 
860  function get_maximum_top_width(this, idcxs, width) result(r)
861  ! modules
863  ! dummy
864  class(swfcxstype) :: this
865  integer(I4B), intent(in) :: idcxs !< cross section id
866  real(dp), intent(in) :: width !< width in reach
867  ! local
868  real(dp) :: r !< calculated hydraulic radius
869  integer(I4B) :: i0
870  integer(I4B) :: i1
871  integer(I4B) :: npts
872  integer(I4B) :: icalcmeth
873  call this%get_cross_section_info(idcxs, i0, i1, npts, icalcmeth)
874  if (npts == 0) then
875  r = width
876  else
877  r = get_saturated_topwidth(npts, this%xfraction(i0:i1), width)
878  end if
879  end function get_maximum_top_width
880 
881 end module swfcxsmodule
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter dtwothirds
real constant 2/3
Definition: Constants.f90:70
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b), parameter lenmempath
maximum length of the memory path
Definition: Constants.f90:27
This module defines variable data types.
Definition: kind.f90:8
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
subroutine, public memorystore_remove(component, subcomponent, context)
This module contains the base numerical package type.
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
subroutine, public store_error_filename(filename, terminate)
Store the erroring file name.
Definition: Sim.f90:203
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
character(len=linelength) idm_context
real(dp) function get_wetted_top_width(this, idcxs, width, depth)
Definition: swf-cxs.f90:838
subroutine source_packagedata(this)
Copy options from IDM into package.
Definition: swf-cxs.f90:315
real(dp) function cxs_wetted_perimeter(this, idcxs, width, depth)
Definition: swf-cxs.f90:710
subroutine check_crosssectiondata(this, ifno)
Check crosssectiondata IFNO against PACKAGEDATA section groupings.
Definition: swf-cxs.f90:497
subroutine source_options(this)
Copy options from IDM into package.
Definition: swf-cxs.f90:176
real(dp) function get_roughness(this, idcxs, width, depth, rough, slope)
Definition: swf-cxs.f90:736
subroutine allocate_arrays(this)
allocate memory for arrays
Definition: swf-cxs.f90:278
subroutine log_options(this, found)
Write user options to list file.
Definition: swf-cxs.f90:202
subroutine calc_iacross(nxspoints, iacross)
Calculate index pointer array iacross from nxspoints.
Definition: swf-cxs.f90:361
real(dp) function get_maximum_top_width(this, idcxs, width)
Definition: swf-cxs.f90:861
real(dp) function get_hydraulic_radius(this, idcxs, width, depth, area)
Definition: swf-cxs.f90:806
subroutine log_dimensions(this, found)
Write user options to list file.
Definition: swf-cxs.f90:258
subroutine get_cross_section_info(this, idcxs, i0, i1, npts, icalcmeth)
Definition: swf-cxs.f90:649
subroutine check_packagedata(this)
Check packagedata.
Definition: swf-cxs.f90:373
real(dp) function get_area(this, idcxs, width, depth)
Definition: swf-cxs.f90:685
subroutine source_dimensions(this)
Copy options from IDM into package.
Definition: swf-cxs.f90:218
subroutine cxs_da(this)
deallocate memory
Definition: swf-cxs.f90:619
subroutine source_crosssectiondata(this)
Copy options from IDM into package.
Definition: swf-cxs.f90:428
real(dp) function cxs_conveyance(this, idcxs, width, depth, rough)
Calculate and return conveyance.
Definition: swf-cxs.f90:774
subroutine allocate_scalars(this)
@ brief Allocate scalars
Definition: swf-cxs.f90:157
subroutine log_crosssectiondata(this, found)
Write user packagedata to list file.
Definition: swf-cxs.f90:528
subroutine, public cxs_cr(pobj, name_model, input_mempath, inunit, iout, dis)
create package
Definition: swf-cxs.f90:91
subroutine log_packagedata(this, found)
Write user packagedata to list file.
Definition: swf-cxs.f90:408
subroutine write_cxs_table(this, idcxs, width, slope, rough, unitconv)
Definition: swf-cxs.f90:554
This module contains stateless sfr subroutines and functions.
Definition: SwfCxsUtils.f90:11
real(dp) function, public get_hydraulic_radius_xf(npts, xfraction, heights, width, d)
Calculate the hydraulic radius for a reach.
real(dp) function, public calc_composite_roughness(npts, depth, width, rough, slope, cxs_xf, cxs_h, cxs_rf, linmeth)
real(dp) function, public get_cross_section_area(npts, xfraction, heights, width, d)
Calculate the cross-sectional area for a reach.
real(dp) function, public get_wetted_topwidth(npts, xfraction, heights, width, d)
Calculate the wetted top width for a reach.
real(dp) function, public get_conveyance(npts, xfraction, heights, cxs_rf, width, rough, d)
Calculate conveyance.
real(dp) function, public get_saturated_topwidth(npts, xfraction, width)
Calculate the saturated top width for a reach.
real(dp) function, public get_wetted_perimeter(npts, xfraction, heights, width, d)
Calculate the wetted perimeter for a reach.
flags indicating which CXS crosssectiondata were found in the input
Definition: swf-cxs.f90:38
flags indicating which CXS dimensions were found in the input
Definition: swf-cxs.f90:26
flags indicating which CXS options were found in the input
Definition: swf-cxs.f90:21
flags indicating which CXS packagedata were found in the input
Definition: swf-cxs.f90:32