 MODULE toolmath
!!======================================================================
!!                   ***  MODULE toolmath  ***
!! Mathematical functions
!!======================================================================
!! * Modules used

   IMPLICIT NONE
   PRIVATE

!! * Accessibility
   PUBLIC sumindex                    ! function called by
   PUBLIC f_alpha

!! * Shared or public module variables

!! * Private variables

 CONTAINS
!!======================================================================
 
  FUNCTION sumindex(nstart,nend)
 
!&E---------------------------------------------------------------------
!&E                 ***  FUNCTION sumindex  ***
!&E
!&E ** Purpose : sum of indexes associated to biological variables
!&E              to make sure all indexes are initialized
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2009-10  (V. Garnier)  Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Arguments
   INTEGER, INTENT( in ) :: nstart,nend          ! index of the begining/end of the series
   INTEGER               :: sumindex             ! result of the function

!! * Local declarations
   INTEGER               :: i                    ! loop index

!!----------------------------------------------------------------------
!! * Executable part

   sumindex = nstart
   i = nstart + 1

   DO WHILE (i <= nend)
     sumindex = sumindex + i
     i = i + 1
   END DO

  END FUNCTION sumindex

!!======================================================================

  FUNCTION f_alpha(alpha)

!&E---------------------------------------------------------------------
!&E                 ***  FUNCTION f_alpha  ***
!&E
!&E ** Purpose : Limiting function for splitting vertical velocity into
!&E              explicit and implicit parts
!&E
!&E ** Description : Adaptive implicit vertical advection
!&E  c_min: Courant number threshold value for fully explicit advection
!&E  c_max: Courant number threshold value for fully implicit advection
!&E  between c_min and c_max: mixed advection
!&E
!&E ** Reference :
!&E       Shchepetkin A.F., 2015: An adaptive, Courant-number-dependent
!&E       implicit scheme for vertical advection in oceanic modeling,
!&E       Ocean Modelling, 91, 38-69.
!&E
!&E ** History :
!&E       !  2016-10  (S. Petton)  Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh

!! * Arguments
   REAL(kind=rsh),INTENT(in) :: alpha    ! local courant number
   REAL(kind=rsh)            :: f_alpha  ! result of the function

!! * Local declarations
   REAL(kind=rsh), PARAMETER :: alpha_min = 0.6_rsh,alpha_max=0.95_rsh
!REAL(kind=rsh), PARAMETER :: alpha_cut = 2.*alpha_max-alpha_min
   REAL(kind=rsh), PARAMETER :: alpha_cut = 1.3_rsh
!REAL(kind=rsh), PARAMETER :: alpha_tmp = 4.*alpha_max*(alpha_max-alpha_min)
   REAL(kind=rsh), PARAMETER :: alpha_tmp = 1.33_rsh

!!----------------------------------------------------------------------
!! * Executable part

! ==> f_alpha = 1 / f(alpha,alpha_max) dans le papier
   
   IF (alpha <= alpha_min) THEN
     f_alpha = 1._rsh
   ELSE IF (alpha < alpha_cut) THEN
     f_alpha = alpha_tmp / ( alpha_tmp + (alpha-alpha_min)**2 ) 
   ELSE
     f_alpha = alpha_max / alpha
   END IF

  END FUNCTION f_alpha


 END MODULE
