SmoothingFunctions.f90

Go to the documentation of this file.

module smoothingmodule
  use kindmodule, only: dp, i4b
  use constantsmodule, only: dzero, dhalf, done, dtwo, dthree, dfour, &
                             dsix, dprec, dem2, dem4, dem5, dem6, dem8, dem14
  implicit none
 
contains
 
  !> @ brief SCurve
  !!
  !! Computes the S curve for smooth derivatives between x=0 and x=1
  !! from mfusg smooth subroutine in gwf2wel7u1.f
  !<
  subroutine sscurve(x, range, dydx, y)
    real(DP), intent(in) :: x
    real(DP), intent(in) :: range
    real(DP), intent(inout) :: dydx
    real(DP), intent(inout) :: y
    !--local variables
    real(DP) :: s
    real(DP) :: xs
    ! -- code
    !
    s = range
    if (s < dprec) s = dprec
    xs = x / s
    if (xs < dzero) xs = dzero
    if (xs <= dzero) then
      y = dzero
      dydx = dzero
    elseif (xs < done) then
      y = -dtwo * xs**dthree + dthree * xs**dtwo
      dydx = -dsix * xs**dtwo + dsix * xs
    else
      y = done
      dydx = dzero
    end if
  end subroutine sscurve
 
  !> @ brief sCubicLinear
  !!
  !! Computes the s curve where dy/dx = 0 at x=0; and dy/dx = 1 at x=1.
  !! Smooths from zero to a slope of 1.
  !<
  subroutine scubiclinear(x, range, dydx, y)
    real(DP), intent(in) :: x
    real(DP), intent(in) :: range
    real(DP), intent(inout) :: dydx
    real(DP), intent(inout) :: y
    !--local variables
    real(DP) :: s
    real(DP) :: xs
    ! -- code
    !
    s = range
    if (s < dprec) s = dprec
    xs = x / s
    if (xs < dzero) xs = dzero
    if (xs <= dzero) then
      y = dzero
      dydx = dzero
    elseif (xs < done) then
      y = -done * xs**dthree + dtwo * xs**dtwo
      dydx = -dthree * xs**dtwo + dfour * xs
    else
      y = done
      dydx = dzero
    end if
  end subroutine scubiclinear
 
  !> @ brief sRamp
  !!
  !! Smoothed ramp function; a smooth approximation of max(x, 0) with a
  !! continuous value and slope. The value is zero with zero slope for x <= 0 and
  !! is the exact line x with unit slope for x >= range, joined by a cubic over
  !! (0, range). Unlike sCubicLinear, which levels off at 1, the value continues
  !! with unit slope, so it is exact outside the transition (max(x, 0) for
  !! |x| >= range).
  !<
  subroutine sramp(x, range, dydx, y)
    real(DP), intent(in) :: x
    real(DP), intent(in) :: range
    real(DP), intent(inout) :: dydx
    real(DP), intent(inout) :: y
    !--local variables
    real(DP) :: s
    real(DP) :: xs
    ! -- code
    !
    s = range
    if (s < dprec) s = dprec
    xs = x / s
    if (xs <= dzero) then
      y = dzero
      dydx = dzero
    elseif (xs < done) then
      y = x * xs * (dtwo - xs)
      dydx = xs * (dfour - dthree * xs)
    else
      y = x
      dydx = done
    end if
  end subroutine sramp
 
  !> @ brief sCubic
  !!
  !! Nonlinear smoothing function returns value between 0-1; cubic function
  !<
  subroutine scubic(x, range, dydx, y)
    real(DP), intent(inout) :: x
    real(DP), intent(inout) :: range
    real(DP), intent(inout) :: dydx
    real(DP), intent(inout) :: y
    !--local variables
    real(DP) :: s, aa, bb
    real(DP) :: cof1, cof2, cof3
    ! -- code
    !
    dydx = dzero
    y = dzero
    if (range < dprec) range = dprec
    if (x < dprec) x = dprec
    s = range
    aa = -dsix / (s**dthree)
    bb = -dsix / (s**dtwo)
    cof1 = x**dtwo
    cof2 = -(dtwo * x) / (s**dthree)
    cof3 = dthree / (s**dtwo)
    y = cof1 * (cof2 + cof3)
    dydx = (aa * x**dtwo - bb * x)
    if (x <= dzero) then
      y = dzero
      dydx = dzero
    else if ((x - s) > -dprec) then
      y = done
      dydx = dzero
    end if
  end subroutine scubic
 
  !> @ brief sLinear
  !!
  !! Linear smoothing function returns value between 0-1
  !<
  subroutine slinear(x, range, dydx, y)
    real(DP), intent(inout) :: x
    real(DP), intent(inout) :: range
    real(DP), intent(inout) :: dydx
    real(DP), intent(inout) :: y
    !--local variables
    real(DP) :: s
    ! -- code
    !
    dydx = dzero
    y = dzero
    if (range < dprec) range = dprec
    if (x < dprec) x = dprec
    s = range
    y = done - (s - x) / s
    dydx = done / s
    if (y > done) then
      y = done
      dydx = dzero
    end if
  end subroutine slinear
 
  !> @ brief sQuadratic
  !!
  !! Nonlinear quadratic smoothing function returns value between 0-1
  !<
  subroutine squadratic(x, range, dydx, y)
    real(DP), intent(inout) :: x
    real(DP), intent(inout) :: range
    real(DP), intent(inout) :: dydx
    real(DP), intent(inout) :: y
    !--local variables
    real(DP) :: s
    ! -- code
    !
    dydx = dzero
    y = dzero
    if (range < dprec) range = dprec
    if (x < dprec) x = dprec
    s = range
    y = (x**dtwo) / (s**dtwo)
    dydx = dtwo * x / (s**dtwo)
    if (y > done) then
      y = done
      dydx = dzero
    end if
  end subroutine squadratic
 
  !> @ brief sChSmooth
  !!
  !! Function to smooth channel variables during channel drying
  !<
  subroutine schsmooth(d, smooth, dwdh)
    real(DP), intent(in) :: d
    real(DP), intent(inout) :: smooth
    real(DP), intent(inout) :: dwdh
    !
    ! -- local variables
    real(DP) :: s
    real(DP) :: diff
    real(DP) :: aa
    real(DP) :: ad
    real(DP) :: b
    real(DP) :: x
    real(DP) :: y
    ! -- code
    !
    smooth = dzero
    s = dem5
    x = d
    diff = x - s
    if (diff > dzero) then
      smooth = done
      dwdh = dzero
    else
      aa = -done / (s**dtwo)
      ad = -dtwo / (s**dtwo)
      b = dtwo / s
      y = aa * x**dtwo + b * x
      dwdh = (ad * x + b)
      if (x <= dzero) then
        y = dzero
        dwdh = dzero
      else if (diff > -dem14) then
        y = done
        dwdh = dzero
      end if
      smooth = y
    end if
  end subroutine schsmooth
 
  !> @ brief sLinearSaturation
  !!
  !! Linear saturation function returns value between 0-1
  !<
  function slinearsaturation(top, bot, x) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: top
    real(dp), intent(in) :: bot
    real(dp), intent(in) :: x
    ! -- local
    real(dp) :: b
    ! -- code
    !
    b = top - bot
    if (x < bot) then
      y = dzero
    else if (x > top) then
      y = done
    else
      y = (x - bot) / b
    end if
  end function slinearsaturation
 
  !> @ brief sCubicSaturation
  !!
  !! Nonlinear cubic saturation function returns value between 0-1
  !<
  function scubicsaturation(top, bot, x, eps) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: top
    real(dp), intent(in) :: bot
    real(dp), intent(in) :: x
    real(dp), intent(in), optional :: eps
    ! -- local
    real(dp) :: teps
    real(dp) :: w
    real(dp) :: b
    real(dp) :: s
    real(dp) :: cof1
    real(dp) :: cof2
    ! -- code
    !
    if (present(eps)) then
      teps = eps
    else
      teps = dem2
    end if
    w = x - bot
    b = top - bot
    s = teps * b
    cof1 = done / (s**dtwo)
    cof2 = dtwo / s
    if (w < dzero) then
      y = dzero
    else if (w < s) then
      y = -cof1 * (w**dthree) + cof2 * (w**dtwo)
    else if (w < (b - s)) then
      y = w / b
    else if (w < b) then
      y = done + cof1 * ((b - w)**dthree) - cof2 * ((b - w)**dtwo)
    else
      y = done
    end if
 
  end function scubicsaturation
 
  !> @ brief sQuadraticSaturation
  !!
  !! Nonlinear quadratic saturation function returns value between 0-1
  !<
  function squadraticsaturation(top, bot, x, eps) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: top
    real(dp), intent(in) :: bot
    real(dp), intent(in) :: x
    real(dp), optional, intent(in) :: eps
    ! -- local
    real(dp) :: teps
    real(dp) :: b
    real(dp) :: br
    real(dp) :: bri
    real(dp) :: av
    ! -- code
    !
    if (present(eps)) then
      teps = eps
    else
      teps = dem6
    end if
    b = top - bot
    if (b > dzero) then
      if (x < bot) then
        br = dzero
      else if (x > top) then
        br = done
      else
        br = (x - bot) / b
      end if
      av = done / (done - teps)
      bri = done - br
      if (br < teps) then
        y = av * dhalf * (br * br) / teps
      elseif (br < (done - teps)) then
        y = av * br + dhalf * (done - av)
      elseif (br < done) then
        y = done - ((av * dhalf * (bri * bri)) / teps)
      else
        y = done
      end if
    else
      if (x < bot) then
        y = dzero
      else
        y = done
      end if
    end if
 
  end function squadraticsaturation
 
  !> @ brief sQuadraticSaturation
  !!
  !! van Genuchten saturation function returns value between 0-1
  !<
  function svangenuchtensaturation(top, bot, x, alpha, beta, sr) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: top
    real(dp), intent(in) :: bot
    real(dp), intent(in) :: x
    real(dp), intent(in) :: alpha
    real(dp), intent(in) :: beta
    real(dp), intent(in) :: sr
    ! -- local
    real(dp) :: b
    real(dp) :: pc
    real(dp) :: gamma
    real(dp) :: seff
    ! -- code
    !
    b = top - bot
    pc = (dhalf * b) - x
    if (pc <= dzero) then
      y = dzero
    else
      gamma = done - (done / beta)
      seff = (done + (alpha * pc)**beta)**gamma
      seff = done / seff
      y = seff * (done - sr) + sr
    end if
 
  end function svangenuchtensaturation
 
  !> @ brief Derivative of the quadratic saturation function
  !!
  !! Derivative of nonlinear smoothing function returns value between 0-1;
  !<
  function squadraticsaturationderivative(top, bot, x, eps) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: top
    real(dp), intent(in) :: bot
    real(dp), intent(in) :: x
    real(dp), optional, intent(in) :: eps
    ! -- local
    real(dp) :: teps
    real(dp) :: b
    real(dp) :: br
    real(dp) :: bri
    real(dp) :: av
    ! -- code
    !
    if (present(eps)) then
      teps = eps
    else
      teps = dem6
    end if
    b = top - bot
    if (x < bot) then
      br = dzero
    else if (x > top) then
      br = done
    else
      br = (x - bot) / b
    end if
    av = done / (done - teps)
    bri = done - br
    if (br < teps) then
      y = av * br / teps
    elseif (br < (done - teps)) then
      y = av
    elseif (br < done) then
      y = av * bri / teps
    else
      y = dzero
    end if
    y = y / b
 
  end function squadraticsaturationderivative
 
  !> @ brief sQSaturation
  !!
  !! Nonlinear smoothing function returns value between 0-1
  !<
  function sqsaturation(top, bot, x, c1, c2) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: top
    real(dp), intent(in) :: bot
    real(dp), intent(in) :: x
    real(dp), intent(in), optional :: c1
    real(dp), intent(in), optional :: c2
    ! -- local
    real(dp) :: w
    real(dp) :: b
    real(dp) :: s
    real(dp) :: cof1
    real(dp) :: cof2
    ! -- code
    !
    ! -- process optional variables
    if (present(c1)) then
      cof1 = c1
    else
      cof1 = -dtwo
    end if
    if (present(c2)) then
      cof2 = c2
    else
      cof2 = dthree
    end if
    !
    ! -- calculate head difference from bottom (w),
    !    calculate range (b), and
    !    calculate normalized head difference from bottom (s)
    w = x - bot
    b = top - bot
    s = w / b
    !
    ! -- divide cof1 and cof2 by range to the power 3 and 2, respectively
    cof1 = cof1 / b**dthree
    cof2 = cof2 / b**dtwo
    !
    ! -- calculate fraction
    if (s < dzero) then
      y = dzero
    else if (s < done) then
      y = cof1 * w**dthree + cof2 * w**dtwo
    else
      y = done
    end if
  end function sqsaturation
 
  !> @ brief sQSaturationDerivative
  !!
  !! Nonlinear smoothing function returns value between 0-1
  !<
  function sqsaturationderivative(top, bot, x, c1, c2) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: top
    real(dp), intent(in) :: bot
    real(dp), intent(in) :: x
    real(dp), intent(in), optional :: c1
    real(dp), intent(in), optional :: c2
    ! -- local
    real(dp) :: w
    real(dp) :: b
    real(dp) :: s
    real(dp) :: cof1
    real(dp) :: cof2
    ! -- code
    !
    !
    ! -- process optional variables
    if (present(c1)) then
      cof1 = c1
    else
      cof1 = -dtwo
    end if
    if (present(c2)) then
      cof2 = c2
    else
      cof2 = dthree
    end if
    !
    ! -- calculate head difference from bottom (w),
    !    calculate range (b), and
    !    calculate normalized head difference from bottom (s)
    w = x - bot
    b = top - bot
    s = w / b
    !
    ! -- multiply cof1 and cof2 by 3 and 2, respectively, and then
    !    divide by range to the power 3 and 2, respectively
    cof1 = cof1 * dthree / b**dthree
    cof2 = cof2 * dtwo / b**dtwo
    !
    ! -- calculate derivative of fraction with respect to x
    if (s < dzero) then
      y = dzero
    else if (s < done) then
      y = cof1 * w**dtwo + cof2 * w
    else
      y = dzero
    end if
  end function sqsaturationderivative
 
  !> @ brief sSlope
  !!
  !! Nonlinear smoothing function returns a smoothed value of y whose limbs are
  !! the lines yi + (sm * dx) for x-values less than xi and yi + (sp * dx) for
  !! x-values greater than xi, where dx = x - xi. The slope change at xi is
  !! smoothed with the hyperbolic smoothed min and max of Chen and Mangasarian
  !! (1996), as applied to hydrologic-model thresholds by Kavetski and Kuczera
  !! (2007, eq. 18; https://dx.doi.org/10.1029/2006WR005195); here the smoothing
  !! constant is b**2 - a**2 with b = a / (sqrt(2) - 1). Currently unused.
  !<
  function sslope(x, xi, yi, sm, sp, ta) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: x
    real(dp), intent(in) :: xi
    real(dp), intent(in) :: yi
    real(dp), intent(in) :: sm
    real(dp), intent(in) :: sp
    real(dp), optional, intent(in) :: ta
    ! -- local
    real(dp) :: a
    real(dp) :: b
    real(dp) :: dx
    real(dp) :: xm
    real(dp) :: xp
    real(dp) :: ym
    real(dp) :: yp
    !
    ! -- set smoothing variable a
    if (present(ta)) then
      a = ta
    else
      a = dem8
    end if
    !
    ! -- calculate b from smoothing variable a
    b = a / (sqrt(dtwo) - done)
    !
    ! -- calculate contributions to y. xm and xp are smoothed values of
    !    min(x, xi) and max(x, xi), so the limbs are yi + sm*dx (x < xi) and
    !    yi + sp*dx (x > xi).
    dx = x - xi
    xm = dhalf * (x + xi - sqrt(dx**dtwo + b**dtwo - a**dtwo))
    xp = dhalf * (x + xi + sqrt(dx**dtwo + b**dtwo - a**dtwo))
    ym = sm * (xm - xi)
    yp = sp * (xp - xi)
    !
    ! -- calculate y from ym and yp contributions
    y = yi + ym + yp
  end function sslope
 
  !> @ brief sSlopeDerivative
  !!
  !! Derivative of sSlope (see that routine for the provenance): smoothly
  !! transitions from slope sm for x-values less than xi to slope sp for x-values
  !! greater than xi, where dx = x - xi.
  !<
  function sslopederivative(x, xi, sm, sp, ta) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: x
    real(dp), intent(in) :: xi
    real(dp), intent(in) :: sm
    real(dp), intent(in) :: sp
    real(dp), optional, intent(in) :: ta
    ! -- local
    real(dp) :: a
    real(dp) :: b
    real(dp) :: dx
    real(dp) :: mu
    real(dp) :: rho
    !
    ! -- set smoothing variable a
    if (present(ta)) then
      a = ta
    else
      a = dem8
    end if
    !
    ! -- calculate b from smoothing variable a
    b = a / (sqrt(dtwo) - done)
    !
    ! -- calculate contributions to derivative
    dx = x - xi
    mu = sqrt(dx**dtwo + b**dtwo - a**dtwo)
    rho = dx / mu
    !
    ! -- calculate derivative from individual contributions
    y = dhalf * (sm + sp) - dhalf * rho * (sm - sp)
  end function sslopederivative
 
  !> @ brief sQuadratic0sp
  !!
  !! Nonlinear smoothing function returns a smoothed value of y that uses a
  !! quadratic to smooth x over range of xi - epsilon to xi + epsilon.
  !! Simplification of sQuadraticSlope with sm = 0, sp = 1, and yi = 0.
  !! From Panday et al. (2013) - eq. 35 - https://dx.doi.org/10.5066/F7R20ZFJ
  !<
  function squadratic0sp(x, xi, tomega) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: x
    real(dp), intent(in) :: xi
    real(dp), optional, intent(in) :: tomega
    ! -- local
    real(dp) :: omega
    real(dp) :: epsilon
    real(dp) :: dx
    !
    ! -- set smoothing interval
    if (present(tomega)) then
      omega = tomega
    else
      omega = dem6
    end if
    !
    ! -- set smoothing interval
    epsilon = dhalf * omega
    !
    ! -- calculate distance from xi
    dx = x - xi
    !
    ! -- evaluate smoothing function
    if (dx < -epsilon) then
      y = xi
    else if (dx < epsilon) then
      y = (dx**dtwo / (dfour * epsilon)) + dhalf * dx + (epsilon / dfour) + xi
    else
      y = x
    end if
  end function squadratic0sp
 
  !> @ brief sQuadratic0spDerivative
  !!
  !! Derivative of nonlinear smoothing function returns a smoothed value of y
  !! that uses a quadratic to smooth x over range of xi - epsilon to xi + epsilon.
  !! Simplification of sQuadraticSlope with sm = 0, sp = 1, and yi = 0.
  !! From Panday et al. (2013) - eq. 35 - https://dx.doi.org/10.5066/F7R20ZFJ
  !<
  function squadratic0spderivative(x, xi, tomega) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: x
    real(dp), intent(in) :: xi
    real(dp), optional, intent(in) :: tomega
    ! -- local
    real(dp) :: omega
    real(dp) :: epsilon
    real(dp) :: dx
    !
    ! -- set smoothing interval
    if (present(tomega)) then
      omega = tomega
    else
      omega = dem6
    end if
    !
    ! -- set smoothing interval
    epsilon = dhalf * omega
    !
    ! -- calculate distance from xi
    dx = x - xi
    !
    ! -- evaluate smoothing function
    if (dx < -epsilon) then
      y = 0
    else if (dx < epsilon) then
      y = (dx / omega) + dhalf
    else
      y = 1
    end if
  end function squadratic0spderivative
 
  !> @ brief sQuadraticSlope
  !!
  !! Quadratic smoothing function returns a smoothed value of y whose limbs are
  !! the lines yi + (sm * dx) for x-values less than xi and yi + (sp * dx) for
  !! x-values greater than xi, where dx = x - xi. The limbs are blended by a
  !! quadratic over xi +/- epsilon, so the value equals yi at xi only in the
  !! limit of zero smoothing; otherwise it differs from yi by a term
  !! proportional to the smoothing interval.
  !<
  function squadraticslope(x, xi, yi, sm, sp, tomega) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: x
    real(dp), intent(in) :: xi
    real(dp), intent(in) :: yi
    real(dp), intent(in) :: sm
    real(dp), intent(in) :: sp
    real(dp), optional, intent(in) :: tomega
    ! -- local
    real(dp) :: omega
    real(dp) :: epsilon
    real(dp) :: dx
    real(dp) :: c
    !
    ! -- set smoothing interval
    if (present(tomega)) then
      omega = tomega
    else
      omega = dem6
    end if
    !
    ! -- set smoothing interval
    epsilon = dhalf * omega
    !
    ! -- calculate distance from xi
    dx = x - xi
    !
    ! -- evaluate smoothing function
    if (dx < -epsilon) then
      y = sm * dx
    else if (dx < epsilon) then
      c = dx / epsilon
      y = dhalf * epsilon * (dhalf * (sp - sm) * (done + c**dtwo) + (sm + sp) * c)
    else
      y = sp * dx
    end if
    !
    ! -- add value at xi
    y = y + yi
  end function squadraticslope
 
  !> @ brief sQuadraticSlopeDerivative
  !!
  !! Derivative of sQuadraticSlope: smoothly transitions from the slope sm for
  !! x-values less than xi to the slope sp for x-values greater than xi, blended
  !! by a line over xi +/- epsilon, where dx = x - xi.
  !<
  function squadraticslopederivative(x, xi, sm, sp, tomega) result(y)
    ! -- return
    real(dp) :: y
    ! -- dummy variables
    real(dp), intent(in) :: x
    real(dp), intent(in) :: xi
    real(dp), intent(in) :: sm
    real(dp), intent(in) :: sp
    real(dp), optional, intent(in) :: tomega
    ! -- local
    real(dp) :: omega
    real(dp) :: epsilon
    real(dp) :: dx
    real(dp) :: c
    !
    ! -- set smoothing interval
    if (present(tomega)) then
      omega = tomega
    else
      omega = dem6
    end if
    !
    ! -- set smoothing interval
    epsilon = dhalf * omega
    !
    ! -- calculate distance from xi
    dx = x - xi
    !
    ! -- evaluate smoothing function
    if (dx < -epsilon) then
      y = sm
    else if (dx < epsilon) then
      c = dx / epsilon
      y = dhalf * ((sp - sm) * c + (sm + sp))
    else
      y = sp
    end if
  end function squadraticslopederivative
 
end module smoothingmodule