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