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 SUBROUTINE JWINDLL( PSHUP, KTRUNC, PSTART, PBUILD, PINTVL, KLUNIT, 12 X KLATO,KLONO,PLEG,PTRIGS,KMFAX,PZFA,KRET) 13C 14C----> 15C**** JWINDLL 16C 17C Purpose 18C ------- 19C This routine converts spectral input fields to standard 20C lat/long grid fields for wind components U and V. 21C 22C 23C Interface 24C --------- 25C CALL JWINDLL( PSHUP, KTRUNC, PSTART, PBUILD, PINTVL, KLUNIT, 26C X KLATO,KLONO,PLEG,PTRIGS,KMFAX,PZFA,KRET) 27C 28C 29C Input parameters 30C ---------------- 31C PSHUP - Spherical harmonics field, unpacked 32C KTRUNC - Truncation number of spherical harmonics field 33C PSTART - Start latitude (northernmost) for output field 34C (must be positive - see comments below) 35C PBUILD - Grid interval used to build the legendre coefficients file 36C PINTVL - Grid latitude interval in degrees 37C KLUNIT - stream number of the legendre function file 38C KLATO - Number of latitude points in output field 39C KLONO - Number of longitude points in output field 40C PLEG - Array used to hold legendre functions 41C PTRIGS - Initialized array of trig.functions (setup by JJSET99) 42C KMFAX - Initialized array of prime factors (setup by JJSET99) 43C 44C 45C Output parameters 46C ----------------- 47C PZFA - Output grid point field; contains upto 32 each of 48C North and South latitude rows symmetrically. 49C KRET - Return status code 50C 0 = OK 51C 52C 53C Common block usage 54C ------------------ 55C JDCNDBG 56C 57C 58C Method 59C ------ 60C None. 61C 62C 63C Externals 64C --------- 65C JREADLL - Reads the legendre functions for a latitude 66C JSPPOLE - Applies correction at North or South pole 67C FFT99 - Carries out FFT 68C INTLOG - Output log message 69C INTLOGR - Output log message (with real value) 70C GETENV - Pick up contents of an environment variable 71C JUVPOLE - Use gaussian grid to calculate U and V at the poles 72C 73C 74C Reference 75C --------- 76C E.C.M.W.F. Research Department technical memorandum no. 56 77C "The forecast and analysis post-processing package" 78C May 1982. J.Haseler. 79C 80C 81C Comments 82C -------- 83C This is a redesign, based on SPECGP.F 84C 85C It handles transformation to a regular lat/long grid. 86C The generated grid is symmetrical about the equator, so 87C PSTART must be positive. 88C 89C It is only for U and V fields (correction is applied at the 90C poles and a scale factor is applied according to latitude). 91C 92C 93C Author 94C ------ 95C J.D.Chambers *ECMWF* Nov 1993 96C 97C 98C Modifications 99C ------------- 100C J.D.Chambers ECMWF Feb 1997 101C Allow for 64-bit pointers 102C 103C J.D.Chambers ECMWF October 2002 104C Add handling for U/V spectral (as opposed to Ucos(theta)/Vcos(theta)) 105C 106C 107C----< 108C -----------------------------------------------------------------| 109C 110C 111 IMPLICIT NONE 112#include "jparams.h" 113#include "parim.h" 114#include "nifld.common" 115C 116C Subroutine arguments 117C 118 COMPLEX PSHUP 119 DIMENSION PSHUP(*) 120 INTEGER KTRUNC 121 REAL PSTART, PBUILD, PINTVL 122 INTEGER KLUNIT, KLATO, KLONO, KMFAX, KRET 123 REAL PLEG, PTRIGS, PZFA 124 DIMENSION PZFA(JPLONO + 2, 64) 125 DIMENSION KMFAX(*), PLEG(*), PTRIGS(*) 126C 127C Parameters 128C 129 INTEGER JPROUTINE 130 PARAMETER ( JPROUTINE = 31400 ) 131C 132C Local variables 133C 134 REAL ZPIBY2, ZDEGR, ZLAT 135 INTEGER ILIM, IMLIM, ILN 136 INTEGER ITAL, ITALA, ITALS, IMN, IMP 137 INTEGER INORTH, ISOUTH 138 INTEGER JM, J245, JNEXTLAT 139 INTEGER NERR 140 INTEGER*8 IOFF 141 INTEGER*8 JDCLOOP 142 REAL ZF 143 REAL ZCOSI 144C 145 DIMENSION ZF(2*JPFFT) 146C 147#ifndef _CRAYFTN 148#ifdef POINTER_64 149 INTEGER*8 IWORK 150#endif 151#endif 152 REAL WORK 153 DIMENSION WORK(1) 154 POINTER ( IWORK, WORK ) 155 COMPLEX ZDUM(JPTRNC + 1) 156 COMPLEX ZSUMS(JPTRNC + 1), ZSUMA(JPTRNC + 1) 157 COMPLEX*16 CHOLD 158 INTEGER LOOP 159C 160 LOGICAL LXFIRST, LPLAINU 161 INTEGER ISIZE, IBLANK 162 CHARACTER*10 PLAINUV 163C 164 DATA LXFIRST/.TRUE./, LPLAINU/.FALSE./ 165 DATA ISIZE/0/ 166 SAVE ISIZE, IWORK, LXFIRST, LPLAINU 167C 168C Statement function 169C 170 LOGICAL GPOLE 171 REAL ANGLE 172 GPOLE(ANGLE) = ( ABS(90.0 - ABS(ANGLE) ) .LT. 1.0E-3 ) 173C = .TRUE. if LAT is 90.0 or -90.0 174C 175C -----------------------------------------------------------------| 176C* Section 1. Initialization. 177C -----------------------------------------------------------------| 178C 179 100 CONTINUE 180 IOFF = 0 181C 182C Check environment variable to see if wind components are plain 183C U and V, or U*cos(theta) and V*cos(theta) (the default). 184C 185 IF( LXFIRST ) THEN 186C 187 LXFIRST = .FALSE. 188 CALL GETENV('INTERP_PLAIN_UV', PLAINUV) 189 IBLANK = INDEX(PLAINUV, ' ') 190C 191 IF( (IBLANK.EQ.3).AND.(PLAINUV(1:2).EQ.'ON') ) LPLAINU =.TRUE. 192 IF( (IBLANK.EQ.3).AND.(PLAINUV(1:2).EQ.'on') ) LPLAINU =.TRUE. 193 IF( (IBLANK.EQ.4).AND.(PLAINUV(1:3).EQ.'YES') ) LPLAINU =.TRUE. 194 IF( (IBLANK.EQ.4).AND.(PLAINUV(1:3).EQ.'yes') ) LPLAINU =.TRUE. 195C 196 IF( (IBLANK.EQ.3).AND.(PLAINUV(1:2).EQ.'NO') ) LPLAINU =.FALSE. 197 IF( (IBLANK.EQ.3).AND.(PLAINUV(1:2).EQ.'no') ) LPLAINU =.FALSE. 198 IF( (IBLANK.EQ.4).AND.(PLAINUV(1:3).EQ.'OFF') ) LPLAINU =.FALSE. 199 IF( (IBLANK.EQ.4).AND.(PLAINUV(1:3).EQ.'off') ) LPLAINU =.FALSE. 200C 201C First time through, dynamically allocate memory for workspace 202C 203 ISIZE = 2*JPFFT*64 204 CALL JMEMHAN( 9, IWORK, ISIZE, 1, KRET) 205 IF( KRET.NE.0 ) THEN 206 CALL INTLOG(JP_ERROR,'JWINDLL: memory allocation error',IWORK) 207 KRET = JPROUTINE + 1 208 GOTO 990 209 ENDIF 210 ENDIF 211C 212 IF( NDBG.GT.1 ) THEN 213 CALL INTLOG(JP_DEBUG, 214 X 'JWINDLL: Spherical harmonic coeffs(first 20):',JPQUIET) 215 DO NDBGLP = 1, 20 216 CALL INTLOGR(JP_DEBUG,' ',PSHUP( NDBGLP )) 217 ENDDO 218 CALL INTLOG(JP_DEBUG, 219 X 'JWINDLL: Input parameters:',JPQUIET) 220 CALL INTLOG(JP_DEBUG, 221 X 'JWINDLL: Spherical harmonic truncation = ', KTRUNC) 222 CALL INTLOGR(JP_DEBUG, 223 X 'JWINDLL: Start latitude(northernmost) = ', PSTART) 224 CALL INTLOGR(JP_DEBUG, 225 X 'JWINDLL: Grid lat. interval(leg. file) = ', PBUILD) 226 CALL INTLOGR(JP_DEBUG, 227 X 'JWINDLL: Grid lat. interval(degrees) = ', PINTVL) 228 CALL INTLOG(JP_DEBUG, 229 X 'JWINDLL: Stream number of leg. file = ', KLUNIT) 230 CALL INTLOG(JP_DEBUG, 231 X 'JWINDLL: Number of lat. rows in output = ', KLATO) 232 CALL INTLOG(JP_DEBUG, 233 X 'JWINDLL: Number of long. pts per row = ', KLONO) 234 CALL INTLOG(JP_DEBUG, 235 X 'JWINDLL: Trig.functions (setup by JJSET99):',JPQUIET) 236 DO NDBGLP = 1, 10 237 CALL INTLOGR(JP_DEBUG,' ',PTRIGS( NDBGLP )) 238 ENDDO 239 CALL INTLOG(JP_DEBUG, 240 X 'JWINDLL: Prime factors (setup by JJSET99):',JPQUIET) 241 DO NDBGLP = 1, 10 242 CALL INTLOG(JP_DEBUG,' ',KMFAX( NDBGLP )) 243 ENDDO 244 ENDIF 245C 246 ZPIBY2 = PPI / 2.0 247 ZDEGR = PPI / 180.0 248 ILIM = KTRUNC + 1 249 IMLIM = KTRUNC + 1 250 INORTH = -1 251 ILN = KLONO + 2 252C 253C -----------------------------------------------------------------| 254C* Section 2. Main loop through latitude rows to 255C* calculate fourier coefficients 256C -----------------------------------------------------------------| 257C 258 200 CONTINUE 259C 260C For each latitude, the north and corresponding south latitude row 261C are calculated at the same time from the same legendre functions. 262C 263 DO 280 JNEXTLAT = 1, KLATO 264 ZLAT = ( PSTART - (PINTVL * REAL(JNEXTLAT - 1)) ) 265C 266 IF( NDBG.GT.1 ) 267 X CALL INTLOGR(JP_DEBUG, 'JWINDLL: Next latitude = ', ZLAT) 268C 269C If required, generate the coefficients 'on the fly' 270C 271 IF( LON_FLY ) THEN 272 CALL NMAKLL( KTRUNC, PBUILD, ZLAT, 1, PLEG, NERR) 273 IOFF = 0 274 ELSE IF( LFILEIO ) THEN 275 CALL JREADLL( KLUNIT, KTRUNC, PBUILD, ZLAT, PLEG, NERR) 276 IF ( NERR .NE. 0 ) THEN 277 CALL INTLOG(JP_ERROR,'JWINDLL: JREADLL error',NERR) 278 KRET = JPROUTINE + 2 279 GOTO 990 280 ENDIF 281 ELSE 282 IOFF = NINT( (90.0 - ZLAT)/PBUILD ) 283 IOFF = IOFF *(KTRUNC+1) 284 IOFF = IOFF *(KTRUNC+4) 285 IOFF = IOFF/2 286 ENDIF 287C 288C Clear unused slots in array. 289C Note there are two slots in the array - one for north latitude 290C and one for the corresponding south latitude. 291 INORTH = INORTH + 2 292 ISOUTH = INORTH + 1 293 DO LOOP = 1, ILN 294 ZF( LOOP) = 0.0 295 PZFA( LOOP, INORTH) = 0.0 296 PZFA( LOOP, ISOUTH) = 0.0 297 ENDDO 298C 299C For spectral Ucos(theta), Vcos(theta), treat the poles as a 300C special case. 301C 302 IF( GPOLE(ZLAT).AND.(.NOT.LPLAINU) ) THEN 303C 304 IF( NDBG.GT.1 ) CALL INTLOG(JP_DEBUG, 305 X 'JWINDLL: Pole special for spectral Ucos(t), Vcos(t)', 306 X JPQUIET) 307C 308 CALL JSPPOLE( PSHUP, 1, KTRUNC, .TRUE., ZF) 309 DO LOOP = 1, ILN 310 PZFA( LOOP, INORTH) = ZF( LOOP) 311 ENDDO 312C 313 CALL JSPPOLE( PSHUP, 0, KTRUNC, .TRUE., ZF) 314 DO LOOP = 1, ILN 315 PZFA( LOOP, ISOUTH) = ZF( LOOP) 316 ENDDO 317C 318 ELSE 319C 320C Otherwise not at pole, or spectral U and V 321C 322C Fill slots which are used 323 IMN = 0 324 IMP = 0 325C 326 DO 247 JM = 1, IMLIM 327 ITAL = ILIM - JM + 1 328 DO 245 J245 = 1, ITAL 329 IF( LFILEIO ) THEN 330 ZDUM(J245) = PLEG(IMP + J245)*PSHUP(IMN + J245) 331 ELSE 332 JDCLOOP = IMP + J245 + IOFF 333 ZDUM(J245) = PLEG(JDCLOOP)*PSHUP(IMN + J245) 334 ENDIF 335 245 CONTINUE 336 IMP = IMP + ITAL + 1 337 IMN = IMN + ITAL 338 ITALS = (ITAL + 1)/2 339 ITALA = ITAL/2 340 CHOLD = (0.0D0, 0.0D0) 341 DO LOOP = 1, 2*ITALS, 2 342 CHOLD = CHOLD + ZDUM(LOOP) 343 ENDDO 344 ZSUMS(JM) = CHOLD 345 CHOLD = (0.0D0, 0.0D0) 346 DO LOOP = 2, 2*ITALA, 2 347 CHOLD = CHOLD + ZDUM(LOOP) 348 ENDDO 349 ZSUMA(JM) = CHOLD 350 247 CONTINUE 351C 352C For the southern hemisphere row, the legendre functions are 353C the complex conjugates of the corresponding northern row - 354C hence the juggling with the signs in the next loop. 355C 356C Note that PZFA is REAL, but the coefficients being calculated 357C are COMPLEX. There are pairs of values for each coefficient 358C (real and imaginary parts) and pairs of values for each 359C latitude (north and south). 360C 361 DO JM = 1, IMLIM 362 PZFA(2*JM -1, INORTH) = REAL(ZSUMS(JM)) + REAL(ZSUMA(JM)) 363 PZFA(2*JM , INORTH) = AIMAG(ZSUMS(JM)) + AIMAG(ZSUMA(JM)) 364 PZFA(2*JM -1, ISOUTH) = REAL(ZSUMS(JM)) - REAL(ZSUMA(JM)) 365 PZFA(2*JM , ISOUTH) = AIMAG(ZSUMS(JM)) - AIMAG(ZSUMA(JM)) 366 ENDDO 367 ENDIF 368C 369C* End of main loop through latitude rows. 370C 371 280 CONTINUE 372C 373C -----------------------------------------------------------------| 374C* Section 3. Fast fourier transform 375C -----------------------------------------------------------------| 376C 377 300 CONTINUE 378C 379 IF( NDBG.GT.1 ) CALL INTLOG(JP_DEBUG, 380 X 'JWINDLL: FFT, no.of rows (N and S) = ', ISOUTH) 381C 382 CALL FFT99(PZFA,WORK,PTRIGS,KMFAX,1,JPLONO+2,KLONO,ISOUTH,1) 383C 384 IF( NDBG.GT.1 ) THEN 385 CALL INTLOG(JP_DEBUG, 386 X 'JWINDLL: Values calculated by FFT:',JPQUIET) 387 DO NDBGLP = 1, 20 388 CALL INTLOGR(JP_DEBUG,' ',PZFA( 1, NDBGLP )) 389 CALL INTLOGR(JP_DEBUG,' ',PZFA( 2, NDBGLP )) 390 ENDDO 391 ENDIF 392C 393C For spectral U, V generate the poles from a reduced gaussian line 394C of latitude. 395C 396 IF( GPOLE(PSTART).AND.(LPLAINU) ) THEN 397C 398 IF( NDBG.GT.1 ) CALL INTLOG(JP_DEBUG, 399 X 'JWINDLL: Poles generated from gaussian for spectral U, V', 400 X JPQUIET) 401C 402 CALL JUVPOLE(PSHUP, KTRUNC, PZFA, KLONO, KRET) 403 IF( KRET.NE.0 ) THEN 404 CALL INTLOG(JP_ERROR,'JWINDLL: pole wind create error',IWORK) 405 KRET = JPROUTINE + 1 406 GOTO 990 407 ENDIF 408C 409 ENDIF 410C 411C -----------------------------------------------------------------| 412C* Section 4. For spectral Ucos(theta), Vcos(theta), apply scale 413C factor to all latitudes except poles. 414C -----------------------------------------------------------------| 415C 416 400 CONTINUE 417C 418 IF( NDBG.GT.0 ) THEN 419 IF( LPLAINU ) THEN 420 CALL INTLOG(JP_DEBUG, 421 X 'JWINDLL: Do not apply scale to latitudes',JPQUIET) 422 ELSE 423 CALL INTLOG(JP_DEBUG, 424 X 'JWINDLL: Apply scale to latitudes',JPQUIET) 425 ENDIF 426 ENDIF 427C 428 IF( .NOT.LPLAINU ) THEN 429 INORTH = -1 430 DO JNEXTLAT = 1, KLATO 431 ZLAT = ( PSTART - (PINTVL * REAL(JNEXTLAT - 1)) ) 432 INORTH = INORTH + 2 433 ISOUTH = INORTH + 1 434 IF( .NOT.GPOLE(ZLAT) ) THEN 435 ZCOSI = 1.0 / COS( ZLAT * ZDEGR ) 436 CALL EMOSLIB_SSCAL( KLONO, ZCOSI, PZFA( 2, INORTH), 1) 437 CALL EMOSLIB_SSCAL( KLONO, ZCOSI, PZFA( 2, ISOUTH), 1) 438 ENDIF 439 ENDDO 440C 441 ENDIF 442C -----------------------------------------------------------------| 443C* Section 9. Return to calling routine. Format statements 444C -----------------------------------------------------------------| 445C 446 900 CONTINUE 447C 448 KRET = 0 449C 450 990 CONTINUE 451 RETURN 452 END 453