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 HGENGG(AREA,POLE,KGNNEW,HTYPE,KNPTS,GNLATS,
12     X                        IGSIZE,RLAT,RLON,NUMPTS)
13C
14C---->
15C**** HGENGG
16C
17C     Purpose
18C     -------
19C
20C     This routine calculates the original latitude and longitude
21C     values (before rotation) for a rotated gaussian grid.
22C
23C
24C     Interface
25C     ---------
26C
27C     IRET = HGENGG(AREA,POLE,KGNNEW,HTYPE,KNPTS,IGSIZE,
28C    X              RLAT,RLON,NUMPTS)
29C
30C
31C     Input parameters
32C     ----------------
33C
34C     AREA     - Limits of area (N/W/S/E)
35C     POLE     - Pole of rotation (lat/long)
36C     KGNNEW   - Gaussian grid number
37C     HTYPE    - Gaussian grid type
38C                = 'R' for reduced ("quasi-regular"),
39C                = 'O' for reduced/octahedral,
40C                = 'F' for full,
41C                = 'U' for a user-defined gaussian grid
42C     KNPTS    - Array containing number of points at each latitude
43C     GNLATS   - Array containing list of gaussian latitudes
44C     IGSIZE   - The size of the array to fill with the gaussian field
45C     NUMPTWE  - Calculate number of grid pts in range from west to east
46C
47C
48C     Output parameters
49C     -----------------
50C
51C     RLAT  - Vector of orginal latitude values for the points in
52C             the rotated grid.
53C     RLON  - Vector of orginal longitude values for the points in
54C             the rotated grid.
55C     NUMPTS - Number of points in new field
56C
57C     Returns 0 if function successful, non-zero otherwise.
58C
59C
60C     Common block usage
61C     ------------------
62C
63C     None.
64C
65C
66C     Method
67C     ------
68C
69C     The vector of points runs from West to East in rows, the rows
70C     run from North to South.
71C
72C
73C     Externals
74C     ---------
75C
76C     INTLOG  - Log error message.
77C     JDEBUG  - Tests if debug output required.
78C     HLL2XYZ - Converts a latitude/longitude position to (x,y,z)
79C               wrt axes through the centre of the globe. The z-axis
80C               runs from the south to north pole. The x- and y-axes
81C               are in the plane of the equator with the x-axis
82C               pointing out through lat/long (0,0).
83C     YROTATE - Rotates the globe about the new y-axis.
84C     XYZ2LL  - Converts an (x,y,z) position to a latitude/longitude.
85C     JMALLOC - Dynamically allocate memory
86C     JFREE   - Free dynamically allocated memory
87C
88C
89C     Reference
90C     ---------
91C
92C     None.
93C
94C
95C     Comments
96C     --------
97C
98C     None.
99C
100C
101C     Author
102C     ------
103C
104C     J.D.Chambers      ECMWF      February 2001
105C
106C
107C     Modifications
108C     -------------
109C
110C     None.
111C
112C----<
113C     -----------------------------------------------------------------|
114C*    Section 0. Definition of variables.
115C     -----------------------------------------------------------------|
116C
117      IMPLICIT NONE
118C
119#include "parim.h"
120#include "nofld.common"
121C
122C     Parameters
123C
124      INTEGER JNORTH, JSOUTH, JWEST, JEAST, JW_E, JN_S, JLAT, JLON
125      REAL JSCALE
126      PARAMETER (JNORTH = 1 )
127      PARAMETER (JWEST  = 2 )
128      PARAMETER (JSOUTH = 3 )
129      PARAMETER (JEAST  = 4 )
130      PARAMETER (JW_E  = 1 )
131      PARAMETER (JN_S  = 2 )
132      PARAMETER (JLAT  = 1 )
133      PARAMETER (JLON  = 2 )
134      PARAMETER (JSCALE  = 1000.0 )
135C
136C     Function arguments
137C
138      CHARACTER*1 HTYPE
139      INTEGER KGNNEW,KNPTS(KGNNEW*2),IGSIZE,NUMPTS
140      REAL AREA(4),POLE(2),RLAT(*),RLON(*), GNLATS(KGNNEW*2)
141C
142C     Local variables
143C
144      INTEGER LOOP, LOOPI, LOOPO, IGRIDNI, NLON
145      INTEGER ISIZE, ISIZOLD, NBYTES
146      REAL GRID, X(1), Y(1), Z(1)
147      REAL RX(1), RY(1), RZ(1)
148      POINTER( IPX, X )
149      POINTER( IPY, Y )
150      POINTER( IPZ, Z )
151      POINTER( IPRX, RX )
152      POINTER( IPRY, RY )
153      POINTER( IPRZ, RZ )
154C
155      DATA ISIZOLD/-1/
156      SAVE ISIZOLD,IPX,IPY,IPZ,IPRX,IPRY,IPRZ
157C
158C     Externals
159C
160      INTEGER NUMPTWE
161#ifdef POINTER_64
162      INTEGER*8 JMALLOC
163#else
164      INTEGER JMALLOC
165#endif
166C
167C     Statement functions
168C
169      REAL A, B
170      LOGICAL SOUTHOF
171      SOUTHOF(A,B) = ((A) - (B)).GT.-1E-4
172C
173C     -----------------------------------------------------------------|
174C     Section 1.  Initialise.
175C     -----------------------------------------------------------------|
176C
177  100 CONTINUE
178C
179      HGENGG = 0
180C
181      CALL JDEBUG()
182C
183C     -----------------------------------------------------------------|
184C     Section 2.  Calculate current grid latitudes and longitudes
185C     -----------------------------------------------------------------|
186C
187  200 CONTINUE
188C
189c     EMOS-199: adjusted for reduced_gg/octahedral
190c     GRID = 360.0 / REAL(KGNNEW*4)
191      GRID = 360.0 / REAL(KNPTS(KGNNEW))
192      IGRIDNI = NINT(GRID*JSCALE)
193C
194      NUMPTS = 0
195      NOPCNT = 0
196      NONS   = 0
197      NOWE   = 0
198      NO1NS  = 0
199C
200      DO LOOPO = 1, KGNNEW*2
201C
202C       Generate points inside the area only.
203C
204        IF( SOUTHOF(AREA(JNORTH),GNLATS(LOOPO)).AND.
205     X      SOUTHOF(GNLATS(LOOPO),AREA(JSOUTH)) ) THEN
206C
207          IF( NO1NS.EQ.0 ) NO1NS = LOOPO
208          NONS = NONS + 1
209C
210C         Grid step varies for "quasi-regular"/octahedral reduced Gaussian grids
211C
212          IF( HTYPE.EQ.'R' .OR. HTYPE.EQ.'r' .OR.
213     X        HTYPE.EQ.'O' .OR. HTYPE.EQ.'o' .OR.
214     X        HTYPE.EQ.'U' .OR. HTYPE.EQ.'u' ) THEN
215            GRID = 360.0 / REAL(KNPTS(LOOPO))
216            IGRIDNI = NINT(GRID*JSCALE)
217          ENDIF
218C
219          NLON = NUMPTWE(AREA(JWEST),AREA(JEAST),GRID)
220          NOWE = NLON
221          DO LOOPI = 1, NLON
222            NUMPTS = NUMPTS + 1
223            RLAT(NUMPTS) = GNLATS(LOOPO)
224            RLON(NUMPTS) = AREA(JWEST)+REAL((LOOPI-1)*IGRIDNI)/JSCALE
225          ENDDO
226        ENDIF
227      ENDDO
228C
229C     -----------------------------------------------------------------|
230C     Section 3.  Get some space for rotating the lat/longs
231C     -----------------------------------------------------------------|
232C
233  300 CONTINUE
234C
235      NOPCNT = NUMPTS
236      ISIZE  = NUMPTS
237      IF( ISIZE.GT.ISIZOLD ) THEN
238C
239        IF( ISIZOLD.GT.0 ) CALL JFREE(IPX)
240C
241        NBYTES = 6*ISIZE*JPRLEN
242C
243        IPX = JMALLOC(NBYTES)
244#ifdef hpR64
245        IPX = IPX/(1024*1024*1024*4)
246#endif
247        IF( IPX.EQ.0 ) THEN
248          CALL INTLOG(JP_WARN,'HGENGG: Memory allocate fail',JPQUIET)
249          HGENGG = 1
250          GOTO 900
251        ENDIF
252C
253        IPY  = IPX  + (ISIZE*JPRLEN)
254        IPZ  = IPY  + (ISIZE*JPRLEN)
255        IPRX = IPZ  + (ISIZE*JPRLEN)
256        IPRY = IPRX + (ISIZE*JPRLEN)
257        IPRZ = IPRY + (ISIZE*JPRLEN)
258C
259        ISIZOLD = ISIZE
260C
261      ENDIF
262C
263C     -----------------------------------------------------------------|
264C     Section 4.  Calculate the lat/longs before rotation
265C     -----------------------------------------------------------------|
266C
267  400 CONTINUE
268C
269C     Convert the rotated row points lat/longs to (x,y,z) coordinates
270C
271      CALL HLL2XYZ(RLAT,RLON,X,Y,Z,ISIZE)
272C
273C     Rotate the rotated row points back through the original latitude
274C     rotation
275C
276      CALL YROTATE(-(90.0+POLE(JLAT)),X,Y,Z,RX,RY,RZ,ISIZE)
277C
278C     Convert the rotated row points adjusted (x,y,z) coordinates to
279C     lat/long in the original grid (but still containing the
280C     longitude rotation)
281C
282      CALL XYZ2LL(RX,RY,RZ,RLAT,RLON,ISIZE)
283C
284C     Adjust the rotated line longitudes to remove the longitude
285C     rotation
286C
287      DO LOOP = 1, ISIZE
288        RLON(LOOP) = RLON(LOOP) + POLE(JLON)
289        IF( RLON(LOOP).LT.0.0 )   RLON(LOOP) = RLON(LOOP) + 360.0
290        IF( RLON(LOOP).GE.360.0 ) RLON(LOOP) = RLON(LOOP) - 360.0
291      ENDDO
292C
293C     -----------------------------------------------------------------|
294C     Section 9.  Return
295C     -----------------------------------------------------------------|
296C
297  900 CONTINUE
298C
299      RETURN
300      END
301