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 HNTFAUS(FLDIN,INLEN) 12C 13C----> 14C**** HNTFAUS 15C 16C Purpose 17C ------- 18C 19C Prepare to interpolate input field... 20C 21C 22C Interface 23C --------- 24C 25C IRET = HNTFAUS( FLDIN,INLEN) 26C 27C Input 28C ----- 29C 30C FLDIN - Input field (unpacked). 31C INLEN - Input field length (words). 32C 33C 34C Output 35C ------ 36C 37C Field unpacked values are in ZNFELDI, rotated if necessary. 38C 39C Returns: 0, if OK. Otherwise, an error occured in rotation. 40C 41C 42C Method 43C ------ 44C 45C Rotate the field values if necessary. 46C 47C Some rotation coefficients cannot be created because the 48C algorithm does not converge for rotations which are smaller 49C than a limit which depends on the spectral truncation. 50C 51C If the rotation is too small, it has to be split into three 52C steps (two forward and one backward). 53C 54C Current know limits are: 55C 56C Truncation Max South pole lat Minimum rotation 57C ---------- ------------------ ---------------- 58C 63 -89.5 0.5 59C 106 -85 5.0 60C 159 -77 13.0 61C 213 -68 22.0 62C 319 -51 39.0 63C 511 -30.5 59.5 64C 65C 66C Externals 67C --------- 68C 69C INTLOG - Log error message. 70C INTLOGR - Log error message. 71C RPHI - Rotate a spectral field by longitude. 72C JACOBI - Rotates spectral coefficients by latitude. 73C 74C 75C Author 76C ------ 77C 78C J.D.Chambers ECMWF February 2001 79C 80C----< 81C 82C -----------------------------------------------------------------| 83C* Section 0. Variables 84C -----------------------------------------------------------------| 85C 86 IMPLICIT NONE 87C 88C Function arguments 89C 90 INTEGER INLEN 91 REAL FLDIN(*) 92C 93#include "parim.h" 94#include "nifld.common" 95#include "nofld.common" 96#include "grfixed.h" 97#include "intf.h" 98C 99C Parameters 100C 101 INTEGER JPROUTINE 102 PARAMETER (JPROUTINE = 40140 ) 103 INTEGER JPLEN, JPNM 104 PARAMETER ( JPNM = JPSTRUNC ) 105 PARAMETER ( JPLEN = (JPNM+1)*(JPNM+6) ) 106C 107C Local variables 108C 109 INTEGER IRET 110 INTEGER ISIZE 111 INTEGER ITRUNC 112 INTEGER LOOP 113 INTEGER NBYTES 114 INTEGER NUMPTS 115 LOGICAL LOK 116 LOGICAL LSPLIT 117 REAL POLELAT 118 REAL SWORK(1) 119 REAL*8 DATA(1) 120 REAL*8 DLAT 121 REAL*8 DLON 122 REAL*8 DWORK(1) 123 REAL*8 LATLIM 124#ifndef _CRAYFTN 125#ifdef POINTER_64 126 INTEGER*8 IDATA 127 INTEGER*8 IDWORK 128 INTEGER*8 ISWORK 129#endif 130#endif 131 POINTER (IDATA, DATA) 132 POINTER (IDWORK, DWORK) 133 POINTER (ISWORK, SWORK) 134C 135 SAVE IDATA 136 SAVE IDWORK 137 SAVE ISIZE 138 SAVE ISWORK 139C 140 DATA ISIZE/0/ 141 DATA NBYTES/4/ 142C 143C Externals 144 LOGICAL JACOBI 145 EXTERNAL JACOBI 146C 147C -----------------------------------------------------------------| 148C* Section 1. Initialise 149C -----------------------------------------------------------------| 150C 151 100 CONTINUE 152C 153 HNTFAUS = 0 154 IRET = 0 155C 156C -----------------------------------------------------------------| 157C* Section 2. Store input data. 158C -----------------------------------------------------------------| 159C 160 200 CONTINUE 161C 162C Allocate work array ZNFELDI if not already done. 163C 164 IF( IZNJDCI.NE.1952999238 ) THEN 165 CALL JMEMHAN( 19, IZNFLDI, JPEXPAND, 1, IRET) 166 IF( IRET.NE.0 ) THEN 167 CALL INTLOG(JP_WARN,'HNTFAUS: ZNFELDI allocate fail',JPQUIET) 168 HNTFAUS = IRET 169 GOTO 900 170 ENDIF 171 IZNJDCI = 1952999238 172 ENDIF 173C 174C Move unpacked values in from user array 175C 176 DO LOOP = 1, INLEN 177 ZNFELDI( LOOP ) = FLDIN( LOOP ) 178 ENDDO 179C 180C -----------------------------------------------------------------| 181C* Section 3. Handle rotation, if necessary. 182C -----------------------------------------------------------------| 183C 184 300 CONTINUE 185C 186C Is a rotation of the SH required? 187C 188 IF( LNOROTA ) THEN 189C 190C Avoid rotating U or V spectral coefficients 191C 192 IF( (NIPARAM.EQ.JP_U).OR.(NIPARAM.EQ.JP_V) ) THEN 193 CALL INTLOG(JP_ERROR, 194 X 'HNTFAUS: Unable to rotate spectral U or V:',NIPARAM) 195 HNTFAUS = JPROUTINE + 3 196 GOTO 900 197 ENDIF 198C 199C Dynamically allocate memory for workspace. 200C 201C (Note the multiplication by 2 for REAL*8). 202C 203 NUMPTS = INLEN 204 ISIZE = NUMPTS*2 205 CALL JMEMHAN( 11, IDWORK, ISIZE, 1, IRET) 206 IF( IRET.NE.0 ) THEN 207 CALL INTLOG(JP_ERROR,'HNTFAUS: memory allocation fail',IDWORK) 208 HNTFAUS = JPROUTINE + 3 209 GOTO 900 210 ENDIF 211C 212C (Note the multiplication by 2 for REAL*8). 213C 214 ISIZE = JPLEN*2 215 CALL JMEMHAN( 15, IDATA, ISIZE, 1, IRET) 216 IF( IRET.NE.0 ) THEN 217 CALL INTLOG(JP_ERROR,'HNTFAUS: memory allocation fail',IDATA) 218 HNTFAUS = JPROUTINE + 3 219 GOTO 900 220 ENDIF 221C 222 CALL INTLOG(JP_DEBUG,'HNTFAUS: Rotate SH field.',JPQUIET) 223 CALL INTLOG(JP_DEBUG,'HNTFAUS: South pole latitude ',NOROTA(1)) 224 CALL INTLOG(JP_DEBUG,'HNTFAUS: South pole longitude ',NOROTA(2)) 225 IF ( NOREPR.EQ.JPSPHERE ) NOREPR = JPSPHROT 226 IF ( NOREPR.EQ.JPREGULAR ) NOREPR = JPREGROT 227 ITRUNC = NIRESO 228C 229C Some rotation coefficients cannot be created because the 230C algorithm does not converge for rotations which are smaller 231C than a limit which depends on the spectral truncation. 232C 233 IF( ITRUNC.LE.63 ) THEN 234 LATLIM = -89.5 235 ELSE IF( ITRUNC.LE.106 ) THEN 236 LATLIM = -85.0 237 ELSE IF( ITRUNC.LE.159 ) THEN 238 LATLIM = -77.0 239 ELSE IF( ITRUNC.LE.213 ) THEN 240 LATLIM = -68.0 241 ELSE IF( ITRUNC.LE.319 ) THEN 242 LATLIM = -51.0 243 ELSE IF( ITRUNC.LE.511 ) THEN 244 LATLIM = -29.5 245 ELSE 246 HNTFAUS = JPROUTINE + 3 247 GOTO 900 248 ENDIF 249 CALL INTLOGR(JP_DEBUG, 250 X 'HNTFAUS: South pole latitude limit = ',LATLIM) 251C 252 POLELAT = (REAL(NOROTA(1))/JPMULT) 253C 254C If the rotation is too small, it has to be split into three 255C steps (two forward and one backward). 256C 257 LSPLIT = POLELAT.LT.LATLIM 258 IF( LSPLIT ) THEN 259 DLAT = -90.0 - LATLIM 260 CALL INTLOG(JP_DEBUG, 261 X 'HNTFAUS: Rotation has been split into three steps',JPQUIET) 262 ELSE 263 DLAT = -90.0 - POLELAT 264 ENDIF 265#if (defined REAL_8) 266C 267C Double precision REALs 268C 269 DLON = - REAL(NOROTA(2))/JPMULT 270C 271C Rotate the spectral field by longitude. 272 CALL RPHI( ZNFELDI, ITRUNC, DWORK, DLON) 273C 274C Rotate the spectral field by latitude. 275C 276 IF( LSPLIT ) THEN 277C 278C Two rotations forward .. 279C 280 CALL INTLOGR(JP_DEBUG, 281 X 'HNTFAUS: Twice forward through angle ', DLAT) 282 LOK = JACOBI(ZNFELDI, ITRUNC, DWORK, DLAT) 283 IF(.NOT.LOK) THEN 284 CALL INTLOG(JP_ERROR, 285 X 'HNTFAUS: Spectral rotation failed',JPQUIET) 286 HNTFAUS = JPROUTINE + 3 287 GOTO 900 288 ENDIF 289 LOK = JACOBI(ZNFELDI, ITRUNC, DWORK, DLAT) 290 IF(.NOT.LOK) THEN 291 CALL INTLOG(JP_ERROR, 292 X 'HNTFAUS: Spectral rotation failed',JPQUIET) 293 HNTFAUS = JPROUTINE + 3 294 GOTO 900 295 ENDIF 296C 297C .. and one backward 298C 299 DLAT = -90.0 - ( (LATLIM*2.0) - POLELAT) 300 DLAT= -DLAT 301 CALL INTLOGR(JP_DEBUG, 302 X 'HNTFAUS: Once backward through angle ', DLAT) 303 ENDIF 304C 305 LOK = JACOBI( ZNFELDI, ITRUNC, DWORK, DLAT) 306 IF(.NOT.LOK) THEN 307 CALL INTLOG(JP_ERROR, 308 X 'HNTFAUS: Spectral rotation failed',JPQUIET) 309 HNTFAUS = JPROUTINE + 3 310 GOTO 900 311 ENDIF 312#else 313C 314C Single precision REALs; so convert to double precision. 315C 316 DLAT = -90.0 - DBLE(NOROTA(1))/JPMULT 317 DLON = - DBLE(NOROTA(2))/JPMULT 318C 319C Expand spectral coefficients to REAL*8 320C 321 DO LOOP = 1, NUMPTS 322 DATA(LOOP) = DBLE(ZNFELDI(LOOP)) 323 ENDDO 324C 325C Rotate the spectral field by longitude. 326C 327 CALL RPHI( DATA, ITRUNC, DWORK, DLON) 328C 329C Rotate the spectral field by latitude. 330C 331 IF( LSPLIT ) THEN 332C 333C Two rotations forward .. 334C 335 CALL INTLOGR(JP_DEBUG, 336 X 'HNTFAUS: Twice forward through angle ', DLAT) 337 LOK = JACOBI(ZNFELDI, ITRUNC, DWORK, DLAT) 338 IF(.NOT.LOK) THEN 339 CALL INTLOG(JP_ERROR, 340 X 'HNTFAUS: Spectral rotation failed',JPQUIET) 341 HNTFAUS = JPROUTINE + 3 342 GOTO 900 343 ENDIF 344 LOK = JACOBI(ZNFELDI, ITRUNC, DWORK, DLAT) 345 IF(.NOT.LOK) THEN 346 CALL INTLOG(JP_ERROR, 347 X 'HNTFAUS: Spectral rotation failed',JPQUIET) 348 HNTFAUS = JPROUTINE + 3 349 GOTO 900 350 ENDIF 351C 352C .. and one backward 353C 354 DLAT = -90.0 - ( (LATLIM*2.0) - POLELAT) 355 DLAT= -DLAT 356 CALL INTLOGR(JP_DEBUG, 357 X 'HNTFAUS: Once backward through angle ', DLAT) 358 ENDIF 359C 360 LOK = JACOBI( DATA, ITRUNC, DWORK, DLAT) 361#endif 362C 363 IF( .NOT.LOK ) THEN 364 CALL INTLOG(JP_ERROR,'HNTFAUS: Rotation failed',JPQUIET) 365 HNTFAUS = JPROUTINE + 3 366 GOTO 900 367 ENDIF 368#if (!defined REAL_8) 369C 370C Single precision REALs; repack spectral coefficients to REAL*4. 371C 372 DO LOOP = 1, NUMPTS 373 ZNFELDI(LOOP) = SNGL(DATA(LOOP)) 374 ENDDO 375#endif 376C 377 ENDIF 378C 379C -----------------------------------------------------------------| 380C* Section 9. Closedown. 381C -----------------------------------------------------------------| 382C 383 900 CONTINUE 384C 385 RETURN 386 END 387