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
|
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
INTEGER FUNCTION SPROTAT(INGRIB,INLEN,SLON,SLAT,OUTGRIB,OUTLEN)
C
C---->
C**** SPROTAT
C
C Purpose
C -------
C
C Rotates a spherical harmonic field.
C
C
C Interface
C ---------
C
C IRET = SPROTAT(INGRIB,INLEN,SLON,SLAT,OUTGRIB,OUTLEN)
C
C Input
C -----
C
C INGRIB - Input GRIB field of spherical harmonics.
C INLEN - Length in bytes of the input GRIB.
C SLON - Longitude rotation required (degrees)
C SLAT - Latitude rotation required (degrees)
C
C
C Output
C ------
C
C OUTGRIB - Output GRIB field of rotated spherical harmonics.
C OUTLEN - Length in bytes of the output GRIB.
C
C Function returns 0 if all OK, otherwise the rotation failed.
C
C
C Method
C ------
C
C Externals
C ---------
C
C GRIBEX - Decode/encode a GRIB product.
C RPHI - Rotates spectral coefficients by longitude.
#ifdef __uxp__
C JACOBIF - Rotates spectral coefficients by latitude.
#else
C JACOBI - Rotates spectral coefficients by latitude.
#endif
C INTLOG - Logs messages.
C
C Author
C ------
C
C J.D.Chambers ECMWF October, 1995.
C
C
C---------------------------------------------------------------------
C----<
C
IMPLICIT NONE
#include "parim.h"
C
C Parameters
INTEGER JPROUTINE
PARAMETER (JPROUTINE = 28300)
INTEGER JPLEN, JPNM
PARAMETER ( JPNM = JPSTRUNC )
PARAMETER ( JPLEN = (JPNM+1)*(JPNM+2) )
C
C Function arguments.
INTEGER INGRIB, INLEN, OUTGRIB, OUTLEN
DIMENSION INGRIB(*), OUTGRIB(*)
REAL SLON, SLAT
C
C Local variables
INTEGER KSEC0, KSEC1, KSEC2, KSEC3, KSEC4, KWORD, KLENP, KLENG
DIMENSION KSEC0(2),KSEC1(512),KSEC2(512),KSEC3(2),KSEC4(512)
REAL PSEC4, PSEC2, PSEC3
DIMENSION PSEC2(512),PSEC3(2),PSEC4(JPLEN)
REAL*8 DLON, DLAT
REAL*8 WORK
DIMENSION WORK(2*(JPNM+1)*(JPNM+6))
INTEGER ITRUNC, IRET, NBYTES, LOOP
LOGICAL LOK
REAL*8 DATA
DIMENSION DATA(JPLEN)
C
#if (defined CRAY)
DATA NBYTES/8/
#else
DATA NBYTES/4/
#endif
C
C Externals
#ifdef __uxp__
LOGICAL JACOBIF
EXTERNAL JACOBIF
#else
LOGICAL JACOBI
EXTERNAL JACOBI
#endif
C
C _______________________________________________________
C
C* Section 1. Initialise.
C _______________________________________________________
C
100 CONTINUE
C
SPROTAT = 0
IRET = 0
#if (defined CRAY) || (defined REAL_8)
DLON = SLON
DLAT = -90.0 - SLAT
#else
DLON = DBLE(SLON)
DLAT = -90.0 - DBLE(SLAT)
#endif
C
C Decode the field.
KLENP = JPLEN
KLENG = (INLEN+NBYTES-1)/NBYTES
IRET = 1
CALL GRIBEX(KSEC0, KSEC1, KSEC2, PSEC2, KSEC3, PSEC3, KSEC4,
X PSEC4, KLENP, INGRIB, KLENG, KWORD, 'D' , IRET )
IF(IRET.NE.0)THEN
CALL INTLOG(JP_FATAL,'SPROTAT: Error decoding GRIB,IRET=', IRET)
SPROTAT = JPROUTINE + 1
GOTO 900
ENDIF
C
C Check that it really is a spherical harmonic field.
IF( (KSEC2(1).NE.50) .AND. (KSEC2(1).NE.60) ) THEN
CALL INTLOG(JP_FATAL,'SPROTAT: GRIB field not spectral',JPQUIET)
CALL INTLOG(JP_FATAL,'SPROTAT: GRIB KSEC2(1) = ', KSEC2(1))
SPROTAT = JPROUTINE + 2
GOTO 900
ENDIF
C
ITRUNC = KSEC2(2)
C
C _______________________________________________________
C
C* Section 2. Rotate the spectral coefficients.
C _______________________________________________________
C
200 CONTINUE
C
#if (defined CRAY) || (defined REAL_8)
C*********************************************************************
C
C (REALs are already double precision)
C
C Rotate the spectral field by longitude.
C Positive DLON => frame rotated from west to east.
C
CALL RPHI( PSEC4, ITRUNC, WORK, DLON)
C
C Rotate the spectral field by latitude.
C Negative DLAT => rotate counter-clockwise about new polar axis.
C
#if (defined __uxp__)
LOK = JACOBIF( PSEC4, ITRUNC, WORK, DLAT)
#else
LOK = JACOBI ( PSEC4, ITRUNC, WORK, DLAT)
#endif
#else
C*********************************************************************
C
C (REALs are single precision)
C
C Expand spectral coefficients to REAL*8
DO 210 LOOP = 1, KSEC4(1)
DATA(LOOP) = DBLE(PSEC4(LOOP))
210 CONTINUE
C
C Rotate the spectral field by longitude.
C Positive DLON => frame rotated from west to east.
C
CALL RPHI( DATA, ITRUNC, WORK, DLON)
C
C Rotate the spectral field by latitude.
C Negative DLAT => rotate counter-clockwise about new polar axis.
C
#if (defined __uxp__)
LOK = JACOBIF( DATA, ITRUNC, WORK, DLAT)
#else
LOK = JACOBI ( DATA, ITRUNC, WORK, DLAT)
#endif
IF(.NOT.LOK) THEN
CALL INTLOG(JP_FATAL,'SPROTAT: JACOBI failed.', JPQUIET)
SPROTAT = JPROUTINE + 3
GOTO 900
ENDIF
C
C Repack spectral coefficients to REAL*4.
DO 220 LOOP = 1, KSEC4(1)
PSEC4(LOOP) = SNGL(DATA(LOOP))
220 CONTINUE
C*********************************************************************
#endif
C
C _______________________________________________________
C
C* Section 3. Repack the spectral coefficients.
C _______________________________________________________
C
300 CONTINUE
C
C Put in details of southern pole of rotation
KSEC2(13) = NINT(SLAT*1000.0)
KSEC2(14) = NINT(SLON*1000.0)
C
C Put in indicator to show field has been rotated ..
KSEC2(1) = 60
C
C .. unless south pole of rotation is (0,-90).
IF( (KSEC2(14).EQ.0) .AND. (KSEC2(13).EQ.-90000) ) KSEC2(1) = 50
C
C Encode the field.
KLENP = KSEC4(1)
KLENG = (OUTLEN+NBYTES-1)/NBYTES
IRET = 1
CALL GRIBEX(KSEC0, KSEC1, KSEC2, PSEC2, KSEC3, PSEC3, KSEC4,
X PSEC4, KLENP, OUTGRIB, KLENG, KWORD, 'C' , IRET )
IF(IRET.NE.0)THEN
CALL INTLOG(JP_FATAL,'SPROTAT: Error encoding GRIB,IRET=', IRET)
SPROTAT = JPROUTINE + 4
GOTO 900
ENDIF
OUTLEN = KSEC0(1)
C
C _______________________________________________________
C
C* Section 9. Return.
C _______________________________________________________
C
900 CONTINUE
C
RETURN
END
|
