FreundlichIsotherm.f90

Go to the documentation of this file.

Module freundlichisothermmodule
 
  use kindmodule, only: dp, i4b
  use isotherminterfacemodule, only: isothermtype
  use constantsmodule, only: dzero
 
  Implicit None
  Private
  Public :: freundlichisothermtype
 
  !> @brief Freundlich isotherm implementation of `IsothermType`.
  !>
  !> Sorbed concentration is cs = Kf*c^a.
  !>
  !> However, this expression has a singularity at c = 0 when a < 1,
  !> leading to infinite derivative. To avoid this, the Freundlich isotherm
  !> is modified as follows.
  !>
  !> Modified Sorbed concentration is cs = Kf*(c + eps)^a - Kf*eps^a
  !> where eps = (K/(a*Kf))^(1/(a-1)) and K is a large constant (default 10).
  !> This ensures that the derivative at c = 0 is below K.
  !<
  type, extends(isothermtype) :: freundlichisothermtype
    real(dp), pointer, dimension(:) :: kf => null() !< Freundlich constant
    real(dp), pointer, dimension(:) :: a => null() !< Freundlich exponent
  contains
    procedure :: value
    procedure :: derivative
  end type freundlichisothermtype
 
  interface freundlichisothermtype
    module procedure constructor
  end interface freundlichisothermtype
 
contains
 
  !> @brief Constructor for Freundlich isotherm
  !<
  function constructor(Kf, a) Result(isotherm)
    type(freundlichisothermtype) :: isotherm
    ! -- dummy
    real(dp), pointer, dimension(:), intent(in) :: kf
    real(dp), pointer, dimension(:), intent(in) :: a
    ! -- local
    isotherm%Kf => kf
    isotherm%a => a
 
  end function constructor
 
  !> @brief Evaluate the isotherm at a given node
  !<
  function value(this, c, n) result(val)
    ! -- return
    real(dp) :: val !< isotherm value
    ! -- dummy
    class(freundlichisothermtype), intent(in) :: this
    real(dp), dimension(:), intent(in) :: c !< concentration array
    integer(I4B), intent(in) :: n !< node index
    real(dp), parameter :: k = 10.0_dp !< constant to limit derivative at c=0
    real(dp) :: eps !< small concentration offset
 
    eps = (k / (this%a(n) * this%Kf(n)))**(1.0_dp / (this%a(n) - 1.0_dp))
 
    if (c(n) > dzero) then
      val = this%Kf(n) * (c(n) + eps)**this%a(n) - this%Kf(n) * eps**this%a(n)
    else
      val = 0.0_dp
    end if
  end function value
 
  !> @brief Evaluate derivative of the isotherm at a given node
  !<
  function derivative(this, c, n) result(derv)
    ! -- return
    real(dp) :: derv !< derivative d(value)/dc evaluated at c
    ! -- dummy
    class(freundlichisothermtype), intent(in) :: this
    real(dp), dimension(:), intent(in) :: c !< concentration array
    integer(I4B), intent(in) :: n !< node index
    real(dp), parameter :: k = 10.0_dp !< constant to limit derivative at c=0
    real(dp) :: eps !< small concentration offset
 
    eps = (k / (this%a(n) * this%Kf(n)))**(1.0_dp / (this%a(n) - 1.0_dp))
 
    if (c(n) > dzero) then
      derv = this%a(n) * this%Kf(n) * ((c(n) + eps)**(this%a(n) - 1.0_dp))
    else
      derv = 0.0_dp
    end if
  end function derivative
 
end module freundlichisothermmodule