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