! (C) Copyright 1987- ECMWF.
! (C) Copyright 1987- Meteo-France.
! 
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!

MODULE SUPOL_MOD
CONTAINS
SUBROUTINE SUPOL(KNSMAX,PDDMU,PFN,PDDPOL)

!**** *SUPOL * - Routine to compute the Legendre polynomials

!     Purpose.
!     --------
!           For a given value of mu, computes the Legendre polynomials.

!**   Interface.
!     ----------
!        *CALL* *SUPOL(...)

!        Explicit arguments :
!        --------------------
!        KNSMAX    :  Truncation  (triangular)                            [in]
!        PDDMU     :  Abscissa at which the polynomials are computed (mu) [in]
!        PFN       :  Fourier coefficients of series expansion 
!                     for the ordinary Legendre polynomials               [in]
!        PDDPOL    :  Polynomials (the first index is m and the second n) [out]

!        Implicit arguments :   None
!        --------------------

!     Method.
!     -------
!        See documentation about spectral transforms
!         (doc (IDTS) by K. Yessad, appendix 3, or doc (NTA30) by M. Rochas)

!     Externals.
!     ----------

!     Reference.
!     ----------
!        ECMWF Research Department documentation of the IFS

!     Author.
!     -------
!        Mats Hamrud and Philippe Courtier  *ECMWF*

!     Modifications.
!     --------------
!        Original : 87-10-15
!        K. YESSAD (MAY 1998): modification to avoid underflow.
!        M.Hamrud      01-Oct-2003 CY28 Cleaning
!        R. El Khatib 11-Apr-2007 Emulation of vectorized quadruple precision
!                                 on NEC
!        K. YESSAD (NOV 2008): make consistent arp/SUPOLA and tfl/SUPOL.
!        Nils Wedi + Mats Hamrud, 2009-02-05 revised following Swarztrauber, 2002
!        R. El Khatib 30-Apr-2013 Open-MP parallelization
!      F. Vana  05-Mar-2015  Support for single precision
!     ------------------------------------------------------------------

USE EC_PARKIND ,ONLY : JPRD, JPIM
USE TPM_POL    ,ONLY : DDI, DDA, DDH, DDE, DDC, DDD

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)  :: KNSMAX
REAL(KIND=JPRD)   ,INTENT(IN)  :: PDDMU
REAL(KIND=JPRD)   ,INTENT(IN)  :: PFN(0:KNSMAX,0:KNSMAX)

REAL(KIND=JPRD)   ,INTENT(OUT) :: PDDPOL(0:KNSMAX,0:KNSMAX)

REAL(KIND=JPRD) :: ZDLX,ZDLX1,ZDLSITA,ZDL1SITA,ZDLS,ZDLK,ZDLLDN

INTEGER(KIND=JPIM) :: JM, JN, JK
REAL(KIND=JPRD) :: Z

!     ------------------------------------------------------------------

!*       1. First two columns.
!           ------------------

ZDLX=PDDMU
ZDLX1=ACOS(ZDLX)
ZDLSITA=SQRT(1.0_JPRD-ZDLX*ZDLX)

PDDPOL(0,0)=1._JPRD
ZDLLDN = 0.0_JPRD

! IF WE ARE LESS THAN 1Meter FROM THE POLE,
IF(ABS(REAL(ZDLSITA,KIND(Z))) <= SQRT(EPSILON(Z)))THEN
  ZDLX=1._JPRD
  ZDLSITA=0._JPRD
  ZDL1SITA=0._JPRD
ELSE
  ZDL1SITA=1.0_JPRD/ZDLSITA
ENDIF

!*          ordinary Legendre polynomials from series expansion
!           ---------------------------------------------------

! even N
!$OMP PARALLEL DO PRIVATE(JN,ZDLK,ZDLLDN,JK)
DO JN=2,KNSMAX,2
  ZDLK = 0.5_JPRD*PFN(JN,0)
  ZDLLDN = 0.0_JPRD
  ! represented by only even k
  DO JK=2,JN,2
    ! normalised ordinary Legendre polynomial == \overbar{P_n}^0
    ZDLK = ZDLK + PFN(JN,JK)*COS(DDI(JK)*ZDLX1)
    ! normalised associated Legendre polynomial == \overbar{P_n}^1
    ZDLLDN = ZDLLDN + DDA(JN)*PFN(JN,JK)*DDI(JK)*SIN(DDI(JK)*ZDLX1)
  ENDDO
  PDDPOL(0,JN) = ZDLK
  PDDPOL(1,JN) = ZDLLDN
ENDDO
!$OMP END PARALLEL DO
! odd N
!$OMP PARALLEL DO PRIVATE(JN,ZDLK,ZDLLDN,JK)
DO JN=1,KNSMAX,2
  ZDLK = 0.0_JPRD
  ZDLLDN = 0.0_JPRD
  ! represented by only odd k
  DO JK=1,JN,2
    ! normalised ordinary Legendre polynomial == \overbar{P_n}^0
    ZDLK = ZDLK + PFN(JN,JK)*COS(DDI(JK)*ZDLX1)
    ! normalised associated Legendre polynomial == \overbar{P_n}^1
    ZDLLDN = ZDLLDN + DDA(JN)*PFN(JN,JK)*DDI(JK)*SIN(DDI(JK)*ZDLX1)
  ENDDO
  PDDPOL(0,JN) = ZDLK
  PDDPOL(1,JN) = ZDLLDN
ENDDO
!$OMP END PARALLEL DO

!     ------------------------------------------------------------------

!*       2. Diagonal (the terms 0,0 and 1,1 have already been computed)
!           Belousov, equation (23)
!           -----------------------------------------------------------

ZDLS=ZDL1SITA*TINY(ZDLS)

#ifdef VPP
!OCL SCALAR
#endif
DO JN=2,KNSMAX
  PDDPOL(JN,JN)=PDDPOL(JN-1,JN-1)*ZDLSITA*DDH(JN)
  IF ( ABS(PDDPOL(JN,JN)) < ZDLS ) PDDPOL(JN,JN)=0.0_JPRD
ENDDO

!     ------------------------------------------------------------------

!*       3. General recurrence (Belousov, equation 17)
!           -----------------------------------------

DO JN=3,KNSMAX
!DIR$ IVDEP
!OCL NOVREC
  DO JM=2,JN-1
    PDDPOL(JM,JN)=DDC(JM,JN)*PDDPOL(JM-2,JN-2)&
     &-DDD(JM,JN)*PDDPOL(JM-2,JN-1)*ZDLX &
     &+DDE(JM,JN)*PDDPOL(JM  ,JN-1)*ZDLX
  ENDDO
ENDDO

!     ------------------------------------------------------------------

END SUBROUTINE SUPOL
END MODULE SUPOL_MOD
