splegend.f

Go to the documentation of this file.

C> @file
C>
C> Compute Legendre polynomials
C> @author IREDELL @date 92-10-31
 
C> Evaluates the orthonormal associated Legendre polynomials in the
C> spectral domain at a given latitude. Subprogram splegend should
C> be called already. If l is the zonal wavenumber, N is the total
C> wavenumber, and EPS(L,N)=SQRT((N**2-L**2)/(4*N**2-1)) then the
C> following bootstrapping formulas are used:
C>
C> <pre>
C> PLN(0,0)=SQRT(0.5)
C> PLN(L,L)=PLN(L-1,L-1)*CLAT*SQRT(FLOAT(2*L+1)/FLOAT(2*L))
C> PLN(L,N)=(SLAT*PLN(L,N-1)-EPS(L,N-1)*PLN(L,N-2))/EPS(L,N)
C> </pre>
C>
C> Synthesis at the pole needs only two zonal wavenumbers. Scalar
C> fields are synthesized with zonal wavenumber 0 while vector
C> fields are synthesized with zonal wavenumber 1. (Thus polar
C> vector fields are implicitly divided by clat.) The following
C> bootstrapping formulas are used at the pole:
C>
C> <pre>
C> PLN(0,0)=SQRT(0.5)
C> PLN(1,1)=SQRT(0.75)
C> PLN(L,N)=(PLN(L,N-1)-EPS(L,N-1)*PLN(L,N-2))/EPS(L,N)
C> </pre>
C>
C> PROGRAM HISTORY LOG:
C> -  91-10-31  MARK IREDELL
C> - 98-06-10  MARK IREDELL  GENERALIZE PRECISION
C>
C> @param I        - INTEGER SPECTRAL DOMAIN SHAPE
C>                (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
C> @param M        - INTEGER SPECTRAL TRUNCATION
C> @param SLAT     - REAL SINE OF LATITUDE
C> @param CLAT     - REAL COSINE OF LATITUDE
C> @param EPS      - REAL ((M+1)*((I+1)*M+2)/2) SQRT((N**2-L**2)/(4*N**2-1))
C> @param EPSTOP   - REAL (M+1) SQRT((N**2-L**2)/(4*N**2-1)) OVER TOP
C> @param[out] PLN - REAL ((M+1)*((I+1)*M+2)/2) LEGENDRE POLYNOMIAL
C> @param[out] PLNTOP - REAL (M+1) LEGENDRE POLYNOMIAL OVER TOP
C>
      SUBROUTINE splegend(I,M,SLAT,CLAT,EPS,EPSTOP,PLN,PLNTOP)
 
CFPP$ NOCONCUR R
      REAL EPS((M+1)*((I+1)*M+2)/2),EPSTOP(M+1)
      REAL PLN((M+1)*((I+1)*M+2)/2),PLNTOP(M+1)
      REAL(KIND=selected_real_kind(15,45)):: dln((m+1)*((i+1)*m+2)/2)
      REAL :: TINYREAL=tiny(1.0), rdln1, rdln2
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  ITERATIVELY COMPUTE PLN WITHIN SPECTRAL DOMAIN AT POLE
      m1=m+1
      m2=2*m+i+1
      mx=(m+1)*((i+1)*m+2)/2
      IF(abs(clat).LT.tinyreal) THEN
        dln(1)=sqrt(0.5)
        IF(m.GT.0) THEN
          dln(m1+1)=sqrt(0.75)
          dln(2)=slat*dln(1)/eps(2)
        ENDIF
        IF(m.GT.1) THEN
          dln(m1+2)=slat*dln(m1+1)/eps(m1+2)
          dln(3)=(slat*dln(2)-eps(2)*dln(1))/eps(3)
          DO n=3,m
            k=1+n
            dln(k)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/eps(k)
            k=m1+n
            dln(k)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/eps(k)
          ENDDO
          IF(i.EQ.1) THEN
            k=m2
            dln(k)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/eps(k)
          ENDIF
          DO k=m2+1,mx
            dln(k)=0.
          ENDDO
        ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COMPUTE POLYNOMIALS OVER TOP OF SPECTRAL DOMAIN
        k=m1+1
        rdln1=real(dln(k-1))
        rdln2=real(dln(k-2))
        plntop(1)=(slat*rdln1-eps(k-1)*rdln2)/epstop(1)
        IF(m.GT.0) THEN
          k=m2+1
          rdln1=real(dln(k-1))
          rdln2=real(dln(k-2))
          plntop(2)=(slat*rdln1-eps(k-1)*rdln2)/epstop(2)
          DO l=2,m
            plntop(l+1)=0.
          ENDDO
        ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  ITERATIVELY COMPUTE PLN(L,L) (BOTTOM HYPOTENUSE OF DOMAIN)
      ELSE
        nml=0
        k=1
        dln(k)=sqrt(0.5)
        DO l=1,m+(i-1)*nml
          kp=k
          k=l*(2*m+(i-1)*(l-1))/2+l+nml+1
          dln(k)=dln(kp)*clat*sqrt(float(2*l+1)/float(2*l))
        ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COMPUTE PLN(L,L+1) (DIAGONAL NEXT TO BOTTOM HYPOTENUSE OF DOMAIN)
        nml=1
CDIR$ IVDEP
        DO l=0,m+(i-1)*nml
          k=l*(2*m+(i-1)*(l-1))/2+l+nml+1
          dln(k)=slat*dln(k-1)/eps(k)
        ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COMPUTE REMAINING PLN IN SPECTRAL DOMAIN
        DO nml=2,m
CDIR$ IVDEP
          DO l=0,m+(i-1)*nml
            k=l*(2*m+(i-1)*(l-1))/2+l+nml+1
            dln(k)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/eps(k)
          ENDDO
        ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  COMPUTE POLYNOMIALS OVER TOP OF SPECTRAL DOMAIN
        DO l=0,m
          nml=m+1+(i-1)*l
          k=l*(2*m+(i-1)*(l-1))/2+l+nml+1
          rdln1=real(dln(k-1))
          rdln2=real(dln(k-2))
          plntop(l+1)=(slat*rdln1-eps(k-1)*rdln2)/epstop(l+1)
        ENDDO
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  RETURN VALUES
      DO k=1,mx
        pln(k)=real(dln(k))
      ENDDO
      RETURN
      END