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
|
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 JSLLGP( PSHUP, KTRUNC, PNORTH, PSOUTH, PWEST, PEAST,
X PLATINC, PLONINC, POUTF, KRET)
C
C---->
C**** JSLLGP
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 JSLLGP( PSHUP, KTRUNC, PNORTH, PSOUTH, PWEST, PEAST,
C X PLATINC, PLONINC, POUTF, KRET)
C
C Input parameters
C ________________
C
C PSHUP - Spherical harmonics field, unpacked
C KTRUNC - Truncation number of spherical harmonics field
C PNORTH - Northernmost latitude for output field (degrees)
C PSOUTH - Southernmost latitude for output field (degrees)
C PWEST - Westernmost longitude for output field (degrees)
C PEAST - Easternmost longitude for output field (degrees)
C PLATINC - Grid interval between latitudes in degrees
C PLONINC - Grid interval between longitudes in degrees
C INTLOG - Output log message
C INTLOGR - Output log message (with real value)
C
C Output parameters
C ________________
C
C POUTF - Output grid point field
C KRET - Return status code
C 0 = OK
C
C Common block usage
C __________________
C
C JDCNDBG
C
C Method
C ______
C
C Generates one latitude at a time.
C
C Externals
C _________
C
C JJSET99 - Sets up FFT
C JSTRLL - Generate latitude line from spherical harmonics
C JSTRWLL - Generate latitude line from spherical harmonics (wind)
C JDEBUG - Checks environment variable to switch on/off debug
C INTLOG - Logs output messages
C INTLOGR - Logs output messages
C
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
C Comments
C ________
C
C None
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* Section 0. Definition of variables.
C _______________________________________________________
C
IMPLICIT NONE
#include "jparams.h"
#include "parim.h"
#include "nifld.common"
C
C Parameters
C
INTEGER JPROUTINE
PARAMETER ( JPROUTINE = 30000 )
C
C Subroutine arguments
C
COMPLEX PSHUP
DIMENSION PSHUP(*)
INTEGER KTRUNC
REAL PNORTH, PSOUTH, PWEST, PEAST, PLATINC, PLONINC
REAL POUTF
DIMENSION POUTF(*)
INTEGER KRET
C
C Local variables
C
REAL PSTART, PNORMX, PDIFF
REAL DEG2RAD, ZBUILD
INTEGER ILON, IOLON, ITOUP1, IJPWR, NEXT, NROWS, NUMLON
INTEGER IOFFN, LONGIT, N360, NSTART, INEXT, LOOPCL
C
INTEGER IMFAX
REAL ATRIGS
DIMENSION IMFAX(10),ATRIGS(JPLONO*2)
#ifdef POINTER_64
INTEGER*8 IALEG, IZFA
#endif
#ifdef CRAY
REAL ALEG
DIMENSION ALEG(JPK)
REAL ZFA
DIMENSION ZFA(JPLONO+2)
#else
REAL ALEG
DIMENSION ALEG(1)
POINTER ( IALEG, ALEG )
REAL ZFA
DIMENSION ZFA(1)
POINTER ( IZFA, ZFA )
#endif
C
INTEGER ISIZE, ISIZE2
DATA ISIZE/0/, ISIZE2/0/
SAVE ISIZE, ISIZE2, IZFA
C
SAVE IMFAX, ATRIGS, IALEG
C
C _______________________________________________________
C
C* Section 1. Initialization.
C _______________________________________________________
C
100 CONTINUE
C
CALL JDEBUG( )
C
C First time through, dynamically allocate memory for workspace
C
IF( ISIZE.EQ.0 ) THEN
ISIZE = (JPLONO+2) * 2
CALL JMEMHAN( 10, IZFA, ISIZE, 1, KRET)
IF( KRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'JSLLGP: memory allocation error.',IZFA)
KRET = JPROUTINE + 1
GOTO 900
ENDIF
ENDIF
C
IF( NDBG.GT.0 ) THEN
CALL INTLOG(JP_DEBUG,'JSLLGP: Input parameters:', JPQUIET)
CALL INTLOG(JP_DEBUG,
X 'JSLLGP: Spherical harmonic coeffs(first 20):', JPQUIET)
DO 101 NDBGLP = 1, 20
CALL INTLOGR(JP_DEBUG,' ', PSHUP( NDBGLP ))
101 CONTINUE
CALL INTLOG(JP_DEBUG,
X 'JSLLGP: Spherical harmonic truncation = ', KTRUNC)
CALL INTLOGR(JP_DEBUG,
X 'JSLLGP: Northern latitude for output = ', PNORTH)
CALL INTLOGR(JP_DEBUG,
X 'JSLLGP: Southern latitude for output = ', PSOUTH)
CALL INTLOGR(JP_DEBUG,
X 'JSLLGP: Western longitude for output = ', PWEST)
CALL INTLOGR(JP_DEBUG,
X 'JSLLGP: Eastern longitude for output = ', PEAST)
CALL INTLOGR(JP_DEBUG,
X 'JSLLGP: Latitude grid interval (deg) = ', PLATINC)
CALL INTLOGR(JP_DEBUG,
X 'JSLLGP: Longitude grid interval (deg) = ', PLONINC)
ENDIF
C
C Coefficients are to generated 'on the fly',
C
ISIZE2 = (KTRUNC+1)*(KTRUNC+2)
CALL JMEMHAN( 14, IALEG, ISIZE2, 1, KRET)
IF( KRET.NE.0 ) THEN
CALL INTLOG(JP_ERROR,'JSLLGP: memory allocate error.',IALEG)
KRET = JPROUTINE + 2
GOTO 900
ENDIF
C
DEG2RAD = PPI / 180.0
IOLON = NINT( 360.0 / PLATINC )
ILON = NINT( 360.0 / PLONINC )
ITOUP1 = KTRUNC
C
C For calculation purposes, ILON, the number of longitude points,
C has to be greater than 2*(output truncation) to ensure that the
C fourier transform is exact (see Reference, page 10).
C The following loop increases ILON*IJPWR until ILON=>2*ITOUP1.
C Later, when filling the output array POUTF, the longitude points
C have to be taken selectively to avoid the intermediate generated
C points, picking up values only at the required longitudes.
C
IJPWR = 1
110 CONTINUE
IF( ILON.GT.2*ITOUP1 ) GOTO 120
IJPWR = IJPWR*2
ILON = ILON*2
GOTO 110
120 CONTINUE
C
C Initialise FFT constants
C
IF( NDBG.GT.0 ) THEN
CALL INTLOG(JP_DEBUG,
X 'JSLLGP: Initialise FFT constants with ', ILON)
CALL INTLOG(JP_DEBUG,
X 'JSLLGP: longitude points.', JPQUIET)
ENDIF
C
CALL JJSET99( ATRIGS, IMFAX, ILON, KRET)
IF( KRET.NE.0 ) GOTO 900
C
C _______________________________________________________
C
C* Section 2. Work out geometry of the request
C _______________________________________________________
C
200 CONTINUE
C
PNORMX = PNORTH
PDIFF = PNORTH - PSOUTH
NROWS = NINT( PDIFF/PLATINC ) + 1
IF( ((NROWS-1)*PLATINC).GT.PDIFF ) NROWS = NROWS - 1
C
NUMLON = NINT( (PEAST - PWEST)/PLONINC ) + 1
ZBUILD = PLATINC
C
IF( NDBG.GT.0 ) THEN
CALL INTLOG(JP_DEBUG,'JSLLGP: Number of longitudes = ', NUMLON)
CALL INTLOGR(JP_DEBUG,
X 'JSLLGP: North latitude used for calc = ', PNORMX)
CALL INTLOG(JP_DEBUG,'JSLLGP: No.lats used in calc = ', NROWS)
ENDIF
C _______________________________________________________
C
C* Section 4. Process one row at a time.
C _______________________________________________________
C
400 CONTINUE
C
DO NEXT = 1, NROWS
PSTART = PNORMX - REAL(NEXT-1)*PLATINC
C
IF( NDBG.GT.0 )
X CALL INTLOGR(JP_DEBUG, 'JSLLGP: Next latitude = ', PSTART)
C
C Generate the next row
C
IF( LWIND ) THEN
CALL JSTRWLL(PSHUP, KTRUNC, PSTART, ZBUILD, ILON,
X ALEG, ATRIGS, IMFAX, ZFA, KRET)
ELSE
CALL JSTRLL(PSHUP, KTRUNC, PSTART, ZBUILD, ILON,
X ALEG, ATRIGS, IMFAX, ZFA, KRET)
ENDIF
IF( KRET.NE.0 ) GOTO 900
C
C Find start offset in storage array
C
IOFFN = NINT( (PNORTH - PSTART)/PLATINC ) + 1
LONGIT = NINT( (360.0 + PWEST)/PLONINC )
N360 = NINT( 360.0/PLONINC )
NSTART = MOD( LONGIT, N360)*IJPWR
C
DO LOOPCL = 1, NUMLON
INEXT = 2 + MOD( NSTART + (LOOPCL-1)*IJPWR , ILON)
POUTF( LOOPCL + (IOFFN-1)*NUMLON ) = ZFA( INEXT)
ENDDO
C
ENDDO
C
KRET = 0
C
C _______________________________________________________
C
C* Section 9. Return to calling routine. Format statements
C _______________________________________________________
C
900 CONTINUE
C
RETURN
END
|
