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
|
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 WV2DI32(KOLDNUM,NUMOLD,OLDLATS,OLDWEST,OLDEAST,
X KNEWNUM,NUMNEW,NEWLATS,NEWWAVE)
C
C---->
C*****WV2DI32*
C
C PURPOSE
C -------
C
C Determines which nearest-neighbour values of an input global wave
C 2D-spectra grid field to use for an output global wave 2D-spectra
C grid field.
C (32-bit REALs version)
C
C
C INTERFACE
C ---------
C
C IRET = WV2DI32(KOLDNUM,NUMOLD,OLDLATS,OLDWEST,OLDEAST,
C X KNEWNUM,NUMNEW,NEWLATS,NEWWAVE)
C
C Input arguments
C ---------------
C
C KOLDNUM - No. of meridians from North to South pole (input field)
C NUMOLD - Array giving number of points along each latitude
C (empty latitudes have entry 0)
C OLDLATS - input field latitudes
C OLDWEST - western longitude of the input field (degrees)
C OLDEAST - eastern longitude of the input field (degrees)
C
C KNEWNUM - No. of meridians from North to South pole (output field)
C NUMNEW - Array giving number of points along each latitude
C (empty latitudes have entry 0)
C NEWLATS - output field latitudes
C
C Output arguments
C ----------------
C
C NEWWAVE - Indices of points to use
C
C Function returns 0 if the interpolation was OK.
C
C
C METHOD
C ------
C
C The index of the nearest neighbouring grid point value of
C the input field is assigned to each point of the output
C grid.
C The index is zero if the output grid point is 'missing'.
C
C The input field can be regular or quasi-regular, and can be
C global or a subarea.
C The output field can be regular or quasi-regular.
C
C
C EXTERNALS
C ---------
C
C INTLOG - Log error message.
C
C
C REFERENCE
C ---------
C
C None
C
C
C Author.
C -------
C
C J.D.Chambers ECMWF October 1997
C
C
C----<
C
IMPLICIT NONE
C
#include "parim.h"
C
C Parameters
C
INTEGER JPROUTINE
PARAMETER ( JPROUTINE = 19410 )
INTEGER JPLLMAX
PARAMETER ( JPLLMAX = 721 )
C `--> allow upto 0.25 degree resolution
INTEGER JPNMOUT
PARAMETER ( JPNMOUT = 720 )
C `--> allow upto 0.25 degree resolution
C
C Function arguments
C
INTEGER KOLDNUM, NUMOLD, KNEWNUM, NUMNEW
DIMENSION NUMOLD(*), NUMNEW(*)
REAL*4 OLDWEST, OLDEAST, OLDLATS, NEWLATS
DIMENSION OLDLATS(*), NEWLATS(*)
INTEGER NEWWAVE
DIMENSION NEWWAVE(*)
C
C Local arguments
C
INTEGER K, NEWCOL, NEWROW, INCOL, LOOP, NEXT
INTEGER I_NW, I_SW, I_N, I_S, I_NE, I_SE
REAL*8 DELONGN, DELONGS, DIFF, RESOL, DINC
REAL*8 DIST_NW, DIST_SW, DIST_NE, DIST_SE, DIST_N, DIST_S
REAL*4 DISNW, DISNE, DISSW, DISSE
REAL*8 ZXIN, ZXIS, LON, LAT, OWEST, OEAST
DIMENSION LON(JPNMOUT*2)
LOGICAL LINGNS, LINGWE
INTEGER INDEXI
DIMENSION INDEXI(JPLLMAX)
C
C ---------------------------------------------------------------------
C* Section 1. Initalisation.
C ---------------------------------------------------------------------
C
100 CONTINUE
C
WV2DI32 = 0
CALL INTLOG(JP_DEBUG,
X 'WV2DI32: Wave interpolation requested.',JPQUIET)
C
C Check latitude/longitude grid specification
IF( KOLDNUM.GT.JPLLMAX ) THEN
CALL INTLOG(JP_ERROR,
X 'WV2DI32: Number of latitudes in input grid = ',KOLDNUM)
CALL INTLOG(JP_ERROR,
X 'WV2DI32: And is greater than allowed maximum = ',JPLLMAX)
WV2DI32 = JPROUTINE + 1
GOTO 900
ENDIF
C
C Set up INDEXI for latitude lines in input lat/lon array
C
INDEXI(1) = 0
DO LOOP = 2, KOLDNUM
INDEXI(LOOP) = INDEXI(LOOP-1) + NUMOLD(LOOP-1)
ENDDO
C
OEAST = OLDEAST
OWEST = OLDWEST
C
C ---------------------------------------------------------------------
C* Section 2. Interpolation.
C ---------------------------------------------------------------------
C
200 CONTINUE
C
NEXT = 1
C
DO 220 NEWROW = 1, KNEWNUM
C
C Find old latitude rows to north and south of new latitude and
C set up the distance between new row and the old rows to N and S.
C
DO LOOP = 1, KOLDNUM -1
IF( (NEWLATS(NEWROW).LE.OLDLATS(LOOP)) .AND.
X (NEWLATS(NEWROW).GT.OLDLATS(LOOP+1)) ) THEN
I_N = LOOP
I_S = MIN(LOOP+1,KOLDNUM)
GOTO 205
ENDIF
ENDDO
C
IF( NEWLATS(NEWROW).EQ.OLDLATS(KOLDNUM) ) THEN
I_N = KOLDNUM
I_S = KOLDNUM
ELSE
I_N = -1
I_S = -1
ENDIF
C
205 CONTINUE
IF( I_N.NE.I_S) THEN
DIFF = OLDLATS(I_N) - OLDLATS(I_S)
DIST_N = (OLDLATS(I_N) - NEWLATS(NEWROW)) / DIFF
ELSE
DIST_N = 1.0
ENDIF
DIST_S = 1.0 - DIST_N
C
C Check if the new interpolated row lies between 2 old rows which
C have data points
C
LINGNS = .FALSE.
IF( I_N.GT.0 ) THEN
IF( (NUMOLD(I_N).GT.0).AND.(NUMOLD(I_S).GT.0) ) THEN
C
C Yes, so set up the grid increments.
C
DINC = 360.0/DBLE(NUMOLD(I_N))
DELONGN = (DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))
X +DINC) / DBLE(NUMOLD(I_N))
DINC = 360.0/DBLE(NUMOLD(I_S))
DELONGS = (DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))
X +DINC) / DBLE(NUMOLD(I_S))
LINGNS = .TRUE.
ENDIF
C
C The equator is given special treatment so that the northern
C and southern hemispheres for the parameter 250 can (later)
C be stitched together.
C
C If the input field finishs at the equator and the output
C field has a line at the equator, use the input equator for
C interpolation.
C
IF( NEWLATS(I_N).EQ.0.0 .AND.
x ((NUMOLD(I_N).GT.0).AND.(NUMOLD(I_S).EQ.0)) ) THEN
C
C Yes, so set up the grid increments.
C
DINC = 360.0/DBLE(NUMOLD(I_N))
DELONGN = (DMOD((OEAST-OWEST+DBLE(360.0)), DBLE(360.0))
X +DINC) / DBLE(NUMOLD(I_N))
DELONGS = DELONGN
DIST_N = 0.0
DIST_S = 1.0
LINGNS = .TRUE.
ENDIF
ENDIF
C
C Setup longitudes for the current row
C
INCOL = NUMNEW(NEWROW)
IF( INCOL.NE.0 ) THEN
RESOL = 360.0/DBLE(INCOL)
DO LOOP = 1, INCOL
LON(LOOP) = DBLE(LOOP-1)*RESOL
ENDDO
C
DO 210 NEWCOL = 1, INCOL
C
LINGWE = (LON(NEWCOL).GT.OWEST).OR.(LON(NEWCOL).LT.OEAST)
C
IF( LINGNS .AND. LINGWE ) THEN
C
C If point to be interpolated is inside the input field ,
C work out distances from new point to four neighbours.
C
ZXIN = DMOD( (LON(NEWCOL) + (360.0-OWEST) + DBLE(720.0)),
X DBLE(360.0))
ZXIS = ZXIN
ZXIN = ZXIN/DELONGN + 1.0
ZXIS = ZXIS/DELONGS + 1.0
I_NW = INT(ZXIN)
I_NE = I_NW + 1
I_SW = INT(ZXIS)
I_SE = I_SW + 1
C
C Allow wrap-around
C
C
I_NW = MOD(I_NW,NUMOLD(I_N))
I_NE = MOD(I_NE,NUMOLD(I_N))
I_SW = MOD(I_SW,NUMOLD(I_S))
I_SE = MOD(I_SE,NUMOLD(I_S))
IF( I_NW.EQ.0 ) I_NW = NUMOLD(I_N)
IF( I_SW.EQ.0 ) I_SW = NUMOLD(I_S)
IF( I_NE.EQ.0 ) I_NE = NUMOLD(I_N)
IF( I_SE.EQ.0 ) I_SE = NUMOLD(I_S)
C
C Calculate distance from interpolated point to its neighbours
C
DIST_NW = ZXIN - REAL(I_NW)
DIST_NE = 1.0 - DIST_NW
DIST_SW = ZXIS - REAL(I_SW)
DIST_SE = 1.0 - DIST_SW
C
C Take nearest grid point value.
C
DISNW = DIST_NW*DIST_NW + DIST_N*DIST_N
DISNE = DIST_NE*DIST_NE + DIST_N*DIST_N
DISSW = DIST_SW*DIST_SW + DIST_S*DIST_S
DISSE = DIST_SE*DIST_SE + DIST_S*DIST_S
C
IF( (DISNW.LE.DISNE).AND.(DISNW.LE.DISSW).AND.
X (DISNW.LE.DISSE) ) THEN
C
C Use north-west neighbour
C
NEWWAVE(NEXT) = INDEXI(I_N)+I_NW
C
ELSE IF( (DISNE.LE.DISSW).AND.(DISNE.LE.DISSE) ) THEN
C
C Use north-east neighbour
C
NEWWAVE(NEXT) = INDEXI(I_N)+I_NE
C
ELSE IF( DISSW.LE.DISSE ) THEN
C
C Use south-west neighbour
C
NEWWAVE(NEXT) = INDEXI(I_S)+I_SW
ELSE
C
C Use south-east neighbour
C
NEWWAVE(NEXT) = INDEXI(I_S)+I_SE
ENDIF
C
C If the point to be interpolated is outside the input field,
C set it to 'land'.
C
ELSE
NEWWAVE(NEXT) = 0
ENDIF
C
NEXT = NEXT + 1
C
210 CONTINUE
C
ENDIF
C
220 CONTINUE
C
C ---------------------------------------------------------------------
C* Section 9. Closedown.
C ---------------------------------------------------------------------
C
900 CONTINUE
C
IF( WV2DI32.EQ.0 ) THEN
CALL INTLOG(JP_DEBUG, 'WV2DI32: successful.',JPQUIET)
ELSE
CALL INTLOG(JP_ERROR, 'WV2DI32: failed.',JPQUIET)
ENDIF
C
RETURN
END
|
