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 IGTOG (PIFELD, KIWE, KINS, KOWE, KONS, KWEIND,
12     1   KNSIND, PWFACT, POFELD, KPR, KERR)
13C
14C---->
15C**** *IGTOG*
16C
17C     Purpose
18C     -------
19C
20C     Perform basic interpolation between regular input and output
21C     fields.
22C
23C     Interface
24C     ---------
25C
26C     IERR = IGTOG (PIFELD, KIWE, KINS, KOWE, KONS, KWEIND, KNSIND,
27C    1   PWFACT, POFELD, KPR, KERR)
28C
29C     Input parameters
30C     ----------------
31C
32C     PIFELD     - The input field provided by the calling routine.
33C
34C     KIWE       - The number of points in the West-East direction in
35C                  the input field.
36C
37C     KINS       - The number of points in the North-South direction
38C                  in the input field.
39C
40C     KOWE       - The number of points in the West-East direction in
41C                  the output field.
42C
43C     KONS       - The number of points in the North-South direction
44C                  in the output field.
45C
46C     KWEIND     - This array contains the array offsets of the West
47C                  and East points in the input array required for
48C                  interpolation.
49C
50C     KNSIND     - This array contains the array offsets of the North
51C                  and South points in the input array required for
52C                  interpolation.
53C
54C     PWFACT     - The array of interpolating weights to the four
55C                  neighbouring points for every output point.
56C
57C     KPR        - The debug print switch.
58C                  0  , No debugging output.
59C                  1  , Produce debugging output.
60C
61C     KERR       - The error control flag.
62C                  -ve, No error message. Return error code.
63C                  0  , Hard failure with error message.
64C                  +ve, Print error message. Return error code.
65C
66C     Output parameters
67C     -----------------
68C
69C     POFELD     - The output field returned to the calling routine.
70C
71C     Return value
72C     ------------
73C
74C     The error indicator (INTEGER).
75C
76C     Common block usage
77C     ------------------
78C
79C     None
80C
81C     Externals
82C     ---------
83C
84C     INTLOG - Logs messages
85C     FORCED_NEAREST_NEIGHBOUR - check forced interpolation method
86C
87C
88C     Method
89C     ------
90C
91C     This routine performs basic linear interpolation using the four
92C     neighbouring points in the input array to generate the output
93C     array.
94C
95C     Reference
96C     ---------
97C
98C     None
99C
100C     Comments
101C     --------
102C
103C     None
104C
105C     Author
106C     ------
107C
108C     K. Fielding      *ECMWF*      Oct 1993
109C
110C     Modifications
111C     -------------
112C
113C     Allow for missing data values
114C     J.D.Chambers      ECMWF       August 2000
115C
116C     Force nearest neighbour processing with env variable or
117C     INTOUT parameter
118C     S.Curic           ECMWF       September 2005
119C
120C----<
121C     -----------------------------------------------------------------|
122C*    Section 0. Definition of variables.
123C     -----------------------------------------------------------------|
124C
125      IMPLICIT NONE
126C
127#include "parim.h"
128#include "nifld.common"
129#include "nofld.common"
130C
131C     Function arguments
132C
133      INTEGER KIWE, KINS, KOWE, KONS, KPR, KERR
134      INTEGER KWEIND (2, KOWE), KNSIND (2, KONS)
135      REAL PIFELD (KIWE, KINS), POFELD (KOWE, KONS)
136      REAL PWFACT (4, KOWE, KONS)
137C
138C     Local variables
139C
140      INTEGER ILATN, ILATS, JOLAT, JOLON, COUNT
141      REAL NEAREST
142      LOGICAL LVEGGY
143C
144C     Externals
145C
146      LOGICAL FORCED_NEAREST_NEIGHBOUR
147C
148C     Statement function
149C
150      REAL A, B
151      LOGICAL NOTEQ
152      NOTEQ(A,B) = (ABS((A)-(B)).GT.(ABS(A)*1E-3))
153C
154C     -----------------------------------------------------------------|
155C*    Section 1. Initialisation
156C     -----------------------------------------------------------------|
157C
158  100 CONTINUE
159C
160      IF( KPR.GE.1 ) CALL INTLOG(JP_DEBUG,'IGTOG: Section 1.',JPQUIET)
161C
162      IGTOG = 0
163C
164      IF( KPR.GE.1 ) THEN
165        CALL INTLOG(JP_DEBUG,'IGTOG: Input parameters.',JPQUIET)
166        CALL INTLOG(JP_DEBUG,'IGTOG: No.input fld longitudes = ',KIWE)
167        CALL INTLOG(JP_DEBUG,'IGTOG: No.input fld latitudes  = ',KINS)
168        CALL INTLOG(JP_DEBUG,'IGTOG: No.output fld longitudes = ',KOWE)
169        CALL INTLOG(JP_DEBUG,'IGTOG: No.output fld latitudes  = ',KONS)
170      ENDIF
171C
172C     Use nearest neighbour if required
173      LVEGGY = FORCED_NEAREST_NEIGHBOUR(LMETHOD,NITABLE,NIPARAM)
174      IF( LVEGGY ) CALL INTLOG(JP_DEBUG,
175     X  'IGTOG: nearest neighbour processing',JPQUIET)
176C
177C     -----------------------------------------------------------------|
178C*    Section 2. Basic interpolation
179C     -----------------------------------------------------------------|
180C
181  200 CONTINUE
182C
183      IF( KPR.GE.1 ) CALL INTLOG(JP_DEBUG,'IGTOG: Section 2.',JPQUIET)
184C
185      DO JOLAT = 1, KONS
186C
187        ILATN = KNSIND(JP_I_N,JOLAT)
188        ILATS = KNSIND(JP_I_S,JOLAT)
189C
190        DO JOLON = 1, KOWE
191C
192C         Count non-missing data values
193C
194          IF( LIMISSV ) THEN
195            COUNT = 0
196            IF( NOTEQ(PIFELD(KWEIND(JP_I_W,JOLON),ILATN), RMISSGV) )
197     X        COUNT = COUNT + 1
198            IF( NOTEQ(PIFELD(KWEIND(JP_I_W,JOLON),ILATS), RMISSGV) )
199     X        COUNT = COUNT + 1
200            IF( NOTEQ(PIFELD(KWEIND(JP_I_E,JOLON),ILATN), RMISSGV) )
201     X        COUNT = COUNT + 1
202            IF( NOTEQ(PIFELD(KWEIND(JP_I_E,JOLON),ILATS), RMISSGV) )
203     X        COUNT = COUNT + 1
204          ELSE
205            COUNT = 4
206          ENDIF
207C
208C         Interpolate using four neighbours if none are missing
209C
210          IF( (COUNT.EQ.4).AND.(.NOT.LVEGGY) ) THEN
211            POFELD(JOLON,JOLAT) =
212     X        PIFELD(KWEIND(JP_I_W,JOLON),ILATN) *
213     X          PWFACT(JP_I_NW,JOLON,JOLAT) +
214     X        PIFELD(KWEIND(JP_I_W,JOLON),ILATS) *
215     X          PWFACT(JP_I_SW,JOLON,JOLAT) +
216     X        PIFELD(KWEIND(JP_I_E,JOLON),ILATN) *
217     X          PWFACT(JP_I_NE,JOLON,JOLAT) +
218     X        PIFELD(KWEIND(JP_I_E,JOLON),ILATS) *
219     X          PWFACT(JP_I_SE,JOLON,JOLAT)
220C
221C         Set missing if all neighbours are missing
222C
223          ELSE IF( COUNT.EQ.0 ) THEN
224            POFELD(JOLON,JOLAT) = RMISSGV
225C
226C         Otherwise, use the nearest neighbour
227C
228          ELSE
229            NEAREST = PWFACT(JP_I_NW,JOLON,JOLAT)
230            POFELD(JOLON,JOLAT) =
231     X        PIFELD(KWEIND(JP_I_W,JOLON),ILATN)
232C
233            IF( PWFACT(JP_I_NE,JOLON,JOLAT).GT.NEAREST ) THEN
234              NEAREST = PWFACT(JP_I_NE,JOLON,JOLAT)
235              POFELD(JOLON,JOLAT) =
236     X          PIFELD(KWEIND(JP_I_E,JOLON),ILATN)
237            ENDIF
238C
239            IF( PWFACT(JP_I_SW,JOLON,JOLAT).GT.NEAREST ) THEN
240              NEAREST = PWFACT(JP_I_SW,JOLON,JOLAT)
241              POFELD(JOLON,JOLAT) =
242     X          PIFELD(KWEIND(JP_I_W,JOLON),ILATS)
243            ENDIF
244C
245            IF( PWFACT(JP_I_SE,JOLON,JOLAT).GT.NEAREST ) THEN
246              NEAREST = PWFACT(JP_I_SE,JOLON,JOLAT)
247              POFELD(JOLON,JOLAT) =
248     X          PIFELD(KWEIND(JP_I_E,JOLON),ILATS)
249            ENDIF
250C
251          ENDIF
252C
253        ENDDO
254C
255      ENDDO
256C
257C     -----------------------------------------------------------------|
258C*    Section 9. Return to calling routine. Format statements
259C     -----------------------------------------------------------------|
260C
261  900 CONTINUE
262C
263      IF( KPR.GE.1 ) CALL INTLOG(JP_DEBUG,'IGTOG: Section 9.',JPQUIET)
264C
265      RETURN
266      END
267
268