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