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 FIXAREA( )
12C
13C---->
14C**** FIXAREA
15C
16C     Purpose
17C     -------
18C
19C     Fixup input/output field area definitions.
20C
21C
22C     Interface
23C     ---------
24C
25C     IRET = FIXAREA( )
26C
27C     Input
28C     -----
29C
30C     None.
31C
32C
33C     Output
34C     ------
35C
36C     None.
37C
38C
39C     Method
40C     ------
41C
42C     If default (0/0/0/0) selected for input, input area is set
43C       - to global for spherical harmonics
44C       - to global for lat/long grid
45C       - to global for gaussian grid
46C
47C     If default (0/0/0/0) selected for output, output area is set
48C       - to same as input area for lat/long grid
49C       - to same as input area for regular gaussian grid
50C       - to global for reduced gaussian grid
51C
52C     (Currently, subareas are not supported for reduced gaussian fields.
53C      Should work OK; but need to correct setup of values in GRIB
54C      product for number of points in each latitude row)
55C
56C     Output area is adjusted to fit the given grid step.
57C
58C
59C     Externals
60C     ---------
61C
62C     AREACHK - Match input/output field area definitions according to
63C               grid specification.
64C     GETENV   - Get value of an environment variable
65C     INTLOG   - Logs output messages
66C
67C
68C     Author
69C     ------
70C
71C     J.D.Chambers     ECMWF     Jan 1995
72C
73C----<
74C
75      IMPLICIT NONE
76C
77C     Function arguments
78C
79#include "parim.h"
80#include "nifld.common"
81#include "nofld.common"
82#include "grfixed.h"
83#include "intf.h"
84#include "current.h"
85C
86C     Parameters
87C
88      INTEGER JPROUTINE
89      PARAMETER (JPROUTINE = 19200 )
90C
91C     Local variables
92C
93      INTEGER INORTH, ISOUTH, IEAST
94      INTEGER ITEMP, IRET
95      REAL EW, NS, NORTH, SOUTH, EAST, WEST, FACTOR3, FACTOR4
96      LOGICAL LDEFIN, LDEFOUT, LGLOBAL, LOVERDE
97      CHARACTER*20 OVERIDE
98      INTEGER IBLANK
99C
100C     Externals
101C
102      INTEGER AREACHK
103C
104C ------------------------------------------------------------------
105C*    Section 1.   Initialise
106C ------------------------------------------------------------------
107C
108  100 CONTINUE
109      FIXAREA = 0
110      FACTOR3 = 0.
111      FACTOR4 = 0.
112      LDEFIN  = .FALSE.
113      LGLOBAL = .FALSE.
114C
115C ------------------------------------------------------------------
116C*    Section 2.   Fixup input area if default (0/0/0/0).
117C ------------------------------------------------------------------
118C
119  200 CONTINUE
120C
121C     Fixup input ocean area definition
122C
123      IF( NILOCAL.EQ.4 ) THEN
124        FACTOR3 = 0.0
125        FACTOR4 = 0.0
126C
127        IF( NIFORM.EQ.1 ) THEN
128          IF( ISEC1(60).EQ.1 ) FACTOR3 = 1.0
129          IF( ISEC1(60).EQ.2 ) THEN
130            IF( ISEC1(47).EQ.160 ) FACTOR3 = 1000.0
131          ENDIF
132          IF( ISEC1(60).EQ.3 ) FACTOR3 = 1000000.0
133          IF( ISEC1(60).EQ.4 ) FACTOR3 = 1000000.0
134C
135          IF( ISEC1(61).EQ.1 ) FACTOR4 = 1.0
136          IF( ISEC1(61).EQ.2 ) THEN
137            IF( ISEC1(47).EQ.160 ) FACTOR4 = 1000.0
138          ENDIF
139          IF( ISEC1(61).EQ.3 ) FACTOR4 = 1000000.0
140          IF( ISEC1(61).EQ.4 ) FACTOR4 = 1000000.0
141C
142          IF( (FACTOR3.EQ.0).OR.(FACTOR4.EQ.0) ) THEN
143            FIXAREA = 2
144            GOTO 900
145          ENDIF
146          NORTH = REAL(ISEC1(62))/FACTOR4
147          WEST  = REAL(ISEC1(63))/FACTOR3
148          SOUTH = REAL(ISEC1(64))/FACTOR4
149          EAST  = REAL(ISEC1(65))/FACTOR3
150        ELSE
151          IF( NIOCO3.EQ.1 ) FACTOR3 = 1.0
152          IF( NIOCO3.EQ.2 ) THEN
153            IF( NIVCDEF.EQ.160 ) FACTOR3 = 1000.0
154          ENDIF
155          IF( NIOCO3.EQ.3 ) FACTOR3 = 1000000.0
156          IF( NIOCO3.EQ.4 ) FACTOR3 = 1000000.0
157C
158          IF( NIOCO4.EQ.1 ) FACTOR4 = 1.0
159          IF( NIOCO4.EQ.2 ) THEN
160            IF( NIVCDEF.EQ.160 ) FACTOR4 = 1000.0
161          ENDIF
162          IF( NIOCO4.EQ.3 ) FACTOR4 = 1000000.0
163          IF( NIOCO4.EQ.4 ) FACTOR4 = 1000000.0
164C
165          IF( (FACTOR3.EQ.0).OR.(FACTOR4.EQ.0) ) THEN
166            FIXAREA = 2
167            GOTO 900
168          ENDIF
169          NORTH = REAL(NIOCO4F)/FACTOR4
170          WEST  = REAL(NIOCO3F)/FACTOR3
171          SOUTH = REAL(NIOCO4L)/FACTOR4
172          EAST  = REAL(NIOCO3L)/FACTOR3
173        ENDIF
174        NIAREA(1) = NINT( NORTH * PPMULT + 0.1)
175        NIAREA(2) = NINT( WEST  * PPMULT + 0.1)
176        NIAREA(3) = NINT( SOUTH * PPMULT + 0.1)
177        NIAREA(4) = NINT( EAST  * PPMULT + 0.1)
178        GOTO 300
179      ENDIF
180C
181C Sinisa bug fix for grid global fields
182C       If input is lat/lon, check if it is global to within
183C       a tolerance of 0.1 degrees
184       IF ( NIREPR .EQ. JPREGULAR ) THEN
185          IF ( NIAREA(4).GT.0 ) THEN
186            IEAST = NIAREA(4)
187          ELSE
188            IEAST = JP360 + NIAREA(4)
189          ENDIF
190         IF(NIAREA(1).EQ.JP90.AND.NIAREA(3).EQ.-JP90.AND.
191     X      NIAREA(2).EQ.0.AND.
192     X      IABS(JP360 - IEAST - NIGRID(1)).LT.1000) THEN
193              NIAREA(1) = 0
194              NIAREA(2) = 0
195              NIAREA(3) = 0
196              NIAREA(4) = 0
197csinisa
198cs          ELSEIF(IABS(JP360 - IEAST - NIGRID(1)).LT.1000) THEN
199cs              NOAREA(1) = NIAREA(1)
200cs              NOAREA(2) = 0
201cs              NOAREA(3) = NIAREA(3)
202cs              NOAREA(4) = 36000000
203          ENDIF
204        ENDIF
205
206      LDEFIN = ( (NIAREA(1) .EQ. 0) .AND. (NIAREA(2) .EQ. 0) .AND.
207     X           (NIAREA(3) .EQ. 0) .AND. (NIAREA(4) .EQ. 0) )
208
209C
210      IF( LDEFIN ) THEN
211C
212        IF ( (NIREPR.EQ.JPSPHERE) .OR. (NIREPR.EQ.JPSPHROT) ) THEN
213C
214C         Spectral input ..
215          EW = 0.0
216          NS = 0.0
217C
218        ELSE IF ( NIREPR .EQ. JPREGULAR ) THEN
219C
220C         Regular lat/long grid ..
221          EW = FLOAT( NIGRID(1) ) / PPMULT
222          NS = FLOAT( NIGRID(2) ) / PPMULT
223C
224        ELSE
225C
226C         Gaussian grid ..
227          EW = FLOAT( NIGAUSS )
228          NS = 0.0
229        ENDIF
230C
231        NORTH = 0.0
232        WEST  = 0.0
233        SOUTH = 0.0
234        EAST  = 0.0
235        IRET = AREACHK( EW, NS, NORTH, WEST, SOUTH, EAST )
236        IF( IRET.NE.0 ) THEN
237          FIXAREA = IRET
238          GOTO 900
239        ENDIF
240        NIAREA(1) = NINT( NORTH * PPMULT + 0.1)
241        NIAREA(2) = NINT( WEST  * PPMULT + 0.1)
242        NIAREA(3) = NINT( SOUTH * PPMULT + 0.1)
243        NIAREA(4) = NINT( EAST  * PPMULT + 0.1)
244C
245      ENDIF
246C
247C ------------------------------------------------------------------
248C*    Section 3.   Fixup output area if default (0/0/0/0).
249C ------------------------------------------------------------------
250C
251  300 CONTINUE
252C
253C     See if environment variable has been specified to override
254C     output area specification
255C
256      LOVERDE = .FALSE.
257      CALL GETENV('OVERRIDE_OUTPUT_AREA', OVERIDE)
258      IBLANK = INDEX(OVERIDE, ' ')
259      IF( IBLANK.GT.1 ) THEN
260        IF( OVERIDE(1:2).EQ.'ON' ) THEN
261          LOVERDE = .TRUE.
262          CALL INTLOG(JP_DEBUG,
263     X      'FIXAREA: OVERRIDE_OUTPUT_AREA is ON',JPQUIET)
264        ENDIF
265      ENDIF
266C
267      LDEFOUT = ( (NOAREA(1) .EQ. 0) .AND. (NOAREA(2) .EQ. 0) .AND.
268     X            (NOAREA(3) .EQ. 0) .AND. (NOAREA(4) .EQ. 0) )
269     X          .OR. LOVERDE
270C
271C     Fixup output ocean area definition
272C
273      IF( NILOCAL.EQ.4 ) THEN
274        IF( LDEFOUT ) THEN
275          IF( NIFORM.EQ.1 ) THEN
276              SOUTH  = REAL(ISEC1(64))/FACTOR4
277              NORTH  = REAL(ISEC1(62))/FACTOR4
278              WEST   = REAL(ISEC1(63))/FACTOR3
279              EAST   = REAL(ISEC1(65))/FACTOR3
280
281              EW = FLOAT( NOGRID(1) ) / PPMULT
282              NS = FLOAT( NOGRID(2) ) / PPMULT
283
284            IF( ISEC1(60).EQ.3 ) THEN
285              WEST  = 0.0
286              EAST  = 360.0
287              IRET = AREACHK( EW, NS, NORTH, WEST, SOUTH, EAST )
288              IF ( IRET .NE. 0 ) THEN
289                 FIXAREA = IRET
290                 GOTO 900
291              ENDIF
292            ELSEIF( ISEC1(60).EQ.4 ) THEN
293              WEST  = -90.0
294              EAST  =  90.0
295            ENDIF
296
297            IF( ISEC1(61).EQ.3 ) THEN
298              SOUTH  = 0.0
299              NORTH  = 360.0
300            ELSEIF( ISEC1(61).EQ.4 ) THEN
301              SOUTH  = -90.0
302              NORTH  =  90.0
303            ENDIF
304          ELSE
305              NORTH = REAL(NIOCO4F)/FACTOR4
306              WEST  = REAL(NIOCO3F)/FACTOR3
307              SOUTH = REAL(NIOCO4L)/FACTOR4
308              EAST  = REAL(NIOCO3L)/FACTOR3
309              EW = FLOAT( NOGRID(1) ) / PPMULT
310              NS = FLOAT( NOGRID(2) ) / PPMULT
311            IF( NIOCO3.EQ.3 ) THEN
312              WEST  = 0.0
313              EAST  = 360.0
314              IRET = AREACHK( EW, NS, NORTH, WEST, SOUTH, EAST )
315              IF ( IRET .NE. 0 ) THEN
316                 FIXAREA = IRET
317                 GOTO 900
318              ENDIF
319            ELSEIF( NIOCO3.EQ.4 ) THEN
320              WEST  = -90.0
321              EAST  =  90.0
322            ENDIF
323
324            IF( NIOCO4.EQ.3 ) THEN
325              SOUTH  = 0.0
326              NORTH  = 360.0
327            ELSEIF( NIOCO4.EQ.4 ) THEN
328              SOUTH  = -90.0
329              NORTH  =  90.0
330            ENDIF
331          ENDIF
332
333          NOAREA(1) = NINT( NORTH * PPMULT + 0.1)
334          NOAREA(2) = NINT( WEST  * PPMULT + 0.1)
335          NOAREA(3) = NINT( SOUTH * PPMULT + 0.1)
336          NOAREA(4) = NINT( EAST  * PPMULT + 0.1)
337        ELSE
338          NORTH = REAL(NOAREA(1)) / PPMULT
339          WEST  = REAL(NOAREA(2)) / PPMULT
340          SOUTH = REAL(NOAREA(3)) / PPMULT
341          EAST  = REAL(NOAREA(4)) / PPMULT
342        ENDIF
343        GOTO 900
344      ENDIF
345C
346      IF( LDEFOUT ) THEN
347C
348C       If input is gaussian, check if it is global to within
349C       a tolerance of 0.1 degrees.
350C
351        IF ( (NIREPR.EQ.JPGAUSSIAN) .OR. (NIREPR.EQ.JPQUASI) ) THEN
352          INORTH = NINT( RIGAUSS(1) * PPMULT )
353          ISOUTH = -INORTH
354c         EMOS-199: adjusted for reduced_gg/octahedral
355c         IEAST  = JP360 - (JP90/NIGAUSS)
356          IEAST  = JP360 - (JP360/MILLEN(NIGAUSS))
357          LGLOBAL = (NIAREA(2).EQ.0)
358     X      .AND. (IABS(NIAREA(1)-INORTH).LE.1000)
359     X      .AND. (IABS(NIAREA(3)-ISOUTH).LE.1000)
360     X      .AND. (IABS(NIAREA(4)-IEAST ).LE.1000 .OR.
361     X        NIREPR.EQ.JPGAUSSIAN .OR. .TRUE. )
362C         NOTE: a similar (but useless) test was done here before adding
363C         .TRUE. be replicate the same (useless) behaviour
364        ENDIF
365C
366C       Spectral output ..
367C
368C       .. skip section fitting output area to grid
369C
370        IF ( (NOREPR.EQ.JPSPHERE).OR.(NOREPR.EQ.JPSPHROT) ) GOTO 900
371C
372C       Regular lat/long grid ..
373C
374        IF ( (NOREPR.EQ.JPREGULAR).OR.(NOREPR.EQ.JPREGROT) ) THEN
375          IF( LDEFIN .OR. LGLOBAL.OR.(NIREPR.EQ.JPREDLL) ) THEN
376            EW = FLOAT( NOGRID(1) ) / PPMULT
377            NS = FLOAT( NOGRID(2) ) / PPMULT
378            NORTH = 90.0
379            WEST  = 0.0
380            SOUTH = -NORTH
381            EAST  = 360.0
382C
383C           Use GRIB header values for reduced lat/long grids
384C           (maybe 'Mediterranean' sub-area).
385C
386            IF( NIREPR.EQ.JPREDLL ) THEN
387              NORTH = REAL(NIAREA(1))/PPMULT
388              WEST  = REAL(NIAREA(2))/PPMULT
389              SOUTH = REAL(NIAREA(3))/PPMULT
390              EAST  = REAL(NIAREA(4))/PPMULT
391            ENDIF
392            IRET = AREACHK( EW, NS, NORTH, WEST, SOUTH, EAST )
393            IF ( IRET .NE. 0 ) THEN
394               FIXAREA = IRET
395               GOTO 900
396            ENDIF
397            NOAREA(1) = NINT( NORTH * PPMULT + 0.1)
398            NOAREA(2) = NINT( WEST  * PPMULT + 0.1)
399            NOAREA(3) = NINT( SOUTH * PPMULT + 0.1)
400            NOAREA(4) = NINT( EAST  * PPMULT + 0.1)
401C
402          ELSE
403            NOAREA(1) = NIAREA(1)
404            NOAREA(2) = NIAREA(2)
405            NOAREA(3) = NIAREA(3)
406            NOAREA(4) = NIAREA(4)
407          ENDIF
408C
409        ENDIF
410C
411C       Gaussian ..
412C
413        IF ( (NOREPR.EQ.JPGAUSSIAN) .OR. (NOREPR.EQ.JPQUASI) ) THEN
414C
415C         Is output grid specification the same as the input?
416          IF( NIGAUSS.EQ.NOGAUSS) THEN
417            IF( (NOREPR.EQ.JPGAUSSIAN) .AND. LGLOBAL ) THEN
418              NOAREA(1) = JP90
419              NOAREA(2) = 0
420              NOAREA(3) = -JP90
421c     EMOS-199: adjusted for reduced_gg/octahedral
422c             NOAREA(4) = JP360 - (JP90/NOGAUSS)
423              NOAREA(4) = JP360 - (JP360/NOLPTS(NOGAUSS))
424            ELSE
425              NOAREA(1) = NIAREA(1)
426              NOAREA(2) = NIAREA(2)
427              NOAREA(3) = NIAREA(3)
428              NOAREA(4) = NIAREA(4)
429           ENDIF
430C
431C           Skip section fitting output area to grid
432            GOTO 900
433C
434          ELSE
435C
436C           Different grid resolutions
437            EW = FLOAT( NOGAUSS )
438            NS = 0.0
439            IF( (NOREPR.EQ.JPGAUSSIAN) .AND. LGLOBAL ) THEN
440              NORTH =  90.0
441              WEST  =   0.0
442              SOUTH = -90.0
443c     EMOS-199: adjusted for reduced_gg/octahedral
444c             EAST  = 360.0 - (360.0/(EW*4.0))
445              EAST  = 360.0 - (360.0/FLOAT(NOLPTS(NOGAUSS)))
446            ELSE
447              NORTH = FLOAT( NIAREA(1) ) / PPMULT
448              WEST  = FLOAT( NIAREA(2) ) / PPMULT
449              SOUTH = FLOAT( NIAREA(3) ) / PPMULT
450              EAST  = FLOAT( NIAREA(4) ) / PPMULT
451            ENDIF
452            IRET = AREACHK( EW, NS, NORTH, WEST, SOUTH, EAST )
453            IF ( IRET .NE. 0 ) THEN
454               FIXAREA = IRET
455               GOTO 900
456            ENDIF
457            NOAREA(1) = NINT( NORTH * PPMULT + 0.1)
458            NOAREA(2) = NINT( WEST  * PPMULT + 0.1)
459            NOAREA(3) = NINT( SOUTH * PPMULT + 0.1)
460            NOAREA(4) = NINT( EAST  * PPMULT + 0.1)
461C
462C           If regular lat/long input, check if output west longitude
463C           reachs to full globe
464            ITEMP = NIAREA(4)+NIGRID(1)
465            IF( (NIREPR.EQ.JPREGULAR .AND. ITEMP.EQ.JP360) .OR.
466     X          (NIREPR.EQ.JPSPHERE) .OR.
467     X          (NIREPR.EQ.JPSPHROT) ) THEN
468c     EMOS-199: adjusted for reduced_gg/octahedral
469c             NOAREA(4) = JP360 - (JP90/NOGAUSS)
470              NOAREA(4) = JP360 - (JP360/NOLPTS(NOGAUSS))
471            ENDIF
472          ENDIF
473        ENDIF
474C
475      ENDIF
476C
477C ------------------------------------------------------------------
478C*    Section 4.   Now fixup output areas to correspond to the grid.
479C ------------------------------------------------------------------
480C
481  400 CONTINUE
482C
483      IF ( (NOREPR.EQ.JPREGULAR).OR.(NOREPR.EQ.JPREGROT) ) THEN
484C
485C       Regular lat/long grid ..
486        EW = FLOAT( NOGRID(1) ) / PPMULT
487        NS = FLOAT( NOGRID(2) ) / PPMULT
488C
489      ELSE
490C
491C       Gaussian grid ..
492        EW = FLOAT( NOGAUSS )
493        NS = 0.0
494      ENDIF
495C
496      NORTH = FLOAT( NOAREA(1) ) / PPMULT
497      WEST  = FLOAT( NOAREA(2) ) / PPMULT
498      SOUTH = FLOAT( NOAREA(3) ) / PPMULT
499      EAST  = FLOAT( NOAREA(4) ) / PPMULT
500      IRET = AREACHK( EW, NS, NORTH, WEST, SOUTH, EAST )
501      IF ( IRET .NE. 0 ) THEN
502         FIXAREA = IRET
503         GOTO 900
504      ENDIF
505      NOAREA(1) = NINT( NORTH * PPMULT + 0.1)
506      NOAREA(2) = NINT( WEST  * PPMULT + 0.1)
507      NOAREA(3) = NINT( SOUTH * PPMULT + 0.1)
508      NOAREA(4) = NINT( EAST  * PPMULT + 0.1)
509
510C EMOS-199: [F|N] grids assume max(pl)=4*N, however O grids have
511C max(pl)=4*+16 and AREACHK doesn't know of the difference
512C The code below only checks if the adjusted EAST matches what
513C would happen on a [F|N] global grid, and corrects it for an O grid.
514      IF ( (NOREPR.EQ.JPQUASI)
515     X  .AND. (HOGAUST.EQ.'O'
516     X  .AND. NOAREA(2).EQ.0
517     X  .AND. NOAREA(4).EQ.(JP360-JP90/NOGAUSS) ))
518     X  NOAREA(4) = JP360 - (JP360/NOLPTS(NOGAUSS))
519
520C
521C ------------------------------------------------------------------
522C*    Section 9.   Closedown.
523C ------------------------------------------------------------------
524C
525  900 CONTINUE
526C
527      RETURN
528      END
529