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