1C Copyright 1981-2016 ECMWF.
2C
3C This software is licensed under the terms of the Apache Licence
4C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5C
6C In applying this licence, ECMWF does not waive the privileges and immunities
7C granted to it by virtue of its status as an intergovernmental organisation
8C nor does it submit to any jurisdiction.
9C
10
11      SUBROUTINE JSPLEG1( PLEG, PLAT, KTRUNC)
12      IMPLICIT NONE
13C
14C---->
15C**** JSPLEG1 - Routine to calculate legendre functions
16C
17C     Purpose
18C     --------
19C
20C     This routine calculates the legendre functions for one latitude.
21C     (but not their derivatives)
22C
23C
24C     Interface
25C     ----------
26C
27C     CALL JSPLEG1( PLEG, PLAT, KTRUNC)
28C
29C
30C     Input parameters
31C     ----------------
32C
33C     PLAT      - Latitude in radians
34C     KTRUNC    - Spectral truncation
35C
36C
37C     Output parameters
38C     -----------------
39C
40C     PLEG      - Array of legendre functions for one latitude.
41C                 The array must be at least (KTRUNC+1)*(KTRUNC+4)/2
42C                 words long.
43C
44C
45C     Common block usage
46C     ------------------
47C
48C     None
49C
50C
51C     Method
52C     ------
53C
54C     Recurrence relation with explicit relations for P(m,m) and
55C     P(m,m+1)
56C
57C
58C     Externals
59C     ---------
60C
61C     None
62C
63C
64C     Reference
65C     ---------
66C
67C     None
68C
69C
70C     Comments
71C     --------
72C
73C     Rewritten from SPLEG1 to avoid using common blocks.
74C     Work arrays ZHLPx currently dimensioned for T213 using
75C     parameter JPTRP1.
76C
77C
78C     AUTHOR
79C     ------
80C
81C     J.D.Chambers         ECMWF        9 November 1993
82C
83C
84C     Modifications
85C     -------------
86C
87C     None
88C
89C----<
90C     _______________________________________________________
91C
92C*    Section 0. Definition of variables.
93C     _______________________________________________________
94C
95C*    Prefix conventions for variable names
96C
97C     Logical      L (but not LP), global or common.
98C                  O, dummy argument
99C                  G, local variable
100C                  LP, parameter.
101C     Character    C, global or common.
102C                  H, dummy argument
103C                  Y (but not YP), local variable
104C                  YP, parameter.
105C     Integer      M and N, global or common.
106C                  K, dummy argument
107C                  I, local variable
108C                  J (but not JP), loop control
109C                  JP, parameter.
110C     Real         A to F and Q to X, global or common.
111C                  P (but not PP), dummy argument
112C                  Z, local variable
113C                  PP, parameter.
114#include "jparams.h"
115C
116C     Subroutine arguments
117C
118      REAL PLEG, PLAT
119      DIMENSION PLEG(*)
120      INTEGER KTRUNC
121C
122C     Local variables
123C
124      REAL ZHLP1, ZHLP2, ZHLP3
125      DIMENSION ZHLP1(JPTRP1), ZHLP2(JPTRP1), ZHLP3(JPTRP1)
126      INTEGER ITOUT1, I1M, ILM, JM, JCN, IM2
127      REAL ZSIN, ZCOS, ZF1M, ZRE1, ZF2M, ZM, ZN, ZE1, ZE2
128C
129C     _______________________________________________________
130C
131C
132C*    Section 1. Initialization.
133C     _______________________________________________________
134C
135  100 CONTINUE
136C
137      ITOUT1 = KTRUNC+1
138      ZSIN   = SIN(PLAT)
139      ZCOS   = SQRT(1.-ZSIN*ZSIN)
140C
141C     Step 1.        M = 0, N = 0 and N = 1
142C
143      ILM     = 2
144      PLEG(1) = 1.0
145      ZF1M    = SQRT(3.0)
146      PLEG(2) = ZF1M*ZSIN
147C     _______________________________________________________
148C
149C
150C     Step 2.       Sum for M = 0 to T (T = truncation)
151C     _______________________________________________________
152C
153 200  CONTINUE
154C
155      DO 205 JM = 2,ITOUT1
156        ZM        = JM-1
157        ZHLP1(JM) = SQRT(2.*ZM+3.)
158        ZHLP2(JM) = 1./SQRT(2.*ZM)
159 205  CONTINUE
160      ZHLP1(1) = SQRT(3.)
161C
162      DO 570 JM = 1,ITOUT1
163C
164        I1M  = JM-1
165        ZM   = I1M
166        ZRE1 = ZHLP1(JM)
167        ZE1  = 1./ZRE1
168C     _______________________________________________________
169C
170C
171C     Step 3.       M > 0 only
172C     _______________________________________________________
173C
174 300    CONTINUE
175C
176        IF (I1M.NE.0) THEN
177          ZF2M = ZF1M*ZCOS*ZHLP2(JM)
178          ZF1M = ZF2M*ZRE1
179C     _______________________________________________________
180C
181C
182C     Step 4.       N = M and N = M+1
183C     _______________________________________________________
184C
185 400    CONTINUE
186C
187          ILM       = ILM+1
188          PLEG(ILM) = ZF2M
189          ILM       = ILM+1
190          PLEG(ILM) = ZF1M*ZSIN
191C
192C         When output truncation is reached, return to calling program
193          IF ( JM .EQ. ITOUT1 ) GOTO 570
194C
195        ENDIF
196C     _______________________________________________________
197C
198C
199C     Step 5.       Sum for N = M+2 to T+1
200C     _______________________________________________________
201C
202 500    CONTINUE
203C
204        IM2 = I1M+2
205C
206        DO 520 JCN = IM2,ITOUT1
207          ZN         = JCN
208          ZHLP3(JCN) = SQRT((4.*ZN*ZN-1.)/(ZN*ZN-ZM*ZM))
209 520    CONTINUE
210C
211        DO 560 JCN = IM2,ITOUT1
212C
213          ZE2       = ZHLP3(JCN)
214          ILM       = ILM+1
215          PLEG(ILM) = ZE2*(ZSIN*PLEG(ILM-1)-ZE1*PLEG(ILM-2))
216          ZE2       = 1./ZE2
217          ZE1       = ZE2
218C
219 560    CONTINUE
220C
221 570  CONTINUE
222C
223C     _______________________________________________________
224C
225C
226C*    Section 9. Return to calling routine.
227C     _______________________________________________________
228C
229  900 CONTINUE
230C
231      RETURN
232      END
233