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