Skip to content

Commit

Permalink
* update
Browse files Browse the repository at this point in the history
  • Loading branch information
@jdhughes-usgs
jdhughes-usgs committed Jul 9, 2024
1 parent 5c4c0fe commit 3404680
Show file tree
Hide file tree
Showing 3 changed files with 165 additions and 54 deletions.
64 changes: 64 additions & 0 deletions src/Model/GroundWaterFlow/gwf-sfr.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,7 @@ module SfrModule
procedure, private :: sfr_check_reaches
procedure, private :: sfr_check_connections
procedure, private :: sfr_check_diversions
procedure, private :: sfr_check_initialstages
procedure, private :: sfr_check_ustrf
! -- budget
procedure, private :: sfr_setup_budobj
Expand Down Expand Up @@ -927,6 +928,11 @@ subroutine sfr_ar(this)
if (this%idiversions /= 0) then
call this%sfr_check_diversions()
end if

! -- check the diversion data
if (this%istorage == 1) then
call this%sfr_check_initialstages()
end if
!
! -- terminate if errors were detected in any of the static sfr data
ierr = count_errors()
Expand Down Expand Up @@ -1878,6 +1884,7 @@ subroutine sfr_read_initial_stages(this)

rval = this%parser%GetDouble()
this%stage(n) = rval
this%depth(n) = rval - this%strtop(n)

if (rval < this%strtop(n)) then
write (errmsg, '(a,g0,a,1x,i0,1x,a,g0,a)') &
Expand Down Expand Up @@ -4924,6 +4931,63 @@ subroutine sfr_check_diversions(this)
return
end subroutine sfr_check_diversions

!> @brief Check initial stage data
!!
!! Method to check initial data for a SFR package and calculates
!! the initial upstream and downstream flows for the reach based
!! on the initial staalso creates the tables used to print input
!! data, if this option in enabled in the SFR package.
!!
!<
subroutine sfr_check_initialstages(this)
class(SfrType) :: this !< SfrType object

character(len=LINELENGTH) :: title
character(len=LINELENGTH) :: text
character(len=5) :: crch
integer(I4B) :: n
real(DP) :: qman

! skip check if storage is not activated
if (this%istorage == 0) return

! write header
if (this%iprpak /= 0) then
!
! -- reset the input table object
title = trim(adjustl(this%text))//' PACKAGE ('// &
trim(adjustl(this%packName))//') REACH INITIAL STAGE DATA'
call table_cr(this%inputtab, this%packName, title)
call this%inputtab%table_df(this%maxbound, 4, this%iout)
text = 'REACH'
call this%inputtab%initialize_column(text, 10, alignment=TABCENTER)
text = 'INITIAL STAGE'
call this%inputtab%initialize_column(text, 10, alignment=TABCENTER)
text = 'INITIAL DEPTH'
call this%inputtab%initialize_column(text, 10, alignment=TABCENTER)
text = 'INITIAL FLOW'
call this%inputtab%initialize_column(text, 10, alignment=TABCENTER)
end if
!
! -- check that data are correct
do n = 1, this%maxbound
write (crch, '(i5)') n

! calculate the initial flows
call this%sfr_calc_qman(n, this%depth(n), qman)
this%usinflow(n) = qman
this%dsflow(n) = qman

! add terms to the table
if (this%iprpak /= 0) then
call this%inputtab%add_term(n)
call this%inputtab%add_term(this%stage(n))
call this%inputtab%add_term(this%depth(n))
call this%inputtab%add_term(qman)
end if
end do
end subroutine sfr_check_initialstages

!> @brief Check upstream fraction data
!!
!! Method to check upstream fraction data for a SFR package.
Expand Down
83 changes: 51 additions & 32 deletions src/Model/GroundWaterFlow/submodules/gwf-sfr-transient.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,13 +12,13 @@

integer(I4B) :: igwfconn
integer(I4B) :: number_picard
integer(I4B) :: iconverged
integer(I4B) :: i
integer(I4B) :: j
real(DP) :: kinematic_residual
real(DP) :: kinematic_storage
real(DP) :: weight
! real(DP) :: weightinv
real(DP) :: celerity
real(DP) :: courant
real(DP) :: dq
real(DP) :: qtol
real(DP) :: qsrc
Expand All @@ -33,14 +33,16 @@
real(DP) :: xsa_a
real(DP) :: xsa_b
real(DP) :: xsa_c
! real(DP) :: xsa_d
real(DP) :: q
real(DP) :: q2
real(DP) :: d
real(DP) :: d2
real(DP) :: a
real(DP) :: a2
real(DP) :: d1old
real(DP) :: qd2
real(DP) :: ad
real(DP) :: ad2
real(DP) :: qgwf_pul
! real(DP) :: f11
! real(DP) :: f12
real(DP) :: residual
real(DP) :: residual2
real(DP) :: residual_final
Expand All @@ -53,38 +55,60 @@
dq = this%deps !DEM6 !DEM4 !this%deps
qtol = dq * DTWO !1e-9 !dq * DTWO

celerity = DZERO
qgwf = DZERO
qgwf_pul = DZERO

! calculate the flow at end of the reach
! excluding groundwater leakage
qsrc = qu + qi + qr - qe + qro + qfrommvr

qlat = (qr + qro - qe) / this%length(n)

this%usinflow(n) = qu + qi + qfrommvr

qa = this%usinflowold(n)
qb = this%dsflowold(n)
qc = this%usinflow(n)
call this%sfr_calc_reach_depth(n, qa, da)
call this%sfr_calc_reach_depth(n, qb, db)

qc = this%usinflow(n)
call this%sfr_calc_reach_depth(n, qc, dc)

xsa_a = this%calc_area_wet(n, da)
xsa_b = this%calc_area_wet(n, db)
xsa_c = this%calc_area_wet(n, dc)

! estimate qd
qd = (qc + qb) * DHALF
qd = this%dsflow(n)
if (qd == DZERO) then
qd = (qc + qb) * DHALF
end if
call this%sfr_calc_reach_depth(n, qd, dd)
ad = this%calc_area_wet(n, dd)

! -- estimate the depth at the midpoint
! estimate the depth at the midpoint
d1 = (dc + dd) * DHALF
d1old = d1

! estimate qgwf
igwfconn = this%sfr_gwf_conn(n)
if (igwfconn == 1) then
q = qu + qi + qr - qe + qro + qfrommvr
call this%sfr_calc_qgwf(n, d1, hgwf, qgwf)
qgwf = -qgwf
if (qgwf > q) then
qgwf = q
end if
end if

! calculate maximum wave speed and courant number
q = qc + qlat - qgwf
call this%sfr_calc_reach_depth(n, q, d)
a = this%calc_area_wet(n, d)
if (d > DZERO) then
q2 = q + dq
call this%sfr_calc_reach_depth(n, q2, d2)
a2 = this%calc_area_wet(n, d2)
celerity = (q2 - q) / (a2 - a)
courant = celerity * delt / this%length(n)
end if

number_picard = this%maxsfrpicard
if (igwfconn == 1) then
number_picard = this%maxsfrpicard
Expand All @@ -94,9 +118,11 @@

kinematicpicard: do i = 1, number_picard
if (igwfconn == 1) then
q = qu + qi + qr - qe + qro + qfrommvr
call this%sfr_calc_qgwf(n, d1, hgwf, qgwf)
if (qgwf > qsrc) then
qgwf = qsrc
qgwf = -qgwf
if (qgwf > q) then
qgwf = q
end if
end if

Expand All @@ -109,11 +135,13 @@

residual = kinematic_residual(qa, qb, qc, qd, &
xsa_a, xsa_b, xsa_c, ad, &
qsrc, this%length(n), weight, delt)
qsrc, this%length(n), weight, delt, &
courant)

residual2 = kinematic_residual(qa, qb, qc, qd2, &
xsa_a, xsa_b, xsa_c, ad2, &
qsrc, this%length(n), weight, delt)
qsrc, this%length(n), weight, delt, &
courant)
qderv = (residual2 - residual) / dq
if (qderv > DZERO) then
delq = -residual / qderv
Expand All @@ -131,7 +159,8 @@
ad = this%calc_area_wet(n, dd)
residual_final = kinematic_residual(qa, qb, qc, qd, &
xsa_a, xsa_b, xsa_c, ad, &
qsrc, this%length(n), weight, delt)
qsrc, this%length(n), weight, delt, &
courant)

if (abs(delq) < qtol) then ! .and. abs(residual_final) < qtol) then
exit newton
Expand All @@ -141,27 +170,17 @@

qd = max(qd, DZERO)
d1 = (dc + dd) * DHALF
! if (qd == DZERO) then
! d1 = DZERO
! else
! d1 = (dc + dd) * DHALF
! end if
delh = (d1 - d1old)

iconverged = 0
if (i == 1 .and. qd == DZERO) then
iconverged = 1
end if
if (i > 1 .and. abs(delh) < this%dmaxchg) then
iconverged = 1
end if
if (iconverged == 1) then
exit kinematicpicard
end if

end do kinematicpicard

this%storage(n) = kinematic_storage(qc, qd)
this%storage(n) = kinematic_storage(xsa_a, xsa_b, xsa_c, ad, &
this%length(n), delt, &
courant)

end procedure sfr_calc_transient

Expand Down
72 changes: 50 additions & 22 deletions src/Model/ModelUtilities/GwfSfrStatic.f90
View file
Original file line number Diff line number Diff line change
@@ -1,42 +1,70 @@
!> @brief Kinematic routing equation residual
!!
!! Method to calculate the kinematic-wave routing
!! residual.
!!
!<
function kinematic_residual(qa, qb, qc, qd, &
aa, ab, ac, ad, &
qsrc, length, weight, delt)
qsrc, length, weight, delt, &
courant)
use KindModule, only: DP
use ConstantsModule, only: DZERO, DHALF, DONE, DTWO
use ConstantsModule, only: DZERO, DHALF, DONE
! -- return variable
real(DP) :: kinematic_residual !< kinematic-wave residual
! -- dummy variables
real(DP), intent(in) :: qa
real(DP), intent(in) :: qb
real(DP), intent(in) :: qc
real(DP), intent(in) :: qd
real(DP), intent(in) :: aa
real(DP), intent(in) :: ab
real(DP), intent(in) :: ac
real(DP), intent(in) :: ad
real(DP), intent(in) :: qsrc
real(DP), intent(in) :: length
real(DP), intent(in) :: weight
real(DP), intent(in) :: delt
real(DP), intent(in) :: qa !< upstream flow at previous time
real(DP), intent(in) :: qb !< downstream flow at previous time
real(DP), intent(in) :: qc !< upstream flow
real(DP), intent(in) :: qd !< downstream flow
real(DP), intent(in) :: aa !< upstream area at previous time
real(DP), intent(in) :: ab !< downstream area at previous time
real(DP), intent(in) :: ac !< upstream area
real(DP), intent(in) :: ad !< downstream area
real(DP), intent(in) :: qsrc !< lateral flow term (L3T-1L-1)
real(DP), intent(in) :: length !< reach length
real(DP), intent(in) :: weight !< termporal weight
real(DP), intent(in) :: delt !< time step length
real(DP), intent(in) :: courant !< courant number
! --local variables
real(DP) :: f11
real(DP) :: f12

f11 = (weight * (qd - qc) + (DONE - weight) * (qb - qa)) / length
f12 = DHALF * ((ad - ab) + (ac - aa)) / delt
if (courant > DONE) then
f11 = (qd - qc) / length
f12 = (ac - aa) / delt
else
f11 = (weight * (qd - qc) + (DONE - weight) * (qb - qa)) / length
f12 = DHALF * ((ad - ab) + (ac - aa)) / delt
end if
kinematic_residual = f11 + f12 - qsrc

end function kinematic_residual

function kinematic_storage(qc, qd)
!> @brief Kinematic routing equation storage term
!!
!! Method to calculate the kinematic-wave routing
!! storage term.
!!
!<
function kinematic_storage(aa, ab, ac, ad, length, delt, courant)
use KindModule, only: DP
use ConstantsModule, only: DHALF
! -- return variable

real(DP) :: kinematic_storage !< kinematic-wave storage change
! -- dummy variables
real(DP), intent(in) :: qc
real(DP), intent(in) :: qd

kinematic_storage = (qd - qc)
real(DP), intent(in) :: aa !< upstream area at previous time
real(DP), intent(in) :: ab !< downstream area at previous time
real(DP), intent(in) :: ac !< upstream area
real(DP), intent(in) :: ad !< downstream area
real(DP), intent(in) :: length !< reach length
real(DP), intent(in) :: delt !< time step length
real(DP), intent(in) :: courant !< courant number

if (courant > DONE) then
kinematic_storage = (aa - ac) * length / delt
else
kinematic_storage = DHALF * ((ac + ad) - (aa + ab)) * length / delt
end if

end function kinematic_storage

0 comments on commit 3404680

Please sign in to comment.