!MNH_LIC Copyright 2016-2019 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! Modifications:
!  P.Wautelet 22/01/2019: replace double precision declarations by real(kind(0.0d0)) (to allow compilation by NAG compiler)
!
!-----------------------------------------------------------------
!#################################
        MODULE MODI_LIMA_FUNCTIONS
!#################################
!
INTERFACE
!
FUNCTION COUNTJV(LTAB,I1,I2,I3) RESULT(IC)
  LOGICAL, DIMENSION(:,:,:),   INTENT(IN)           :: LTAB 
  INTEGER, DIMENSION(:),       INTENT(INOUT)        :: I1,I2,I3 
  INTEGER                                           :: IC
END FUNCTION COUNTJV
!
FUNCTION MOMG (PALPHA,PNU,PP) RESULT (PMOMG)
  REAL,                        INTENT(IN)           :: PALPHA 
  REAL,                        INTENT(IN)           :: PNU    
  REAL,                        INTENT(IN)           :: PP     
  REAL                                              :: PMOMG  
END FUNCTION MOMG
!
FUNCTION RECT(PA,PB,PX,PX1,PX2)  RESULT(PRECT)
  REAL,                        INTENT(IN)           :: PA
  REAL,                        INTENT(IN)           :: PB
  REAL, DIMENSION(:),          INTENT(IN)           :: PX
  REAL,                        INTENT(IN)           :: PX1
  REAL,                        INTENT(IN)           :: PX2
  REAL, DIMENSION(SIZE(PX,1))                       :: PRECT
END FUNCTION RECT
!
FUNCTION DELTA(PA,PB,PX,PX1,PX2)  RESULT(PDELTA)
  REAL,                        INTENT(IN)           :: PA
  REAL,                        INTENT(IN)           :: PB
  REAL, DIMENSION(:),          INTENT(IN)           :: PX
  REAL,                        INTENT(IN)           :: PX1
  REAL,                        INTENT(IN)           :: PX2
  REAL, DIMENSION(SIZE(PX,1))                       :: PDELTA
END FUNCTION DELTA
!
FUNCTION DELTA_VEC(PA,PB,PX,PX1,PX2)  RESULT(PDELTA_VEC)
  REAL,                        INTENT(IN)           :: PA
  REAL,                        INTENT(IN)           :: PB
  REAL, DIMENSION(:),          INTENT(IN)           :: PX
  REAL, DIMENSION(:),          INTENT(IN)           :: PX1
  REAL, DIMENSION(:),          INTENT(IN)           :: PX2
  REAL, DIMENSION(SIZE(PX,1))                       :: PDELTA_VEC
END FUNCTION DELTA_VEC
!
SUBROUTINE GAULAG(x,w,n,alf)
  INTEGER,                     INTENT(IN)           :: n
  REAL,                        INTENT(IN)           :: alf
  REAL, DIMENSION(n),          INTENT(INOUT)        :: w, x
END SUBROUTINE GAULAG
!
SUBROUTINE GAUHER(x,w,n)
  INTEGER,                     INTENT(IN)           :: n
  REAL, DIMENSION(n),          INTENT(INOUT)        :: w, x
END SUBROUTINE GAUHER
!
END INTERFACE
!
END MODULE MODI_LIMA_FUNCTIONS
!
!------------------------------------------------------------------------------
!
!#########################################
FUNCTION COUNTJV(LTAB,I1,I2,I3) RESULT(IC)
!#########################################
!
  IMPLICIT NONE
!
  LOGICAL, DIMENSION(:,:,:) :: LTAB ! Mask
  INTEGER, DIMENSION(:) :: I1,I2,I3 ! Used to replace the COUNT and PACK
  INTEGER :: JI,JJ,JK,IC
!
  IC = 0
  DO JK = 1,SIZE(LTAB,3)
     DO JJ = 1,SIZE(LTAB,2)
        DO JI = 1,SIZE(LTAB,1)
           IF( LTAB(JI,JJ,JK) ) THEN
              IC = IC +1
              I1(IC) = JI
              I2(IC) = JJ
              I3(IC) = JK
           END IF
        END DO
     END DO
  END DO
!
END FUNCTION COUNTJV
!
!------------------------------------------------------------------------------
!
!###########################################
FUNCTION MOMG (PALPHA,PNU,PP) RESULT (PMOMG)
!###########################################
!
! auxiliary routine used to compute the Pth moment order of the generalized
! gamma law
!
  USE MODI_GAMMA
!
  IMPLICIT NONE
!
  REAL     :: PALPHA ! first shape parameter of the dimensionnal distribution
  REAL     :: PNU    ! second shape parameter of the dimensionnal distribution
  REAL     :: PP     ! order of the moment
  REAL     :: PMOMG  ! result: moment of order ZP
!
  PMOMG = GAMMA_X0D(PNU+PP/PALPHA)/GAMMA_X0D(PNU)
!
END FUNCTION MOMG
!
!------------------------------------------------------------------------------
!
!#############################################
FUNCTION RECT(PA,PB,PX,PX1,PX2)  RESULT(PRECT)
!#############################################
!
!     PRECT takes the value PA if PX1<=PX<PX2, and PB outside the [PX1;PX2[ interval
!
  IMPLICIT NONE
!
  REAL,                        INTENT(IN)           :: PA
  REAL,                        INTENT(IN)           :: PB
  REAL, DIMENSION(:),          INTENT(IN)           :: PX
  REAL,                        INTENT(IN)           :: PX1
  REAL,                        INTENT(IN)           :: PX2
  REAL, DIMENSION(SIZE(PX,1))                       :: PRECT
!
  PRECT(:) = PB
  WHERE (PX(:).GE.PX1 .AND. PX(:).LT.PX2)
     PRECT(:) = PA
  END WHERE
  RETURN
!
END FUNCTION RECT
!
!-------------------------------------------------------------------------------
!
!###############################################
FUNCTION DELTA(PA,PB,PX,PX1,PX2)  RESULT(PDELTA)
!###############################################
!
!     PDELTA takes the value PA if PX<PX1, and PB if PX>=PX2
!     PDELTA is a cubic interpolation between PA and PB for PX between PX1 and PX2 
!
  IMPLICIT NONE
!
  REAL,                        INTENT(IN)           :: PA
  REAL,                        INTENT(IN)           :: PB
  REAL, DIMENSION(:),          INTENT(IN)           :: PX
  REAL,                        INTENT(IN)           :: PX1
  REAL,                        INTENT(IN)           :: PX2
  REAL, DIMENSION(SIZE(PX,1))                       :: PDELTA
!
!*       local variable
!
  REAL                                              :: ZA
!
  ZA = 6.0*(PA-PB)/(PX2-PX1)**3
  WHERE     (PX(:).LT.PX1)
     PDELTA(:) = PA
  ELSEWHERE (PX(:).GE.PX2)
     PDELTA(:) = PB
  ELSEWHERE
     PDELTA(:) =   PA + ZA*PX1**2*(PX1/6.0 - 0.5*PX2)              &
          + ZA*PX1*PX2*                          (PX(:))    &
          - (0.5*ZA*(PX1+PX2))*                  (PX(:)**2) & 
          + (ZA/3.0)*                            (PX(:)**3)
  END WHERE
  RETURN
!
END FUNCTION DELTA
!
!-------------------------------------------------------------------------------
!
!#######################################################
FUNCTION DELTA_VEC(PA,PB,PX,PX1,PX2)  RESULT(PDELTA_VEC)
!#######################################################
!
!    Same as DELTA for vectorized PX1 and PX2 arguments
!
  IMPLICIT NONE
!
  REAL,                        INTENT(IN)           :: PA
  REAL,                        INTENT(IN)           :: PB
  REAL, DIMENSION(:),          INTENT(IN)           :: PX
  REAL, DIMENSION(:),          INTENT(IN)           :: PX1
  REAL, DIMENSION(:),          INTENT(IN)           :: PX2
  REAL, DIMENSION(SIZE(PX,1))                       :: PDELTA_VEC
!
!*       local variable
!
  REAL, DIMENSION(SIZE(PX,1))                       :: ZA
!
  ZA(:) = 0.0
  wHERE     (PX(:)<=PX1(:))
     PDELTA_VEC(:) = PA
  ELSEWHERE (PX(:)>=PX2(:))
     PDELTA_VEC(:) = PB
  ELSEWHERE
     ZA(:)         = 6.0*(PA-PB)/(PX2(:)-PX1(:))**3
     PDELTA_VEC(:) =   PA + ZA(:)*PX1(:)**2*(PX1(:)/6.0 - 0.5*PX2(:))           &
          + ZA(:)*PX1(:)*PX2(:)*                          (PX(:))    &
          - (0.5*ZA(:)*(PX1(:)+PX2(:)))*                  (PX(:)**2) & 
          + (ZA(:)/3.0)*                                  (PX(:)**3)
  END WHERE
  RETURN
!
END FUNCTION DELTA_VEC
!
!-------------------------------------------------------------------------------
!
!###########################
SUBROUTINE gaulag(x,w,n,alf)
!###########################
  INTEGER n,MAXIT
  REAL alf,w(n),x(n)
  REAL(kind=8) :: EPS
  PARAMETER (EPS=3.D-14,MAXIT=10)
  INTEGER i,its,j
  REAL ai
  REAL(kind=8) :: p1,p2,p3,pp,z,z1
!
  REAL SUMW
!
  do 13 i=1,n
     if(i.eq.1)then
        z=(1.+alf)*(3.+.92*alf)/(1.+2.4*n+1.8*alf)
     else if(i.eq.2)then
        z=z+(15.+6.25*alf)/(1.+.9*alf+2.5*n)
     else
        ai=i-2
        z=z+((1.+2.55*ai)/(1.9*ai)+1.26*ai*alf/(1.+3.5*ai))* &
             (z-x(i-2))/(1.+.3*alf)
     endif
     do 12 its=1,MAXIT
        p1=1.d0
        p2=0.d0
        do 11 j=1,n
           p3=p2
           p2=p1
           p1=((2*j-1+alf-z)*p2-(j-1+alf)*p3)/j
11      continue
        pp=(n*p1-(n+alf)*p2)/z
        z1=z
        z=z1-p1/pp
        if(abs(z-z1).le.EPS)goto 1
12   continue
1    x(i)=z
     w(i)=-exp(gammln(alf+n)-gammln(float(n)))/(pp*n*p2)
13 continue
!
! NORMALISATION
!
  SUMW = 0.0
  DO 14 I=1,N
     SUMW = SUMW + W(I)
14 CONTINUE
  DO 15 I=1,N
     W(I) = W(I)/SUMW
15 CONTINUE
!
  return
END SUBROUTINE gaulag
!
!------------------------------------------------------------------------------
!
!##########################################
SUBROUTINE gauher(x,w,n)
!##########################################
  INTEGER n,MAXIT
  REAL w(n),x(n)
  REAL(kind=8) :: EPS,PIM4
  PARAMETER (EPS=3.D-14,PIM4=.7511255444649425D0,MAXIT=10)
  INTEGER i,its,j,m
  REAL(kind=8) :: p1,p2,p3,pp,z,z1
!
  REAL SUMW
!
  m=(n+1)/2
  do 13 i=1,m
     if(i.eq.1)then
        z=sqrt(float(2*n+1))-1.85575*(2*n+1)**(-.16667)
     else if(i.eq.2)then
        z=z-1.14*n**.426/z
     else if (i.eq.3)then
        z=1.86*z-.86*x(1)
     else if (i.eq.4)then
        z=1.91*z-.91*x(2)
     else
        z=2.*z-x(i-2)
     endif
     do 12 its=1,MAXIT
        p1=PIM4
        p2=0.d0
        do 11 j=1,n
           p3=p2
           p2=p1
           p1=z*sqrt(2.d0/j)*p2-sqrt(dble(j-1)/dble(j))*p3
11      continue
        pp=sqrt(2.d0*n)*p2
        z1=z
        z=z1-p1/pp
        if(abs(z-z1).le.EPS)goto 1
12   continue
1    x(i)=z
     x(n+1-i)=-z
     pp=pp/PIM4 ! NORMALIZATION
     w(i)=2.0/(pp*pp)
     w(n+1-i)=w(i)
13 continue
!
! NORMALISATION
!
  SUMW = 0.0
  DO 14 I=1,N
     SUMW = SUMW + W(I)
14 CONTINUE
  DO 15 I=1,N
     W(I) = W(I)/SUMW
15 CONTINUE
!
  return
END SUBROUTINE gauher
!
!------------------------------------------------------------------------------
