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 JSLLGP( PSHUP, KTRUNC, PNORTH, PSOUTH, PWEST, PEAST, 12 X PLATINC, PLONINC, POUTF, KRET) 13C 14C----> 15C**** JSLLGP 16C 17C PURPOSE 18C _______ 19C 20C This routine converts stretched spectral input fields to 21C standard lat/long grid fields. 22C 23C INTERFACE 24C _________ 25C 26C CALL JSLLGP( PSHUP, KTRUNC, PNORTH, PSOUTH, PWEST, PEAST, 27C X PLATINC, PLONINC, POUTF, KRET) 28C 29C Input parameters 30C ________________ 31C 32C PSHUP - Spherical harmonics field, unpacked 33C KTRUNC - Truncation number of spherical harmonics field 34C PNORTH - Northernmost latitude for output field (degrees) 35C PSOUTH - Southernmost latitude for output field (degrees) 36C PWEST - Westernmost longitude for output field (degrees) 37C PEAST - Easternmost longitude for output field (degrees) 38C PLATINC - Grid interval between latitudes in degrees 39C PLONINC - Grid interval between longitudes in degrees 40C INTLOG - Output log message 41C INTLOGR - Output log message (with real value) 42C 43C Output parameters 44C ________________ 45C 46C POUTF - Output grid point field 47C KRET - Return status code 48C 0 = OK 49C 50C Common block usage 51C __________________ 52C 53C JDCNDBG 54C 55C Method 56C ______ 57C 58C Generates one latitude at a time. 59C 60C Externals 61C _________ 62C 63C JJSET99 - Sets up FFT 64C JSTRLL - Generate latitude line from spherical harmonics 65C JSTRWLL - Generate latitude line from spherical harmonics (wind) 66C JDEBUG - Checks environment variable to switch on/off debug 67C INTLOG - Logs output messages 68C INTLOGR - Logs output messages 69C 70C 71C Reference 72C _________ 73C 74C E.C.M.W.F. Research Department technical memorandum no. 56 75C "The forecast and analysis post-processing package" 76C May 1982. J.Haseler. 77C 78C 79C Comments 80C ________ 81C 82C None 83C 84C 85C AUTHOR 86C ______ 87C 88C J.D.Chambers ECMWF July 1999 89C 90C MODIFICATIONS 91C _____________ 92C 93C None 94C 95C----< 96C _______________________________________________________ 97C 98C* Section 0. Definition of variables. 99C _______________________________________________________ 100C 101 IMPLICIT NONE 102#include "jparams.h" 103#include "parim.h" 104#include "nifld.common" 105C 106C Parameters 107C 108 INTEGER JPROUTINE 109 PARAMETER ( JPROUTINE = 30000 ) 110C 111C Subroutine arguments 112C 113 COMPLEX PSHUP 114 DIMENSION PSHUP(*) 115 INTEGER KTRUNC 116 REAL PNORTH, PSOUTH, PWEST, PEAST, PLATINC, PLONINC 117 REAL POUTF 118 DIMENSION POUTF(*) 119 INTEGER KRET 120C 121C Local variables 122C 123 REAL PSTART, PNORMX, PDIFF 124 REAL DEG2RAD, ZBUILD 125 INTEGER ILON, IOLON, ITOUP1, IJPWR, NEXT, NROWS, NUMLON 126 INTEGER IOFFN, LONGIT, N360, NSTART, INEXT, LOOPCL 127C 128 INTEGER IMFAX 129 REAL ATRIGS 130 DIMENSION IMFAX(10),ATRIGS(JPLONO*2) 131#ifndef _CRAYFTN 132#ifdef POINTER_64 133 INTEGER*8 IALEG, IZFA 134#endif 135#endif 136 REAL ALEG 137 DIMENSION ALEG(1) 138 POINTER ( IALEG, ALEG ) 139 REAL ZFA 140 DIMENSION ZFA(1) 141 POINTER ( IZFA, ZFA ) 142C 143 INTEGER ISIZE, ISIZE2 144 DATA ISIZE/0/, ISIZE2/0/ 145 SAVE ISIZE, ISIZE2, IZFA 146C 147 SAVE IMFAX, ATRIGS, IALEG 148C 149C _______________________________________________________ 150C 151C* Section 1. Initialization. 152C _______________________________________________________ 153C 154 100 CONTINUE 155C 156 CALL JDEBUG( ) 157C 158C First time through, dynamically allocate memory for workspace 159C 160 IF( ISIZE.EQ.0 ) THEN 161 ISIZE = (JPLONO+2) * 2 162 CALL JMEMHAN( 10, IZFA, ISIZE, 1, KRET) 163 IF( KRET.NE.0 ) THEN 164 CALL INTLOG(JP_ERROR,'JSLLGP: memory allocation error.',IZFA) 165 KRET = JPROUTINE + 1 166 GOTO 900 167 ENDIF 168 ENDIF 169C 170 IF( NDBG.GT.0 ) THEN 171 CALL INTLOG(JP_DEBUG,'JSLLGP: Input parameters:', JPQUIET) 172 CALL INTLOG(JP_DEBUG, 173 X 'JSLLGP: Spherical harmonic coeffs(first 20):', JPQUIET) 174 DO 101 NDBGLP = 1, 20 175 CALL INTLOGR(JP_DEBUG,' ', PSHUP( NDBGLP )) 176 101 CONTINUE 177 CALL INTLOG(JP_DEBUG, 178 X 'JSLLGP: Spherical harmonic truncation = ', KTRUNC) 179 CALL INTLOGR(JP_DEBUG, 180 X 'JSLLGP: Northern latitude for output = ', PNORTH) 181 CALL INTLOGR(JP_DEBUG, 182 X 'JSLLGP: Southern latitude for output = ', PSOUTH) 183 CALL INTLOGR(JP_DEBUG, 184 X 'JSLLGP: Western longitude for output = ', PWEST) 185 CALL INTLOGR(JP_DEBUG, 186 X 'JSLLGP: Eastern longitude for output = ', PEAST) 187 CALL INTLOGR(JP_DEBUG, 188 X 'JSLLGP: Latitude grid interval (deg) = ', PLATINC) 189 CALL INTLOGR(JP_DEBUG, 190 X 'JSLLGP: Longitude grid interval (deg) = ', PLONINC) 191 ENDIF 192C 193C Coefficients are to generated 'on the fly', 194C 195 ISIZE2 = (KTRUNC+1)*(KTRUNC+2) 196 CALL JMEMHAN( 14, IALEG, ISIZE2, 1, KRET) 197 IF( KRET.NE.0 ) THEN 198 CALL INTLOG(JP_ERROR,'JSLLGP: memory allocate error.',IALEG) 199 KRET = JPROUTINE + 2 200 GOTO 900 201 ENDIF 202C 203 DEG2RAD = PPI / 180.0 204 IOLON = NINT( 360.0 / PLATINC ) 205 ILON = NINT( 360.0 / PLONINC ) 206 ITOUP1 = KTRUNC 207C 208C For calculation purposes, ILON, the number of longitude points, 209C has to be greater than 2*(output truncation) to ensure that the 210C fourier transform is exact (see Reference, page 10). 211C The following loop increases ILON*IJPWR until ILON=>2*ITOUP1. 212C Later, when filling the output array POUTF, the longitude points 213C have to be taken selectively to avoid the intermediate generated 214C points, picking up values only at the required longitudes. 215C 216 IJPWR = 1 217 110 CONTINUE 218 IF( ILON.GT.2*ITOUP1 ) GOTO 120 219 IJPWR = IJPWR*2 220 ILON = ILON*2 221 GOTO 110 222 120 CONTINUE 223C 224C Initialise FFT constants 225C 226 IF( NDBG.GT.0 ) THEN 227 CALL INTLOG(JP_DEBUG, 228 X 'JSLLGP: Initialise FFT constants with ', ILON) 229 CALL INTLOG(JP_DEBUG, 230 X 'JSLLGP: longitude points.', JPQUIET) 231 ENDIF 232C 233 CALL JJSET99( ATRIGS, IMFAX, ILON, KRET) 234 IF( KRET.NE.0 ) GOTO 900 235C 236C _______________________________________________________ 237C 238C* Section 2. Work out geometry of the request 239C _______________________________________________________ 240C 241 200 CONTINUE 242C 243 PNORMX = PNORTH 244 PDIFF = PNORTH - PSOUTH 245 NROWS = NINT( PDIFF/PLATINC ) + 1 246 IF( ((NROWS-1)*PLATINC).GT.PDIFF ) NROWS = NROWS - 1 247C 248 NUMLON = NINT( (PEAST - PWEST)/PLONINC ) + 1 249 ZBUILD = PLATINC 250C 251 IF( NDBG.GT.0 ) THEN 252 CALL INTLOG(JP_DEBUG,'JSLLGP: Number of longitudes = ', NUMLON) 253 CALL INTLOGR(JP_DEBUG, 254 X 'JSLLGP: North latitude used for calc = ', PNORMX) 255 CALL INTLOG(JP_DEBUG,'JSLLGP: No.lats used in calc = ', NROWS) 256 ENDIF 257C _______________________________________________________ 258C 259C* Section 4. Process one row at a time. 260C _______________________________________________________ 261C 262 400 CONTINUE 263C 264 DO NEXT = 1, NROWS 265 PSTART = PNORMX - REAL(NEXT-1)*PLATINC 266C 267 IF( NDBG.GT.0 ) 268 X CALL INTLOGR(JP_DEBUG, 'JSLLGP: Next latitude = ', PSTART) 269C 270C Generate the next row 271C 272 IF( LWIND ) THEN 273 CALL JSTRWLL(PSHUP, KTRUNC, PSTART, ZBUILD, ILON, 274 X ALEG, ATRIGS, IMFAX, ZFA, KRET) 275 ELSE 276 CALL JSTRLL(PSHUP, KTRUNC, PSTART, ZBUILD, ILON, 277 X ALEG, ATRIGS, IMFAX, ZFA, KRET) 278 ENDIF 279 IF( KRET.NE.0 ) GOTO 900 280C 281C Find start offset in storage array 282C 283 IOFFN = NINT( (PNORTH - PSTART)/PLATINC ) + 1 284 LONGIT = NINT( (360.0 + PWEST)/PLONINC ) 285 N360 = NINT( 360.0/PLONINC ) 286 NSTART = MOD( LONGIT, N360)*IJPWR 287C 288 DO LOOPCL = 1, NUMLON 289 INEXT = 2 + MOD( NSTART + (LOOPCL-1)*IJPWR , ILON) 290 POUTF( LOOPCL + (IOFFN-1)*NUMLON ) = ZFA( INEXT) 291 ENDDO 292C 293 ENDDO 294C 295 KRET = 0 296C 297C _______________________________________________________ 298C 299C* Section 9. Return to calling routine. Format statements 300C _______________________________________________________ 301C 302 900 CONTINUE 303C 304 RETURN 305 END 306