Project

General

Profile

transforms4py.F90

Update Gauss spectral transforms for getting derivatives (v1.1.3) - alexandre mary, 09/21/2016 05:29 PM

 
1
SUBROUTINE W_SPEC_SETUP(KRETURNCODE, KSIZEI, KSIZEJ, KPHYSICALSIZEI, KPHYSICALSIZEJ, &
2
                       &KTRUNCX, KTRUNCY, KNUMMAXRESOL, KLOEN, LDLAM, &
3
                       &KSIZEKLOEN, PDELTAX, PDELTAY, &
4
                       &KIDENTRESOL, LDSTOP)
5
! ** PURPOSE
6
!    Setup spectral transform for LAM and global
7
!
8
! ** DUMMY ARGUMENTS
9
!    KRETURNCODE: error code
10
!    KSIZEI, KSIZEJ: size of grid-point field (with extension zone for LAM), put max size for KSIZEI in global
11
!    KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field for LAM (put 0 for global)
12
!    KTRUNCX, KTRUNCY: troncatures for LAM (only KTRUNCX is used for global, put 0 for KTRUNCY)
13
!    KNUMMAXRESOL: maximum number of troncatures handled
14
!    KLOEN: number of points on each latitude row
15
!    KSIZEKLOEN: size of KLOEN array
16
!    PDELTAX: x resolution
17
!    PDELTAY: y resolution
18
!    LDLAM: LAM (.TRUE.) or global (.FALSE.)
19
!    KIDENTRESOL: identification of resolution
20
!    LDSTOP: exception raised?
21
!
22
! ** AUTHOR
23
!    9 April 2014, S. Riette
24
!
25
! ** MODIFICATIONS
26
!    6 Jan 2016, S. Riette: PDELTAX and PDELTAY added
27
!
28
! I. Dummy arguments declaration
29
IMPLICIT NONE
30
INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
31
INTEGER, INTENT(IN) :: KSIZEI, KSIZEJ
32
INTEGER, INTENT(IN) :: KPHYSICALSIZEI, KPHYSICALSIZEJ
33
INTEGER, INTENT(IN) :: KTRUNCX, KTRUNCY
34
INTEGER, INTENT(IN) :: KNUMMAXRESOL
35
INTEGER, DIMENSION(KSIZEKLOEN), INTENT(IN) :: KLOEN
36
LOGICAL, INTENT(IN) :: LDLAM
37
INTEGER, INTENT(IN) :: KSIZEKLOEN
38
REAL(KIND=8), INTENT(IN) :: PDELTAX
39
REAL(KIND=8), INTENT(IN) :: PDELTAY
40
INTEGER, INTENT(OUT) :: KIDENTRESOL
41
LOGICAL, INTENT(OUT) :: LDSTOP
42
!
43
! II. Local variables declaration
44
INTEGER, DIMENSION(2*KSIZEJ) :: ILOEN
45
INTEGER :: JI
46
LOGICAL, SAVE :: LLFIRSTCALL=.TRUE.
47
LOGICAL :: LLNEWRESOL
48
INTEGER, SAVE :: INBRESOL=0
49
INTEGER(KIND=8) :: ICODEILOEN
50
INTEGER, SAVE :: INUMMAXRESOLSAVE=-1
51
INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ITRUNCXSAVE, ITRUNCYSAVE, &
52
                                            IPHYSICALSIZEISAVE, &
53
                                            IPHYSICALSIZEJSAVE, &
54
                                            ISIZEISAVE, ISIZEJSAVE, &
55
                                            IIDENTRESOLSAVE
56
INTEGER(KIND=8), DIMENSION(:), ALLOCATABLE, SAVE :: ILOENSAVE
57
REAL(KIND=8), DIMENSION(:), ALLOCATABLE, SAVE :: ZDELTAXSAVE, &
58
                                                 ZDELTAYSAVE
59
REAL(KIND=8) :: ZEXWN, ZEYWN
60

    
61
#include "setup_trans0.h"
62
#include "esetup_trans.h"
63
#include "setup_trans.h"
64

    
65
KRETURNCODE=0
66
LDSTOP=.FALSE.
67
! III. Setup
68

    
69
! III.a Setup LAM and global spectral transform - all resolutions
70
! Maximum number of resolution is set now and cannot be change anymore
71
IF (LLFIRSTCALL) THEN
72
  !This code is called only once, whatever is the number of resolutions
73
  CALL SETUP_TRANS0(KPRINTLEV=0, LDMPOFF=.TRUE., KMAX_RESOL=KNUMMAXRESOL)
74
  ALLOCATE(ITRUNCXSAVE(KNUMMAXRESOL))
75
  ALLOCATE(ITRUNCYSAVE(KNUMMAXRESOL))
76
  ALLOCATE(IPHYSICALSIZEISAVE(KNUMMAXRESOL))
77
  ALLOCATE(IPHYSICALSIZEJSAVE(KNUMMAXRESOL))
78
  ALLOCATE(ISIZEJSAVE(KNUMMAXRESOL))
79
  ALLOCATE(ISIZEISAVE(KNUMMAXRESOL))
80
  ALLOCATE(ILOENSAVE(KNUMMAXRESOL))
81
  ALLOCATE(IIDENTRESOLSAVE(KNUMMAXRESOL))
82
  ALLOCATE(ZDELTAXSAVE(KNUMMAXRESOL))
83
  ALLOCATE(ZDELTAYSAVE(KNUMMAXRESOL))
84
  ITRUNCXSAVE=-1
85
  ITRUNCYSAVE=-1
86
  IPHYSICALSIZEISAVE=-1
87
  IPHYSICALSIZEJSAVE=-1
88
  ISIZEJSAVE=-1
89
  ISIZEISAVE=-1
90
  ILOENSAVE=-1
91
  IIDENTRESOLSAVE=-1
92
  ZDELTAXSAVE=-1.
93
  ZDELTAXSAVE=-1.
94
  LLFIRSTCALL=.FALSE.
95
  INUMMAXRESOLSAVE=KNUMMAXRESOL
96
ENDIF
97
!
98
! III.b Is-it a new resolution?
99
LLNEWRESOL=.TRUE.
100
IF(LDLAM) THEN
101
  ILOEN(:)=KSIZEI
102
ELSE
103
  ILOEN(:)=0
104
  ILOEN(1:MIN(SIZE(ILOEN),SIZE(KLOEN)))=KLOEN(1:MIN(SIZE(ILOEN),SIZE(KLOEN)))
105
ENDIF
106
ICODEILOEN=0
107
DO JI=1, SIZE(ILOEN)
108
  ICODEILOEN=ICODEILOEN+ILOEN(JI)*JI**4
109
ENDDO
110
DO JI=1, INBRESOL
111
  IF (KTRUNCX==ITRUNCXSAVE(JI) .AND. KTRUNCY==ITRUNCYSAVE(JI) .AND. &
112
   &KPHYSICALSIZEI==IPHYSICALSIZEISAVE(JI) .AND. &
113
   &KPHYSICALSIZEJ==IPHYSICALSIZEJSAVE(JI) .AND. &
114
   &KSIZEJ==ISIZEJSAVE(JI) .AND. KSIZEI==ISIZEISAVE(JI) .AND. &
115
   &ICODEILOEN==ILOENSAVE(JI) .AND. &
116
   &PDELTAX==ZDELTAXSAVE(JI) .AND. PDELTAY==ZDELTAYSAVE(JI)) THEN
117
    KIDENTRESOL=IIDENTRESOLSAVE(JI)
118
    LLNEWRESOL=.FALSE.
119
  ENDIF
120
ENDDO
121
IF(LLNEWRESOL) THEN
122
  INBRESOL=INBRESOL+1
123
  IF(INBRESOL>INUMMAXRESOLSAVE) THEN
124
    PRINT*, "Error in W_SPEC_SETUP: Maximum number of resolution is exceeded."
125
    KRETURNCODE=-999
126
    LDSTOP=.TRUE.
127
  ENDIF
128
ENDIF
129
!
130
! III.c Setup LAM or global spectral transform - once by resolution
131
IF(LLNEWRESOL .AND. .NOT. LDSTOP) THEN
132
  ! The following code is exectuded once for each resolution
133
  ITRUNCXSAVE(INBRESOL)=KTRUNCX
134
  ITRUNCYSAVE(INBRESOL)=KTRUNCY
135
  IPHYSICALSIZEISAVE(INBRESOL)=KPHYSICALSIZEI
136
  IPHYSICALSIZEJSAVE(INBRESOL)=KPHYSICALSIZEJ
137
  ISIZEISAVE(INBRESOL)=KSIZEI
138
  ISIZEJSAVE(INBRESOL)=KSIZEJ
139
  ILOENSAVE(INBRESOL)=ICODEILOEN
140
  ZDELTAXSAVE(INBRESOL)=PDELTAX
141
  ZDELTAYSAVE(INBRESOL)=PDELTAY
142
  IF(LDLAM) THEN
143
    ZEXWN=2*3.141592653589797/(KSIZEI*PDELTAX)
144
    ZEYWN=2*3.141592653589797/(KSIZEJ*PDELTAY)
145
    CALL ESETUP_TRANS(KMSMAX=ITRUNCXSAVE(INBRESOL), KSMAX=ITRUNCYSAVE(INBRESOL), &
146
                     &KDGUX=IPHYSICALSIZEJSAVE(INBRESOL), &
147
                     &KDGL=ISIZEJSAVE(INBRESOL), KLOEN=ILOEN(:), KRESOL=IIDENTRESOLSAVE(INBRESOL), &
148
                     &PEXWN=ZEXWN, PEYWN=ZEYWN)
149
  ELSE
150
    PRINT*, "Setup spectral transform"
151
    CALL SETUP_TRANS(KSMAX=ITRUNCXSAVE(INBRESOL), KDGL=ISIZEJSAVE(INBRESOL), &
152
                    &KLOEN=ILOEN(1:ISIZEJSAVE(INBRESOL)), KRESOL=IIDENTRESOLSAVE(INBRESOL))
153
    PRINT*, "End Setup spectral transform"
154
  ENDIF
155
  KIDENTRESOL=IIDENTRESOLSAVE(INBRESOL)
156
ENDIF
157
END SUBROUTINE W_SPEC_SETUP
158

    
159
!__________________________________________________________________________
160

    
161
SUBROUTINE W_TRANS_INQ(KRETURNCODE, KSIZEJ, KTRUNC, KSLOEN, KLOEN, KNUMMAXRESOL, &
162
                      &KGPTOT, KSPEC, KNMENG)
163
! ** PURPOSE
164
!    Simplified wrapper to TRANS_INQ
165
!
166
! ** DUMMY ARGUMENTS
167
!    KSIZEJ: number of latitudes in grid-point space
168
!    KTRUNC: troncature
169
!    KSLOEN: Size of KLOEN
170
!    KLOEN: number of points on each latitude row
171
!    KNUMMAXRESOL: maximum number of troncatures handled
172
!    KGPTOT: number of gridpoints
173
!    KSPEC: number of spectral coefficients
174
!    KNMENG: cut-off zonal wavenumber
175
!
176
! ** AUTHOR
177
!    9 April 2014, S. Riette
178
!
179
! ** MODIFICATIONS
180
!    6 Jan., S. Riette: w_spec_setup interfaced modified
181
!
182
! I. Dummy arguments declaration
183
IMPLICIT NONE
184
INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
185
INTEGER(KIND=8), INTENT(IN) :: KSIZEJ
186
INTEGER(KIND=8), INTENT(IN) :: KTRUNC
187
INTEGER(KIND=8), INTENT(IN) :: KSLOEN
188
INTEGER(KIND=8), DIMENSION(KSLOEN), INTENT(IN) :: KLOEN
189
INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
190
INTEGER(KIND=8), INTENT(OUT) :: KGPTOT
191
INTEGER(KIND=8), INTENT(OUT) :: KSPEC
192
INTEGER(KIND=8), DIMENSION(KSLOEN), INTENT(OUT) :: KNMENG
193
!
194
! II. Local variables declaration
195
INTEGER, DIMENSION(SIZE(KLOEN)) :: ILOEN
196
INTEGER :: ISIZEI, ISIZEJ, &
197
         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
198
         & ITRUNCX, ITRUNCY, &
199
         & INUMMAXRESOL
200
LOGICAL :: LLSTOP
201
INTEGER :: IIDENTRESOL
202
INTEGER :: IGPTOT, ISPEC
203
INTEGER, DIMENSION(SIZE(KLOEN)) :: INMENG
204
REAL(KIND=8) :: ZDELTAX, ZDELTAY
205
#include "trans_inq.h"
206

    
207
ILOEN(:)=KLOEN(:)
208
ISIZEI=0
209
ISIZEJ=KSIZEJ
210
IPHYSICALSIZEI=0
211
IPHYSICALSIZEJ=0
212
ITRUNCX=KTRUNC
213
ITRUNCY=0
214
INUMMAXRESOL=KNUMMAXRESOL
215
INMENG(:)=KNMENG(:)
216
!
217
! III. Setup
218
ZDELTAX=0.
219
ZDELTAY=0.
220
CALL W_SPEC_SETUP(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
221
                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .FALSE., SIZE(ILOEN), &
222
                  &ZDELTAX, ZDELTAY, IIDENTRESOL, LLSTOP)
223
IF (.NOT. LLSTOP) THEN
224
  CALL TRANS_INQ(KRESOL=IIDENTRESOL, KGPTOT=IGPTOT, KSPEC=ISPEC, KNMENG=INMENG)
225
  KGPTOT=IGPTOT
226
  KSPEC=ISPEC
227
  KNMENG=INMENG
228
ENDIF
229
!
230
END SUBROUTINE W_TRANS_INQ
231

    
232
!___________________________________________________________________________________________
233

    
234
SUBROUTINE W_ETRANS_INQ(KRETURNCODE, KSIZEI, KSIZEJ, KPHYSICALSIZEI, KPHYSICALSIZEJ, &
235
                       &KTRUNCX, KTRUNCY, KNUMMAXRESOL, PDELTAX, PDELTAY, &
236
                       &KGPTOT, KSPEC)
237
! ** PURPOSE
238
!    Simplified wrapper to ETRANS_INQ
239
!
240
! ** DUMMY ARGUMENTS
241
!    KSIZEI, KSIZEJ: size of grid-point field (with extension zone)
242
!    KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field
243
!    KTRUNCX, KTRUNCY: troncatures
244
!    KNUMMAXRESOL: maximum number of troncatures handled
245
!    PDELTAX: x resolution
246
!    PDELTAY: y resolution
247
!    KGPTOT: number of gridpoints
248
!    KSPEC: number of spectral coefficients
249
!
250
! ** AUTHOR
251
!    9 April 2014, S. Riette
252
!
253
! ** MODIFICATIONS
254
!    6 Jan., S. Riette: PDELTAX and PDELTAY added
255
!
256
! I. Dummy arguments declaration
257
IMPLICIT NONE
258
INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
259
INTEGER(KIND=8), INTENT(IN) :: KSIZEI, KSIZEJ
260
INTEGER(KIND=8), INTENT(IN) :: KPHYSICALSIZEI, KPHYSICALSIZEJ
261
INTEGER(KIND=8), INTENT(IN) :: KTRUNCX, KTRUNCY
262
INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
263
REAL(KIND=8),    INTENT(IN) :: PDELTAX
264
REAL(KIND=8),    INTENT(IN) :: PDELTAY
265
INTEGER(KIND=8), INTENT(OUT) :: KGPTOT
266
INTEGER(KIND=8), INTENT(OUT) :: KSPEC
267
!
268
! II. Local variables declaration
269
INTEGER, DIMENSION(0:KTRUNCX) :: IESM0
270
INTEGER :: ISIZEI, ISIZEJ, &
271
         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
272
         & ITRUNCX, ITRUNCY, &
273
         & INUMMAXRESOL
274
LOGICAL :: LLSTOP
275
INTEGER :: IIDENTRESOL
276
INTEGER, DIMENSION(1) :: ILOEN
277
INTEGER :: IGPTOT, ISPEC
278

    
279
#include "etrans_inq.h"
280

    
281
ISIZEI=KSIZEI
282
ISIZEJ=KSIZEJ
283
IPHYSICALSIZEI=KPHYSICALSIZEI
284
IPHYSICALSIZEJ=KPHYSICALSIZEJ
285
ITRUNCX=KTRUNCX
286
ITRUNCY=KTRUNCY
287
INUMMAXRESOL=KNUMMAXRESOL
288

    
289
! III. Setup
290
CALL W_SPEC_SETUP(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
291
                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .TRUE., 1, &
292
                  &PDELTAX, PDELTAY, IIDENTRESOL, LLSTOP)
293
IF (.NOT. LLSTOP) THEN
294
  CALL ETRANS_INQ(KRESOL=IIDENTRESOL, KGPTOT=IGPTOT, KSPEC=ISPEC, KESM0=IESM0)
295
  KGPTOT=IGPTOT
296
  KSPEC=ISPEC
297
ENDIF
298
!
299
END SUBROUTINE W_ETRANS_INQ
300

    
301
!______________________________________________________________________
302

    
303
SUBROUTINE W_SPEC2GPT_LAM(KRETURNCODE, KSIZEI, KSIZEJ, KPHYSICALSIZEI, KPHYSICALSIZEJ, &
304
                         &KTRUNCX, KTRUNCY, KNUMMAXRESOL, KSIZE, LGRADIENT, LREORDER, PDELTAX, PDELTAY, &
305
                         &PSPEC, PGPT, PGPTM, PGPTL)
306
! ** PURPOSE
307
!    Transform spectral coefficients into grid-point values
308
!
309
! ** DUMMY ARGUMENTS
310
!    KRETURNCODE: error code
311
!    KSIZEI, KSIZEJ: size of grid-point field (with extension zone)
312
!    KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field
313
!    KTRUNCX, KTRUNCY: troncatures
314
!    KNUMMAXRESOL: maximum number of troncatures handled
315
!    KSIZE: size of PSPEC
316
!    LREORDER: switch to reorder spectral coefficients or not
317
!    LGRADIENT: switch to compute or not gradient
318
!    PDELTAX: x resolution
319
!    PDELTAY: y resolution
320
!    PSPEC: spectral coefficient array
321
!    PGPT: grid-point field
322
!    PGPTM: N-S derivative if LGRADIENT
323
!    PGPTL: E-W derivative if LGRADIENT
324
!
325
! ** AUTHOR
326
!    9 April 2014, S. Riette
327
!
328
! ** MODIFICATIONS
329
!    5 Jan., S. Riette: PDELTAX, PDELTAY, LGRADIENT, PGPTM and PGPTL added
330
!    March, 2016, A.Mary: LREORDER
331
!
332
! I. Dummy arguments declaration
333
IMPLICIT NONE
334
INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
335
INTEGER(KIND=8), INTENT(IN) :: KSIZEI, KSIZEJ
336
INTEGER(KIND=8), INTENT(IN) :: KPHYSICALSIZEI, KPHYSICALSIZEJ
337
INTEGER(KIND=8), INTENT(IN) :: KTRUNCX, KTRUNCY
338
INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
339
INTEGER(KIND=8), INTENT(IN) :: KSIZE
340
LOGICAL, INTENT(IN) :: LGRADIENT
341
LOGICAL, INTENT(IN) :: LREORDER
342
REAL(KIND=8), INTENT(IN) :: PDELTAX
343
REAL(KIND=8), INTENT(IN) :: PDELTAY
344
REAL(KIND=8), DIMENSION(KSIZE), INTENT(IN) :: PSPEC
345
REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ), INTENT(OUT) :: PGPT
346
REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ), INTENT(OUT) :: PGPTM
347
REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ), INTENT(OUT) :: PGPTL
348
!
349
! II. Local variables declaration
350
INTEGER, DIMENSION(0:KTRUNCX) :: IESM0
351
INTEGER :: IGPTOT, ISPEC
352
INTEGER, DIMENSION(0:KTRUNCY) :: ISPECINI, ISPECEND
353
REAL(KIND=8), DIMENSION(1, KSIZE) :: ZSPBUF
354
REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ, 3, 1) :: ZGPBUF
355
INTEGER :: JI, JM, JN, IINDEX, IIDENTRESOL
356
LOGICAL :: LLSTOP
357
INTEGER :: ISIZEI, ISIZEJ, &
358
         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
359
         & ITRUNCX, ITRUNCY, &
360
         & INUMMAXRESOL
361
INTEGER, DIMENSION(1) :: ILOEN
362

    
363
#include "einv_trans.h"
364
#include "etrans_inq.h"
365

    
366
KRETURNCODE=0
367
LLSTOP=.FALSE.
368
ISIZEI=KSIZEI
369
ISIZEJ=KSIZEJ
370
IPHYSICALSIZEI=KPHYSICALSIZEI
371
IPHYSICALSIZEJ=KPHYSICALSIZEJ
372
ITRUNCX=KTRUNCX
373
ITRUNCY=KTRUNCY
374
INUMMAXRESOL=KNUMMAXRESOL
375
ILOEN(:)=0
376

    
377
! III. Setup
378
CALL W_SPEC_SETUP(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
379
                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .TRUE., 1, &
380
                  &PDELTAX, PDELTAY, IIDENTRESOL, LLSTOP)
381

    
382
! IV. Transformation
383

    
384
! IV.a Shape of coefficient array
385
!IGPTOT is the total number of points in grid-point space
386
!ISPEC is the number of spectral coefficients
387
!IESM0(m) is the index of spectral coefficient (m,0) in model
388
!ISPECINI(n) is the index of the first of the 4 spectral coefficient (0,n) in FA file
389
!ISPECEND(n) is the index of the last of the last 4 spectral coefficients (:,n) in FA file
390
IF (.NOT. LLSTOP) THEN
391
  CALL ETRANS_INQ(KRESOL=IIDENTRESOL, KGPTOT=IGPTOT, KSPEC=ISPEC, KESM0=IESM0)
392
  JI=1
393
  DO JN=0, ITRUNCY
394
    ISPECINI(JN)=(JI-1)*4+1
395
    JI=JI+COUNT(IESM0(1:ITRUNCX)-IESM0(0:ITRUNCX-1)>JN*4)
396
    IF (ISPEC-IESM0(ITRUNCX)>JN*4) JI=JI+1
397
    ISPECEND(JN)=(JI-1)*4
398
  ENDDO
399
ENDIF
400

    
401
! III.b Reordering
402
! reorder Aladin :  file ordering = coeffs per blocks of m, 4 reals per coeff
403
!           Aladin array ordering = coeffs per blocks of n, 4 reals per coeff
404
IF (LREORDER) THEN
405
  IF (.NOT. LLSTOP) THEN
406
    ZSPBUF(:,:)=0.
407
    JI=1
408
    DO JM=0,ITRUNCX+1
409
      DO JN=0,ITRUNCY
410
        IF (ISPECINI(JN)+JM*4+3<=ISPECEND(JN)) THEN
411
          DO IINDEX=ISPECINI(JN)+JM*4, ISPECINI(JN)+JM*4+3
412
            ZSPBUF(1,JI)=PSPEC(IINDEX)
413
            JI=JI+1
414
          ENDDO
415
        ENDIF
416
      ENDDO
417
    ENDDO
418
    IF (JI/=ISPEC+1) THEN
419
      PRINT*, "Internal error in W_SPEC2GPT_LAM (spectral reordering)"
420
      KRETURNCODE=-999
421
      LLSTOP=.TRUE.
422
    ENDIF
423
  ENDIF
424
ELSE
425
  ZSPBUF(1,:) = PSPEC(:)
426
ENDIF
427

    
428
! III.c Inverse transform
429
IF (.NOT. LLSTOP) THEN
430
  IF (.NOT. LGRADIENT) THEN
431
      CALL EINV_TRANS(PSPSCALAR=ZSPBUF(:,:), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL)
432
      PGPT(:)=ZGPBUF(:,1,1)
433
  ELSE
434
      CALL EINV_TRANS(PSPSCALAR=ZSPBUF(:,:), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL, LDSCDERS=.TRUE.)
435
      PGPT(:)=ZGPBUF(:,1,1)
436
      PGPTM(:)=ZGPBUF(:,2,1)
437
      PGPTL(:)=ZGPBUF(:,3,1)
438
  ENDIF
439
ENDIF
440

    
441
END SUBROUTINE W_SPEC2GPT_LAM
442

    
443
!____________________________________________________________________________
444

    
445
SUBROUTINE W_SPEC2GPT_GAUSS(KRETURNCODE, KSIZEJ, KTRUNC, KNUMMAXRESOL, KGPTOT, KSLOEN, KLOEN, KSIZE, &
446
                          & LGRADIENT, LREORDER, PSPEC, PGPT, PGPTM, PGPTL)
447
! ** PURPOSE
448
!    Transform spectral coefficients into grid-point values
449
!
450
! ** DUMMY ARGUMENTS
451
!    KSIZEJ: Number of latitudes
452
!    KTRUNC: troncature
453
!    KNUMMAXRESOL: maximum number of troncatures handled
454
!    KGPTOT: number of grid-points
455
!    KSLOEN: Size of KLOEN
456
!    KLOEN:
457
!    KSIZE: Size of PSPEC
458
!    LREORDER: switch to reorder spectral coefficients or not
459
!    LGRADIENT: switch to compute or not gradient
460
!    PSPEC: spectral coefficient array
461
!    PGPT: grid-point field
462
!    PGPTM: N-S derivative if LGRADIENT
463
!    PGPTL: E-W derivative if LGRADIENT
464
!
465
! ** AUTHOR
466
!    9 April 2014, S. Riette
467
!
468
! ** MODIFICATIONS
469
!    6 Jan., S. Riette: w_spec_setup interface modified
470
!    March, 2016, A.Mary: LREORDER
471
!    Sept., 2016, A.Mary: LGRADIENT
472
!
473
! I. Dummy arguments declaration
474
IMPLICIT NONE
475
INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
476
INTEGER(KIND=8), INTENT(IN) :: KSIZEJ
477
INTEGER(KIND=8), INTENT(IN) :: KTRUNC
478
INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
479
INTEGER(KIND=8), INTENT(IN) :: KGPTOT
480
INTEGER(KIND=8), INTENT(IN) :: KSLOEN
481
INTEGER(KIND=8), DIMENSION(KSLOEN), INTENT(IN) :: KLOEN
482
INTEGER(KIND=8), INTENT(IN) :: KSIZE
483
LOGICAL, INTENT(IN) :: LGRADIENT
484
LOGICAL, INTENT(IN) :: LREORDER
485
REAL(KIND=8), DIMENSION(KSIZE), INTENT(IN) :: PSPEC
486
REAL(KIND=8), DIMENSION(KGPTOT), INTENT(OUT) :: PGPT
487
REAL(KIND=8), DIMENSION(KGPTOT), INTENT(OUT) :: PGPTM
488
REAL(KIND=8), DIMENSION(KGPTOT), INTENT(OUT) :: PGPTL
489
!
490
! II. Local variables declaration
491
INTEGER, DIMENSION(SIZE(KLOEN)) :: ILOEN
492
INTEGER :: ISIZEI, ISIZEJ, &
493
         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
494
         & ITRUNCX, ITRUNCY, &
495
         & INUMMAXRESOL
496
LOGICAL :: LLSTOP
497
INTEGER :: IIDENTRESOL
498
INTEGER :: JI, JM, JN
499
INTEGER, DIMENSION(0:KTRUNC) :: NASM0
500
REAL(KIND=8), DIMENSION(1, KSIZE) :: ZSPBUF
501
REAL(KIND=8), DIMENSION(KGPTOT, 3, 1) :: ZGPBUF
502
REAL(KIND=8) :: ZDELTAX, ZDELTAY
503
#include "trans_inq.h"
504
#include "inv_trans.h"
505

    
506
ILOEN(:)=KLOEN(:)
507
ISIZEI=0
508
ISIZEJ=KSIZEJ
509
IPHYSICALSIZEI=0
510
IPHYSICALSIZEJ=0
511
ITRUNCX=KTRUNC
512
ITRUNCY=0
513
INUMMAXRESOL=KNUMMAXRESOL
514
!
515
! III. Setup
516
ZDELTAX=0.
517
ZDELTAY=0.
518
CALL W_SPEC_SETUP(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
519
                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .FALSE., SIZE(ILOEN), &
520
                  &ZDELTAX, ZDELTAY, IIDENTRESOL, LLSTOP)
521
!
522
! IV. Transformation
523
IF (LREORDER) THEN
524
  ! IV.a Shape of coefficient array
525
  IF (.NOT. LLSTOP) THEN
526
    JI=1
527
    DO JN=0, KTRUNC
528
      NASM0(JN)=JI
529
      JI=JI+1+JN+(JN+1)
530
    ENDDO
531
  ENDIF
532

    
533
  ! IV.b Reordering
534
  IF(.NOT. LLSTOP) THEN
535
    ZSPBUF(1,:)=0.
536
    JI=1
537
    DO JM=0, KTRUNC
538
      DO JN=JM, KTRUNC
539
        ZSPBUF(1,JI)=PSPEC(NASM0(JN)+JM)
540
        JI=JI+1
541
        IF(JM==0) THEN
542
          ZSPBUF(1,JI)=0
543
        ELSE
544
          ZSPBUF(1,JI)=PSPEC(NASM0(JN)-JM)
545
        ENDIF
546
        JI=JI+1
547
      ENDDO
548
    ENDDO
549
  ENDIF
550
ELSE
551
  ZSPBUF(1,:) = PSPEC(:)
552
ENDIF
553

    
554
! IV.c Inverse transform
555
IF (.NOT. LLSTOP) THEN
556
  IF (.NOT. LGRADIENT) THEN
557
    CALL INV_TRANS(PSPSCALAR=ZSPBUF(:,:), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL)
558
    PGPT(:)=ZGPBUF(:,1,1)
559
  ELSE
560
    CALL INV_TRANS(PSPSCALAR=ZSPBUF(:,:), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL, LDSCDERS=.TRUE.)
561
    PGPT(:)=ZGPBUF(:,1,1)
562
    PGPTM(:)=ZGPBUF(:,2,1)
563
    PGPTL(:)=ZGPBUF(:,3,1)
564
  ENDIF
565
ENDIF
566
END SUBROUTINE W_SPEC2GPT_GAUSS
567

    
568
!______________________________________________________________________
569

    
570
SUBROUTINE W_GPT2SPEC_LAM(KRETURNCODE, KSIZE, KSIZEI, KSIZEJ, KPHYSICALSIZEI, KPHYSICALSIZEJ, &
571
                         &KTRUNCX, KTRUNCY, KNUMMAXRESOL, PDELTAX, PDELTAY, LREORDER, PGPT, PSPEC)
572
! ** PURPOSE
573
!    Transform grid point values into spectral coefficients
574
!
575
! ** DUMMY ARGUMENTS
576
!    KRETURNCODE: error code
577
!    KSIZE: size of spectral field
578
!    KSIZEI, KSIZEJ: size of grid-point field (with extension zone)
579
!    KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field
580
!    KTRUNCX, KTRUNCY: troncatures
581
!    KNUMMAXRESOL: maximum number of troncatures handled
582
!    PDELTAX: x resolution
583
!    PDELTAY: y resolution
584
!    LREORDER: switch to reorder spectral coefficients or not
585
!    PGPT: grid-point field
586
!    PSPEC: spectral coefficient array
587
!
588
! ** AUTHOR
589
!    9 April 2014, S. Riette
590
!
591
! ** MODIFICATIONS
592
!    6 Jan., S. Riette: PDELTAX and PDELTAY added
593
!    March, 2016, A.Mary: LREORDER
594
!
595
! I. Dummy arguments declaration
596
IMPLICIT NONE
597
INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
598
INTEGER(KIND=8), INTENT(IN) :: KSIZE, KSIZEI, KSIZEJ
599
INTEGER(KIND=8), INTENT(IN) :: KPHYSICALSIZEI, KPHYSICALSIZEJ
600
INTEGER(KIND=8), INTENT(IN) :: KTRUNCX, KTRUNCY
601
INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
602
REAL(KIND=8), INTENT(IN) :: PDELTAX
603
REAL(KIND=8), INTENT(IN) :: PDELTAY
604
LOGICAL, INTENT(IN) :: LREORDER
605
REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ), INTENT(IN) :: PGPT
606
REAL(KIND=8), DIMENSION(KSIZE), INTENT(OUT) :: PSPEC
607
!
608
! II. Local variables declaration
609
INTEGER, DIMENSION(0:KTRUNCX) :: IESM0
610
INTEGER :: IGPTOT, ISPEC
611
INTEGER, DIMENSION(0:KTRUNCY) :: ISPECINI, ISPECEND
612
REAL(KIND=8), DIMENSION(1, KSIZEI*KSIZEJ) :: ZSPBUF !size over-evaluated
613
REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ, 1, 1) :: ZGPBUF
614
INTEGER :: JI, JM, JN, IIDENTRESOL
615
LOGICAL :: LLSTOP
616
INTEGER :: ISIZEI, ISIZEJ, &
617
         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
618
         & ITRUNCX, ITRUNCY, &
619
         & INUMMAXRESOL
620
INTEGER, DIMENSION(1) :: ILOEN
621

    
622
#include "edir_trans.h"
623
#include "etrans_inq.h"
624

    
625
KRETURNCODE=0
626
LLSTOP=.FALSE.
627
ISIZEI=KSIZEI
628
ISIZEJ=KSIZEJ
629
IPHYSICALSIZEI=KPHYSICALSIZEI
630
IPHYSICALSIZEJ=KPHYSICALSIZEJ
631
ITRUNCX=KTRUNCX
632
ITRUNCY=KTRUNCY
633
INUMMAXRESOL=KNUMMAXRESOL
634

    
635
! III. Setup
636
CALL W_SPEC_SETUP(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
637
                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .TRUE., 1, &
638
                  &PDELTAX, PDELTAY, IIDENTRESOL, LLSTOP)
639

    
640
! IV. Transformation
641

    
642
! IV.a Shape of coefficient array
643
!IGPTOT is the total number of points in grid-point space
644
!ISPEC is the number of spectral coefficients
645
!IESM0(m) is the index of spectral coefficient (m,0) in model
646
!ISPECINI(n) is the index of the first of the 4 spectral coefficient (0,n) in FA file
647
!ISPECEND(n) is the index of the last of the last 4 spectral coefficients (:,n) in FA file
648
IF (.NOT. LLSTOP) THEN
649
  CALL ETRANS_INQ(KRESOL=IIDENTRESOL, KGPTOT=IGPTOT, KSPEC=ISPEC, KESM0=IESM0)
650
  JI=1
651
  DO JN=0, ITRUNCY
652
    ISPECINI(JN)=(JI-1)*4+1
653
    JI=JI+COUNT(IESM0(1:ITRUNCX)-IESM0(0:ITRUNCX-1)>JN*4)
654
    IF (ISPEC-IESM0(ITRUNCX)>JN*4) JI=JI+1
655
    ISPECEND(JN)=(JI-1)*4
656
  ENDDO
657
ENDIF
658

    
659
! III.b transform
660
IF (.NOT. LLSTOP) THEN
661
  ZGPBUF(:,1,1)=PGPT(:)
662
  CALL EDIR_TRANS(PSPSCALAR=ZSPBUF(:,:), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL)
663
ENDIF
664

    
665
! III.c Reordering
666
! reorder Aladin :  file ordering = coeffs per blocks of m, 4 reals per coeff
667
!           Aladin array ordering = coeffs per blocks of n, 4 reals per coeff
668
IF (LREORDER) THEN
669
  IF (.NOT. LLSTOP) THEN
670
    JI=1
671
    PSPEC(:)=0.
672
    DO JM=0,ITRUNCX*4+4,4
673
      DO JN=0,ITRUNCY
674
        IF (ISPECINI(JN)+JM+3<=ISPECEND(JN)) THEN
675
          PSPEC(ISPECINI(JN)+JM:ISPECINI(JN)+JM+3) = ZSPBUF(1,JI:JI+3)
676
          JI=JI+4
677
        ENDIF
678
      ENDDO
679
    ENDDO
680
    IF(JI/=ISPEC+1) THEN
681
      PRINT*, "Internal error in W_GPT2SPEC_LAM (spectral reordering)"
682
      KRETURNCODE=-999
683
    ENDIF
684
  ENDIF
685
ELSE
686
  PSPEC(1:KSIZE) = ZSPBUF(1,1:KSIZE)
687
ENDIF
688

    
689
END SUBROUTINE W_GPT2SPEC_LAM
690

    
691
!____________________________________________________________________________
692

    
693
SUBROUTINE W_GPT2SPEC_GAUSS(KRETURNCODE, KSPEC, KSIZEJ, KTRUNC, KNUMMAXRESOL, KSLOEN, KLOEN, KSIZE, LREORDER, PGPT, PSPEC)
694
! ** PURPOSE
695
!    Transform spectral coefficients into grid-point values
696
!
697
! ** DUMMY ARGUMENTS
698
!    KRETURNCODE: error code
699
!    KSPEC: size of spectral coefficients array
700
!    KSIZEJ: Number of latitudes
701
!    KTRUNC: troncature
702
!    KNUMMAXRESOL: maximum number of troncatures handled
703
!    KSLOEN: Size ok KLOEN
704
!    KLOEN
705
!    KSIZE: Size of PGPT
706
!    LREORDER: switch to reorder spectral coefficients or not
707
!    PGPT: grid-point field
708
!    PSPEC: spectral coefficient array
709
!
710
! ** AUTHOR
711
!    9 April 2014, S. Riette
712
!
713
! ** MODIFICATIONS
714
!    6 Jan. 2016, S. Riette: w_spec_setup interface modified
715
!    March, 2016, A.Mary: LREORDER
716
!
717
! I. Dummy arguments declaration
718
IMPLICIT NONE
719
INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
720
INTEGER(KIND=8), INTENT(IN) :: KSPEC
721
INTEGER(KIND=8), INTENT(IN) :: KSIZEJ
722
INTEGER(KIND=8), INTENT(IN) :: KTRUNC
723
INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
724
INTEGER(KIND=8), INTENT(IN) :: KSLOEN
725
INTEGER(KIND=8), DIMENSION(KSLOEN), INTENT(IN) :: KLOEN
726
INTEGER(KIND=8), INTENT(IN) :: KSIZE
727
LOGICAL, INTENT(IN) :: LREORDER
728
REAL(KIND=8), DIMENSION(KSIZE), INTENT(IN) :: PGPT
729
REAL(KIND=8), DIMENSION(KSPEC), INTENT(OUT) :: PSPEC
730
!
731
! II. Local variables declaration
732
INTEGER, DIMENSION(SIZE(KLOEN)) :: ILOEN
733
INTEGER :: ISIZEI, ISIZEJ, &
734
         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
735
         & ITRUNCX, ITRUNCY, &
736
         & INUMMAXRESOL
737
LOGICAL :: LLSTOP
738
INTEGER :: IIDENTRESOL
739
INTEGER :: JI, JM, JN
740
INTEGER, DIMENSION(0:KTRUNC) :: NASM0
741
REAL(KIND=8), DIMENSION(1, SIZE(PGPT)) :: ZSPBUF !size over-evaluated
742
REAL(KIND=8), DIMENSION(SIZE(PGPT), 1, 1) :: ZGPBUF
743
REAL(KIND=8) :: ZDELTAX, ZDELTAY
744

    
745
#include "trans_inq.h"
746
#include "dir_trans.h"
747
KRETURNCODE=0
748
ILOEN(:)=KLOEN(:)
749
ISIZEI=0
750
ISIZEJ=KSIZEJ
751
IPHYSICALSIZEI=0
752
IPHYSICALSIZEJ=0
753
ITRUNCX=KTRUNC
754
ITRUNCY=0
755
INUMMAXRESOL=KNUMMAXRESOL
756
!
757
! III. Setup
758
ZDELTAX=0.
759
ZDELTAY=0.
760
CALL W_SPEC_SETUP(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
761
                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .FALSE., SIZE(ILOEN), &
762
                  &ZDELTAX, ZDELTAY, IIDENTRESOL, LLSTOP)
763
!
764
! IV. Transformation
765
! IV.a Shape of coefficient array
766
IF (.NOT. LLSTOP) THEN
767
  JI=1
768
  DO JN=0, KTRUNC
769
    NASM0(JN)=JI
770
    JI=JI+1+JN+(JN+1)
771
  ENDDO
772
ENDIF
773

    
774
! IV.b Direct transform
775
IF (.NOT. LLSTOP) THEN
776
  ZGPBUF(:,1,1)=PGPT(:)
777
  CALL DIR_TRANS(PSPSCALAR=ZSPBUF(:,:), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL)
778
ENDIF
779

    
780
! IV.c Reordering
781
IF (LREORDER) THEN
782
  IF(.NOT. LLSTOP) THEN
783
    PSPEC(:)=0.
784
    JI=1
785
    DO JM=0, KTRUNC
786
      DO JN=JM, KTRUNC
787
        PSPEC(NASM0(JN)+JM)=ZSPBUF(1,JI)
788
        JI=JI+1
789
        IF(JM/=0) THEN
790
          PSPEC(NASM0(JN)-JM)=ZSPBUF(1,JI)
791
        ENDIF
792
        JI=JI+1
793
      ENDDO
794
    ENDDO
795
    IF(JI-1/=KSPEC) THEN
796
      PRINT*, "Internal error in W_GPT2SPEC_GAUSS (spectral reordering)"
797
      KRETURNCODE=-999
798
    ENDIF
799
  ENDIF
800
ELSE
801
  PSPEC(1:KSPEC) = ZSPBUF(1,1:KSPEC)
802
ENDIF
803

    
804
END SUBROUTINE W_GPT2SPEC_GAUSS
805

    
806
SUBROUTINE W_SPEC2GPT_FFT1D(KSIZES, KTRUNC, PSPEC, KSIZEG, PGPT)
807
! ** PURPOSE
808
!    Transform spectral coefficients into grid-point values,
809
!    for a 1D array (vertical section academic model)
810
!
811
! ** DUMMY ARGUMENTS
812
!    KSIZES size of PSPEC
813
!    KTRUNC: troncature
814
!    PSPEC: spectral coefficient array
815
!    KSIZEG: size of grid-point field (with extension zone)
816
!    PGPT: grid-point field
817
!
818
! ** AUTHOR
819
!    26 March 2015, A. Mary, from utilities/pinuts/module/fa_datas_mod.F90
820
!
821
! ** MODIFICATIONS
822
!
823
! I. Dummy arguments declaration
824
IMPLICIT NONE
825

    
826
INTEGER(KIND=8), INTENT(IN) :: KSIZES
827
INTEGER(KIND=8), INTENT(IN) :: KTRUNC
828
REAL(KIND=8), DIMENSION(KSIZES), INTENT(IN) :: PSPEC
829
INTEGER(KIND=8), INTENT(IN) :: KSIZEG
830
REAL(KIND=8), DIMENSION(KSIZEG), INTENT(OUT) :: PGPT
831

    
832
INTEGER(KIND=8)                             :: NSEFRE, NFTM, NDGLSUR
833
REAL(KIND=8), DIMENSION(:), ALLOCATABLE     :: SP2, TRIGSE
834
INTEGER(KIND=8), DIMENSION(:), ALLOCATABLE  :: NFAXE
835
INTEGER(KIND=8), PARAMETER                  :: NZERO=0
836

    
837
NDGLSUR = KSIZEG+MOD(KSIZEG,2)+2
838
NFTM    = 2*(KTRUNC+1)
839
ALLOCATE(SP2(NDGLSUR*NFTM))
840
SP2     = 0.0
841
SP2     = CONVRT2FFT(PSPEC,NZERO,KTRUNC,NDGLSUR)
842
ALLOCATE(NFAXE(1:10))
843
ALLOCATE(TRIGSE(1:KSIZEG)) 
844
CALL SET99(TRIGSE,NFAXE,KSIZEG)
845
CALL FFT992(SP2(1:KSIZEG), TRIGSE, NFAXE, 1, NDGLSUR, KSIZEG, 1, 1)
846
DEALLOCATE(TRIGSE)
847
DEALLOCATE(NFAXE)
848
PGPT(:) = SP2(1:KSIZEG)
849

    
850
CONTAINS
851

    
852
! from utilities/pinuts/module/fa_datas_mod.F90
853
! and utilities/pinuts/module/array_lib_mod.F90
854

    
855
FUNCTION CONVRT2FFT(IN,X,Y,N) RESULT(OU)
856
REAL(KIND=8),DIMENSION(:),INTENT(IN)             :: IN
857
INTEGER(KIND=8),INTENT(IN)                       :: X, Y, N
858
REAL(KIND=8),DIMENSION(N*2*(X+1))                :: OU
859

    
860
INTEGER(KIND=8),DIMENSION(2*(X+1),(N/2))         :: MINQ 
861
INTEGER(KIND=8),DIMENSION((N/2),2*(X+1))         :: TMINQ 
862
REAL(KIND=8),DIMENSION(2*(X+1),(N/2))            :: OMINQ, EMINQ
863
REAL(KIND=8),DIMENSION((N/2),2*(X+1))            :: TOMINQ, TEMINQ   
864
REAL(KIND=8),DIMENSION(N*(X+1))                  :: OINI, EINI
865
REAL(KIND=8), PARAMETER                          :: ZZERO=0.0
866

    
867
CALL SPLIT_ODEV(IN,OINI,EINI)
868
MINQ   = MASQ(X,Y,N)
869
OMINQ  = UNPACK(OINI,MINQ == 1,ZZERO)
870
TOMINQ = TRANSPOSE(OMINQ)
871
EMINQ  = UNPACK(EINI,MINQ == 1,ZZERO)
872
TEMINQ = TRANSPOSE(EMINQ)
873
TMINQ  = 1
874
OINI   = PACK(TOMINQ,TMINQ > 0)
875
EINI   = PACK(TEMINQ,TMINQ > 0)
876
OU     = MIX_ODEV(OINI,EINI)
877
END FUNCTION CONVRT2FFT
878

    
879
FUNCTION MASQ(X,Y,N) RESULT(T)
880
INTEGER(KIND=8),INTENT(IN)                       :: X, Y, N
881
INTEGER(KIND=8),DIMENSION(1:2*(X+1),1:(N/2))     :: T
882

    
883
INTEGER(KIND=8)                                  :: I, J
884
INTEGER(KIND=8),DIMENSION(0:X)                   :: KM
885
INTEGER(KIND=8),DIMENSION(0:Y)                   :: KN
886
CALL ELLIPS64(X,Y,KN,KM) 
887
T = 0
888
DO I=0,Y
889
  DO J=0,2*KN(I)+1
890
    T(J+1,I+1)=1
891
  END DO
892
END DO
893
END FUNCTION MASQ
894

    
895
FUNCTION MIX_ODEV(TO,TE) RESULT(T)
896
REAL(KIND=8),DIMENSION(:),INTENT(IN)                  :: TO,TE
897
REAL(KIND=8),DIMENSION(SIZE(TO)+SIZE(TE))             :: T
898

    
899
INTEGER(KIND=8)                             :: I
900

    
901
DO I=1,(SIZE(TO)+SIZE(TE))/2
902
  T((2*I)-1)=TE(I)
903
  T(2*I)=TO(I)
904
END DO
905
END FUNCTION MIX_ODEV
906

    
907
SUBROUTINE SPLIT_ODEV(T,TO,TE)
908
REAL(KIND=8),DIMENSION(:),INTENT(IN)                  :: T
909
REAL(KIND=8),DIMENSION(SIZE(T)/2),INTENT(OUT)    :: TO,TE
910

    
911
INTEGER(KIND=8)                             :: I
912

    
913
DO I=1,SIZE(T)/2
914
  TO(I)=T(2*I)
915
  TE(I)=T((2*I)-1)
916
END DO
917
END SUBROUTINE SPLIT_ODEV
918

    
919
END SUBROUTINE W_SPEC2GPT_FFT1D