gwf-sfr-transient.f90

Go to the documentation of this file.

submodule(sfrmodule) sfrmoduletransient
contains
 
  !> @brief Method for solving for transient reaches
    !!
    !! Method to solve the continuity equation for a SFR package
    !! reach using the kinematic wave approximation.
    !!
  !<
  module procedure sfr_calc_transient
  use tdismodule, only: delt
 
  integer(I4B) :: igwfconn
  integer(I4B) :: number_picard
  integer(I4B) :: i
  integer(I4B) :: j
  real(DP) :: kinematic_residual
  real(DP) :: kinematic_storage
  real(DP) :: weight
  real(DP) :: celerity
  real(DP) :: courant
  real(DP) :: dq
  real(DP) :: qtol
  real(DP) :: qsrc
  real(DP) :: qlat
  real(DP) :: da
  real(DP) :: db
  real(DP) :: dc
  real(DP) :: dd
  real(DP) :: qa
  real(DP) :: qb
  real(DP) :: qc
  real(DP) :: xsa_a
  real(DP) :: xsa_b
  real(DP) :: xsa_c
  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) :: residual
  real(DP) :: residual2
  real(DP) :: residual_final
  real(DP) :: qderv
  real(DP) :: delq
  real(DP) :: delh
  real(DP) :: dd2
  real(DP) :: en1
  real(DP) :: en2
  real(DP) :: dmid
  logical :: lbisect
 
  weight = this%storage_weight
  dq = this%deps
  qtol = dq * dtwo
 
  celerity = dzero
  qgwf = dzero
  lbisect = .false.
 
  qlat = qr + qro - qe
 
  this%usinflow(n) = qu + qi + qfrommvr
 
  qa = this%usinflowold(n)
  qb = this%dsflowold(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 = 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
  d1 = (dc + dd) * dhalf
  d1old = d1
 
  ! estimate qgwf
  igwfconn = this%sfr_gwf_conn(n)
  if (igwfconn == 1) then
    ! -- a dry gaining reach with no inflow can get stuck at zero depth under
    !    Picard iteration: with a small time step zero depth is a stable fixed
    !    point because conductance and saturation vanish there. Only in that
    !    case (aquifer head above the streambed top and effectively no routed
    !    inflow) solve the depth by bisection instead, bracketed between the
    !    wet-streambed depth (DEM5) and the aquifer depth above the streambed
    !    top, so the reach rewets at any dt. A reach with meaningful inflow
    !    (dc >= DEM5) is left on the untouched Picard path.
    en2 = hgwf - this%strtop(n)
    if (en2 > dtwo * dem5 .and. dc < dem5) then
      lbisect = .true.
      en1 = dem5
      d1 = dhalf * (en1 + en2)
    end if
    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)
  else
    celerity = dzero
  end if
  courant = celerity * delt / this%length(n)
  if (courant > dzero) then
    if (courant < this%crmin(n)) this%crmin(n) = courant
    if (courant > this%crmax(n)) this%crmax(n) = courant
    this%crsum(n) = this%crsum(n) + courant
    this%crcnt(n) = this%crcnt(n) + 1
  end if
 
  qlat = qlat / this%length(n)
 
  number_picard = this%maxsfrpicard
  if (igwfconn == 1) then
    number_picard = this%maxsfrpicard
  else
    number_picard = 1
  end if
 
  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)
      qgwf = -qgwf
      if (qgwf > q) then
        qgwf = q
      end if
    end if
 
    qsrc = qlat - qgwf / this%length(n)
 
    newton: do j = 1, this%maxsfrit
      qd2 = qd + dq
      call this%sfr_calc_reach_depth(n, qd2, dd2)
      ad2 = this%calc_area_wet(n, dd2)
 
      residual = kinematic_residual(qa, qb, qc, qd, &
                                    xsa_a, xsa_b, xsa_c, ad, &
                                    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, &
                                     courant)
      qderv = (residual2 - residual) / dq
      if (qderv > dzero) then
        delq = -residual / qderv
      else
        delq = dzero
      end if
 
      if (qd + delq < dem30) then
        delq = -qd
      end if
 
      qd = qd + delq
 
      call this%sfr_calc_reach_depth(n, qd, dd)
      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, &
                                          courant)
 
      if (abs(delq) < qtol .and. abs(residual_final) < qtol) then
        exit newton
      end if
 
    end do newton
 
    qd = max(qd, dzero)
    if (lbisect) then
      ! -- bisection step. dmid is the depth implied by the routed flow at the
      !    current depth d1; the root is where the two are equal. Move the
      !    bracket end that keeps the root inside it.
      dmid = (dc + dd) * dhalf
      if (dmid > d1) then
        en1 = d1
      else
        en2 = d1
      end if
      d1 = dhalf * (en1 + en2)
      delh = en2 - en1
    else
      d1 = (dc + dd) * dhalf
      delh = (d1 - d1old)
    end if
 
    if (i > 1 .and. abs(delh) < this%dmaxchg) then
      exit kinematicpicard
    end if
 
  end do kinematicpicard
 
  this%storage(n) = kinematic_storage(xsa_a, xsa_b, xsa_c, ad, &
                                      this%length(n), delt, &
                                      courant)
 
  end procedure sfr_calc_transient
 
end submodule