log_mean Function

public pure function log_mean(a, b) result(am)

Numerically stable logarithmic mean (Ismail-Roe 2009 / Ranocha 2018).

_log = (a - b) / (ln(a) - ln(b)) for a /= b _log = (a + b) / 2 for a == b

Falls back to a Taylor series in u = f^2 with f = (a-b)/(a+b) when a and b are close, to avoid 0/0.

Arguments

TypeIntentOptionalAttributesName
real(kind=prec), intent(in) :: a
real(kind=prec), intent(in) :: b

Return Value real(kind=prec)


Called by

proc~~log_mean~2~~CalledByGraph proc~log_mean~2 log_mean proc~twopointflux3d_esatmo3d_t twopointflux3d_ESAtmo3D_t proc~twopointflux3d_esatmo3d_t->proc~log_mean~2 proc~sourcemethod_esatmo3d_t SourceMethod_ESAtmo3D_t proc~sourcemethod_esatmo3d_t->proc~log_mean~2

Contents

Source Code


Source Code

  pure function log_mean(a,b) result(am)
    !! Numerically stable logarithmic mean (Ismail-Roe 2009 / Ranocha 2018).
    !!
    !!   <a>_log = (a - b) / (ln(a) - ln(b))   for a /= b
    !!   <a>_log = (a + b) / 2                 for a == b
    !!
    !! Falls back to a Taylor series in u = f^2 with f = (a-b)/(a+b)
    !! when a and b are close, to avoid 0/0.
    real(prec),intent(in) :: a,b
    real(prec) :: am
    ! Local
    real(prec) :: zeta,f,u,F_F
    real(prec),parameter :: eps_lm = 1.0e-2_prec

    zeta = a/b
    f = (zeta-1.0_prec)/(zeta+1.0_prec)
    u = f*f
    if(u < eps_lm) then
      F_F = 1.0_prec+u/3.0_prec+u*u/5.0_prec+u*u*u/7.0_prec
    else
      F_F = log(zeta)/(2.0_prec*f)
    endif
    am = (a+b)/(2.0_prec*F_F)
  endfunction log_mean