w3fa12.f

Go to the documentation of this file.

C> @file
C> @brief Computes legendre polynomials at a given latitude.
C> @author Joe Sela @date 1980-10-28
 
C> Subroutine computes legendre polynomials at a given latitude.
C>
C> Program history log:
C> - Joe Sela 1980-10-20
C> - Ralph Jones 1984-06-01 Change to ibm vs fortran.
C> - Ralph Jones 1993-04-12 Changes for cray, double precision to real.
C>
C> @param[out] PLN Real locations contain legendre
C> polynomials, size is (jcap+2)*(jcap+1)
C> @param[in] COLRAD Colatitude in radians of desired point.
C> @param[in] JCAP For rhomboiadal truncation of zonal wave
C> @param[in] EPS Coeff. used in recursion equation.
C> Dimension of eps is (jcap+2)*(jcap+1)
C>
C> @author Joe Sela @date 1980-10-28
       SUBROUTINE w3fa12(PLN,COLRAD,JCAP,EPS)
       REAL A
       REAL B
       REAL COLRAD
       REAL COS2
       REAL EPS(*)
       REAL FL
       REAL PROD
       REAL P1
       REAL P2
       REAL P3
       REAL SINLAT
       REAL PLN(*)
C
       SAVE
C
         sinlat = cos(colrad)
         cos2   = 1.0 - sinlat * sinlat
         prod   = 1.0
         a      = 1.0
         b      = 0.0
         jcap1  = jcap+1
         jcap2  = jcap+2
C
         DO 300 ll = 1,jcap1
           l   = ll - 1
           fl  = l
           jle = l * jcap2
           IF (l.EQ.0) GO TO 100
             a = a + 2.0
             b = b + 2.0
             prod = prod * cos2 * a / b
 100       CONTINUE
           p1 = sqrt(0.5 * prod)
           pln(jle+1) = p1
           p2 = sqrt(2.0 * fl + 3.0) * sinlat * p1
           pln(jle+2) = p2
C
           DO 200 n = 3,jcap2
             lindex = jle + n
             p3 = (sinlat*p2 - eps(lindex-1)*p1)/eps(lindex)
             pln(lindex) = p3
             p1 = p2
             p2 = p3
200        CONTINUE
300      CONTINUE
         RETURN
C
       END