supol_mod.F90

Go to the documentation of this file.

MODULE supol_mod
CONTAINS
SUBROUTINE supol(KNSMAX,DDMU,DDPOL,DDA,DDB,DDC,DDD,DDE,DDF,DDG,DDH,DDI)

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

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

!**   Interface.
!     ----------
!        *CALL* *SUPOL(KNSMAX,DDMU,DDPOL,DDA,DDB,DDC,DDD,DDE
!        ,DDF,DDG,DDH,DDI)

!        Explicit arguments :
!        --------------------
!              KNSMAX   :  Truncation  (triangular)
!              DDMU     :  Abscissa at which the polynomials are computed (mu)
!              DDPOL    :  Polynomials (the first index is m and the second n)


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

!     Method.
!     -------
!        See documentation

!     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.
!     ------------------------------------------------------------------

USE parkind1  ,ONLY : jpim     ,jprb
USE parkind2  ,ONLY : jprh

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)  :: KNSMAX
REAL(KIND=JPRH)   ,INTENT(IN)  :: DDMU
REAL(KIND=JPRH)   ,INTENT(IN)  :: DDC(0:knsmax,0:knsmax)
REAL(KIND=JPRH)   ,INTENT(IN)  :: DDD(0:knsmax,0:knsmax)
REAL(KIND=JPRH)   ,INTENT(IN)  :: DDE(0:knsmax,0:knsmax)
REAL(KIND=JPRH)   ,INTENT(IN)  :: DDA(0:knsmax),DDB(0:knsmax),DDF(0:knsmax)
REAL(KIND=JPRH)   ,INTENT(IN)  :: DDG(0:knsmax),DDH(0:knsmax),DDI(0:knsmax)
REAL(KIND=JPRH)   ,INTENT(OUT) :: DDPOL(0:knsmax,0:knsmax)

REAL(KIND=JPRH) :: DLX,DLSITA,DL1SITA,DLKM2,DLKM1,DLK,DL1,DLS

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

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

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

dlx=ddmu
dlsita=sqrt(1.0_jprb-dlx*dlx)

! IF WE ARE LESS THAN 1Meter FROM THE POLE,
IF(abs(REAL(dlsita,kind(z))) <= sqrt(epsilon(z)))THEN
  dlx=1._jprb
  dlsita=0._jprb
  dl1sita=0._jprb
ELSE
  dl1sita=1.0_jprb/dlsita
ENDIF
dlkm2=1._jprb
dlkm1=dlx
ddpol(0,0)=dlkm2
ddpol(0,1)=dlkm1*dda(1)
ddpol(1,1)=dlsita*ddb(1)
DO jn=2,knsmax
  dlk=ddf(jn)*dlx*dlkm1-ddg(jn)*dlkm2
  dl1=ddi(jn)*(dlkm1-dlx*dlk)*dl1sita
  ddpol(0,jn)=dlk*dda(jn)
  ddpol(1,jn)=dl1*ddb(jn)
  dlkm2=dlkm1
  dlkm1=dlk
ENDDO

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

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

dls=dl1sita*tiny(dls)

!OCL SCALAR
DO jn=2,knsmax
  ddpol(jn,jn)=ddpol(jn-1,jn-1)*dlsita*ddh(jn)
  IF ( abs(ddpol(jn,jn))  <  dls ) ddpol(jn,jn)=0.0_jprb
ENDDO

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

!*       3. General recurrence.
!           -------------------

DO jn=3,knsmax
!DIR$ IVDEP
!OCL NOVREC
  DO jm=2,jn-1
    ddpol(jm,jn)=ddc(jm,jn)*ddpol(jm-2,jn-2)&
     &-ddd(jm,jn)*ddpol(jm-2,jn-1)*dlx &
     &+dde(jm,jn)*ddpol(jm  ,jn-1)*dlx
  ENDDO
ENDDO

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

END SUBROUTINE supol
END MODULE supol_mod