1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
|
C Copyright 1981-2007 ECMWF
C
C Licensed under the GNU Lesser General Public License which
C incorporates the terms and conditions of version 3 of the GNU
C General Public License.
C See LICENSE and gpl-3.0.txt for details.
C
SUBROUTINE JSYMGG( PSHUP, KTRUNC, KSTART, KLUNIT, KLATO, KLONO,
X PLAT, PLEG, PTRIGS, KMFAX, PZFA, KRET)
C
C---->
C**** JSYMGG
C
C PURPOSE
C _______
C
C This routine converts spectral input fields to standard
C lat/long grid fields.
C
C INTERFACE
C _________
C
C CALL JSYMGG( PSHUP, KTRUNC, KSTART, KLUNIT, KLATO, KLONO,
C X PLAT, PLEG, PTRIGS, KMFAX, PZFA, KRET)
C
C Input parameters
C ________________
C
C PSHUP - Spherical harmonics field, unpacked
C KTRUNC - Truncation number of spherical harmonics field
C KSTART - Number of start latitude row (northernmost) for output
C field (must be positive - see comments below)
C KLUNIT - stream number of the legendre function file
C KLATO - Number of latitude rows in output field
C KLONO - Number of longitude points in output field
C PLAT - Array of gaussian latitudes
C PLEG - Array used to hold legendre functions
C PTRIGS - Initialized array of trig.functions (setup by JJSET99)
C KMFAX - Initialized array of prime factors (setup by JJSET99)
C
C Output parameters
C ________________
C
C PZFA - Output grid point field; contains upto 32 each of
C North and South latitude rows symmetrically.
C KRET - Return status code
C 0 = OK
C
C Common block usage
C __________________
C
C JDCNDBG
C
C Method
C ______
C
C None.
C
C Externals
C _________
C
C JREADGG - Reads the legendre functions for a latitude
C FFT99 - Carries out FFT
C INTLOG - Output log message
C INTLOGR - Output log message (with real value)
C NMAKGG - Make interpolation coefficients one latitude at a time
C
C Reference
C _________
C
C E.C.M.W.F. Research Department technical memorandum no. 56
C "The forecast and analysis post-processing package"
C May 1982. J.Haseler.
C
C Comments
C ________
C
C This is a redesign, based on SPECGP.F
C
C It handles transformation to a gaussian grid.
C The generated grid is symmetrical about the equator, so
C KSTART must be positive.
C
C It is not for U and V fields (no correction is applied at the
C poles).
C
C
C AUTHOR
C ______
C
C J.D.Chambers *ECMWF* Jan 1994
C
C MODIFICATIONS
C _____________
C
C J.D.Chambers ECMWF Feb 1997
C Allow for 64-bit pointers
C
C----<
C _______________________________________________________
C
IMPLICIT NONE
#include "jparams.h"
#include "parim.h"
#include "nifld.common"
C
C Subroutine arguments
COMPLEX PSHUP
DIMENSION PSHUP(*)
INTEGER KTRUNC
INTEGER KSTART
INTEGER KLUNIT, KLATO, KLONO, KMFAX, KRET
REAL PLAT, PLEG, PTRIGS, PZFA
DIMENSION PZFA(JPLONO + 2, 64)
DIMENSION KMFAX(*), PLAT(*), PLEG(*), PTRIGS(*)
C
C Parameters
INTEGER JPROUTINE
PARAMETER ( JPROUTINE = 31200 )
C
C Local variables
INTEGER ILIM, IMLIM, ILN
INTEGER ITAL, ITALA, ITALS, IMN, IMP
INTEGER INORTH, ISOUTH
INTEGER JM, J242, JNEXTLAT, JF
INTEGER NERR
INTEGER*8 IOFF
INTEGER*8 JDCLOOP
C
#ifdef POINTER_64
INTEGER*8 IWORK
#endif
REAL WORK
DIMENSION WORK(1)
POINTER ( IWORK, WORK )
COMPLEX ZDUM(JPTRNC + 1)
COMPLEX ZSUMS(JPTRNC + 1), ZSUMA(JPTRNC + 1)
COMPLEX*16 CHOLD
INTEGER*8 LOOP
C
INTEGER ISIZE
DATA ISIZE/0/
SAVE ISIZE, IWORK
C
C _______________________________________________________
C
C* Section 1. Initialization.
C _______________________________________________________
C
100 CONTINUE
C
C First time through, dynamically allocate memory for workspace
C
IF( ISIZE.EQ.0 ) THEN
ISIZE = 2*JPFFT*64
CALL JMEMHAN( 9, IWORK, ISIZE, 1, KRET)
IF( KRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'JSYMGG: memory allocation error.',IWORK)
KRET = JPROUTINE + 1
GOTO 990
ENDIF
ENDIF
C
IF ( NDBG .GT. 1) THEN
CALL INTLOG(JP_DEBUG,
X 'JSYMGG: Spherical harmonic coeffs(first 20):',JPQUIET)
DO 101 NDBGLP = 1, 20
CALL INTLOGR(JP_DEBUG,' ',PSHUP( NDBGLP ))
101 CONTINUE
CALL INTLOG(JP_DEBUG,'JSYMGG: Input parameters:',JPQUIET)
CALL INTLOG(JP_DEBUG,
X 'JSYMGG: Spherical harmonic truncation = ', KTRUNC)
CALL INTLOG(JP_DEBUG,
X 'JSYMGG: Start latitude(northernmost) = ', KSTART)
CALL INTLOG(JP_DEBUG,
X 'JSYMGG: Stream number of leg. file = ', KLUNIT)
CALL INTLOG(JP_DEBUG,
X 'JSYMGG: Number of lat. rows in output = ', KLATO)
CALL INTLOG(JP_DEBUG,
X 'JSYMGG: Number of long. pts per row = ', KLONO)
CALL INTLOG(JP_DEBUG,
X 'JSYMGG: Trig.functions (setup by JJSET99):',JPQUIET)
DO 102 NDBGLP = 1, 10
CALL INTLOGR(JP_DEBUG,' ',PTRIGS( NDBGLP ))
102 CONTINUE
CALL INTLOG(JP_DEBUG,
X 'JSYMGG: Prime factors (setup by JJSET99):',JPQUIET)
DO 103 NDBGLP = 1, 10
CALL INTLOG(JP_DEBUG,' ',KMFAX( NDBGLP ))
103 CONTINUE
ENDIF
C
ILIM = KTRUNC + 1
IMLIM = KTRUNC + 1
INORTH = -1
ILN = KLONO + 2
C
C _______________________________________________________
C
C* Section 2. Main loop through latitude rows to
C* calculate fourier coefficients
C _______________________________________________________
C
200 CONTINUE
C
C For each latitude, the north and corresponding south latitude row
C are calculated at the same time from the same legendre functions.
C
DO 280 JNEXTLAT = KSTART, KSTART+KLATO-1
C
IF ( NDBG .GT. 1)
X CALL INTLOG(JP_DEBUG,'JSYMGG: Next latitude = ', JNEXTLAT)
C
C If required, generate the coefficients 'on the fly'
C
IF( LON_FLY ) THEN
CALL NMAKGG( KTRUNC, JNEXTLAT, PLAT, 1, PLEG, NERR)
IOFF = 0
ELSE IF( LFILEIO ) THEN
CALL JREADGG( KLUNIT, KTRUNC, JNEXTLAT, PLEG, NERR)
IF ( NERR .NE. 0 ) THEN
CALL INTLOG(JP_ERROR,'JSYMGG: JREADGG error',NERR)
KRET = JPROUTINE + 2
GOTO 990
ENDIF
ELSE
IOFF = (JNEXTLAT-1)*(KTRUNC+1)*(KTRUNC+4)/2
ENDIF
C
C Clear unused slots in array.
C and one for the corresponding south latitude.
C
INORTH = INORTH + 2
ISOUTH = INORTH + 1
DO 241 JF = 2*IMLIM + 1, ILN
PZFA(JF, INORTH) = 0.0
PZFA(JF, ISOUTH) = 0.0
241 CONTINUE
C
C Now fill slots which are used
C
IMN = 0
IMP = 0
C
DO 244 JM = 1, IMLIM
ITAL = ILIM - JM + 1
DO 242 J242 = 1, ITAL
#ifndef __uxp__
IF( LFILEIO ) THEN
ZDUM(J242) = PLEG(IMP + J242)*PSHUP(IMN + J242)
ELSE
JDCLOOP = IOFF + IMP + J242
ZDUM(J242) = PLEG(JDCLOOP)*PSHUP(IMN + J242)
ENDIF
#else
JDCLOOP = IOFF + IMP + J242
ZDUM(J242) = PLEG(JDCLOOP)*PSHUP(IMN + J242)
#endif
242 CONTINUE
IMP = IMP + ITAL + 1
IMN = IMN + ITAL
ITALS = (ITAL + 1)/2
ITALA = ITAL/2
#ifndef CRAY
CHOLD = (0.0D0, 0.0D0)
#else
CHOLD = (0.0, 0.0)
#endif
DO LOOP = 1, 2*ITALS, 2
CHOLD = CHOLD + ZDUM(LOOP)
ENDDO
ZSUMS(JM) = CHOLD
#ifndef CRAY
CHOLD = (0.0D0, 0.0D0)
#else
CHOLD = (0.0, 0.0)
#endif
DO LOOP = 2, 2*ITALA, 2
CHOLD = CHOLD + ZDUM(LOOP)
ENDDO
ZSUMA(JM) = CHOLD
244 CONTINUE
C
C For the southern hemisphere row, the legendre functions are
C the complex conjugates of the corresponding northern row -
C hence the juggling with the signs in the next loop.
C
C Note that PZFA is REAL, but the coefficients being calculated
C are COMPLEX. There are pairs of values for each coefficient
C (real and imaginary parts) and pairs of values for each
C latitude (north and south).
C
DO 246 JM = 1, IMLIM
PZFA(2*JM -1, INORTH) = REAL(ZSUMS(JM)) + REAL(ZSUMA(JM))
PZFA(2*JM , INORTH) = AIMAG(ZSUMS(JM)) + AIMAG(ZSUMA(JM))
PZFA(2*JM -1, ISOUTH) = REAL(ZSUMS(JM)) - REAL(ZSUMA(JM))
PZFA(2*JM , ISOUTH) = AIMAG(ZSUMS(JM)) - AIMAG(ZSUMA(JM))
246 CONTINUE
C
C* End of main loop through latitude rows.
C
280 CONTINUE
C
C _______________________________________________________
C
C* Section 3. Fast fourier transform
C _______________________________________________________
C
300 CONTINUE
C
IF ( NDBG .GT. 1) CALL INTLOG(JP_DEBUG,
X 'JSYMGG: FFT, no.of rows (N and S) = ',ISOUTH)
C
CALL FFT99(PZFA,WORK,PTRIGS,KMFAX,1,J2NFFT,KLONO,ISOUTH,1)
C
IF ( NDBG .GT. 1) THEN
CALL INTLOG(JP_DEBUG,
X 'JSYMGG: Values calculated by FFT:',JPQUIET)
DO 301 NDBGLP = 1, 20
CALL INTLOGR(JP_DEBUG,' ',PZFA( 1, NDBGLP ))
CALL INTLOGR(JP_DEBUG,' ',PZFA( 2, NDBGLP ))
301 CONTINUE
ENDIF
C
C _______________________________________________________
C
C* Section 9. Return to calling routine. Format statements
C _______________________________________________________
C
C
900 CONTINUE
C
KRET = 0
C
990 CONTINUE
RETURN
END
|
