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