NCEPLIBS-sp  2.3.3
splegend.f
Go to the documentation of this file.
1 C> @file
2 C>
3 C> Compute legendre polynomials
4 C> @author IREDELL @date 92-10-31
5 
6 C> Evaluates the orthonormal associated legendre polynomials in the
7 C> spectral domain at a given latitude. Subprogram splegend should
8 C> be called already. If l is the zonal wavenumber, N is the total
9 C> wavenumber, and EPS(L,N)=SQRT((N**2-L**2)/(4*N**2-1)) then the
10 C> following bootstrapping formulas are used:
11 C>
12 C> <pre>
13 C> PLN(0,0)=SQRT(0.5)
14 C> PLN(L,L)=PLN(L-1,L-1)*CLAT*SQRT(FLOAT(2*L+1)/FLOAT(2*L))
15 C> PLN(L,N)=(SLAT*PLN(L,N-1)-EPS(L,N-1)*PLN(L,N-2))/EPS(L,N)
16 C> </pre>
17 C>
18 C> Synthesis at the pole needs only two zonal wavenumbers. Scalar
19 C> fields are synthesized with zonal wavenumber 0 while vector
20 C> fields are synthesized with zonal wavenumber 1. (Thus polar
21 C> vector fields are implicitly divided by clat.) The following
22 C> bootstrapping formulas are used at the pole:
23 C>
24 C> <pre>
25 C> PLN(0,0)=SQRT(0.5)
26 C> PLN(1,1)=SQRT(0.75)
27 C> PLN(L,N)=(PLN(L,N-1)-EPS(L,N-1)*PLN(L,N-2))/EPS(L,N)
28 C> </pre>
29 C>
30 C> PROGRAM HISTORY LOG:
31 C> - 91-10-31 MARK IREDELL
32 C> - 98-06-10 MARK IREDELL GENERALIZE PRECISION
33 C>
34 C> @param I - INTEGER SPECTRAL DOMAIN SHAPE
35 C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
36 C> @param M - INTEGER SPECTRAL TRUNCATION
37 C> @param SLAT - REAL SINE OF LATITUDE
38 C> @param CLAT - REAL COSINE OF LATITUDE
39 C> @param EPS - REAL ((M+1)*((I+1)*M+2)/2) SQRT((N**2-L**2)/(4*N**2-1))
40 C> @param EPSTOP - REAL (M+1) SQRT((N**2-L**2)/(4*N**2-1)) OVER TOP
41 C> @param[out] PLN - REAL ((M+1)*((I+1)*M+2)/2) LEGENDRE POLYNOMIAL
42 C> @param[out] PLNTOP - REAL (M+1) LEGENDRE POLYNOMIAL OVER TOP
43 C>
44  SUBROUTINE splegend(I,M,SLAT,CLAT,EPS,EPSTOP,PLN,PLNTOP)
45 
46 CFPP$ NOCONCUR R
47  REAL EPS((M+1)*((I+1)*M+2)/2),EPSTOP(M+1)
48  REAL PLN((M+1)*((I+1)*M+2)/2),PLNTOP(M+1)
49  REAL(KIND=selected_real_kind(15,45)):: dln((m+1)*((i+1)*m+2)/2)
50 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51 C ITERATIVELY COMPUTE PLN WITHIN SPECTRAL DOMAIN AT POLE
52  m1=m+1
53  m2=2*m+i+1
54  mx=(m+1)*((i+1)*m+2)/2
55  IF(clat.EQ.0.) THEN
56  dln(1)=sqrt(0.5)
57  IF(m.GT.0) THEN
58  dln(m1+1)=sqrt(0.75)
59  dln(2)=slat*dln(1)/eps(2)
60  ENDIF
61  IF(m.GT.1) THEN
62  dln(m1+2)=slat*dln(m1+1)/eps(m1+2)
63  dln(3)=(slat*dln(2)-eps(2)*dln(1))/eps(3)
64  DO n=3,m
65  k=1+n
66  dln(k)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/eps(k)
67  k=m1+n
68  dln(k)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/eps(k)
69  ENDDO
70  IF(i.EQ.1) THEN
71  k=m2
72  dln(k)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/eps(k)
73  ENDIF
74  DO k=m2+1,mx
75  dln(k)=0.
76  ENDDO
77  ENDIF
78 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79 C COMPUTE POLYNOMIALS OVER TOP OF SPECTRAL DOMAIN
80  k=m1+1
81  plntop(1)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/epstop(1)
82  IF(m.GT.0) THEN
83  k=m2+1
84  plntop(2)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/epstop(2)
85  DO l=2,m
86  plntop(l+1)=0.
87  ENDDO
88  ENDIF
89 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90 C ITERATIVELY COMPUTE PLN(L,L) (BOTTOM HYPOTENUSE OF DOMAIN)
91  ELSE
92  nml=0
93  k=1
94  dln(k)=sqrt(0.5)
95  DO l=1,m+(i-1)*nml
96  kp=k
97  k=l*(2*m+(i-1)*(l-1))/2+l+nml+1
98  dln(k)=dln(kp)*clat*sqrt(float(2*l+1)/float(2*l))
99  ENDDO
100 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101 C COMPUTE PLN(L,L+1) (DIAGONAL NEXT TO BOTTOM HYPOTENUSE OF DOMAIN)
102  nml=1
103 CDIR$ IVDEP
104  DO l=0,m+(i-1)*nml
105  k=l*(2*m+(i-1)*(l-1))/2+l+nml+1
106  dln(k)=slat*dln(k-1)/eps(k)
107  ENDDO
108 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109 C COMPUTE REMAINING PLN IN SPECTRAL DOMAIN
110  DO nml=2,m
111 CDIR$ IVDEP
112  DO l=0,m+(i-1)*nml
113  k=l*(2*m+(i-1)*(l-1))/2+l+nml+1
114  dln(k)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/eps(k)
115  ENDDO
116  ENDDO
117 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118 C COMPUTE POLYNOMIALS OVER TOP OF SPECTRAL DOMAIN
119  DO l=0,m
120  nml=m+1+(i-1)*l
121  k=l*(2*m+(i-1)*(l-1))/2+l+nml+1
122  plntop(l+1)=(slat*dln(k-1)-eps(k-1)*dln(k-2))/epstop(l+1)
123  ENDDO
124  ENDIF
125 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126 C RETURN VALUES
127  DO k=1,mx
128  pln(k)=dln(k)
129  ENDDO
130  RETURN
131  END
splegend
subroutine splegend(I, M, SLAT, CLAT, EPS, EPSTOP, PLN, PLNTOP)
Evaluates the orthonormal associated legendre polynomials in the spectral domain at a given latitude.
Definition: splegend.f:45