GwfConductanceUtils.f90

Go to the documentation of this file.

!> @brief This module contains stateless conductance functions
!!
!! This module contains the functions to calculate the horizontal
!! and vertical conductance between two cells that are used in the
!! the node property flow (NPF) package. It also contains functions
!! to calculate the wetted cell fraction. This module does not
!! depend on the NPF package.
!!
!<
module gwfconductanceutilsmodule
  use kindmodule, only: dp, i4b
  use constantsmodule, only: dzero, dhalf, done, &
                             dlnlow, dlnhigh, &
                             c3d_staggered
 
  implicit none
  private
  public :: hcond
  public :: vcond
  public :: condmean
  public :: thksatnm
  public :: staggered_thkfrac
  public :: ccond_hmean
 
  !> @brief enumerator that defines the conductance options
  !<
  ENUM, BIND(C)
    ENUMERATOR :: ccond_hmean = 0 !< Harmonic mean
    ENUMERATOR :: ccond_lmean = 1 !< Logarithmic mean
    ENUMERATOR :: ccond_amtlmk = 2 !< Arithmetic-mean thickness and logarithmic-mean hydraulic conductivity
    ENUMERATOR :: ccond_amthmk = 3 !< Arithmetic-mean thickness and harmonic-mean hydraulic conductivity
  END ENUM
 
contains
 
  !> @brief Horizontal conductance between two cells
  !!
  !! This function uses a weighted transmissivity in the harmonic mean
  !! conductance calculations. This differs from the MODFLOW-NWT and
  !! MODFLOW-USG conductance calculations for the Newton-Raphson formulation
  !! which use a weighted hydraulic conductivity.
  !<
  function hcond(ibdn, ibdm, ictn, ictm, iupstream, ihc, icellavg, &
                 condsat, hn, hm, satn, satm, hkn, hkm, topn, topm, &
                 botn, botm, cln, clm, fawidth) &
    result(condnm)
    ! return variable
    real(dp) :: condnm !< horizontal conductance between two cells
    ! dummy
    integer(I4B), intent(in) :: ibdn !< cell n active flag
    integer(I4B), intent(in) :: ibdm !< cell m active flag
    integer(I4B), intent(in) :: ictn !< cell n convertible cell flag
    integer(I4B), intent(in) :: ictm !< cell m convertible cell flag
    integer(I4B), intent(in) :: iupstream !< flag for upstream weighting
    integer(I4B), intent(in) :: ihc !< connection type
    integer(I4B), intent(in) :: icellavg !< cell averaging option
    real(dp), intent(in) :: condsat !< saturated conductance
    real(dp), intent(in) :: hn !< cell n head
    real(dp), intent(in) :: hm !< cell m head
    real(dp), intent(in) :: satn !< cell n wetted cell fraction
    real(dp), intent(in) :: satm !< cell m wetted cell fraction
    real(dp), intent(in) :: hkn !< horizontal hydraulic conductivity for cell n (in the direction of cell m)
    real(dp), intent(in) :: hkm !< horizontal hydraulic conductivity for cell m (in the direction of cell n)
    real(dp), intent(in) :: topn !< top of cell n
    real(dp), intent(in) :: topm !< top of cell m
    real(dp), intent(in) :: botn !< bottom of cell n
    real(dp), intent(in) :: botm !< bottom of cell m
    real(dp), intent(in) :: cln !< distance from the center of cell n to the shared face with cell m
    real(dp), intent(in) :: clm !< distance from the center of cell m to the shared face with cell n
    real(dp), intent(in) :: fawidth !< width of cell perpendicular to flow
 
    ! n or m is inactive
    if (ibdn == 0 .or. ibdm == 0) then
      condnm = dzero
      ! both cells are non-convertible
    elseif (ictn == 0 .and. ictm == 0) then
      condnm = condsat
    else if (iupstream == 1) then
      condnm = convertible_upstream(hn, hm, satn, satm, condsat)
    else
      condnm = convertible_standard(ihc, icellavg, &
                                    satn, satm, hkn, hkm, &
                                    topn, topm, botn, botm, &
                                    cln, clm, fawidth)
    end if
  end function hcond
 
  !> @brief Convertible cell(s) with upstream weighted horizontal conductance
  !<
  function convertible_upstream(hn, hm, satn, satm, condsat) &
    result(condnm)
    ! return variable
    real(dp) :: condnm !< horizontal conductance between two cells
    ! dummy
    real(dp), intent(in) :: condsat !< saturated conductance
    real(dp), intent(in) :: hn !< cell n head
    real(dp), intent(in) :: hm !< cell m head
    real(dp), intent(in) :: satn !< cell n wetted cell fraction
    real(dp), intent(in) :: satm !< cell m wetted cell fraction
    !  local
    real(dp) :: sat_up
 
    if (hn > hm) then
      sat_up = satn
    else
      sat_up = satm
    end if
    condnm = sat_up * condsat
  end function convertible_upstream
 
  !> @brief Convertible cell(s) with standard weighted horizontal conductance
  !<
  function convertible_standard(ihc, icellavg, satn, satm, hkn, hkm, &
                                topn, topm, botn, botm, cln, clm, fawidth) &
    result(condnm)
    ! return variable
    real(dp) :: condnm !< horizontal conductance between two cells
    ! dummy
    integer(I4B), intent(in) :: ihc !< connection type
    integer(I4B), intent(in) :: icellavg !< cell averaging option
    real(dp), intent(in) :: satn !< cell n wetted cell fraction
    real(dp), intent(in) :: satm !< cell m wetted cell fraction
    real(dp), intent(in) :: hkn !< horizontal hydraulic conductivity for cell n (in the direction of cell m)
    real(dp), intent(in) :: hkm !< horizontal hydraulic conductivity for cell m (in the direction of cell n)
    real(dp), intent(in) :: topn !< top of cell n
    real(dp), intent(in) :: topm !< top of cell m
    real(dp), intent(in) :: botn !< bottom of cell n
    real(dp), intent(in) :: botm !< bottom of cell m
    real(dp), intent(in) :: cln !< distance from the center of cell n to the shared face with cell m
    real(dp), intent(in) :: clm !< distance from the center of cell m to the shared face with cell n
    real(dp), intent(in) :: fawidth !< width of cell perpendicular to flow
    !  local
    real(dp) :: thksatn
    real(dp) :: thksatm
 
    if (ihc == c3d_staggered) then
      thksatn = staggered_thkfrac(topn, botn, satn, topm, botm)
      thksatm = staggered_thkfrac(topm, botm, satm, topn, botn)
    else
      thksatn = satn * (topn - botn)
      thksatm = satm * (topm - botm)
    end if
    condnm = condmean(hkn, hkm, thksatn, thksatm, cln, clm, &
                      fawidth, icellavg)
  end function convertible_standard
 
  !> @brief Vertical conductance between two cells
  !<
  function vcond(ibdn, ibdm, ictn, ictm, inewton, ivarcv, idewatcv, &
                 condsat, hn, hm, vkn, vkm, satn, satm, topn, topm, botn, &
                 botm, flowarea) result(condnm)
    ! return variable
    real(dp) :: condnm
    ! dummy
    integer(I4B), intent(in) :: ibdn !< cell n active flag
    integer(I4B), intent(in) :: ibdm !< cell m active flag
    integer(I4B), intent(in) :: ictn !< cell n convertible cell flag
    integer(I4B), intent(in) :: ictm !< cell m convertible cell flag
    integer(I4B), intent(in) :: inewton !< flag for Newton-Raphson formulation
    integer(I4B), intent(in) :: ivarcv !< variable vertical conductance flag
    integer(I4B), intent(in) :: idewatcv !< dewatered vertical conductance flag
    real(dp), intent(in) :: condsat !< saturated conductance
    real(dp), intent(in) :: hn !< cell n head
    real(dp), intent(in) :: hm !< cell m head
    real(dp), intent(in) :: vkn !< vertical hydraulic conductivity for cell n (in the direction of cell m)
    real(dp), intent(in) :: vkm !< vertical hydraulic conductivity for cell m (in the direction of cell n)
    real(dp), intent(in) :: satn !< cell n wetted cell fraction
    real(dp), intent(in) :: satm !< cell m wetted cell fraction
    real(dp), intent(in) :: topn !< top of cell n
    real(dp), intent(in) :: topm !< top of cell m
    real(dp), intent(in) :: botn !< bottom of cell n
    real(dp), intent(in) :: botm !< bottom of cell m
    real(dp), intent(in) :: flowarea !< flow area between cell n and m
    ! local
    real(dp) :: satntmp
    real(dp) :: satmtmp
    real(dp) :: bovk1
    real(dp) :: bovk2
    real(dp) :: denom
    !
    ! Either n or m is inactive
    if (ibdn == 0 .or. ibdm == 0) then
      condnm = dzero
      !
      ! constantcv
    elseif (ivarcv == 0) then
      condnm = condsat
      !
      ! both cells are non-convertible
    elseif (ictn == 0 .and. ictm == 0) then
      condnm = condsat
      !
      ! both cells are fully saturated
    elseif (hn >= topn .and. hm >= topm) then
      condnm = condsat
      !
      ! At least one cell is partially saturated, so recalculate vertical
      ! conductance for this connection
      ! todo: upstream weighting?
    else
      !
      ! Default is for CV correction (dewatered option); use underlying
      !    saturation of 1.
      satntmp = satn
      satmtmp = satm
      if (idewatcv == 0) then
        if (botn > botm) then
          satmtmp = done
        else
          satntmp = done
        end if
      else
        if (botn > botm) then
          if (satmtmp < done) then
            satmtmp = dzero
          end if
        else
          if (satntmp < done) then
            satntmp = dzero
          end if
        end if
      end if
      bovk1 = satntmp * (topn - botn) * dhalf / vkn
      bovk2 = satmtmp * (topm - botm) * dhalf / vkm
      denom = (bovk1 + bovk2)
      if (denom /= dzero) then
        condnm = flowarea / denom
      else
        condnm = dzero
      end if
    end if
  end function vcond
 
  !> @brief Calculate the conductance between two cells
  !<
  function condmean(k1, k2, thick1, thick2, cl1, cl2, width, iavgmeth)
    ! return variable
    real(dp) :: condmean !< mean conductance between two cells
    ! dummy
    real(dp), intent(in) :: k1 !< hydraulic conductivity for cell n (in the direction of cell m)
    real(dp), intent(in) :: k2 !< hydraulic conductivity for cell m (in the direction of celln)
    real(dp), intent(in) :: thick1 !< saturated thickness for cell 1
    real(dp), intent(in) :: thick2 !< saturated thickness for cell 2
    real(dp), intent(in) :: cl1 !< distance from the center of cell n to the shared face with cell m
    real(dp), intent(in) :: cl2 !< distance from the center of cell m to the shared face with cell n
    real(dp), intent(in) :: width !< width of cell perpendicular to flow
    integer(I4B), intent(in) :: iavgmeth !< averaging method
    ! local
    real(dp) :: t1
    real(dp) :: t2
    real(dp) :: tmean
    real(dp) :: kmean
    real(dp) :: denom
    !
    ! Initialize
    t1 = k1 * thick1
    t2 = k2 * thick2
 
    ! Averaging method
    select case (iavgmeth)
 
    case (ccond_hmean)
      if (t1 * t2 > dzero) then
        condmean = width * t1 * t2 / (t1 * cl2 + t2 * cl1)
      else
        condmean = dzero
      end if
 
    case (ccond_lmean)
      if (t1 * t2 > dzero) then
        tmean = logmean(t1, t2)
      else
        tmean = dzero
      end if
      condmean = tmean * width / (cl1 + cl2)
 
    case (ccond_amtlmk)
      if (k1 * k2 > dzero) then
        kmean = logmean(k1, k2)
      else
        kmean = dzero
      end if
      condmean = kmean * dhalf * (thick1 + thick2) * width / (cl1 + cl2)
 
    case (ccond_amthmk)
      denom = (k1 * cl2 + k2 * cl1)
      if (denom > dzero) then
        kmean = k1 * k2 / denom
      else
        kmean = dzero
      end if
      condmean = kmean * dhalf * (thick1 + thick2) * width
    end select
  end function condmean
 
  !> @brief Calculate the the logarithmic mean of two double precision numbers
  !!
  !! Use an approximation if the ratio is near 1
  !<
  function logmean(d1, d2)
    ! return variable
    real(dp) :: logmean !< logarithmic mean for two number
    ! dummy
    real(dp), intent(in) :: d1 !< first number
    real(dp), intent(in) :: d2 !< second number
    ! local
    real(dp) :: drat
 
    drat = d2 / d1
    if (drat <= dlnlow .or. drat >= dlnhigh) then
      logmean = (d2 - d1) / log(drat)
    else
      logmean = dhalf * (d1 + d2)
    end if
  end function logmean
 
  !> @brief Calculate wetted cell thickness at interface between two cells
  !<
  function thksatnm(ibdn, ibdm, ictn, ictm, iupstream, ihc, &
                    hn, hm, satn, satm, topn, topm, botn, botm) result(res)
    ! return variable
    real(dp) :: res !< wetted cell thickness for connection nm
    ! dummy
    integer(I4B), intent(in) :: ibdn !< cell n active flag
    integer(I4B), intent(in) :: ibdm !< cell m active flag
    integer(I4B), intent(in) :: ictn !< cell n convertible cell flag
    integer(I4B), intent(in) :: ictm !< cell m convertible cell flag
    integer(I4B), intent(in) :: iupstream !< flag for upstream weighting
    integer(I4B), intent(in) :: ihc !< connection type
    real(dp), intent(in) :: hn !< cell n head
    real(dp), intent(in) :: hm !< cell m head
    real(dp), intent(in) :: satn !< cell n wetted cell fraction
    real(dp), intent(in) :: satm !< cell m wetted cell fraction
    real(dp), intent(in) :: topn !< top of cell n
    real(dp), intent(in) :: topm !< top of cell m
    real(dp), intent(in) :: botn !< bottom of cell n
    real(dp), intent(in) :: botm !< bottom of cell m
    ! local
    real(dp) :: thksatn
    real(dp) :: thksatm
    real(dp) :: sill_top
    real(dp) :: sill_bot
    !
    ! n or m is inactive
    if (ibdn == 0 .or. ibdm == 0) then
      res = dzero
      !
      ! both cells are non-convertible
    elseif (ictn == 0 .and. ictm == 0) then
      if (ihc == c3d_staggered) then
        sill_top = min(topn, topm)
        sill_bot = max(botn, botm)
 
        thksatn = max(sill_top - sill_bot, dzero)
        thksatm = thksatn
      else
        thksatn = topn - botn
        thksatm = topm - botm
      end if
      res = dhalf * (thksatn + thksatm)
      !
      ! At least one of the cells is convertible
    elseif (iupstream == 1) then
      if (hn > hm) then
        res = satn * (topn - botn)
      else
        res = satm * (topm - botm)
      end if
      !
      ! At least one of the cells is convertible and not upstream weighted
    else
      if (ihc == c3d_staggered) then
        thksatn = staggered_thkfrac(topn, botn, satn, topm, botm)
        thksatm = staggered_thkfrac(topm, botm, satm, topn, botn)
      else
        thksatn = satn * (topn - botn)
        thksatm = satm * (topm - botm)
      end if
      res = dhalf * (thksatn + thksatm)
    end if
  end function thksatnm
 
  !> @brief Calculate the thickness fraction for staggered grids
  !<
  function staggered_thkfrac(top, bot, sat, topc, botc) result(res)
    ! return variable
    real(dp) :: res !< staggered thickness fraction for cell
    ! dummy
    real(dp) :: top !< top of cell
    real(dp) :: bot !< bottom of cell
    real(dp) :: sat !< cell saturation
    real(dp) :: topc !< top of connected cell
    real(dp) :: botc !< bottom of connected cells
    ! local
    real(dp) :: sill_top
    real(dp) :: sill_bot
    real(dp) :: tp
 
    sill_top = min(top, topc)
    sill_bot = max(bot, botc)
    tp = bot + sat * (top - bot)
    res = max(min(tp, sill_top) - sill_bot, dzero)
  end function staggered_thkfrac
 
end module gwfconductanceutilsmodule