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
|
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 JSTRLL( PSHUP, KTRUNC, PNLAT, PBUILD, KLONO,
X PLEG, PTRIGS, KMFAX, PZFA, KRET)
C
C---->
C**** JSTRLL
C
C PURPOSE
C _______
C
C This routine converts stretched spectral input fields to
C standard lat/long grid fields.
C
C INTERFACE
C _________
C
C CALL JSTRLL( PSHUP, KTRUNC, PNLAT, PBUILD, KLONO,
C X 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 PNLAT - Start latitude (northernmost) for output field
C PBUILD - Grid interval used to build the legendre coefficients
C KLONO - Number of longitude points along the line of latitude
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 one latitude row.
C KRET - Return status code, 0 = OK
C
C Common block usage
C __________________
C
C JDCNDBG
C
C Method
C ______
C
C None.
C
C Externals
C _________
C
C FFT99 - Carries out FFT
C INTLOG - Output log message
C INTLOGR - Output log message (with real value)
C NMAKLL - 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 It handles transformation to a regular lat/long grid.
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 July 1999
C
C MODIFICATIONS
C _____________
C
C None.
C
C----<
C _______________________________________________________
C
C
IMPLICIT NONE
#include "jparams.h"
#include "parim.h"
#include "nifld.common"
C
C Parameters
C
INTEGER JPROUTINE
PARAMETER ( JPROUTINE = 31300 )
INTEGER JPN, JPS
PARAMETER ( JPN = 1 )
PARAMETER ( JPS = 2 )
C
C Subroutine arguments
C
COMPLEX PSHUP(*)
INTEGER KTRUNC
INTEGER KLONO, KRET
INTEGER KMFAX(*)
REAL PNLAT, PBUILD
REAL PZFA(JPLONO+2), PLEG(*), PTRIGS(*)
C
C Local variables
C
INTEGER ILIM, IMLIM, ILN
INTEGER ITAL, ITALA, ITALS, IMN, IMP
INTEGER INORTH
INTEGER JM, LOOP, JF
INTEGER NERR
C
#ifdef POINTER_64
INTEGER*8 IWORK
#endif
REAL ZNLAT, MFACTOR
REAL WORK
DIMENSION WORK(1)
POINTER ( IWORK, WORK )
COMPLEX ZDUM(JPTRNC + 1)
COMPLEX ZSUMS(JPTRNC + 1), ZSUMA(JPTRNC + 1)
COMPLEX*16 CHOLD
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
CALL JMEMHAN( 9, IWORK, ISIZE, 1, KRET)
IF( KRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'JSTRLL: memory allocation error.',IWORK)
KRET = JPROUTINE + 1
GOTO 990
ENDIF
ENDIF
C
IF( NDBG.GT.1 ) THEN
CALL INTLOG(JP_DEBUG,
X 'JSTRLL: Spherical harmonic coeffs(first 20):',JPQUIET)
DO NDBGLP = 1, 20
CALL INTLOGR(JP_DEBUG,' ',PSHUP( NDBGLP ))
ENDDO
CALL INTLOG(JP_DEBUG,'JSTRLL: Input parameters:',JPQUIET)
CALL INTLOG(JP_DEBUG,
X 'JSTRLL: Spherical harmonic truncation = ', KTRUNC)
CALL INTLOGR(JP_DEBUG,
X 'JSTRLL: Current latitude = ', PNLAT)
CALL INTLOGR(JP_DEBUG,
X 'JSTRLL: Grid lat. interval(leg. file) = ', PBUILD)
CALL INTLOG(JP_DEBUG,
X 'JSTRLL: Number of long. pts per row = ', KLONO)
CALL INTLOG(JP_DEBUG,
X 'JSTRLL: Trig.functions (setup by JJSET99):',JPQUIET)
DO NDBGLP = 1, 10
CALL INTLOGR(JP_DEBUG,' ',PTRIGS( NDBGLP ))
ENDDO
CALL INTLOG(JP_DEBUG,
X 'JSTRLL: Prime factors (setup by JJSET99):',JPQUIET)
DO NDBGLP = 1, 10
CALL INTLOG(JP_DEBUG,' ',KMFAX( NDBGLP ))
ENDDO
ENDIF
C
C If it is a stretched field, calculate the true latitude and
C the map factor.
C
IF( RISTRET.NE.0 ) THEN
CALL STRLAT(PNLAT, RISTRET, 1, ZNLAT, MFACTOR)
ELSE
ZNLAT = PNLAT
ENDIF
C
ILIM = KTRUNC + 1
IMLIM = KTRUNC + 1
INORTH = -1
ILN = KLONO + 2
C
C _______________________________________________________
C
C* Section 2. Calculate fourier coefficients
C _______________________________________________________
C
200 CONTINUE
C
C Generate the legendre coefficients 'on the fly'
C
IF( ZNLAT.GE.0 ) THEN
CALL NMAKLL( KTRUNC, PBUILD, ZNLAT, 1, PLEG, NERR)
ELSE
CALL NMAKLL( KTRUNC, PBUILD, -ZNLAT, 1, PLEG, NERR)
ENDIF
C
C Clear array.
C
DO JF = 1, JPLONO+2
PZFA(JF) = 0.0
ENDDO
C
C Fill slots which are used
C
IMN = 0
IMP = 0
C
DO JM = 1, IMLIM
ITAL = ILIM - JM + 1
DO LOOP = 1, ITAL
ZDUM(LOOP) = PLEG(IMP + LOOP)*PSHUP(IMN + LOOP)
ENDDO
C
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
C
ENDDO
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
IF( ZNLAT.GE.0 ) THEN
DO JM = 1, IMLIM
PZFA(2*JM -1) = REAL(ZSUMS(JM)) + REAL(ZSUMA(JM))
PZFA(2*JM ) = AIMAG(ZSUMS(JM)) + AIMAG(ZSUMA(JM))
ENDDO
ELSE
DO JM = 1, IMLIM
PZFA(2*JM -1) = REAL(ZSUMS(JM)) - REAL(ZSUMA(JM))
PZFA(2*JM ) = AIMAG(ZSUMS(JM)) - AIMAG(ZSUMA(JM))
ENDDO
ENDIF
C
C _______________________________________________________
C
C* Section 3. Fast fourier transform
C _______________________________________________________
C
300 CONTINUE
C
CALL FFT99(PZFA,WORK,PTRIGS,KMFAX,1,J2NFFT,KLONO,1,1)
C
C Apply map factor (squared) to vorticity or divergence
C
IF( (NIPARAM.EQ.138).OR.(NIPARAM.EQ.155) ) THEN
MFACTOR = MFACTOR*MFACTOR
DO LOOP = 1, KLONO
PZFA(LOOP) = MFACTOR*PZFA(LOOP)
ENDDO
ENDIF
C
IF( NDBG.GT.1 ) THEN
CALL INTLOG(JP_DEBUG,
X 'JSTRLL: Values calculated by FFT:',JPQUIET)
DO NDBGLP = 1, 20
CALL INTLOGR(JP_DEBUG,' ', PZFA( 1 + (NDBGLP-1)*2 ))
CALL INTLOGR(JP_DEBUG,' ', PZFA( NDBGLP*2 ))
ENDDO
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
|
