MODFLOW 6  version 6.5.0.dev2
MODFLOW 6 Code Documentation
sfrcrosssectionmanager Module Reference

Data Types

type  sfrcrosssectiontype
 
type  sfrcrosssection
 

Functions/Subroutines

subroutine, public cross_section_cr (this, iout, iprpak, nreaches)
 Create a cross-section object. More...
 
subroutine initialize (this, ncrossptstot, ncrosspts, iacross, station, height, roughfraction)
 Initialize a cross-section object. More...
 
subroutine read_table (this, irch, width, filename)
 Read a cross-section table. More...
 
subroutine validate (this, irch)
 Validate cross-section tables. More...
 
subroutine output (this, widths, roughs, kstp, kper)
 Write cross-section tables. More...
 
integer(i4b) function get_ncrossptstot (this)
 Get the total number of cross-section points. More...
 
subroutine pack (this, ncrossptstot, ncrosspts, iacross, station, height, roughfraction)
 Pack the cross-section object. More...
 
subroutine destroy (this)
 Deallocate the cross-section object. More...
 

Function/Subroutine Documentation

◆ cross_section_cr()

subroutine, public sfrcrosssectionmanager::cross_section_cr ( type(sfrcrosssection), pointer  this,
integer(i4b), intent(in), pointer  iout,
integer(i4b), intent(in), pointer  iprpak,
integer(i4b), intent(in), pointer  nreaches 
)

Subroutine to calculate the maximum top width for a reach using the cross-section station data.

Parameters
thisSfrCrossSection object
[in]ioutmodel listing file
[in]iprpakflag for printing table input data
[in]nreachesnumber of reaches

Definition at line 57 of file SfrCrossSectionManager.f90.

58  ! -- dummy variables
59  type(SfrCrossSection), pointer :: this !< SfrCrossSection object
60  integer(I4B), pointer, intent(in) :: iout !< model listing file
61  integer(I4B), pointer, intent(in) :: iprpak !< flag for printing table input data
62  integer(I4B), pointer, intent(in) :: nreaches !< number of reaches
63  !
64  ! -- check if table already associated and reset if necessary
65  if (associated(this)) then
66  call this%destroy()
67  deallocate (this)
68  nullify (this)
69  end if
70  !
71  ! -- Create the object
72  allocate (this)
73  !
74  ! -- initialize scalars
75  this%iout => iout
76  this%iprpak => iprpak
77  this%nreaches => nreaches
78  !
79  ! -- Return
80  return
Here is the caller graph for this function:

◆ destroy()

subroutine sfrcrosssectionmanager::destroy ( class(sfrcrosssection this)

Subroutine to deallocate the cross-section object.

Parameters
thisSfrCrossSection object

Definition at line 702 of file SfrCrossSectionManager.f90.

703  ! -- dummy variables
704  class(SfrCrossSection) :: this !< SfrCrossSection object
705  ! -- local variables
706  integer(I4B) :: n
707  !
708  ! -- deallocate and nullify pointers
709  deallocate (this%npoints)
710  nullify (this%npoints)
711  do n = 1, this%nreaches
712  deallocate (this%cross_sections(n)%npoints)
713  nullify (this%cross_sections(n)%npoints)
714  deallocate (this%cross_sections(n)%station)
715  nullify (this%cross_sections(n)%station)
716  deallocate (this%cross_sections(n)%height)
717  nullify (this%cross_sections(n)%height)
718  deallocate (this%cross_sections(n)%roughfraction)
719  nullify (this%cross_sections(n)%roughfraction)
720  deallocate (this%cross_sections(n)%valid)
721  nullify (this%cross_sections(n)%valid)
722  end do
723  deallocate (this%cross_sections)
724  nullify (this%cross_sections)
725  !
726  ! -- input table
727  if (associated(this%inputtab)) then
728  call this%inputtab%table_da()
729  deallocate (this%inputtab)
730  nullify (this%inputtab)
731  end if
732  !
733  ! -- deallocate and nullify class scalars
734  deallocate (this%invalid)
735  nullify (this%invalid)
736  !
737  ! -- nullify scalars that are pointers to external variables
738  nullify (this%iout)
739  nullify (this%iprpak)
740  nullify (this%nreaches)
741  !
742  ! -- return
743  return

◆ get_ncrossptstot()

integer(i4b) function sfrcrosssectionmanager::get_ncrossptstot ( class(sfrcrosssection this)

Function to get the total number of cross-section points to get the new size of the station xsheight data for all reaches.

Definition at line 640 of file SfrCrossSectionManager.f90.

641  ! -- dummy variables
642  class(SfrCrossSection) :: this
643  ! -- local variables
644  integer(I4B) :: nptstot
645  integer(I4B) :: n
646  !
647  !
648  nptstot = 0
649  do n = 1, this%nreaches
650  nptstot = nptstot + this%npoints(n)
651  end do
652  !
653  ! -- return
654  return

◆ initialize()

subroutine sfrcrosssectionmanager::initialize ( class(sfrcrosssection this,
integer(i4b), intent(in)  ncrossptstot,
integer(i4b), dimension(this%nreaches), intent(in)  ncrosspts,
integer(i4b), dimension(this%nreaches + 1), intent(in)  iacross,
real(dp), dimension(ncrossptstot), intent(in)  station,
real(dp), dimension(ncrossptstot), intent(in)  height,
real(dp), dimension(ncrossptstot), intent(in)  roughfraction 
)

Subroutine to initialize the cross-section object with the current data.

Parameters
thisSfrCrossSection object
[in]ncrossptstottotal number of cross-section points
[in]ncrossptspointers to cross-section data in data vector
[in]iacrosspointers to cross-section data in data vector
[in]stationcross-section station data
[in]heightcross-section height data
[in]roughfractioncross-section roughness fraction data

Definition at line 89 of file SfrCrossSectionManager.f90.

91  ! -- dummy variables
92  class(SfrCrossSection) :: this !< SfrCrossSection object
93  integer(I4B), intent(in) :: ncrossptstot !< total number of cross-section points
94  integer(I4B), dimension(this%nreaches), intent(in) :: ncrosspts !< pointers to cross-section data in data vector
95  integer(I4B), dimension(this%nreaches + 1), intent(in) :: iacross !< pointers to cross-section data in data vector
96  real(DP), dimension(ncrossptstot), intent(in) :: station !< cross-section station data
97  real(DP), dimension(ncrossptstot), intent(in) :: height !< cross-section height data
98  real(DP), dimension(ncrossptstot), intent(in) :: roughfraction !< cross-section roughness fraction data
99  ! -- local variables
100  integer(I4B) :: i
101  integer(I4B) :: n
102  integer(I4B) :: npoints
103  integer(I4B) :: i0
104  integer(I4B) :: i1
105  integer(I4B) :: ipos
106  !
107  ! -- allocate scalars
108  allocate (this%invalid)
109  !
110  ! -- initialize scalars
111  this%invalid = 0
112  !
113  ! -- create cross-section container
114  allocate (this%filenames(this%nreaches))
115  allocate (this%npoints(this%nreaches))
116  allocate (this%cross_sections(this%nreaches))
117  do n = 1, this%nreaches
118  npoints = ncrosspts(n)
119  allocate (this%cross_sections(n)%npoints)
120  allocate (this%cross_sections(n)%station(npoints))
121  allocate (this%cross_sections(n)%height(npoints))
122  allocate (this%cross_sections(n)%roughfraction(npoints))
123  allocate (this%cross_sections(n)%valid(npoints))
124  end do
125  !
126  ! -- fill cross-section container with current values
127  do n = 1, this%nreaches
128  this%filenames(n) = 'NONE'
129  this%cross_sections(n)%npoints = ncrosspts(n)
130  this%npoints(n) = ncrosspts(n)
131  i0 = iacross(n)
132  i1 = iacross(n + 1) - 1
133  ipos = 1
134  do i = i0, i1
135  this%cross_sections(n)%station(ipos) = station(i)
136  this%cross_sections(n)%height(ipos) = height(i)
137  this%cross_sections(n)%roughfraction(ipos) = roughfraction(i)
138  this%cross_sections(n)%valid(ipos) = .true.
139  ipos = ipos + 1
140  end do
141  end do
142  !
143  ! -- return
144  return

◆ output()

subroutine sfrcrosssectionmanager::output ( class(sfrcrosssection this,
real(dp), dimension(this%nreaches), intent(in)  widths,
real(dp), dimension(this%nreaches), intent(in)  roughs,
integer(i4b), intent(in), optional  kstp,
integer(i4b), intent(in), optional  kper 
)

Subroutine to write the cross-section tables to the model listing file.

Parameters
thisSfrCrossSection object
[in]widthsreach widths
[in]roughsreach Manning's roughness coefficients
[in]kstptime step
[in]kperstress period

Definition at line 475 of file SfrCrossSectionManager.f90.

476  use constantsmodule, only: tableft
477  use simmodule, only: store_warning
478  ! -- dummy variables
479  class(SfrCrossSection) :: this !< SfrCrossSection object
480  real(DP), dimension(this%nreaches), intent(in) :: widths !< reach widths
481  real(DP), dimension(this%nreaches), intent(in) :: roughs !< reach Manning's roughness coefficients
482  integer(I4B), intent(in), optional :: kstp !< time step
483  integer(I4B), intent(in), optional :: kper !< stress period
484  ! -- local variables
485  character(len=LINELENGTH) :: title
486  character(len=LINELENGTH) :: text
487  character(len=LINELENGTH) :: filename
488  character(len=10) :: cvalid
489  logical(LGP) :: transient
490  integer(I4B) :: kkstp
491  integer(I4B) :: kkper
492  integer(I4B) :: irch
493  integer(I4B) :: n
494  integer(I4B) :: ntabcols
495  integer(I4B) :: ninvalid_reaches
496  real(DP) :: width
497  real(DP) :: xfraction
498  real(DP) :: rough
499  real(DP) :: r
500  integer(I4B), dimension(this%nreaches) :: reach_fail
501  !
502  ! -- initialize local variables
503  kkstp = 0
504  kkper = 0
505  !
506  ! -- process optional parameters
507  if (present(kstp)) then
508  kkstp = kstp
509  end if
510  if (present(kper)) then
511  kkper = kper
512  end if
513  !
514  ! -- set transient boolean
515  if (kkstp > 0 .and. kkper > 0) then
516  transient = .true.
517  else
518  transient = .false.
519  end if
520  !
521  ! -- set reach failure
522  do irch = 1, this%nreaches
523  filename = this%filenames(irch)
524  reach_fail(irch) = 0
525  !
526  ! -- output cross-section data read from a file
527  if (trim(adjustl(filename)) /= 'NONE') then
528  do n = 1, this%npoints(irch)
529  if (.not. this%cross_sections(irch)%valid(n)) then
530  reach_fail(irch) = reach_fail(irch) + 1
531  end if
532  end do
533  end if
534  end do
535  !
536  ! -- iterate over each reach
537  do irch = 1, this%nreaches
538  filename = this%filenames(irch)
539  !
540  ! -- output cross-section data read from a file
541  if (trim(adjustl(filename)) /= 'NONE') then
542  !
543  ! -- build and write the table for a reach, if required, or
544  ! the cross-section is invalid
545  if (this%iprpak > 0 .or. reach_fail(irch) > 0) then
546  !
547  ! -- calculate the number of table columns
548  if (reach_fail(irch) > 0) then
549  ntabcols = 6
550  else
551  ntabcols = 5
552  end if
553  !
554  ! -- reset the input table object
555  write (title, '(a,1x,i0,1x,3a)') &
556  'CROSS_SECTION DATA FOR REACH', irch, "FROM TAB6 FILE ('", &
557  trim(adjustl(filename)), "')"
558  call table_cr(this%inputtab, trim(adjustl(filename)), title)
559  call this%inputtab%table_df(this%npoints(irch), ntabcols, &
560  this%iout, finalize=.false., &
561  transient=transient)
562  if (transient) then
563  call this%inputtab%set_kstpkper(kkstp, kkper)
564  end if
565  text = 'XFRACTION'
566  call this%inputtab%initialize_column(text, 20, alignment=tableft)
567  text = 'STATION'
568  call this%inputtab%initialize_column(text, 20, alignment=tableft)
569  text = 'HEIGHT'
570  call this%inputtab%initialize_column(text, 20, alignment=tableft)
571  text = "MANFRACTION"
572  call this%inputtab%initialize_column(text, 20, alignment=tableft)
573  text = "MANNING'S ROUGHNESS COEFFICIENT"
574  call this%inputtab%initialize_column(text, 20, alignment=tableft)
575  if (reach_fail(irch) > 0) then
576  text = 'NEEDS ADJUSTMENT'
577  call this%inputtab%initialize_column(text, 10, alignment=tableft)
578  end if
579  !
580  ! -- set the width and roughness for the reach
581  width = widths(irch)
582  rough = roughs(irch)
583  !
584  ! -- fill the table
585  do n = 1, this%npoints(irch)
586  xfraction = this%cross_sections(irch)%station(n) / width
587  r = this%cross_sections(irch)%roughfraction(n) * rough
588  call this%inputtab%add_term(xfraction)
589  call this%inputtab%add_term(this%cross_sections(irch)%station(n))
590  call this%inputtab%add_term(this%cross_sections(irch)%height(n))
591  call this%inputtab%add_term(&
592  &this%cross_sections(irch)%roughfraction(n))
593  call this%inputtab%add_term(r)
594  if (reach_fail(irch) > 0) then
595  if (this%cross_sections(irch)%valid(n)) then
596  cvalid = ''
597  else
598  cvalid = 'TRUE'
599  end if
600  call this%inputtab%add_term(cvalid)
601  end if
602  end do
603  !
604  ! -- finalize the table
605  call this%inputtab%finalize_table()
606  end if
607  end if
608  end do
609  !
610  ! -- save warning message and write summary information to the listing file
611  if (this%invalid > 0) then
612  ninvalid_reaches = 0
613  do irch = 1, this%nreaches
614  if (reach_fail(irch) > 0) then
615  ninvalid_reaches = ninvalid_reaches + 1
616  end if
617  end do
618  write (warnmsg, '(a,1x,i0,7(1x,a))') &
619  'Cross-section data for', ninvalid_reaches, &
620  'reaches include one or more points that result in a', &
621  'non-unique depth-conveyance relation. This occurs when', &
622  'there are horizontal sections at non-zero heights', &
623  '(for example, flat overbank sections). This can usually', &
624  'be resolved by adding a small slope to these flat', &
625  'sections. See the cross-section tables in the model', &
626  'listing file for more information.'
627  call store_warning(warnmsg)
628  end if
629  !
630  ! -- return
631  return
This module contains simulation constants.
Definition: Constants.f90:9
@ tableft
left justified table column
Definition: Constants.f90:170
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_warning(msg, substring)
Store warning message.
Definition: Sim.f90:236
Here is the call graph for this function:

◆ pack()

subroutine sfrcrosssectionmanager::pack ( class(sfrcrosssection this,
integer(i4b), intent(in)  ncrossptstot,
integer(i4b), dimension(this%nreaches), intent(inout)  ncrosspts,
integer(i4b), dimension(this%nreaches + 1), intent(inout)  iacross,
real(dp), dimension(ncrossptstot), intent(inout)  station,
real(dp), dimension(ncrossptstot), intent(inout)  height,
real(dp), dimension(ncrossptstot), intent(inout)  roughfraction 
)

Subroutine to pack the cross-section object into vectors.

Parameters
thisSfrCrossSection object
[in]ncrossptstottotal number of cross-section points
[in,out]ncrossptspointers to cross-section data in data vector
[in,out]iacrosspointers to cross-section data in data vector
[in,out]stationcross-section station data
[in,out]heightcross-section height data
[in,out]roughfractioncross-section roughness fraction data

Definition at line 662 of file SfrCrossSectionManager.f90.

664  ! -- dummy variables
665  class(SfrCrossSection) :: this !< SfrCrossSection object
666  integer(I4B), intent(in) :: ncrossptstot !< total number of cross-section points
667  integer(I4B), dimension(this%nreaches), intent(inout) :: ncrosspts !< pointers to cross-section data in data vector
668  integer(I4B), dimension(this%nreaches + 1), intent(inout) :: iacross !< pointers to cross-section data in data vector
669  real(DP), dimension(ncrossptstot), intent(inout) :: station !< cross-section station data
670  real(DP), dimension(ncrossptstot), intent(inout) :: height !< cross-section height data
671  real(DP), dimension(ncrossptstot), intent(inout) :: roughfraction !< cross-section roughness fraction data
672  ! -- local variables
673  integer(I4B) :: i
674  integer(I4B) :: n
675  integer(I4B) :: npoints
676  integer(I4B) :: ipos
677  !
678  ! -- pack the data
679  ipos = 1
680  iacross(1) = ipos
681  do n = 1, this%nreaches
682  npoints = this%npoints(n)
683  ncrosspts(n) = npoints
684  do i = 1, npoints
685  station(ipos) = this%cross_sections(n)%station(i)
686  height(ipos) = this%cross_sections(n)%height(i)
687  roughfraction(ipos) = this%cross_sections(n)%roughfraction(i)
688  ipos = ipos + 1
689  end do
690  iacross(n + 1) = ipos
691  end do
692  !
693  ! -- return
694  return

◆ read_table()

subroutine sfrcrosssectionmanager::read_table ( class(sfrcrosssection this,
integer(i4b), intent(in)  irch,
real(dp), intent(in)  width,
character(len=*), intent(in)  filename 
)

Subroutine to read a cross-section table file for a reach.

Parameters
thisSfrCrossSection object
[in]irchcurrent reach
[in]widthreach width
[in]filenametable file with station height data

Definition at line 152 of file SfrCrossSectionManager.f90.

153  use constantsmodule, only: iuoc
154  use inputoutputmodule, only: openfile
155  use simmodule, only: store_error
157  ! -- dummy variables
158  class(SfrCrossSection) :: this !< SfrCrossSection object
159  integer(I4B), intent(in) :: irch !< current reach
160  real(DP), intent(in) :: width !< reach width
161  character(len=*), intent(in) :: filename !< table file with station height data
162  ! -- local variables
163  character(len=LINELENGTH) :: tag
164  character(len=LINELENGTH) :: keyword
165  integer(I4B) :: ierr
166  logical :: isfound
167  logical :: endOfBlock
168  integer(I4B) :: iu
169  integer(I4B) :: n
170  integer(I4B) :: j
171  integer(I4B) :: ipos
172  integer(I4B) :: jmin
173  type(BlockParserType) :: parser
174  !
175  ! -- initialize local variables
176  j = 0
177  n = 0
178  jmin = 2
179  !
180  ! -- create a tag with the file name and reach number
181  write (tag, "('Reach',1x,i0,1x,'(',a, ')')") &
182  irch, trim(adjustl(filename))
183  !
184  ! -- open the table file
185  iu = iuoc
186  call openfile(iu, this%iout, filename, 'SFR TABLE')
187  call parser%Initialize(iu, this%iout)
188  !
189  ! -- get dimensions block
190  call parser%GetBlock('DIMENSIONS', isfound, ierr, supportopenclose=.true.)
191  !
192  ! -- parse table dimensions block if detected
193  if (isfound) then
194  !
195  ! -- process the table dimension data
196  if (this%iprpak /= 0) then
197  write (this%iout, '(/1x,a)') &
198  'PROCESSING '//trim(adjustl(tag))//' DIMENSIONS'
199  end if
200  readdims: do
201  call parser%GetNextLine(endofblock)
202  if (endofblock) exit
203  call parser%GetStringCaps(keyword)
204  select case (keyword)
205  case ('NROW')
206  n = parser%GetInteger()
207  if (n < 1) then
208  write (errmsg, '(a)') 'Table NROW must be > 0'
209  call store_error(errmsg)
210  end if
211  case ('NCOL')
212  j = parser%GetInteger()
213  jmin = 2
214  if (j < jmin) then
215  write (errmsg, '(a,1x,i0)') 'Table NCOL must be >= ', jmin
216  call store_error(errmsg)
217  end if
218  case default
219  write (errmsg, '(a,a)') &
220  'UNKNOWN '//trim(adjustl(tag))//' DIMENSIONS keyword: ', &
221  trim(keyword)
222  call store_error(errmsg)
223  end select
224  end do readdims
225  if (this%iprpak /= 0) then
226  write (this%iout, '(1x,a)') &
227  'END OF '//trim(adjustl(tag))//' DIMENSIONS'
228  end if
229  else
230  call store_error('Required DIMENSIONS block not found.')
231  end if
232  !
233  ! -- check that ncol and nrow have been specified
234  if (n < 1) then
235  write (errmsg, '(a)') &
236  'NROW not specified in the table DIMENSIONS block'
237  call store_error(errmsg)
238  end if
239  if (j < 1) then
240  write (errmsg, '(a)') &
241  'NCOL not specified in the table DIMENSIONS block'
242  call store_error(errmsg)
243  end if
244  !
245  ! -- only read the table data if n and j are specified to be greater
246  ! than zero - an error condition exists if n * j = 0
247  if (n * j > 0) then
248  !
249  ! -- set the filename and reset the number of points
250  this%filenames(irch) = filename
251  this%npoints(irch) = n
252  !
253  ! -- deallocate
254  deallocate (this%cross_sections(irch)%npoints)
255  deallocate (this%cross_sections(irch)%station)
256  deallocate (this%cross_sections(irch)%height)
257  deallocate (this%cross_sections(irch)%roughfraction)
258  deallocate (this%cross_sections(irch)%valid)
259  !
260  ! -- reallocate
261  allocate (this%cross_sections(irch)%npoints)
262  allocate (this%cross_sections(irch)%station(n))
263  allocate (this%cross_sections(irch)%height(n))
264  allocate (this%cross_sections(irch)%roughfraction(n))
265  allocate (this%cross_sections(irch)%valid(n))
266  !
267  ! -- initialize
268  this%cross_sections(irch)%npoints = n
269  !
270  ! -- get table block
271  call parser%GetBlock('TABLE', isfound, ierr, supportopenclose=.true.)
272  !
273  ! -- parse well_connections block if detected
274  if (isfound) then
275 
276  ! -- process the table data
277  if (this%iprpak /= 0) then
278  write (this%iout, '(/1x,a)') &
279  'PROCESSING '//trim(adjustl(tag))//' TABLE'
280  end if
281  ipos = 0
282  readtabledata: do
283  call parser%GetNextLine(endofblock)
284  if (endofblock) exit
285  ipos = ipos + 1
286  if (ipos > this%npoints(irch)) then
287  cycle readtabledata
288  end if
289  this%cross_sections(irch)%station(ipos) = parser%GetDouble() * width
290  this%cross_sections(irch)%height(ipos) = parser%GetDouble()
291  if (j > 2) then
292  this%cross_sections(irch)%roughfraction(ipos) = parser%GetDouble()
293  else
294  this%cross_sections(irch)%roughfraction(ipos) = done
295  end if
296  this%cross_sections(irch)%valid(ipos) = .true.
297  end do readtabledata
298 
299  if (this%iprpak /= 0) then
300  write (this%iout, '(1x,a)') &
301  'END OF '//trim(adjustl(tag))//' TABLE'
302  end if
303  else
304  call store_error('Required TABLE block not found.')
305  end if
306  !
307  ! -- error condition if number of rows read are not equal to nrow
308  if (ipos /= this%npoints(irch)) then
309  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
310  'NROW set to', this%npoints(irch), 'but', ipos, 'rows were read'
311  call store_error(errmsg)
312  end if
313  end if
314  !
315  ! -- close the open table file
316  if (iu > 0) then
317  close (iu)
318  end if
319  !
320  ! -- validate the table
321  call this%validate(irch)
322  !
323  ! -- return
324  return
This module contains block parser methods.
Definition: BlockParser.f90:7
integer(i4b), parameter iuoc
open/close file unit number
Definition: Constants.f90:55
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
Here is the call graph for this function:

◆ validate()

subroutine sfrcrosssectionmanager::validate ( class(sfrcrosssection this,
integer(i4b), intent(in)  irch 
)

Subroutine to validate a cross-section table.

Parameters
thisSfrCrossSection object
[in]irchcurrent reach

Definition at line 332 of file SfrCrossSectionManager.f90.

333  use constantsmodule, only: dem6, dtwothirds
334  use simmodule, only: store_error
335  use sortmodule, only: unique_values
338  ! -- dummy variables
339  class(SfrCrossSection) :: this !< SfrCrossSection object
340  integer(I4B), intent(in) :: irch !< current reach
341  ! -- local variables
342  logical(LGP) :: station_error
343  logical(LGP) :: height_error
344  logical(LGP) :: height_zero_error
345  logical(LGP) :: roughness_error
346  character(len=LINELENGTH) :: filename
347  integer(I4B) :: npts
348  integer(I4B) :: n
349  integer(I4B) :: i
350  integer(I4B) :: ipos
351  real(DP) :: station
352  real(DP) :: height
353  real(DP) :: roughfraction
354  real(DP) :: aw
355  real(DP) :: rh
356  real(DP) :: dc0
357  real(DP) :: dc1
358  real(DP), dimension(:), allocatable :: heights
359  real(DP), dimension(:), allocatable :: unique_heights
360  real(DP), dimension(3) :: factor
361  !
362  ! -- initialize local variables
363  station_error = .false.
364  height_error = .false.
365  height_zero_error = .true.
366  roughness_error = .false.
367  npts = this%npoints(irch)
368  !
369  ! -- validate the station and height data
370  do n = 1, npts
371  station = this%cross_sections(irch)%station(n)
372  if (station < dzero) then
373  station_error = .true.
374  end if
375  height = this%cross_sections(irch)%height(n)
376  if (height < dzero) then
377  height_error = .true.
378  else if (height == dzero) then
379  height_zero_error = .false.
380  end if
381  roughfraction = this%cross_sections(irch)%roughfraction(n)
382  if (roughfraction <= dzero) then
383  roughness_error = .true.
384  end if
385  if (station_error .and. height_error .and. &
386  roughness_error) then
387  exit
388  end if
389  end do
390  !
391  ! -- write error messages
392  if (station_error .or. height_error .or. &
393  height_zero_error .or. roughness_error) then
394  filename = this%filenames(irch)
395  if (station_error) then
396  write (errmsg, '(3a,1x,i0,1x,a)') &
397  "All xfraction data in '", trim(adjustl(filename)), &
398  "' for reach", irch, 'must be greater than or equal to zero.'
399  call store_error(errmsg)
400  end if
401  if (height_error) then
402  write (errmsg, '(3a,1x,i0,1x,a)') &
403  "All height data in '", trim(adjustl(filename)), &
404  "' for reach", irch, 'must be greater than or equal to zero.'
405  call store_error(errmsg)
406  end if
407  if (height_zero_error) then
408  write (errmsg, '(3a,1x,i0,1x,a)') &
409  "At least one height data value in '", trim(adjustl(filename)), &
410  "' for reach", irch, 'must be equal to zero.'
411  call store_error(errmsg)
412  end if
413  if (roughness_error) then
414  write (errmsg, '(3a,1x,i0,1x,a)') &
415  "All manfraction data in '", trim(adjustl(filename)), &
416  "' for reach", irch, 'must be greater than zero.'
417  call store_error(errmsg)
418  end if
419  end if
420  !
421  ! -- initialize and fill heights
422  allocate (heights(npts))
423  do n = 1, npts
424  heights(n) = this%cross_sections(irch)%height(n)
425  end do
426  !
427  ! -- get unique heights
428  call unique_values(heights, unique_heights)
429  !
430  ! -- calculate the product of the area and the hydraulic radius to
431  ! the 2/3 power
432  do n = 1, size(unique_heights)
433  if (unique_heights(n) <= dzero) cycle
434  ipos = 1
435  do i = -1, 1, 1
436  height = unique_heights(n) + real(i, dp) * dem6
437  aw = get_cross_section_area(npts, this%cross_sections(irch)%station, &
438  this%cross_sections(irch)%height, height)
439  rh = get_hydraulic_radius(npts, this%cross_sections(irch)%station, &
440  this%cross_sections(irch)%height, height)
441  factor(ipos) = aw * rh**dtwothirds
442  ipos = ipos + 1
443  end do
444  !
445  ! -- calculate the derivative
446  dc0 = (factor(2) - factor(1)) / dem6
447  dc1 = (factor(3) - factor(2)) / dem6
448  !
449  ! -- evaluate the difference
450  if (dc0 < dzero .or. dc1 < dzero) then
451  this%invalid = this%invalid + 1
452  height = unique_heights(n)
453  do i = 1, npts
454  if (this%cross_sections(irch)%height(i) == height) then
455  this%cross_sections(irch)%valid(i) = .false.
456  end if
457  end do
458  end if
459  end do
460  !
461  ! -- deallocate local storage
462  deallocate (heights)
463  deallocate (unique_heights)
464  !
465  ! -- return
466  return
real(dp), parameter dtwothirds
real constant 2/3
Definition: Constants.f90:69
real(dp), parameter dem6
real constant 1e-6
Definition: Constants.f90:108
This module contains stateless sfr subroutines and functions.
real(dp) function, public get_cross_section_area(npts, stations, heights, d)
Calculate the cross-sectional area for a reach.
real(dp) function, public get_hydraulic_radius(npts, stations, heights, d)
Calculate the hydraulic radius for a reach.
Here is the call graph for this function: