
!===============================================================================
! AUTHOR : Sergey Budaev | http://bio.uib.no/te/sb/ | sergey.budaev@uib.no
! Repackaged Visual Range Fortran module with some accessory parts
! Copied from the AHA Model code SVN r6724.
!
! Contents:
! - REAL_PRECISION -- Portable precision constants;
! - REAL_FUNCTS -- Float point approximate equality;
! - VISRANGE -- Visual range calculation routines.
!
! Tested with GNU gfortran and Intel Fortran compilers.
!
! Version information:
! $Id: visrange.f90 7502 2018-11-09 18:53:09Z sbu062 $
!===============================================================================

!===============================================================================
!> @name Portable precision constants
module REAL_PRECISION

  implicit none

  !> Precision definition constants.
  !> SRP is the standard precision, it can normally be either 'single'
  !! or 'double'.
  !!  - 'single' is defined by `selected_real_kind( 6,   37)`
  !!  - 'double' is defined by `selected_real_kind(15,  307)`
  integer, parameter, public :: SRP = selected_real_kind( 15,   307)

  !> HRP defines a 'high' precision model. It is defined such that it has
  !! a minimum of 33 digits of precision and allows a range of values
  !! -10^4931 to -10^-4931 and 10^-4931 to 10^4931.
  !! HRP must be used in the VISRANGE module, anything smaller can produce
  !! FPU overflows.
  integer, parameter, public :: HRP = selected_real_kind(33, 4931)

  !> Some parameters should never be zero or below. In such cases they could
  !! be set to some smallest distinguishable non-zero value. Here set as
  !! the Fortran intrinsic `epsilon` function, a value that is almost
  !! negligible compared to one, i.e. the smallest `real` number *E* such
  !! that @f$ 1 + E > 1 @f$. In some cases it is also reasonable to set the
  !! tolerance limit to this parameter (see @ref intro_computation_real
  !! "Float point computations").
  !! @note The value of this parameter computed on a x86_64 Linux platform is
  !!       `1.19209290E-07`.
  real(SRP), parameter, public :: ZERO = epsilon(0.0_SRP)

  !> The smallest positive number in the commondata::srp standard real model.
  !! @note This parameter is used for definition of numerical tolerance.
  !! See @ref intro_computation_real "Float point computations".
  !! @note The value of this parameter computed on a x86_64 Linux platform is
  !!       `1.17549435E-38`.
  real(SRP), parameter, public :: TINY_SRP=tiny(1.0_SRP)

  !> The smallest positive number in the commondata::hrp high precision real
  !! model. See @ref intro_computation_real "Float point computations".
  !! @note This parameter is used for definition of numerical tolerance.
  real(HRP), parameter, public :: TINY_HRP=tiny(1.0_HRP)

  !> Default value of *low* tolerance (*high precision*). This is the
  !! *standard* commondata::srp precision. See @ref intro_computation_real
  !! "Float point computations".
  !! @note The value of this parameter computed on a x86_64 Linux platform is
  !!       `5.87747175E-38`.
  real(SRP), parameter, public :: TOLERANCE_LOW_DEF_SRP = TINY_SRP * 5.0_SRP

  !> Default value of *low* tolerance (*high precision*). This is the *high*
  !! commondata::hrp precision. See @ref intro_computation_real
  !! "Float point computations".
  real(HRP), parameter, public :: TOLERANCE_LOW_DEF_HRP = TINY_HRP * 5.0_HRP

  !> Default value of *high* tolerance (*low precision*). This is the
  !! *standard* commondata::srp precision real. See @ref intro_computation_real
  !! "Float point computations".
  !! @note The value of this parameter computed on a x86_64 Linux platform is
  !!       `1.19209290E-04`.
  real(SRP), parameter, public :: TOLERANCE_HIGH_DEF_SRP = ZERO * 1000.0_SRP

  !> Default value of *high* tolerance (*low precision*). This is the
  !! *high* commondata::hrp precision real. See @ref intro_computation_real
  !! "Float point computations".
  !! @note The value of this parameter computed on a x86_64 Linux platform is
  !!       `1.92592994438723585305597794258492732E-0031`.
  real(HRP), parameter, public :: TOLERANCE_HIGH_DEF_HRP =                    &
                                                  epsilon(0.0_HRP) * 1000.0_HRP

end module REAL_PRECISION

!===============================================================================
!> @name Float point approximate equality
!! This module defines an approximate float arithmetic type approximate
!! equality functions and operators.
module REAL_FUNCTS

  use REAL_PRECISION
  implicit none

  !> "Float equality" operator: Check if two real values are *nearly equal*
  !! using the commondata::is_near_zero() function. Thus function can be used
  !! for comparing two real values like the below.
  !! - The exact real comparison (incorrect due to possible rounding):
  !!   @code
  !!     if ( a == b ) ...
  !!   @endcode
  !! - should be substituted by such comparison:
  !!   @code
  !!     if ( a .feq. b) ...
  !!   @endcode
  !! .
  !! See the backend procedures commondata::float_equal_srp_operator() and
  !! commondata::float_equal_hrp_operator() for details and
  !! @ref intro_computation_real "Float point" for an introduction.
  interface operator (.feq.)
    procedure float_equal_srp_operator
    procedure float_equal_hrp_operator
  end interface operator (.feq.)
  ! @note Note that all component procedures are private.
  private :: float_equal_srp_operator, float_equal_hrp_operator

  !> Checks if a real number is near 0.0.
  !! Thus function can be used for comparing two real values like this:
  !! @code
  !!   if ( is_near_zero(a) ) then ...
  !! @endcode
  !! See the backend procedures commondata::is_near_zero_srp() and
  !! commondata::is_near_zero_hrp() for details.
  interface is_near_zero
    module procedure is_near_zero_srp
    module procedure is_near_zero_hrp
  end interface is_near_zero
  ! @note Note that all component procedures are private.
  private :: is_near_zero_srp, is_near_zero_hrp

  contains

  !-----------------------------------------------------------------------------
  !> This is a wrapper for the commondata::float_equal_srp() for building
  !! the user defined operator `.feq.` with default tolerance
  !! (`epsilon` parameter).
  !! The exact real comparison (incorrect due to possible rounding):
  !! @code
  !!   if ( a == b ) ...
  !! @endcode
  !! should be substituted by such comparison:
  !! @code
  !!   if ( a .feq. b) ...
  !! @endcode
  !! @warning This function is not intended for direct use.
  elemental function float_equal_srp_operator(value1, value2) result (equal)
    !> @param[in] value1 The first value for comparison.
    real(SRP), intent(in) :: value1
    !> @param[in] value2 The second value for comparison.
    real(SRP), intent(in) :: value2
    !> @return TRUE if the values `value1` and `value2` are nearly equal.
    logical :: equal
    ! EPSILON_LOC sets the tolerance limit for the default .feq. operator.
    real(SRP), parameter :: EPSILON_LOC = TOLERANCE_LOW_DEF_SRP
    equal = is_near_zero( abs(value1 - value2), EPSILON_LOC )
  end function float_equal_srp_operator

  !-----------------------------------------------------------------------------
  !> This is a wrapper for the commondata::float_equal_hrp() for building
  !! the user defined operator `.feq.` with default tolerance
  !! (`epsilon` parameter).
  !! The exact real comparison (incorrect due to possible rounding):
  !! @code
  !!   if ( a == b ) ...
  !! @endcode
  !! should be substituted by such comparison:
  !! @code
  !!   if ( a .feq. b) ...
  !! @endcode
  !! @warning This function is not intended for direct use.
  elemental function float_equal_hrp_operator(value1, value2) result (equal)
    !> @param[in] value1 The first value for comparison.
    real(HRP), intent(in) :: value1
    !> @param[in] value2 The second value for comparison.
    real(HRP), intent(in) :: value2
    !> @return TRUE if the values `value1` and `value2` are nearly equal.
    logical :: equal
    ! EPSILON_LOC sets the tolerance limit for the default .feq. operator.
    real(HRP), parameter :: EPSILON_LOC = TOLERANCE_LOW_DEF_HRP
    equal = is_near_zero( abs(value1 - value2), EPSILON_LOC )
  end function float_equal_hrp_operator

!-----------------------------------------------------------------------------
  !> Checks if a real number is near 0.0.
  !! Thus function can be used for comparing two real values like the below.
  !! @code
  !!   if ( is_near_zero(a) ) then ...
  !! @endcode
  !! @note Note that commondata::float_equal() can be used for approximate real
  !!       comparisons of two real type values.
  !! @note Modified from `Near0_dp()` function from Clerman & Spector 2012,
  !!       p. 250-251.
  !! @note This is the **standard** precision function (commondata::srp).
  elemental function is_near_zero_srp(test_number, epsilon) result (ret_val)
    !> @param[in] test_number the number to check for being near-zero.
    real(SRP), intent (in) :: test_number
    !> @param[in] epsilon optional (very small) tolerance value.
    real(SRP), intent (in), optional :: epsilon
    !> @return TRUE if the `test_number` is near-zero.
    logical :: ret_val

    real(kind(epsilon)) :: local_epsilon
    local_epsilon = TOLERANCE_LOW_DEF_SRP

    if ( present(epsilon) ) then
      if ( abs(epsilon) >= TINY_SRP ) local_epsilon = abs(epsilon)
    end if

    ret_val = abs(test_number) < local_epsilon

  end function is_near_zero_srp

  !-----------------------------------------------------------------------------
  !> Checks if a real number is near 0.0.
  !! Thus function can be used for comparing two real values like the below.
  !! @code
  !!   if ( is_near_zero(a) ) then ...
  !! @endcode
  !! @note Note that commondata::float_equal() can be used for approximate real
  !!       comparisons of two real type values.
  !! @note Modified from `Near0_dp()` function from Clerman & Spector 2012,
  !!       p. 250-251.
  !! @note This is the **high** precision function (commondata::hrp).
  elemental function is_near_zero_hrp(test_number, epsilon) result (ret_val)
    !> @param[in] test_number the number to check for being near-zero.
    real(HRP), intent (in) :: test_number
    !> @param[in] epsilon optional (very small) tolerance value.
    real(HRP), intent (in), optional :: epsilon
    !> @return TRUE if the `test_number` is near-zero.
    logical :: ret_val

    real(kind(epsilon)) :: local_epsilon
    local_epsilon = TOLERANCE_LOW_DEF_HRP

    if ( present(epsilon) ) then
      if ( abs(epsilon) >= TINY_HRP ) local_epsilon = abs(epsilon)
    end if

    ret_val = abs(test_number) < local_epsilon

  end function is_near_zero_hrp

end module REAL_FUNCTS

!===============================================================================
!> @name Visual range calculation backend.
module VISRANGE

  !> @note Note that the module REAL_PRECISION is called in REAL_FUNCTS and
  !!       therefore its "use" could be omitted.
  use REAL_PRECISION !, HRP_OTHER => HRP
  use REAL_FUNCTS    !, HRP_OTHER => HRP
  implicit none

  !-----------------------------------------------------------------------------
  !> HRP defines a 'high' precision model. It is defined such that it has
  !! a minimum of 33 digits of precision and allows a range of values
  !! -10^4931 to -10^-4931 and 10^-4931 to 10^4931.
  !! HRP must be used in the VISRANGE module, anything smaller can produce
  !! FPU overflows.
  !! @note Here HRP is defined as a private constant. Because HR is also
  !!       defined in the REAL_PRECISION module this upstream module HRP is
  !!       aliased as HRP_OTHER to avoid name conflicts. This is done to
  !!       make the module potentially usable without any upstream
  !!       dependencies. Therefore, a single HRP definition from
  !!       REAL_PRECISION can be used. In such a case delete this HRP and
  !!       the HRP_OTHER => HRP aliasing in the "use" construct above.
  !  integer, parameter, private :: HRP = selected_real_kind(33, 4931)
  !-----------------------------------------------------------------------------

  !> The subroutine ::srgetr() is defined through the generic interface
  !! including both real_precision::hrp and real_precision::srp versions.
  !! Then ::srgetr() can be called with parameters defined in both precision
  !! modes.
  !! Note that only ::srgetr() needs to be referred within the calling code
  !! and the appropriate version will be selected automatically by the
  !! compiler.
  interface srgetr
    module procedure srgetr_hrp
    module procedure srgetr_srp
  end interface srgetr

  contains  ! ..................................................................

  !-----------------------------------------------------------------------------
  !> Obtain visual range by solving the non-linear equation
  !! by means of Newton-Raphson iteration and derivation in
  !! subroutine the_environment::deriv(). Initial value is calculated in
  !! the_environment::easyr(). The calculation is based on the model described
  !! in Aksnes & Utne (1997) Sarsia 83:137-147.
  !! @note Programmed and tested 29 January 2001 Dag L Aksnes.
  !! @note This subroutine is left almost intact with only the most crucial
  !!       changes.
  !!       (a) added commondata::hrp precision specifier (128 bit precision
  !!           model) to real type specifiers and `_HRP` for literal constants;
  !!       (b) restored diagnostic `IER` output from archival Hed11.f90.
  !!       (c) added explicit `intent` and declared the procedures as
  !!           `pure` that is required for being parallel-friendly.
  ! **Input parameters:**
  !> @param[in]  RST  start value of r calculated by the_environment::easyr();
  !! @param[in]  c    beam attenuation coefficient `(m-1)`;
  !! @param[in]  C0   prey inherent contrast;
  !! @param[in]  Ap   prey area `(m^2)`;
  !! @param[in]  Vc   parameter characterising visual capacity (d.l.)
  !! @param[in]       this parameter is denoted E' in Aksnes & Utne;
  !! @param[in]  Ke   saturation parameter `(uE m-2 s-1)`;
  !! @param[in]  Eb   background irradiance at depth `DEPTH`;
  ! **Output parameters:**
  !> @param[out] r    the predator's visual range (when `F1=0`);
  !! @param[out] IER  **1** = No convergence after IEND steps, error return;
  !!                  **2** = Return in case of zero divisor;
  !!                  **3** = r out of allowed range (negative);
  !!                  **0** = valid r returned.
  !!
  !> **Notable variables:**
  !> - **F1**     function value of equation in `deriv`;
  !! - **FDER**   the derivative of the function.
  !! .
  elemental subroutine srgetr_hrp(r, c, C0, Ap, Vc, Ke, Eb, IER)
  ! Input parameters
  !     RST       : start value of r calculated by `easyr`
  !     c         : beam attenuation coefficient (m-1)
  !     C0        : prey inherent contrast
  !     Ap        : prey area (m^2)
  !     Vc        : parameter characterising visual capacity (d.l.)
  !                 this parameter is denoted E' in Aksnes & Utne
  !     Ke        : saturation parameter (uE m-2 s-1)
  !     Eb        : background irradiance at depth DEPTH
  ! Output parameters
  !     F1     : function value of equation in `deriv`
  !     FDER   : the derivative of the function
  !     r      : the predator's visual range (when F1=0)
  !     IER    : = 1, No convergence after IEND steps. Error return.
  !              = 2, Return in case of zero divisor.
  !              = 3, r out of allowed range (negative)
  !              = 0, valid r returned
  real(HRP), intent(in)   :: c, C0, Ap, Vc, Ke, Eb
  real(HRP), intent(out)  :: r
  integer, optional, intent(out) :: IER

  real(HRP) :: AS, EPS, RST, TOL, TOLF, F1, FDER, DX
  integer   :: IEND, I

  !.............................................................................

    ! Initial guess of visual range (RST)
    call easyr(RST,C0,Ap,Vc,Ke,Eb)

    ! Upper boundary of allowed error of visual range.
    EPS = TOLERANCE_HIGH_DEF_HRP
    ! Maximum number of iteration steps
    IEND = 100

    ! Prepare iteration
    r = RST
    TOL = r

    call deriv(r,F1,FDER,c,C0,Ap,Vc,Ke,Eb)
    TOLF = 100.0_HRP * EPS

    ! Start iteration expmt
    ! @warning Cannot probably be converted to `do concurrent` due to `exit`s
    !!         from the loop.
    do 6 I = 1, IEND
      if ( is_near_zero(F1) ) goto 7

      ! Equation is not satisfied by r
      if ( is_near_zero(FDER) ) goto 8

      ! Iteration is possible
      DX = F1/FDER
      r = r-DX

      ! Test on allowed range
      if (r .LT. 0.0_HRP) goto 9

      TOL = r

      call deriv(r,F1,FDER,c,C0,Ap,Vc,Ke,Eb)

      ! Test on satisfactory accuracy
      TOL = EPS
      AS = abs(r)
      if ((AS-1.0_HRP) > 0.0_HRP) TOL = TOL*AS
      TOL = TOL*AS
      if ((abs(DX)-TOL) > 0.0_HRP) goto 6
      if ((abs(F1)-TOLF) .LE. 0.0_HRP) goto 7
  6   continue

      ! No convergence after IEND steps. Error return.
      if (present(IER)) IER = 1
  7   return
      ! Return in case of zero divisor
  8   if (present(IER)) IER = 2
      return
      ! r out of allowed range (negative)
  9   if (present(IER)) IER = 3
      return

  end subroutine srgetr_hrp

  !-----------------------------------------------------------------------------
  !> Obtain visual range by solving the non-linear equation
  !! by means of Newton-Raphson iteration and derivation in
  !! subroutine the_environment::deriv(). Initial value is calculated in
  !! the_environment::easyr(). The calculation is based on the model described
  !! in Aksnes & Utne (1997) Sarsia 83:137-147.
  !!
  !! @note This is a wrapper to ::srgetr_hrp() function converting all
  !!       real parameters to standard precision real_precision::srp from
  !!       real_precision::hrp
  elemental subroutine srgetr_srp(r, c, C0, Ap, Vc, Ke, Eb, IER)

  ! Input parameters are the same as in srgetr_hrp()
  real(SRP), intent(in)   :: c, C0, Ap, Vc, Ke, Eb
  real(SRP), intent(out)  :: r
  integer, optional, intent(out) :: IER

  ! Local copies of output parameters.
  real(HRP) :: r_loc
  integer :: IER_loc

  r_loc = -9999.0_HRP
  IER_loc = 0

  ! In this wrapper, ::srgetr_hrp() function is called with inline
  ! conversion of all input parameters from real_precision::hrp to
  ! real_precision::srp.
  call srgetr_hrp( r_loc,              &
                   real(c,  HRP),      &
                   real(C0, HRP),      &
                   real(Ap, HRP),      &
                   real(Vc, HRP),      &
                   real(Ke, HRP),      &
                   real(Eb, HRP),      &
                   IER_loc             )

  r = r_loc
  if (present(IER)) IER = IER_loc

  end subroutine srgetr_srp

  !-----------------------------------------------------------------------------
  !> Obtain a first estimate of visual range by using a simplified
  !! expression of visual range. See `srgetr()` for more details.
  !! @note This subroutine is left almost intact, only (a) added
  !!       commondata::hrp for real type (`HRP` is *high real precision*
  !!       (128 bit) and is defined in @ref commondata).
  elemental subroutine easyr(r, C0, Ap, Vc, Ke, Eb)
    real(HRP), intent(out) :: r
    real(HRP), intent(in)  :: C0, Ap, Vc, Ke, Eb
    real(HRP) :: R2
    ! See the calling routine the_environment::srgetr() for explanation of
    ! parameters
    R2 = abs(C0)*Ap*Vc*Eb/(Ke+Eb)
    r  = sqrt(R2)
    return
  end subroutine easyr

  !-----------------------------------------------------------------------------
  !> Derivation of equation for visual range of a predator.
  !! See the_environment::srgetr() for more details.
  !! @note This is a high precision version. But higher precision alone is
  !!       not sufficient to prevent numerical exponentiation overflow.
  !! @note This subroutine is left almost intact, only (a) added
  !!       commondata::hrp for literal constants and real type (`HRP` is the
  !!       *high real precision* (128 bit) and is defined in @ref commondata),
  !!       (b) added numerical overflow safeguard code based on `MAX_LOG`
  !!       exponentiability limit; added logging of overflow using `LOG_MSG`.
  ! **Output parameters:**
  !> @param[out]   F1     function value of equation in
  !!                      the_environment::deriv();
  !! @param[out]   FDER   the derivative of the function;
  !! @param[inout] r      the predator's visual range (when `F1=0`).
  !!
  !> **Input parameters:**
  !!
  !> See explanation in calling routine the_environment::srgetr().
  !!
  !> The function and the derivative is calculated on the basis of the
  !! log-transformed expression.
  elemental subroutine deriv(r, F1, FDER, c, C0, Ap, Vc, Ke, Eb)
    ! Input parameters
    !     See explanation in calling routine
    ! Output parameters
    !     F1     : function value of equation in `deriv`
    !     FDER   : the derivative of the function
    !     r      : the predator's visual range (when F1=0)
    !
    !    The function and the derivative is calculated on the basis of the
    !    log-transformed expression
    real(HRP), intent(inout) :: r
    real(HRP), intent(out)   :: F1, FDER
    real(HRP), intent(in)    :: c, C0, Ap, Vc, Ke, Eb

    real(HRP) :: FR1, FR2

    !> ### Implementation notes ###
    !! `MAX_LOG` is a parameter determining the safe limit of `exp` function
    !! overflow in the current float point precision model, well below this.
    !! We cannot calculate precise exponent of a value exceeding this
    !! parameter. The maximum possible exponentiation-able value is set
    !! to the maximum **128-bit** real value kind `Q_PREC_128`, this bottom
    !! line value would set a safe limit for `HRP` calculations.
    !! A benefit of this approach is that it doesn't require IEEE exception
    !! handling that depends on not fully portable IEEE modules.
    !! @note  The real 128 bit limit (Q_PREC_128) is sufficient to calculate
    !!        visual range up to the fish length of approximately 700 cm
    !!        (area 153.9 cm2). At this level, the maximum visual range is
    !!        1343.34 cm. This would be sufficient for non-whales.
    real(HRP), parameter :: HUGE_REAL = huge(0.0_HRP)
    real(HRP), parameter :: MAX_LOG = log(HUGE_REAL)

    ! PROCNAME is the procedure name for logging and debugging (with MODNAME).
    character(len=*), parameter :: PROCNAME = "(deriv)"

    FR2=log(abs(C0)*Ap*Vc)

    ! @note    Here `exp(c*r)` can result in huge FPU overflow (Infinity),
    !          so we now check if the exponent is likely to be so high and
    !          **if yes**, just set the huge ceiling. In such a case one needs
    !          to *rescale* variables (e.g. from meters to tenth or hundreds
    !          of meters).
    ! @warning The resulting calculations are then most probably grossly
    !          wrong but we nonetheless avoid FPU runtime error: incorrect
    !          arithmetic calculations.
    if (c*r < MAX_LOG) then
      FR1=log(((Ke+Eb)/Eb)*r*r*exp(c*r))
    else
      FR1=HUGE_REAL
      ! Subroutine `LOG_MSG` is defined in `LOGGER` module. It reports
      ! numerical overflow errors.
      ! @warning `LOG_MSG` is **disabled** here as it cannot be **pure** and
      !          hampers purifying `deriv` for parallel processing.
      !          A negative side is no error report logging. But the `HRP`
      !          64 bit precision model seems to work okay. There is also
      !          a guarantee against overflow and gross errors with the
      !          `HUGE_REAL` parameter.
      !call LOG_MSG (  "ERROR in " // PROCNAME // ":  FR1 overflow "  //      &
      !                "exceeding " // TOSTR(HUGE_REAL) // " limit :" //      &
      !                " c*r =" // TOSTR(c*r)  )
    end if

    F1 = FR1-FR2
    FDER = c + 2.0_HRP/r
    return

  end subroutine deriv

end module VISRANGE
