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      INTEGER FUNCTION IGGLAT(KLAT, PGAUSS, KPR, KERR)
12C
13C---->
14C**** *IGGLAT*
15C
16C     PURPOSE
17C     _______
18C
19C     Compute Gaussian lines of latitude for a given truncation.
20C
21C
22C     INTERFACE
23C     _________
24C
25C     IERR = IGGLAT(KLAT, PGAUSS, KPR, KERR)
26C
27C
28C     Input parameters
29C     ________________
30C
31C     KLAT - The number of lines of latitude from pole to pole.
32C
33C     KPR  - Debug print switch:
34C            0, No debugging output.
35C            1, Produce debugging output.
36C
37C     KERR - Error control flag.
38C            -ve, No error message. Return error code.
39C            0  , Hard failure with error message.
40C            +ve, Print error message. Return error code.
41C
42C
43C     Output parameters
44C     ________________
45C
46C     PGAUSS - The Gaussian lines of latitude for truncation KLAT/2
47C
48C
49C     Return value
50C     ____________
51C
52C     0     - OK
53C     24701 - The calculation of Gaussian lines of latitude failed.
54C
55C
56C     Common block usage
57C     __________________
58C
59C     None
60C
61C
62C     EXTERNALS
63C     _________
64C
65C     INTLOG(R) - Logs messages.
66C     IGBESS    - This routine is used to generate initial
67C                 approximations to the Gaussian latitudes.
68C
69C
70C     METHOD
71C     ______
72C
73C     IGBESS is used to provide zeros of the Bessel function J0,
74C     which are used as starting approximations to the Gaussian
75C     latitudes. Newton iteration is used to generate the latitudes
76C     from these approximations.
77C
78C     Calculated latitudes are stored internally and reused if a
79C     later call asks for the same calculation.
80C
81C
82C     REFERENCE
83C     _________
84C
85C     None
86C
87C
88C     COMMENTS
89C     ________
90C
91C     This routine is adapted from that in the old Marsint library.
92C     The interface and the variable names have been modified.
93C
94C
95C     AUTHOR
96C     ______
97C
98C     K. Fielding      *ECMWF*      Oct 1993
99C
100C
101C     MODIFICATIONS
102C     _____________
103C
104C     J.Chambers        ECMWF       Aug 1998
105C     Add reuse of calculated latitudes.
106C
107C----<
108C     -------------------------------------------------------
109C*    Section 0. Definition of variables.
110C     -------------------------------------------------------
111C
112      IMPLICIT NONE
113C
114#include "jparams.h"
115#include "parim.h"
116C
117C     Function arguments
118C
119      INTEGER KLAT, KPR, KERR
120      REAL PGAUSS (*)
121C
122C     Parameters
123C
124      INTEGER JPROUTINE
125      PARAMETER (JPROUTINE = 24700)
126C     -------------------------------------------------------
127C     The convergence criteria is machine dependent
128C     -------------------------------------------------------
129C
130      REAL PACRCY
131#ifdef REAL_8
132      PARAMETER (PACRCY = 1.0E-14)
133#else
134      PARAMETER (PACRCY = 1.0E-7)
135#endif
136C
137C     Local variables
138C
139      LOGICAL LDEBUG
140      INTEGER ITRUNC, IERR
141      INTEGER JLAT, JITER, JLN, LOOP, JOLDLAT
142      REAL ZRADDEG, ZCON, ZLAT, ZROOT
143      REAL ZKM1, ZKM2, ZLN, ZFUNC, ZDERIV, ZMOVE
144      REAL ROLDGS
145      DIMENSION ROLDGS(JPGTRUNC*2)
146C
147      SAVE JOLDLAT, ROLDGS
148C
149C     Externals
150C
151      INTEGER IGBESS
152C
153C     -------------------------------------------------------
154C*    Section 1. Set constants and get initial approximation.
155C     -------------------------------------------------------
156C
157  100 CONTINUE
158C
159      ZFUNC  = 0.
160      IGGLAT = 0
161      LDEBUG = ( KPR.GE.1 )
162C
163      IF( LDEBUG ) THEN
164        CALL INTLOG(JP_DEBUG,'IGGLAT: Section 1.',JPQUIET)
165        CALL INTLOG(JP_DEBUG,
166     X  'IGGLAT: No.lines lat from pole to pole = ',KLAT)
167      ENDIF
168C
169C     Using previously calculated values if truncation is the same.
170C
171      IF( KLAT.EQ.JOLDLAT ) THEN
172        IF( LDEBUG ) CALL INTLOG(JP_DEBUG,
173     X    'IGGLAT: Using previously calculated values',JPQUIET)
174        DO LOOP = 1, KLAT
175          PGAUSS(LOOP) = ROLDGS(LOOP)
176        ENDDO
177        GOTO 900
178      ENDIF
179C
180      JOLDLAT = KLAT
181C
182      ZRADDEG = 180.0/PPI
183C
184      ZCON   = (PPONE - (PPTWO / PPI) ** 2) * PPQUART
185C
186      ZLAT   = KLAT
187      ITRUNC = KLAT/2
188C
189      IERR = IGBESS(ITRUNC, PGAUSS, KPR, KERR)
190C
191      IF( IERR.GT.0 ) THEN
192         IGGLAT = IERR
193         GOTO 900
194      ENDIF
195C
196C     -------------------------------------------------------
197C*    Section 2. Compute abscissae
198C     -------------------------------------------------------
199C
200  200 CONTINUE
201C
202      IF( LDEBUG ) CALL INTLOG(JP_DEBUG,'IGGLAT: Section 2.',JPQUIET)
203C
204      DO 240 JLAT = 1, ITRUNC
205C
206C       First approximation for ZROOT
207C
208        ZROOT = COS(PGAUSS(JLAT) / SQRT( (ZLAT+PPHALF)**2 + ZCON) )
209C
210C       Perform loop of Newton iterations
211C
212        DO 220 JITER = 1, JPMAXITER
213C
214          ZKM2 = PPONE
215          ZKM1 = ZROOT
216C
217C         Compute Legendre polynomial
218C
219          DO 210 JLN = 2, KLAT
220C
221            ZLN = JLN
222C
223            ZFUNC = ( (PPTWO * ZLN - PPONE) * ZROOT * ZKM1 -
224     1              (ZLN - PPONE) * ZKM2) / ZLN
225C
226            ZKM2 = ZKM1
227            ZKM1 = ZFUNC
228C
229  210     CONTINUE
230C
231C         Perform Newton iteration
232C
233          ZDERIV = (ZLAT * (ZKM2 - ZROOT * ZFUNC) ) /
234     1             (PPONE - ZROOT ** 2)
235C
236          ZMOVE = ZFUNC / ZDERIV
237          ZROOT = ZROOT - ZMOVE
238C
239C         Leave iteration loop when sufficient accuracy achieved.
240C
241          IF( ABS(ZMOVE).LE.PACRCY ) GOTO 230
242C
243  220   CONTINUE
244C
245C       Routine fails if no convergence after JPMAXITER iterations.
246C
247        IGGLAT = JPROUTINE + 1
248C
249        IF( KERR.GE.0 ) CALL INTLOG(JP_ERROR,
250     X    'IGGLAT: Calculation of Gaussian lats failed.',JPQUIET)
251C
252        IF( KERR.EQ.0 ) CALL INTLOG(JP_FATAL,
253     X    'IGGLAT: interpolation failed.',IGGLAT)
254C
255        GOTO 900
256C
257  230   CONTINUE
258C
259C*      Set North and South values using symmetry.
260C
261        PGAUSS(JLAT) = ASIN(ZROOT) * ZRADDEG
262        PGAUSS(KLAT + 1 - JLAT) = - PGAUSS(JLAT)
263C
264  240 CONTINUE
265C
266      IF( KLAT.NE.(ITRUNC*2) ) PGAUSS(ITRUNC + 1) = PPZERO
267C
268C     Store calculated values for re-use
269C
270      DO LOOP = 1, KLAT
271        ROLDGS(LOOP) = PGAUSS(LOOP)
272      ENDDO
273C
274C     -------------------------------------------------------
275C*    Section 9. Return to calling routine.
276C     -------------------------------------------------------
277C
278  900 CONTINUE
279C
280      IF( LDEBUG ) CALL INTLOG(JP_DEBUG,'IGGLAT: Section 9.',JPQUIET)
281C
282      RETURN
283      END
284