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
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
|
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 JWINDLL( PSHUP, KTRUNC, PSTART, PBUILD, PINTVL, KLUNIT,
X KLATO,KLONO,PLEG,PTRIGS,KMFAX,PZFA,KRET)
C
C---->
C**** JWINDLL
C
C Purpose
C -------
C This routine converts spectral input fields to standard
C lat/long grid fields for wind components U and V.
C
C
C Interface
C ---------
C CALL JWINDLL( PSHUP, KTRUNC, PSTART, PBUILD, PINTVL, KLUNIT,
C X KLATO,KLONO,PLEG,PTRIGS,KMFAX,PZFA,KRET)
C
C
C Input parameters
C ----------------
C PSHUP - Spherical harmonics field, unpacked
C KTRUNC - Truncation number of spherical harmonics field
C PSTART - Start latitude (northernmost) for output field
C (must be positive - see comments below)
C PBUILD - Grid interval used to build the legendre coefficients file
C PINTVL - Grid latitude interval in degrees
C KLUNIT - stream number of the legendre function file
C KLATO - Number of latitude points in output field
C KLONO - Number of longitude points in output field
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
C Output parameters
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
C Common block usage
C ------------------
C JDCNDBG
C
C
C Method
C ------
C None.
C
C
C Externals
C ---------
C JREADLL - Reads the legendre functions for a latitude
C JSPPOLE - Applies correction at North or South pole
C FFT99 - Carries out FFT
C INTLOG - Output log message
C INTLOGR - Output log message (with real value)
C GETENV - Pick up contents of an environment variable
C JUVPOLE - Use gaussian grid to calculate U and V at the poles
C
C
C Reference
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
C Comments
C --------
C This is a redesign, based on SPECGP.F
C
C It handles transformation to a regular lat/long grid.
C The generated grid is symmetrical about the equator, so
C PSTART must be positive.
C
C It is only for U and V fields (correction is applied at the
C poles and a scale factor is applied according to latitude).
C
C
C Author
C ------
C J.D.Chambers *ECMWF* Nov 1993
C
C
C Modifications
C -------------
C J.D.Chambers ECMWF Feb 1997
C Allow for 64-bit pointers
C
C J.D.Chambers ECMWF October 2002
C Add handling for U/V spectral (as opposed to Ucos(theta)/Vcos(theta))
C
C
C----<
C -----------------------------------------------------------------|
C
C
IMPLICIT NONE
#include "jparams.h"
#include "parim.h"
#include "nifld.common"
C
C Subroutine arguments
C
COMPLEX PSHUP
DIMENSION PSHUP(*)
INTEGER KTRUNC
REAL PSTART, PBUILD, PINTVL
INTEGER KLUNIT, KLATO, KLONO, KMFAX, KRET
REAL PLEG, PTRIGS, PZFA
DIMENSION PZFA(JPLONO + 2, 64)
DIMENSION KMFAX(*), PLEG(*), PTRIGS(*)
C
C Parameters
C
INTEGER JPROUTINE
PARAMETER ( JPROUTINE = 31400 )
C
C Local variables
C
REAL ZPIBY2, ZDEGR, ZLAT
INTEGER ILIM, IMLIM, ILN
INTEGER ITAL, ITALA, ITALS, IMN, IMP
INTEGER INORTH, ISOUTH
INTEGER JM, J245, JNEXTLAT
INTEGER NERR
INTEGER IOFF
INTEGER JDCLOOP
REAL ZF
REAL ZCOSI
C
DIMENSION ZF(2*JPFFT)
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 LOOP
C
LOGICAL LXFIRST, LPLAINU
INTEGER ISIZE, IBLANK
CHARACTER*10 PLAINUV
C
DATA LXFIRST/.TRUE./, LPLAINU/.FALSE./
DATA ISIZE/0/
SAVE ISIZE, IWORK, LXFIRST, LPLAINU
C
C Statement function
C
LOGICAL GPOLE
REAL ANGLE
GPOLE(ANGLE) = ( ABS(90.0 - ABS(ANGLE) ) .LT. 1.0E-3 )
C = .TRUE. if LAT is 90.0 or -90.0
C
C -----------------------------------------------------------------|
C* Section 1. Initialization.
C -----------------------------------------------------------------|
C
100 CONTINUE
C
C Check environment variable to see if wind components are plain
C U and V, or U*cos(theta) and V*cos(theta) (the default).
C
IF( LXFIRST ) THEN
C
LXFIRST = .FALSE.
CALL GETENV('INTERP_PLAIN_UV', PLAINUV)
IBLANK = INDEX(PLAINUV, ' ')
C
IF( (IBLANK.EQ.3).AND.(PLAINUV(1:2).EQ.'ON') ) LPLAINU =.TRUE.
IF( (IBLANK.EQ.3).AND.(PLAINUV(1:2).EQ.'on') ) LPLAINU =.TRUE.
IF( (IBLANK.EQ.4).AND.(PLAINUV(1:3).EQ.'YES') ) LPLAINU =.TRUE.
IF( (IBLANK.EQ.4).AND.(PLAINUV(1:3).EQ.'yes') ) LPLAINU =.TRUE.
C
IF( (IBLANK.EQ.3).AND.(PLAINUV(1:2).EQ.'NO') ) LPLAINU =.FALSE.
IF( (IBLANK.EQ.3).AND.(PLAINUV(1:2).EQ.'no') ) LPLAINU =.FALSE.
IF( (IBLANK.EQ.4).AND.(PLAINUV(1:3).EQ.'OFF') ) LPLAINU =.FALSE.
IF( (IBLANK.EQ.4).AND.(PLAINUV(1:3).EQ.'off') ) LPLAINU =.FALSE.
C
C First time through, dynamically allocate memory for workspace
C
ISIZE = 2*JPFFT*64
CALL JMEMHAN( 9, IWORK, ISIZE, 1, KRET)
IF( KRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'JWINDLL: memory allocation error',IWORK)
KRET = JPROUTINE + 1
GOTO 990
ENDIF
ENDIF
C
IF( NDBG.GT.1 ) THEN
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Spherical harmonic coeffs(first 20):',JPQUIET)
DO NDBGLP = 1, 20
CALL INTLOGR(JP_DEBUG,' ',PSHUP( NDBGLP ))
ENDDO
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Input parameters:',JPQUIET)
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Spherical harmonic truncation = ', KTRUNC)
CALL INTLOGR(JP_DEBUG,
X 'JWINDLL: Start latitude(northernmost) = ', PSTART)
CALL INTLOGR(JP_DEBUG,
X 'JWINDLL: Grid lat. interval(leg. file) = ', PBUILD)
CALL INTLOGR(JP_DEBUG,
X 'JWINDLL: Grid lat. interval(degrees) = ', PINTVL)
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Stream number of leg. file = ', KLUNIT)
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Number of lat. rows in output = ', KLATO)
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Number of long. pts per row = ', KLONO)
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Trig.functions (setup by JJSET99):',JPQUIET)
DO NDBGLP = 1, 10
CALL INTLOGR(JP_DEBUG,' ',PTRIGS( NDBGLP ))
ENDDO
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Prime factors (setup by JJSET99):',JPQUIET)
DO NDBGLP = 1, 10
CALL INTLOG(JP_DEBUG,' ',KMFAX( NDBGLP ))
ENDDO
ENDIF
C
ZPIBY2 = PPI / 2.0
ZDEGR = PPI / 180.0
ILIM = KTRUNC + 1
IMLIM = KTRUNC + 1
INORTH = -1
ILN = KLONO + 2
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 = 1, KLATO
ZLAT = ( PSTART - (PINTVL * REAL(JNEXTLAT - 1)) )
C
IF( NDBG.GT.1 )
X CALL INTLOGR(JP_DEBUG, 'JWINDLL: Next latitude = ', ZLAT)
C
C If required, generate the coefficients 'on the fly'
C
IF( LON_FLY ) THEN
CALL NMAKLL( KTRUNC, PBUILD, ZLAT, 1, PLEG, NERR)
IOFF = 0
ELSE IF( LFILEIO ) THEN
CALL JREADLL( KLUNIT, KTRUNC, PBUILD, ZLAT, PLEG, NERR)
IF ( NERR .NE. 0 ) THEN
CALL INTLOG(JP_ERROR,'JWINDLL: JREADLL error',NERR)
KRET = JPROUTINE + 2
GOTO 990
ENDIF
ELSE
IOFF = NINT( (90.0 - ZLAT)/PBUILD )*(KTRUNC+1)*(KTRUNC+4)/2
ENDIF
C
C Clear unused slots in array.
C Note there are two slots in the array - one for north latitude
C and one for the corresponding south latitude.
INORTH = INORTH + 2
ISOUTH = INORTH + 1
DO LOOP = 1, ILN
ZF( LOOP) = 0.0
PZFA( LOOP, INORTH) = 0.0
PZFA( LOOP, ISOUTH) = 0.0
ENDDO
C
C For spectral Ucos(theta), Vcos(theta), treat the poles as a
C special case.
C
IF( GPOLE(ZLAT).AND.(.NOT.LPLAINU) ) THEN
C
IF( NDBG.GT.1 ) CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Pole special for spectral Ucos(t), Vcos(t)',
X JPQUIET)
C
CALL JSPPOLE( PSHUP, 1, KTRUNC, .TRUE., ZF)
DO LOOP = 1, ILN
PZFA( LOOP, INORTH) = ZF( LOOP)
ENDDO
C
CALL JSPPOLE( PSHUP, 0, KTRUNC, .TRUE., ZF)
DO LOOP = 1, ILN
PZFA( LOOP, ISOUTH) = ZF( LOOP)
ENDDO
C
ELSE
C
C Otherwise not at pole, or spectral U and V
C
C Fill slots which are used
IMN = 0
IMP = 0
C
DO 247 JM = 1, IMLIM
ITAL = ILIM - JM + 1
DO 245 J245 = 1, ITAL
#ifndef __uxp__
IF( LFILEIO ) THEN
ZDUM(J245) = PLEG(IMP + J245)*PSHUP(IMN + J245)
ELSE
JDCLOOP = IOFF + IMP + J245
ZDUM(J245) = PLEG(JDCLOOP)*PSHUP(IMN + J245)
ENDIF
#else
JDCLOOP = IOFF + IMP + J245
ZDUM(J245) = PLEG(JDCLOOP)*PSHUP(IMN + J245)
#endif
245 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
247 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 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))
ENDDO
ENDIF
C
C* End of main loop through latitude rows.
C
280 CONTINUE
C
C -----------------------------------------------------------------|
C* Section 3. Fast fourier transform
C -----------------------------------------------------------------|
C
300 CONTINUE
C
IF( NDBG.GT.1 ) CALL INTLOG(JP_DEBUG,
X 'JWINDLL: 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 'JWINDLL: Values calculated by FFT:',JPQUIET)
DO NDBGLP = 1, 20
CALL INTLOGR(JP_DEBUG,' ',PZFA( 1, NDBGLP ))
CALL INTLOGR(JP_DEBUG,' ',PZFA( 2, NDBGLP ))
ENDDO
ENDIF
C
C For spectral U, V generate the poles from a reduced gaussian line
C of latitude.
C
IF( GPOLE(PSTART).AND.(LPLAINU) ) THEN
C
IF( NDBG.GT.1 ) CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Poles generated from gaussian for spectral U, V',
X JPQUIET)
C
CALL JUVPOLE(PSHUP, KTRUNC, PZFA, KLONO, KRET)
IF( KRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'JWINDLL: pole wind create error',IWORK)
KRET = JPROUTINE + 1
GOTO 990
ENDIF
C
ENDIF
C
C -----------------------------------------------------------------|
C* Section 4. For spectral Ucos(theta), Vcos(theta), apply scale
C factor to all latitudes except poles.
C -----------------------------------------------------------------|
C
400 CONTINUE
C
IF( NDBG.GT.0 ) THEN
IF( LPLAINU ) THEN
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Do not apply scale to latitudes',JPQUIET)
ELSE
CALL INTLOG(JP_DEBUG,
X 'JWINDLL: Apply scale to latitudes',JPQUIET)
ENDIF
ENDIF
C
IF( .NOT.LPLAINU ) THEN
INORTH = -1
DO JNEXTLAT = 1, KLATO
ZLAT = ( PSTART - (PINTVL * REAL(JNEXTLAT - 1)) )
INORTH = INORTH + 2
ISOUTH = INORTH + 1
IF( .NOT.GPOLE(ZLAT) ) THEN
ZCOSI = 1.0 / COS( ZLAT * ZDEGR )
CALL SSCAL( KLONO, ZCOSI, PZFA( 2, INORTH), 1)
CALL SSCAL( KLONO, ZCOSI, PZFA( 2, ISOUTH), 1)
ENDIF
ENDDO
C
ENDIF
C -----------------------------------------------------------------|
C* Section 9. Return to calling routine. Format statements
C -----------------------------------------------------------------|
C
900 CONTINUE
C
KRET = 0
C
990 CONTINUE
RETURN
END
|
