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
|
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 INTISL(OLDLSM, OLDFLD, NEWLSM, NEWFLD)
C
C---->
C**** INTISL
C
C Purpose
C -------
C
C Interpolate a field based on an old land-sea mask to a field
C based on a different land-sea mask.
C
C
C Interface
C ---------
C
C IRET = INTISL(OLDLSM, OLDFLD, NEWLSM, NEWFLD)
C
C Input
C -----
C
C OLDLSM - Old land-sea mask in GRIB format
C OLDFLD - Old field in GRIB format
C NEWLSM - New land-sea mask in GRIB format
C
C
C Output
C ------
C
C NEWFLD - New field in GRIB format
C
C Function returns:
C - the size in bytes of the new GRIB product if all is well
C - -1, otherwise.
C
C
C Method
C ------
C
C Unpack input GRIBS.
C Create new field and pack it into GRIB format.
C
C
C Externals
C ---------
C
C GRIBEX - Decode and encode GRIB products.
C NEWISL - Interpolate old field to new using both land-sea masks.
C JMALLOC - Dynamically allocate memory
C INTLOG - Log messages
C
C
C Author
C ------
C
C J.D.Chambers ECMWF August 2000
C
C----<
C
IMPLICIT NONE
C
C Function arguments
C
INTEGER OLDLSM(*), OLDFLD(*), NEWLSM(*), NEWFLD(*)
C
#include "parim.h"
#include "intisl.h"
#include "nifld.common"
#include "nofld.common"
C
C Parameters
C
INTEGER JPINTB, JPREALB
PARAMETER (JPINTB = 4)
#ifdef REAL_8
PARAMETER (JPREALB = 8)
#else
PARAMETER (JPREALB = 8)
#endif
C
C Local variables
C
INTEGER IXSEC4(512)
INTEGER ILSEC0(2),ILSEC1(1024),ILSEC2(1024),ILSEC3(2),ILSEC4(512)
INTEGER IFSEC0(2),IFSEC1(1024),IFSEC2(1024),IFSEC3(2),IFSEC4(512)
INTEGER IMSEC0(2),IMSEC1(1024),IMSEC2(1024),IMSEC3(2),IMSEC4(512)
INTEGER NFSEC0(2),NFSEC1(1024),NFSEC2(1024),NFSEC3(2),NFSEC4(512)
INTEGER KLSEC4, KFSEC4, KMSEC4, KPSEC4, NEWTYPE
#if (!defined __uxp__) && (!defined sgi)
#ifdef POINTER_64
INTEGER*8 IZLSEC4, IZFSEC4, IZMSEC4, IPFSEC4
#else
INTEGER IZLSEC4, IZFSEC4, IZMSEC4, IPFSEC4
#endif
#endif
C
REAL ZDUMMY
REAL ZLSEC2(1500), ZLSEC3(2), ZLSEC4(1)
REAL ZFSEC2(1500), ZFSEC3(2), ZFSEC4(1)
REAL ZMSEC2(1500), ZMSEC3(2), ZMSEC4(1)
REAL PFSEC2(1500), PFSEC3(2), PFSEC4(1)
C
POINTER( IZLSEC4, ZLSEC4 )
POINTER( IZFSEC4, ZFSEC4 )
POINTER( IZMSEC4, ZMSEC4 )
POINTER( IPFSEC4, PFSEC4 )
C
INTEGER IRET, ILENB, IPUNP, IWORD
C
C Externals
C
INTEGER NEWISL
#ifdef POINTER_64
INTEGER*8 JMALLOC
#else
INTEGER JMALLOC
#endif
EXTERNAL NEWISL, JMALLOC
C
DATA KLSEC4/0/, KFSEC4/0/, KMSEC4/0/, KPSEC4/0/
SAVE KLSEC4, KFSEC4, KMSEC4, KPSEC4
SAVE IZLSEC4, IZFSEC4, IZMSEC4, IPFSEC4
C
C -----------------------------------------------------------------|
C* Section 1. Initialise
C -----------------------------------------------------------------|
C
100 CONTINUE
C
INTISL = 0
C
C -----------------------------------------------------------------|
C* Section 2. Dynamically allocate memory for the old land-sea mask
C values and old field values.
C -----------------------------------------------------------------|
C
200 CONTINUE
C
C Old land-sea mask
C
IRET = 1
IPUNP = 1
ILENB = NIPOGRS
CALL GRIBEX (ILSEC0,ILSEC1,ILSEC2,ZLSEC2,ILSEC3,ZLSEC3,IXSEC4,
X ZDUMMY,IPUNP,OLDLSM,ILENB,IWORD,'J',IRET)
IF( IRET.NE.0 ) THEN
WRITE(*,*) 'INTISL: Old lsm gribex return code = ',IRET
INTISL = -1
RETURN
ENDIF
C
IF( KLSEC4.LT.(IXSEC4(1)*JPREALB) ) THEN
IF( KLSEC4.NE.0 ) THEN
CALL JFREE(IZLSEC4)
CALL JFREE(IZFSEC4)
ENDIF
KLSEC4 = IXSEC4(1)*JPREALB
KFSEC4 = KLSEC4
CALL INTLOG(JP_DEBUG,'INTISL: Allocate old memory = ',KLSEC4)
IZLSEC4 = JMALLOC(KLSEC4)
IZFSEC4 = JMALLOC(KFSEC4)
#ifdef hpR64
IZLSEC4 = IZLSEC4/(1024*1024*1024*4)
IZFSEC4 = IZFSEC4/(1024*1024*1024*4)
#endif
IF( (IZLSEC4.EQ.0).OR.(IZFSEC4.EQ.0) ) THEN
CALL INTLOG(JP_ERROR,'INTISL: JMALLOC fail oldies',JPQUIET)
INTISL = -1
GOTO 900
ENDIF
ENDIF
C
IRET = 1
IPUNP = IXSEC4(1)
ILENB = NIPOGRS
CALL GRIBEX (ILSEC0,ILSEC1,ILSEC2,ZLSEC2,ILSEC3,ZLSEC3,ILSEC4,
X ZLSEC4,IPUNP,OLDLSM,ILENB,IWORD,'D',IRET)
IF( IRET.NE.0 ) THEN
WRITE(*,*) 'INTISL: Old lsm gribex return code = ',IRET
INTISL = -1
RETURN
ENDIF
C
C Old field
C
IRET = 1
IPUNP = IXSEC4(1)
ILENB = NIPOGRS
IF( LIMISSV ) THEN
IFSEC3(2) = NINT(RMISSGV)
ZFSEC3(2) = RMISSGV
ENDIF
CALL GRIBEX (IFSEC0,IFSEC1,IFSEC2,ZFSEC2,IFSEC3,ZFSEC3,IFSEC4,
X ZFSEC4,IPUNP,OLDFLD,ILENB,IWORD,'D',IRET)
IF( ((IRET.EQ.-2).OR.(IRET.EQ.-4)).AND.
X (LIMISSV.EQV.(.FALSE.)) ) THEN
WRITE(*,*) 'INTISL: Old field has a bitmap'
WRITE(*,*) 'INTISL: A missing value must be specified'
INTISL = -1
RETURN
ELSE IF( IRET.GT.0 ) THEN
WRITE(*,*) 'INTISL: Old field gribex return code = ',IRET
INTISL = -1
RETURN
ENDIF
C
C -----------------------------------------------------------------|
C* Section 3. Dynamically allocate memory for the master land-sea
C mask values and new field values.
C -----------------------------------------------------------------|
C
300 CONTINUE
C
C New land-sea mask
C
IRET = 1
IPUNP = 1
ILENB = NIPNGRS
CALL GRIBEX (IMSEC0,IMSEC1,IMSEC2,ZMSEC2,IMSEC3,ZMSEC3,IXSEC4,
X ZDUMMY,IPUNP,NEWLSM,ILENB,IWORD,'J',IRET)
IF( IRET.NE.0 ) THEN
WRITE(*,*) 'INTISL: Master lsm gribex return code = ',IRET
INTISL = -1
RETURN
ENDIF
C
C Dynamically allocate memory for the master land-sea mask values
C and new field values
C
IF( KMSEC4.LT.(IXSEC4(1)*JPREALB) ) THEN
IF( KMSEC4.NE.0 ) THEN
CALL JFREE(IZMSEC4)
CALL JFREE(IPFSEC4)
ENDIF
KMSEC4 = IXSEC4(1)*JPREALB
KPSEC4 = KMSEC4
CALL INTLOG(JP_DEBUG,'INTISL: Allocate new memory = ',KMSEC4)
IZMSEC4 = JMALLOC(KMSEC4)
IPFSEC4 = JMALLOC(KPSEC4)
#ifdef hpR64
IZMSEC4 = IZMSEC4/(1024*1024*1024*4)
IPFSEC4 = IPFSEC4/(1024*1024*1024*4)
#endif
IF( (IZMSEC4.EQ.0).OR.(IPFSEC4.EQ.0) ) THEN
CALL INTLOG(JP_ERROR,'INTISL: JMALLOC fail newies',JPQUIET)
INTISL = -1
GOTO 900
ENDIF
ENDIF
C
IRET = 1
IPUNP = IXSEC4(1)
ILENB = NIPNGRS
CALL GRIBEX (IMSEC0,IMSEC1,IMSEC2,ZMSEC2,IMSEC3,ZMSEC3,IMSEC4,
X ZMSEC4,IPUNP,NEWLSM,ILENB,IWORD,'D',IRET)
IF( IRET.NE.0 ) THEN
WRITE(*,*) 'INTISL: Master lsm gribex return code = ',IRET
INTISL = -1
RETURN
ENDIF
C
C Requested field must match the new land-sea mask
C
NEWTYPE = NOREPR
IF( NOREPR.EQ.JPQUASI ) NEWTYPE = JPGAUSSIAN
IF( (NEWTYPE.NE.IMSEC2(1)).OR.
X ((NOREPR.NE.JPQUASI).AND.(IMSEC2(17).NE.0)).OR.
X ((NOREPR.EQ.JPQUASI).AND.(IMSEC2(17).NE.1)) ) THEN
CALL INTLOG(JP_ERROR,
X 'INTISL: Land-sea mask not suitable',JPQUIET)
INTISL = -1
RETURN
ENDIF
C
C -----------------------------------------------------------------|
C* Section 4. Create new field and put in GRIB format
C -----------------------------------------------------------------|
C
400 CONTINUE
C
IRET = NEWISL(ILSEC2, IMSEC2, ZLSEC4, ZFSEC4, ZMSEC4, PFSEC4)
IF( IRET.NE.0 ) THEN
WRITE(*,*) 'INTISL: New field creation failed'
INTISL = -1
RETURN
ENDIF
C
IRET = 1
IPUNP = IXSEC4(1)
ILENB = NIPNGRS
C
C There may be 'missing' values in the new field if the old field
C has 'missing' values.
C
IF( IFSEC1(5).EQ.192 ) THEN
IFSEC1(5) = 192
NFSEC3(1) = 0
NFSEC3(2) = INT(RMISSGV)
PFSEC3(1) = 0.0
PFSEC3(2) = RMISSGV
ENDIF
C
CALL GRIBEX (IFSEC0,IFSEC1,IMSEC2,ZMSEC2,NFSEC3,PFSEC3,IMSEC4,
X PFSEC4,IPUNP,NEWFLD,ILENB,IWORD,'C',IRET)
IF( IRET.NE.0 ) THEN
WRITE(*,*) 'INTISL: New field gribex return code = ',IRET
INTISL = -1
RETURN
ENDIF
C
INTISL = IWORD * JPINTB
C
C -----------------------------------------------------------------|
C* Section 9. Return
C -----------------------------------------------------------------|
C
900 CONTINUE
C
RETURN
END
|
