1C Copyright 1981-2016 ECMWF.
2C
3C This software is licensed under the terms of the Apache Licence
4C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5C
6C In applying this licence, ECMWF does not waive the privileges and immunities
7C granted to it by virtue of its status as an intergovernmental organisation
8C nor does it submit to any jurisdiction.
9C
10
11      INTEGER FUNCTION NEWISL(OLDGEO,NEWGEO,OLDLSM,OLDFLD,MASTER,NEWFLD)
12C
13C---->
14C**** NEWISL
15C
16C     Purpose
17C     -------
18C
19C     Interpolate a field based on an old land-sea mask to a field
20C     based on a different land-sea mask.
21C
22C
23C     Interface
24C     ---------
25C
26C     IRET = NEWISL(OLDGEO,NEWGEO,OLDLSM,OLDFLD,MASTER,NEWFLD)
27C
28C     Input
29C     -----
30C
31C     OLDGEO - GRIB section 2 describing grid of old field.
32C     NEWGEO - GRIB section 2 describing grid of new field.
33C     OLDLSM - Array of land-sea mask values for old field.
34C     OLDFLD - Array of values for old field.
35C     MASTER - Array of land-sea mask values for new field.
36C
37C
38C     Output
39C     ------
40C
41C     NEWFLD - Array of values for new field.
42C
43C     Function returns:
44C       - 0 if all is well
45C       - 1, otherwise.
46C
47C
48C     Method
49C     ------
50C
51C     Build up field offsets from input geometries.
52C     For each point of the new field, find in the old field the four
53C     nearest neighbours' positions, values and types.
54C     Calculate new value from the neighbours.
55C
56C
57C     Externals
58C     ---------
59C
60C     IGGLAT  - Compute gaussian lines of latitude.
61C     JNORSGG - Find gaussian latitudes to north and south of latitude.
62C     ISLPROC - Calculate value of new field point.
63C     INTLOG  - Log messages
64C
65C
66C     Author
67C     ------
68C
69C     J.D.Chambers     ECMWF     August 2000
70C
71C----<
72C
73      IMPLICIT NONE
74C
75C     Function arguments
76C
77      INTEGER OLDGEO(*),NEWGEO(*)
78      REAL OLDLSM(*),OLDFLD(*),MASTER(*),NEWFLD(*)
79C
80#include "intisl.h"
81#include "parim.h"
82C
83C     Parameters
84C
85      INTEGER JPMAXLT, JPGAUSS, JP1000
86      PARAMETER (JPMAXLT=721)
87      PARAMETER (JPGAUSS=4)
88      PARAMETER (JP1000=1000)
89C
90C     Local variables
91C
92      INTEGER TOTAL, NEXT, LOOP, NEWOFF(JPMAXLT), OLDOFF(JPMAXLT)
93      INTEGER LATIT, LONG, NEWTYPE, IRET, NPTS
94      INTEGER LAT(2), LON(4)
95      INTEGER PT(4), TYPE(4)
96      REAL OLAT(2), OLON(4)
97      REAL RLATOLD(JPMAXLT), RLATNEW(JPMAXLT), RLAT, RLON
98C
99C     Externals
100C
101      INTEGER IGGLAT, JNORSGG
102      REAL ISLPROC
103      EXTERNAL IGGLAT, JNORSGG, ISLPROC
104C
105C     -----------------------------------------------------------------|
106C*    Section 1.   Build working values using input geometries.
107C     -----------------------------------------------------------------|
108C
109  100 CONTINUE
110C
111      NEWISL = 0
112C
113C     -----------------------------------------------------------------|
114C*    Section 2.   Calculate number of points in new field and offset
115C                  to start of each latitude in the new grid
116C     -----------------------------------------------------------------|
117C
118  200 CONTINUE
119C
120      IF( NEWGEO(1).EQ.JPGAUSS ) THEN
121C
122C       New field is gaussian
123C
124        CALL INTLOG(JP_DEBUG,'NEWISL: New field is gaussian',JPQUIET)
125C
126        IF( NEWGEO(17).EQ.0 ) THEN
127          CALL INTLOG(JP_DEBUG,'NEWISL: New field is regular',JPQUIET)
128          TOTAL = NEWGEO(2)*NEWGEO(3)
129          NEWOFF(1) = 0
130          DO LOOP = 2, NEWGEO(NJ)
131            NEWOFF(LOOP) = NEWOFF(LOOP-1) + NEWGEO(2)
132          ENDDO
133        ELSE
134          CALL INTLOG(JP_DEBUG,'NEWISL: New field is reduced',JPQUIET)
135          NEWOFF(1) = 0
136          TOTAL = NEWGEO(NPOINTS)
137          DO LOOP = 2, NEWGEO(NJ)
138            NEWOFF(LOOP) = NEWOFF(LOOP-1) + NEWGEO(NPOINTS-2+LOOP)
139            TOTAL = TOTAL + NEWGEO(NPOINTS-1+LOOP)
140          ENDDO
141        ENDIF
142C
143C       Get the gaussian latitudes for the new field
144C
145        IRET = IGGLAT(NEWGEO(NGAUSS)*2,RLATNEW,0,-1)
146        IF( IRET.NE.0 ) THEN
147          WRITE(*,*) 'NEWISL: Problem call igglat for new grid'
148          NEWISL = 1
149          RETURN
150        ENDIF
151C
152      ELSE
153C
154C       New field is lat/long
155C
156        CALL INTLOG(JP_DEBUG,'NEWISL: New field is lat/long',JPQUIET)
157        TOTAL = NEWGEO(2)*NEWGEO(3)
158        DO LOOP = 1, NEWGEO(3)
159          NEWOFF(LOOP) = NEWGEO(2)*(LOOP-1)
160          RLATNEW(LOOP) = 90.0 - (REAL((LOOP-1)*NEWGEO(10))/JP1000)
161        ENDDO
162C
163      ENDIF
164C
165      CALL INTLOG(JP_DEBUG,'NEWISL: No. of pts in new field = ',TOTAL)
166C
167C     -----------------------------------------------------------------|
168C*    Section 3.   Get the gaussian latitudes for the old field and
169C                  setup the offsets to the start of each latitude.
170C     -----------------------------------------------------------------|
171C
172  300 CONTINUE
173C
174      OLDOFF(1) = 0
175      DO LOOP = 2, OLDGEO(NJ)
176        OLDOFF(LOOP) = OLDOFF(LOOP-1) + OLDGEO(NPOINTS-2+LOOP)
177      ENDDO
178C
179      IRET = IGGLAT(OLDGEO(NGAUSS)*2,RLATOLD,0,-1)
180      IF( IRET.NE.0 ) THEN
181        WRITE(*,*) 'NEWISL: Problem call igglat for old grid'
182        NEWISL = 1
183        RETURN
184      ENDIF
185C
186C     -----------------------------------------------------------------|
187C*    Section 4.   Work through the points in the new field.
188C     -----------------------------------------------------------------|
189C
190  400 CONTINUE
191C
192      DO NEXT = 1, TOTAL
193C
194C       Calculate lat/long
195C
196        DO LOOP = 1, NEWGEO(NJ)
197          IF( NEWOFF(LOOP).GE.NEXT ) THEN
198            LATIT = LOOP - 1
199            GOTO 410
200          ENDIF
201        ENDDO
202        LATIT = NEWGEO(NJ)
203C
204  410   CONTINUE
205C
206        LONG = NEXT - NEWOFF(LATIT)
207C
208        RLAT = RLATNEW(LATIT)
209        IF( NEWGEO(1).EQ.JPGAUSS ) THEN
210          IF( NEWGEO(17).EQ.0 ) THEN
211            RLON = REAL((LONG-1)*NEWGEO(9))/JP1000
212          ELSE
213            RLON = (REAL(LONG-1)*360.0)/REAL(NEWGEO(NPOINTS+LATIT-1))
214          ENDIF
215        ELSE
216          RLON = REAL((LONG-1)*NEWGEO(9))/JP1000
217        ENDIF
218C
219C       Find type of point (land or sea)
220C
221        IF( MASTER(NEXT).GT.MASTERTHRESHOLD ) THEN
222          NEWTYPE = LAND
223        ELSE
224          NEWTYPE = SEA
225        ENDIF
226C
227C       Find four neighbours in the old field with their types
228C       (Find NW neighbour and deduce the others).
229C
230        LAT(NORTH) = JNORSGG(RLAT,RLATOLD,OLDGEO(NGAUSS),1)
231        LAT(SOUTH) = JNORSGG(RLAT,RLATOLD,OLDGEO(NGAUSS),0)
232C
233        OLAT(NORTH) = RLATOLD(LAT(NORTH))
234        OLAT(SOUTH) = RLATOLD(LAT(SOUTH))
235C
236        NPTS = OLDGEO(NPOINTS-1+LAT(NORTH))
237        LON(NWEST) = 1 + INT(RLON/(360.0/REAL(NPTS)))
238        LON(NEAST) = LON(NWEST) + 1
239        IF( LON(NEAST).GT.NPTS ) LON(NEAST) = 1
240C
241        OLON(NWEST) = (REAL(LON(NWEST)-1)*360.0)/REAL(NPTS)
242        IF( LON(NEAST).EQ.1 ) THEN
243          OLON(NEAST) = 360.0
244        ELSE
245          OLON(NEAST) = (REAL(LON(NEAST)-1)*360.0)/REAL(NPTS)
246        ENDIF
247C
248        NPTS = OLDGEO(NPOINTS-1+LAT(SOUTH))
249        LON(SWEST) = 1 + INT(RLON/(360.0/REAL(NPTS)))
250        LON(SEAST) = LON(SWEST) + 1
251        IF( LON(SEAST).GT.NPTS ) LON(SEAST) = 1
252C
253        OLON(SWEST) = (REAL(LON(SWEST)-1)*360.0)/REAL(NPTS)
254        IF( LON(SEAST).EQ.1 ) THEN
255          OLON(SEAST) = 360.0
256        ELSE
257          OLON(SEAST) = (REAL(LON(SEAST)-1)*360.0)/REAL(NPTS)
258        ENDIF
259C
260        PT(NWEST) = OLDOFF(LAT(NORTH)) + LON(NWEST)
261        IF( OLDLSM(PT(NWEST)).GT.OLDLSMTHRESHOLD ) THEN
262          TYPE(NWEST) = LAND
263        ELSE
264          TYPE(NWEST) = SEA
265        ENDIF
266C
267        PT(NEAST) = OLDOFF(LAT(NORTH)) + LON(NEAST)
268        IF( OLDLSM(PT(NEAST)).GT.OLDLSMTHRESHOLD ) THEN
269          TYPE(NEAST) = LAND
270        ELSE
271          TYPE(NEAST) = SEA
272        ENDIF
273C
274        PT(SWEST) = OLDOFF(LAT(SOUTH)) + LON(SWEST)
275        IF( OLDLSM(PT(SWEST)).GT.OLDLSMTHRESHOLD ) THEN
276          TYPE(SWEST) = LAND
277        ELSE
278          TYPE(SWEST) = SEA
279        ENDIF
280C
281        PT(SEAST) = OLDOFF(LAT(SOUTH)) + LON(SEAST)
282        IF( OLDLSM(PT(SEAST)).GT.OLDLSMTHRESHOLD ) THEN
283          TYPE(SEAST) = LAND
284        ELSE
285          TYPE(SEAST) = SEA
286        ENDIF
287C
288C       Interpolate the new value
289C
290        NEWFLD(NEXT) =
291     X    ISLPROC(RLAT,RLON,NEWTYPE,OLAT,OLON,TYPE,PT,OLDFLD)
292C
293      ENDDO
294C
295C     -----------------------------------------------------------------|
296C*    Section 9.   Return
297C     -----------------------------------------------------------------|
298C
299  900 CONTINUE
300C
301C
302      RETURN
303      END
304
305