From 4f247989d67e9b1876b8e79ae2b81dd9450d0d47 Mon Sep 17 00:00:00 2001
From: Anne Tipka <anne.tipka@ctbto.org>
Date: Thu, 21 Jul 2022 13:38:30 +0000
Subject: [PATCH] elimination of emoslib in fortran code. See ticket #312
---
Source/Fortran/changelog.txt | 11 +-
Source/Fortran/fft99.f90 | 218 +++++++++
Source/Fortran/jsppole.f90 | 182 ++++++++
Source/Fortran/makefile_cray | 2 +-
Source/Fortran/makefile_debug | 3 +-
Source/Fortran/makefile_ecgate | 2 +-
Source/Fortran/makefile_fast | 3 +-
Source/Fortran/qpassm.f90 | 778 +++++++++++++++++++++++++++++++++
Source/Fortran/rpassm.f90 | 760 ++++++++++++++++++++++++++++++++
Source/Fortran/set99.f90 | 65 +++
10 files changed, 2017 insertions(+), 7 deletions(-)
create mode 100644 Source/Fortran/fft99.f90
create mode 100644 Source/Fortran/jsppole.f90
create mode 100644 Source/Fortran/qpassm.f90
create mode 100644 Source/Fortran/rpassm.f90
create mode 100644 Source/Fortran/set99.f90
diff --git a/Source/Fortran/changelog.txt b/Source/Fortran/changelog.txt
index 18de918..71ade4d 100644
--- a/Source/Fortran/changelog.txt
+++ b/Source/Fortran/changelog.txt
@@ -22,4 +22,13 @@ rwGRIB2.f90 subroutine readlatlon: check alloc status deallocate values
calc_etadot.out -> calc_etadot_fast.out
makefile_local_gfortran -> makefile_fast
-
+2022-07-21 AT
+ solved issue with emoslib!
+ see ticket #312 on flexpart.eu
+
+ identified the source files from emoslib which calc_etadot needs
+ converted them from f77 to f90 code
+ adapted makefile to get rid of emoslib and make use of the fft source files
+
+ new files: set99.f90, fft99.f90, qpassm.f90, rpassm.f90, jsppole.f90
+
diff --git a/Source/Fortran/fft99.f90 b/Source/Fortran/fft99.f90
new file mode 100644
index 0000000..305f350
--- /dev/null
+++ b/Source/Fortran/fft99.f90
@@ -0,0 +1,218 @@
+!c Copyright 1981-2016 ECMWF.
+!c
+!c This software is licensed under the terms of the Apache Licence
+!c Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+!c
+!c In applying this licence, ECMWF does not waive the privileges and immunities
+!c granted to it by virtue of its status as an intergovernmental organisation
+!c nor does it submit to any jurisdiction.
+!c
+
+ SUBROUTINE FFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
+ DIMENSION A(N),WORK(N),TRIGS(N),IFAX(*)
+!C
+!C SUBROUTINE 'FFT99' - MULTIPLE FAST REAL PERIODIC TRANSFORM
+!C SUPERSEDES PREVIOUS ROUTINE 'FFT99'
+!C
+!C REAL TRANSFORM OF LENGTH N PERFORMED BY REMOVING REDUNDANT
+!C OPERATIONS FROM COMPLEX TRANSFORM OF LENGTH N
+!C
+!C A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA
+!C WORK IS AN AREA OF SIZE (N+1)*MIN(LOT,64)
+!C TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
+!C IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N
+!C INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'
+!C (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
+!C JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
+!C N IS THE LENGTH OF THE DATA VECTORS
+!C LOT IS THE NUMBER OF DATA VECTORS
+!C ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
+!C = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
+!C
+!C ORDERING OF COEFFICIENTS:
+!C A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)
+!C WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED
+!C
+!C ORDERING OF DATA:
+!C X(N-1),X(0),X(1),X(2),...,X(N-1),X(0)
+!C I.E. EXPLICIT CYCLIC CONTINUITY; (N+2) LOCATIONS REQUIRED
+!C
+!C VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS
+!C IN PARALLEL
+!C
+!C N MUST BE COMPOSED OF FACTORS 2,3 & 5 BUT DOES NOT HAVE TO BE EVEN
+!C
+!C DEFINITION OF TRANSFORMS:
+!C -------------------------
+!C
+!C ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
+!C WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
+!C
+!C ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N))
+!C B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N))
+!C
+ IF(IFAX(10).NE.N) CALL SET99(TRIGS,IFAX,N)
+ NFAX=IFAX(1)
+ NX=N+1
+ IF (MOD(N,2).EQ.1) NX=N
+ NBLOX=1+(LOT-1)/64
+ NVEX=LOT-(NBLOX-1)*64
+ IF (ISIGN.EQ.-1) GO TO 300
+!C
+!C ISIGN=+1, SPECTRAL TO GRIDPOINT TRANSFORM
+!C -----------------------------------------
+ 100 CONTINUE
+ ISTART=1
+ DO 220 NB=1,NBLOX
+ IA=ISTART
+ I=ISTART
+ DO 110 J=1,NVEX
+ A(I+INC)=0.5*A(I)
+ I=I+JUMP
+ 110 CONTINUE
+ IF (MOD(N,2).EQ.1) GO TO 130
+ I=ISTART+N*INC
+ DO 120 J=1,NVEX
+ A(I)=0.5*A(I)
+ I=I+JUMP
+ 120 CONTINUE
+ 130 CONTINUE
+ IA=ISTART+INC
+ LA=1
+ IGO=+1
+!C
+ DO 160 K=1,NFAX
+ IFAC=IFAX(K+1)
+ IERR=-1
+ IF (IGO.EQ.-1) GO TO 140
+ CALL RPASSM(A(IA),A(IA+LA*INC),WORK(1),WORK(IFAC*LA+1),TRIGS,&
+ INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
+ GO TO 150
+ 140 CONTINUE
+ CALL RPASSM(WORK(1),WORK(LA+1),A(IA),A(IA+IFAC*LA*INC),TRIGS,&
+ 1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
+ 150 CONTINUE
+ IF (IERR.NE.0) GO TO 500
+ LA=IFAC*LA
+ IGO=-IGO
+ IA=ISTART+INC
+ 160 CONTINUE
+!C
+!C IF NECESSARY, COPY RESULTS BACK TO A
+!C ------------------------------------
+ IF (MOD(NFAX,2).EQ.0) GO TO 190
+ IBASE=1
+ JBASE=IA
+ DO 180 JJ=1,NVEX
+ I=IBASE
+ J=JBASE
+ DO 170 II=1,N
+ A(J)=WORK(I)
+ I=I+1
+ J=J+INC
+ 170 CONTINUE
+ IBASE=IBASE+NX
+ JBASE=JBASE+JUMP
+ 180 CONTINUE
+ 190 CONTINUE
+!C
+!C FILL IN CYCLIC BOUNDARY VALUES
+!C ------------------------------
+ IX=ISTART
+ IZ=ISTART+N*INC
+!DIR$ IVDEP
+ DO 200 J=1,NVEX
+ A(IX)=A(IZ)
+ A(IZ+INC)=A(IX+INC)
+ IX=IX+JUMP
+ IZ=IZ+JUMP
+ 200 CONTINUE
+!C
+ ISTART=ISTART+NVEX*JUMP
+ NVEX=64
+ 220 CONTINUE
+ RETURN
+!C
+!C ISIGN=-1, GRIDPOINT TO SPECTRAL TRANSFORM
+!C -----------------------------------------
+ 300 CONTINUE
+ ISTART=1
+ DO 410 NB=1,NBLOX
+ IA=ISTART+INC
+ LA=N
+ IGO=+1
+!C
+ DO 340 K=1,NFAX
+ IFAC=IFAX(NFAX+2-K)
+ LA=LA/IFAC
+ IERR=-1
+ IF (IGO.EQ.-1) GO TO 320
+ CALL QPASSM(A(IA),A(IA+IFAC*LA*INC),WORK(1),WORK(LA+1),TRIGS,&
+ INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
+ GO TO 330
+ 320 CONTINUE
+ CALL QPASSM(WORK(1),WORK(IFAC*LA+1),A(IA),A(IA+LA*INC),TRIGS,&
+ 1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
+ 330 CONTINUE
+ IF (IERR.NE.0) GO TO 500
+ IGO=-IGO
+ IA=ISTART+INC
+ 340 CONTINUE
+!C
+!C IF NECESSARY, COPY RESULTS BACK TO A
+!C ------------------------------------
+ IF (MOD(NFAX,2).EQ.0) GO TO 370
+ IBASE=1
+ JBASE=IA
+ DO 360 JJ=1,NVEX
+ I=IBASE
+ J=JBASE
+ DO 350 II=1,N
+ A(J)=WORK(I)
+ I=I+1
+ J=J+INC
+ 350 CONTINUE
+ IBASE=IBASE+NX
+ JBASE=JBASE+JUMP
+ 360 CONTINUE
+ 370 CONTINUE
+!C
+!C SHIFT A(0) & FILL IN ZERO IMAG PARTS
+!C ------------------------------------
+ IX=ISTART
+ DO 380 J=1,NVEX
+ A(IX)=A(IX+INC)
+ A(IX+INC)=0.0
+ IX=IX+JUMP
+ 380 CONTINUE
+ IF (MOD(N,2).EQ.1) GO TO 400
+ IZ=ISTART+(N+1)*INC
+ DO 390 J=1,NVEX
+ A(IZ)=0.0
+ IZ=IZ+JUMP
+ 390 CONTINUE
+ 400 CONTINUE
+!C
+ ISTART=ISTART+NVEX*JUMP
+ NVEX=64
+ 410 CONTINUE
+ RETURN
+!C
+!C ERROR MESSAGES
+!C --------------
+ 500 CONTINUE
+ GO TO (510,530,550) IERR
+ 510 CONTINUE
+ WRITE(6,520) NVEX
+ 520 FORMAT('1VECTOR LENGTH =',I4,', GREATER THAN 64')
+ GO TO 570
+ 530 CONTINUE
+ WRITE(6,540) IFAC
+ 540 FORMAT( '1FACTOR =',I3,', NOT CATERED FOR')
+ GO TO 570
+ 550 CONTINUE
+ WRITE(6,560) IFAC
+ 560 FORMAT('1FACTOR =',I3,', ONLY CATERED FOR IF LA*IFAC=N')
+ 570 CONTINUE
+ RETURN
+ END
diff --git a/Source/Fortran/jsppole.f90 b/Source/Fortran/jsppole.f90
new file mode 100644
index 0000000..5798ac4
--- /dev/null
+++ b/Source/Fortran/jsppole.f90
@@ -0,0 +1,182 @@
+! Copyright 1981-2016 ECMWF.
+!
+! This software is licensed under the terms of the Apache Licence
+! Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+!
+! In applying this licence, ECMWF does not waive the privileges and immunities
+! granted to it by virtue of its status as an intergovernmental organisation
+! nor does it submit to any jurisdiction.
+!
+
+ SUBROUTINE JSPPOLE(PSHUP,KNUMB,KTRUNC,OMARS,PXF)
+ IMPLICIT NONE
+!
+!---->
+!**** *JSPPOLE* - Calculates fourier coefficient for U or V at pole
+!
+! Purpose
+! -------
+!
+! Calculates fourier coefficient for first harmonic only
+! for U and V wind component at the pole.
+!
+! Interface
+! ---------
+!
+! CALL JSPPOLE(PSHUP,KNUMB,KTRUNC,OMARS,PXF)
+!
+! Input parameters
+! ----------------
+!
+! PSHUP - Unpacked harmonics field, unpacked
+! KNUMB - 1 for North Pole, otherwise South Pole
+! KTRUNC - Number (value) of the trucation
+! OMARS - .TRUE. if data is from MARS
+! PXF - Fourier coefficients (zero on input)
+!
+!
+! Output parameters
+! -----------------
+!
+! PXF(2) - Single fourier coefficient calculated
+!
+!
+! Common block usage
+! -----------------
+!
+! None.
+!
+!
+! Externals
+! ---------
+!
+! None.
+!
+!
+! Author
+! ------
+!
+! J.D.Chambers *ECMWF* Oct 1993
+!
+!
+! Modifications
+! -------------
+!
+! None.
+!
+!
+! Comments
+! --------
+!
+! Created from SPPOLE.
+! Changed to provide all parameters in the call, i.e. no common
+! blocks are used.
+!
+!
+! Method
+! ------
+!
+! None.
+!
+!
+! Reference
+! _________
+!
+! None.
+!
+!----<
+! _______________________________________________________
+!
+!* Section 0. Definition of variables.
+! _______________________________________________________
+!
+!* Prefix conventions for variable names
+!
+! Logical L (but not LP), global or common.
+! O, dummy argument
+! G, local variable
+! LP, parameter.
+! Character C, global or common.
+! H, dummy argument
+! Y (but not YP), local variable
+! YP, parameter.
+! Integer M and N, global or common.
+! K, dummy argument
+! I, local variable
+! J (but not JP), loop control
+! JP, parameter.
+! REAL A to F and Q to X, global or common.
+! P (but not PP), dummy argument
+! Z, local variable
+! PP, parameter.
+!
+! Dummy arguments
+!
+ COMPLEX PSHUP
+ INTEGER KNUMB
+ INTEGER KTRUNC
+ LOGICAL OMARS
+ COMPLEX PXF
+ DIMENSION PSHUP(*)
+ DIMENSION PXF(*)
+!
+! Local variables
+!
+ INTEGER I1, ITIN1, ITOUT1, JN
+ REAL Z1, Z2, ZNORM, ZP1, ZP2, ZPOL
+!
+! -----------------------------------------------------------
+!
+!* 1. Set initial values
+! ------------------
+!
+ 100 CONTINUE
+!
+ ITIN1 = KTRUNC + 1
+ ITOUT1 = KTRUNC
+!
+ ZPOL = 1.
+ IF (KNUMB .NE. 1) ZPOL = -1.0
+!
+ ZP1 = -1.0
+ ZP2 = -3.0 * ZPOL
+ I1 = ITIN1 + 1
+!
+!* 2. Change normalisation (if flagged as necessary)
+! --------------------
+!
+ 200 CONTINUE
+!
+ IF (OMARS) THEN
+ ZNORM = -SQRT(2.0)
+ ELSE
+ ZNORM = 1
+ ENDIF
+!
+!
+!* 3. Calculation
+! -----------
+!
+ 300 CONTINUE
+ PXF(2) = (0.0,0.0)
+!
+! Calculate the fourier coefficient for the first harmonic only.
+ DO 310 JN = 1,ITOUT1,2
+!
+ Z1 = SQRT( (2.0*JN + 1.0)/(2.0*JN*(JN + 1.0)) )
+ Z2 = SQRT( (2.0*(JN + 1.0) +1.0)/(2.0*(JN +1.0)*(JN +2.0)) )
+!
+ IF (JN .EQ. ITOUT1) Z2 = 0.0
+!
+ PXF(2) = PXF(2) +(Z1*ZP1*PSHUP(I1) +Z2*ZP2*PSHUP(I1+1))*ZNORM
+ ZP1 = ZP1 - 2.0*(JN + 1.0) - 1.0
+ ZP2 = ZP2 - (2.0*(JN + 2.0) + 1.0)*ZPOL
+ I1 = I1 + 2
+!
+ 310 CONTINUE
+!
+! -------------------------------------------------------------
+!
+ RETURN
+!
+ END
diff --git a/Source/Fortran/makefile_cray b/Source/Fortran/makefile_cray
index 6993d73..5cc2f7a 100644
--- a/Source/Fortran/makefile_cray
+++ b/Source/Fortran/makefile_cray
@@ -14,7 +14,7 @@
EXE = calc_etadot
-LIB = $(ECCODES_LIB) $(EMOSLIB)
+LIB = $(ECCODES_LIB)
INC = -I. -I$(ECCODES_INCLUDE_DIR)
FC = ftn
diff --git a/Source/Fortran/makefile_debug b/Source/Fortran/makefile_debug
index 93d4a09..f2e0cae 100644
--- a/Source/Fortran/makefile_debug
+++ b/Source/Fortran/makefile_debug
@@ -15,8 +15,7 @@
EXE = calc_etadot_debug.out
ECCODES_LIB= -Bstatic -leccodes_f90 -leccodes -Bdynamic -lm
-EMOSLIB=-lemosR64
-LIB = $(ECCODES_LIB) $(EMOSLIB)
+LIB = $(ECCODES_LIB)
ECCODES_INCLUDE_DIR=/usr/lib/x86_64-linux-gnu/fortran/gfortran-mod-15
INC = -I. -I$(ECCODES_INCLUDE_DIR)
diff --git a/Source/Fortran/makefile_ecgate b/Source/Fortran/makefile_ecgate
index e68960f..ce2be2a 100644
--- a/Source/Fortran/makefile_ecgate
+++ b/Source/Fortran/makefile_ecgate
@@ -14,7 +14,7 @@
EXE = calc_etadot_fast.out
-LIB = $(ECCODES_LIB) $(EMOSLIB)
+LIB = $(ECCODES_LIB)
INC = -I. -I$(ECCODES_INCLUDE_DIR)
FC=gfortran
diff --git a/Source/Fortran/makefile_fast b/Source/Fortran/makefile_fast
index 7b0a215..1537b6f 100644
--- a/Source/Fortran/makefile_fast
+++ b/Source/Fortran/makefile_fast
@@ -15,8 +15,7 @@
EXE = calc_etadot_fast.out
ECCODES_LIB = -Bstatic -leccodes_f90 -leccodes -Bdynamic -lm
-EMOSLIB=-lemosR64
-LIB = $(ECCODES_LIB) $(EMOSLIB)
+LIB = $(ECCODES_LIB)
ECCODES_INCLUDE_DIR=/usr/lib/x86_64-linux-gnu/fortran/gfortran-mod-15
#/usr/local/include/ #oldstable
diff --git a/Source/Fortran/qpassm.f90 b/Source/Fortran/qpassm.f90
new file mode 100644
index 0000000..b6efc6e
--- /dev/null
+++ b/Source/Fortran/qpassm.f90
@@ -0,0 +1,778 @@
+!c Copyright 1981-2016 ECMWF.
+!c
+!c This software is licensed under the terms of the Apache Licence
+!c Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+!c
+!c In applying this licence, ECMWF does not waive the privileges and immunities
+!c granted to it by virtue of its status as an intergovernmental organisation
+!c nor does it submit to any jurisdiction.
+!c
+
+ SUBROUTINE QPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,ILOT,N,IFAC,ILA,IERR)
+ REAL A(*),B(*),C(*),D(*),TRIGS(*)
+!C
+!C SUBROUTINE 'QPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART
+!C OF MULTIPLE REAL FFT (FOURIER ANALYSIS) ROUTINE
+!C
+!C A IS FIRST REAL INPUT VECTOR
+!C EQUIVALENCE B(1) WITH A(IFAC*ILA*INC1+1)
+!C C IS FIRST REAL OUTPUT VECTOR
+!C EQUIVALENCE D(1) WITH C(ILA*INC2+1)
+!C TRIGS IS A PRECALCULATED LIST OF SINES & COSINES
+!C INC1 IS THE ADDRESSING INCREMENT FOR A
+!C INC2 IS THE ADDRESSING INCREMENT FOR C
+!C INC3 IS THE INCREMENT BETWEEN INPUT VECTORS A
+!C INC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS C
+!C ILOT IS THE NUMBER OF VECTORS
+!C N IS THE LENGTH OF THE VECTORS
+!C IFAC IS THE CURRENT FACTOR OF N
+!C ILA = N/(PRODUCT OF FACTORS USED SO FAR)
+!C IERR IS AN ERROR INDICATOR:
+!C 0 - PASS COMPLETED WITHOUT ERROR
+!C 1 - ILOT GREATER THAN 64
+!C 2 - IFAC NOT CATERED FOR
+!C 3 - IFAC ONLY CATERED FOR IF ILA=N/IFAC
+!C
+!C-----------------------------------------------------------------------
+!C
+ SAVE SIN36, SIN72, QRT5, SIN60
+
+ DATA SIN36/0.587785252292473/,SIN72/0.951056516295154/,&
+ QRT5/0.559016994374947/,SIN60/0.866025403784437/
+
+ M=N/IFAC
+ IINK=ILA*INC1
+ JINK=ILA*INC2
+ IJUMP=(IFAC-1)*IINK
+ KSTOP=(N-IFAC)/(2*IFAC)
+
+ IBAD=1
+ IF (ILOT.GT.512) GO TO 910
+ IBASE=0
+ JBASE=0
+ IGO=IFAC-1
+ IF (IGO.EQ.7) IGO=6
+ IBAD=2
+ IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910
+ GO TO (200,300,400,500,600,800),IGO
+
+! CODING FOR FACTOR 2
+! -------------------
+ 200 CONTINUE
+ IA=1
+ IB=IA+IINK
+ JA=1
+ JB=JA+(2*M-ILA)*INC2
+
+ IF (ILA.EQ.M) GO TO 290
+
+ DO 220 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 210 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+A(IB+I)
+ C(JB+J)=A(IA+I)-A(IB+I)
+ I=I+INC3
+ J=J+INC4
+ 210 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 220 CONTINUE
+ JA=JA+JINK
+ JINK=2*JINK
+ JB=JB-JINK
+ IBASE=IBASE+IJUMP
+ IJUMP=2*IJUMP+IINK
+ IF (JA.EQ.JB) GO TO 260
+ DO 250 K=ILA,KSTOP,ILA
+ KB=K+K
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ JBASE=0
+ DO 240 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 230 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+(C1*A(IB+I)+S1*B(IB+I))
+ C(JB+J)=A(IA+I)-(C1*A(IB+I)+S1*B(IB+I))
+ D(JA+J)=(C1*B(IB+I)-S1*A(IB+I))+B(IA+I)
+ D(JB+J)=(C1*B(IB+I)-S1*A(IB+I))-B(IA+I)
+ I=I+INC3
+ J=J+INC4
+ 230 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 240 CONTINUE
+ IBASE=IBASE+IJUMP
+ JA=JA+JINK
+ JB=JB-JINK
+ 250 CONTINUE
+ IF (JA.GT.JB) GO TO 900
+ 260 CONTINUE
+ JBASE=0
+ DO 280 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 270 IJK=1,ILOT
+ C(JA+J)=A(IA+I)
+ D(JA+J)=-A(IB+I)
+ I=I+INC3
+ J=J+INC4
+ 270 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 280 CONTINUE
+ GO TO 900
+
+ 290 CONTINUE
+ Z=1.0/FLOAT(N)
+ DO 294 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 292 IJK=1,ILOT
+ C(JA+J)=Z*(A(IA+I)+A(IB+I))
+ C(JB+J)=Z*(A(IA+I)-A(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 292 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 294 CONTINUE
+ GO TO 900
+
+! CODING FOR FACTOR 3
+! -------------------
+ 300 CONTINUE
+ IA=1
+ IB=IA+IINK
+ IC=IB+IINK
+ JA=1
+ JB=JA+(2*M-ILA)*INC2
+ JC=JB
+
+ IF (ILA.EQ.M) GO TO 390
+
+ DO 320 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 310 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
+ C(JB+J)=A(IA+I)-0.5*(A(IB+I)+A(IC+I))
+ D(JB+J)=SIN60*(A(IC+I)-A(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 310 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 320 CONTINUE
+ JA=JA+JINK
+ JINK=2*JINK
+ JB=JB+JINK
+ JC=JC-JINK
+ IBASE=IBASE+IJUMP
+ IJUMP=2*IJUMP+IINK
+ IF (JA.EQ.JC) GO TO 360
+ DO 350 K=ILA,KSTOP,ILA
+ KB=K+K
+ KC=KB+KB
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ C2=TRIGS(KC+1)
+ S2=TRIGS(KC+2)
+ JBASE=0
+ DO 340 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 330 IJK=1,ILOT
+ A1=(C1*A(IB+I)+S1*B(IB+I))+(C2*A(IC+I)+S2*B(IC+I))
+ B1=(C1*B(IB+I)-S1*A(IB+I))+(C2*B(IC+I)-S2*A(IC+I))
+ A2=A(IA+I)-0.5*A1
+ B2=B(IA+I)-0.5*B1
+ A3=SIN60*((C1*A(IB+I)+S1*B(IB+I))-(C2*A(IC+I)+S2*B(IC+I)))
+ B3=SIN60*((C1*B(IB+I)-S1*A(IB+I))-(C2*B(IC+I)-S2*A(IC+I)))
+ C(JA+J)=A(IA+I)+A1
+ D(JA+J)=B(IA+I)+B1
+ C(JB+J)=A2+B3
+ D(JB+J)=B2-A3
+ C(JC+J)=A2-B3
+ D(JC+J)=-(B2+A3)
+ I=I+INC3
+ J=J+INC4
+ 330 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 340 CONTINUE
+ IBASE=IBASE+IJUMP
+ JA=JA+JINK
+ JB=JB+JINK
+ JC=JC-JINK
+ 350 CONTINUE
+ IF (JA.GT.JC) GO TO 900
+ 360 CONTINUE
+ JBASE=0
+ DO 380 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 370 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+0.5*(A(IB+I)-A(IC+I))
+ D(JA+J)=-SIN60*(A(IB+I)+A(IC+I))
+ C(JB+J)=A(IA+I)-(A(IB+I)-A(IC+I))
+ I=I+INC3
+ J=J+INC4
+ 370 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 380 CONTINUE
+ GO TO 900
+
+ 390 CONTINUE
+ Z=1.0/FLOAT(N)
+ ZSIN60=Z*SIN60
+ DO 394 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 392 IJK=1,ILOT
+ C(JA+J)=Z*(A(IA+I)+(A(IB+I)+A(IC+I)))
+ C(JB+J)=Z*(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))
+ D(JB+J)=ZSIN60*(A(IC+I)-A(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 392 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 394 CONTINUE
+ GO TO 900
+
+! CODING FOR FACTOR 4
+! -------------------
+ 400 CONTINUE
+ IA=1
+ IB=IA+IINK
+ IC=IB+IINK
+ ID=IC+IINK
+ JA=1
+ JB=JA+(2*M-ILA)*INC2
+ JC=JB+2*M*INC2
+ JD=JB
+
+ IF (ILA.EQ.M) GO TO 490
+
+ DO 420 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 410 IJK=1,ILOT
+ C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
+ C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))
+ C(JB+J)=A(IA+I)-A(IC+I)
+ D(JB+J)=A(ID+I)-A(IB+I)
+ I=I+INC3
+ J=J+INC4
+ 410 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 420 CONTINUE
+ JA=JA+JINK
+ JINK=2*JINK
+ JB=JB+JINK
+ JC=JC-JINK
+ JD=JD-JINK
+ IBASE=IBASE+IJUMP
+ IJUMP=2*IJUMP+IINK
+ IF (JB.EQ.JC) GO TO 460
+ DO 450 K=ILA,KSTOP,ILA
+ KB=K+K
+ KC=KB+KB
+ KD=KC+KB
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ C2=TRIGS(KC+1)
+ S2=TRIGS(KC+2)
+ C3=TRIGS(KD+1)
+ S3=TRIGS(KD+2)
+ JBASE=0
+ DO 440 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 430 IJK=1,ILOT
+ A0=A(IA+I)+(C2*A(IC+I)+S2*B(IC+I))
+ A2=A(IA+I)-(C2*A(IC+I)+S2*B(IC+I))
+ A1=(C1*A(IB+I)+S1*B(IB+I))+(C3*A(ID+I)+S3*B(ID+I))
+ A3=(C1*A(IB+I)+S1*B(IB+I))-(C3*A(ID+I)+S3*B(ID+I))
+ B0=B(IA+I)+(C2*B(IC+I)-S2*A(IC+I))
+ B2=B(IA+I)-(C2*B(IC+I)-S2*A(IC+I))
+ B1=(C1*B(IB+I)-S1*A(IB+I))+(C3*B(ID+I)-S3*A(ID+I))
+ B3=(C1*B(IB+I)-S1*A(IB+I))-(C3*B(ID+I)-S3*A(ID+I))
+ C(JA+J)=A0+A1
+ C(JC+J)=A0-A1
+ D(JA+J)=B0+B1
+ D(JC+J)=B1-B0
+ C(JB+J)=A2+B3
+ C(JD+J)=A2-B3
+ D(JB+J)=B2-A3
+ D(JD+J)=-(B2+A3)
+ I=I+INC3
+ J=J+INC4
+ 430 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 440 CONTINUE
+ IBASE=IBASE+IJUMP
+ JA=JA+JINK
+ JB=JB+JINK
+ JC=JC-JINK
+ JD=JD-JINK
+ 450 CONTINUE
+ IF (JB.GT.JC) GO TO 900
+ 460 CONTINUE
+ SIN45=SQRT(0.5)
+ JBASE=0
+ DO 480 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 470 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+SIN45*(A(IB+I)-A(ID+I))
+ C(JB+J)=A(IA+I)-SIN45*(A(IB+I)-A(ID+I))
+ D(JA+J)=-A(IC+I)-SIN45*(A(IB+I)+A(ID+I))
+ D(JB+J)=A(IC+I)-SIN45*(A(IB+I)+A(ID+I))
+ I=I+INC3
+ J=J+INC4
+ 470 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 480 CONTINUE
+ GO TO 900
+
+ 490 CONTINUE
+ Z=1.0/FLOAT(N)
+ DO 494 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 492 IJK=1,ILOT
+ C(JA+J)=Z*((A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)))
+ C(JC+J)=Z*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))
+ C(JB+J)=Z*(A(IA+I)-A(IC+I))
+ D(JB+J)=Z*(A(ID+I)-A(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 492 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 494 CONTINUE
+ GO TO 900
+
+! CODING FOR FACTOR 5
+! -------------------
+ 500 CONTINUE
+ IA=1
+ IB=IA+IINK
+ IC=IB+IINK
+ ID=IC+IINK
+ IE=ID+IINK
+ JA=1
+ JB=JA+(2*M-ILA)*INC2
+ JC=JB+2*M*INC2
+ JD=JC
+ JE=JB
+
+ IF (ILA.EQ.M) GO TO 590
+
+ DO 520 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 510 IJK=1,ILOT
+ A1=A(IB+I)+A(IE+I)
+ A3=A(IB+I)-A(IE+I)
+ A2=A(IC+I)+A(ID+I)
+ A4=A(IC+I)-A(ID+I)
+ A5=A(IA+I)-0.25*(A1+A2)
+ A6=QRT5*(A1-A2)
+ C(JA+J)=A(IA+I)+(A1+A2)
+ C(JB+J)=A5+A6
+ C(JC+J)=A5-A6
+ D(JB+J)=-SIN72*A3-SIN36*A4
+ D(JC+J)=-SIN36*A3+SIN72*A4
+ I=I+INC3
+ J=J+INC4
+ 510 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 520 CONTINUE
+ JA=JA+JINK
+ JINK=2*JINK
+ JB=JB+JINK
+ JC=JC+JINK
+ JD=JD-JINK
+ JE=JE-JINK
+ IBASE=IBASE+IJUMP
+ IJUMP=2*IJUMP+IINK
+ IF (JB.EQ.JD) GO TO 560
+ DO 550 K=ILA,KSTOP,ILA
+ KB=K+K
+ KC=KB+KB
+ KD=KC+KB
+ KE=KD+KB
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ C2=TRIGS(KC+1)
+ S2=TRIGS(KC+2)
+ C3=TRIGS(KD+1)
+ S3=TRIGS(KD+2)
+ C4=TRIGS(KE+1)
+ S4=TRIGS(KE+2)
+ JBASE=0
+ DO 540 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 530 IJK=1,ILOT
+ A1=(C1*A(IB+I)+S1*B(IB+I))+(C4*A(IE+I)+S4*B(IE+I))
+ A3=(C1*A(IB+I)+S1*B(IB+I))-(C4*A(IE+I)+S4*B(IE+I))
+ A2=(C2*A(IC+I)+S2*B(IC+I))+(C3*A(ID+I)+S3*B(ID+I))
+ A4=(C2*A(IC+I)+S2*B(IC+I))-(C3*A(ID+I)+S3*B(ID+I))
+ B1=(C1*B(IB+I)-S1*A(IB+I))+(C4*B(IE+I)-S4*A(IE+I))
+ B3=(C1*B(IB+I)-S1*A(IB+I))-(C4*B(IE+I)-S4*A(IE+I))
+ B2=(C2*B(IC+I)-S2*A(IC+I))+(C3*B(ID+I)-S3*A(ID+I))
+ B4=(C2*B(IC+I)-S2*A(IC+I))-(C3*B(ID+I)-S3*A(ID+I))
+ A5=A(IA+I)-0.25*(A1+A2)
+ A6=QRT5*(A1-A2)
+ B5=B(IA+I)-0.25*(B1+B2)
+ B6=QRT5*(B1-B2)
+ A10=A5+A6
+ A20=A5-A6
+ B10=B5+B6
+ B20=B5-B6
+ A11=SIN72*B3+SIN36*B4
+ A21=SIN36*B3-SIN72*B4
+ B11=SIN72*A3+SIN36*A4
+ B21=SIN36*A3-SIN72*A4
+ C(JA+J)=A(IA+I)+(A1+A2)
+ C(JB+J)=A10+A11
+ C(JE+J)=A10-A11
+ C(JC+J)=A20+A21
+ C(JD+J)=A20-A21
+ D(JA+J)=B(IA+I)+(B1+B2)
+ D(JB+J)=B10-B11
+ D(JE+J)=-(B10+B11)
+ D(JC+J)=B20-B21
+ D(JD+J)=-(B20+B21)
+ I=I+INC3
+ J=J+INC4
+ 530 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 540 CONTINUE
+ IBASE=IBASE+IJUMP
+ JA=JA+JINK
+ JB=JB+JINK
+ JC=JC+JINK
+ JD=JD-JINK
+ JE=JE-JINK
+ 550 CONTINUE
+ IF (JB.GT.JD) GO TO 900
+ 560 CONTINUE
+ JBASE=0
+ DO 580 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 570 IJK=1,ILOT
+ A1=A(IB+I)+A(IE+I)
+ A3=A(IB+I)-A(IE+I)
+ A2=A(IC+I)+A(ID+I)
+ A4=A(IC+I)-A(ID+I)
+ A5=A(IA+I)+0.25*(A3-A4)
+ A6=QRT5*(A3+A4)
+ C(JA+J)=A5+A6
+ C(JB+J)=A5-A6
+ C(JC+J)=A(IA+I)-(A3-A4)
+ D(JA+J)=-SIN36*A1-SIN72*A2
+ D(JB+J)=-SIN72*A1+SIN36*A2
+ I=I+INC3
+ J=J+INC4
+ 570 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 580 CONTINUE
+ GO TO 900
+
+ 590 CONTINUE
+ Z=1.0/FLOAT(N)
+ ZQRT5=Z*QRT5
+ ZSIN36=Z*SIN36
+ ZSIN72=Z*SIN72
+ DO 594 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 592 IJK=1,ILOT
+ A1=A(IB+I)+A(IE+I)
+ A3=A(IB+I)-A(IE+I)
+ A2=A(IC+I)+A(ID+I)
+ A4=A(IC+I)-A(ID+I)
+ A5=Z*(A(IA+I)-0.25*(A1+A2))
+ A6=ZQRT5*(A1-A2)
+ C(JA+J)=Z*(A(IA+I)+(A1+A2))
+ C(JB+J)=A5+A6
+ C(JC+J)=A5-A6
+ D(JB+J)=-ZSIN72*A3-ZSIN36*A4
+ D(JC+J)=-ZSIN36*A3+ZSIN72*A4
+ I=I+INC3
+ J=J+INC4
+ 592 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 594 CONTINUE
+ GO TO 900
+
+! CODING FOR FACTOR 6
+! -------------------
+ 600 CONTINUE
+ IA=1
+ IB=IA+IINK
+ IC=IB+IINK
+ ID=IC+IINK
+ IE=ID+IINK
+ IF=IE+IINK
+ JA=1
+ JB=JA+(2*M-ILA)*INC2
+ JC=JB+2*M*INC2
+ JD=JC+2*M*INC2
+ JE=JC
+ JF=JB
+
+ IF (ILA.EQ.M) GO TO 690
+
+ DO 620 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 610 IJK=1,ILOT
+ A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I))
+ C(JA+J)=(A(IA+I)+A(ID+I))+A11
+ C(JC+J)=(A(IA+I)+A(ID+I)-0.5*A11)
+ D(JC+J)=SIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I)))
+ A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I))
+ C(JB+J)=(A(IA+I)-A(ID+I))-0.5*A11
+ D(JB+J)=SIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I)))
+ C(JD+J)=(A(IA+I)-A(ID+I))+A11
+ I=I+INC3
+ J=J+INC4
+ 610 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 620 CONTINUE
+ JA=JA+JINK
+ JINK=2*JINK
+ JB=JB+JINK
+ JC=JC+JINK
+ JD=JD-JINK
+ JE=JE-JINK
+ JF=JF-JINK
+ IBASE=IBASE+IJUMP
+ IJUMP=2*IJUMP+IINK
+ IF (JC.EQ.JD) GO TO 660
+ DO 650 K=ILA,KSTOP,ILA
+ KB=K+K
+ KC=KB+KB
+ KD=KC+KB
+ KE=KD+KB
+ KF=KE+KB
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ C2=TRIGS(KC+1)
+ S2=TRIGS(KC+2)
+ C3=TRIGS(KD+1)
+ S3=TRIGS(KD+2)
+ C4=TRIGS(KE+1)
+ S4=TRIGS(KE+2)
+ C5=TRIGS(KF+1)
+ S5=TRIGS(KF+2)
+ JBASE=0
+ DO 640 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 630 IJK=1,ILOT
+ A1=C1*A(IB+I)+S1*B(IB+I)
+ B1=C1*B(IB+I)-S1*A(IB+I)
+ A2=C2*A(IC+I)+S2*B(IC+I)
+ B2=C2*B(IC+I)-S2*A(IC+I)
+ A3=C3*A(ID+I)+S3*B(ID+I)
+ B3=C3*B(ID+I)-S3*A(ID+I)
+ A4=C4*A(IE+I)+S4*B(IE+I)
+ B4=C4*B(IE+I)-S4*A(IE+I)
+ A5=C5*A(IF+I)+S5*B(IF+I)
+ B5=C5*B(IF+I)-S5*A(IF+I)
+ A11=(A2+A5)+(A1+A4)
+ A20=(A(IA+I)+A3)-0.5*A11
+ A21=SIN60*((A2+A5)-(A1+A4))
+ B11=(B2+B5)+(B1+B4)
+ B20=(B(IA+I)+B3)-0.5*B11
+ B21=SIN60*((B2+B5)-(B1+B4))
+ C(JA+J)=(A(IA+I)+A3)+A11
+ D(JA+J)=(B(IA+I)+B3)+B11
+ C(JC+J)=A20-B21
+ D(JC+J)=A21+B20
+ C(JE+J)=A20+B21
+ D(JE+J)=A21-B20
+ A11=(A2-A5)+(A4-A1)
+ A20=(A(IA+I)-A3)-0.5*A11
+ A21=SIN60*((A4-A1)-(A2-A5))
+ B11=(B5-B2)-(B4-B1)
+ B20=(B3-B(IA+I))-0.5*B11
+ B21=SIN60*((B5-B2)+(B4-B1))
+ C(JB+J)=A20-B21
+ D(JB+J)=A21-B20
+ C(JD+J)=A11+(A(IA+I)-A3)
+ D(JD+J)=B11+(B3-B(IA+I))
+ C(JF+J)=A20+B21
+ D(JF+J)=A21+B20
+ I=I+INC3
+ J=J+INC4
+ 630 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 640 CONTINUE
+ IBASE=IBASE+IJUMP
+ JA=JA+JINK
+ JB=JB+JINK
+ JC=JC+JINK
+ JD=JD-JINK
+ JE=JE-JINK
+ JF=JF-JINK
+ 650 CONTINUE
+ IF (JC.GT.JD) GO TO 900
+ 660 CONTINUE
+ JBASE=0
+ DO 680 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 670 IJK=1,ILOT
+ C(JA+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))+ SIN60*(A(IB+I)-A(IF+I))
+ D(JA+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))-SIN60*(A(IC+I)+A(IE+I))
+ C(JB+J)=A(IA+I)-(A(IC+I)-A(IE+I))
+ D(JB+J)=A(ID+I)-(A(IB+I)+A(IF+I))
+ C(JC+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))-SIN60*(A(IB+I)-A(IF+I))
+ D(JC+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))+SIN60*(A(IC+I)+A(IE+I))
+ I=I+INC3
+ J=J+INC4
+ 670 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 680 CONTINUE
+ GO TO 900
+
+ 690 CONTINUE
+ Z=1.0/FLOAT(N)
+ ZSIN60=Z*SIN60
+ DO 694 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 692 IJK=1,ILOT
+ A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I))
+ C(JA+J)=Z*((A(IA+I)+A(ID+I))+A11)
+ C(JC+J)=Z*((A(IA+I)+A(ID+I))-0.5*A11)
+ D(JC+J)=ZSIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I)))
+ A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I))
+ C(JB+J)=Z*((A(IA+I)-A(ID+I))-0.5*A11)
+ D(JB+J)=ZSIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I)))
+ C(JD+J)=Z*((A(IA+I)-A(ID+I))+A11)
+ I=I+INC3
+ J=J+INC4
+ 692 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 694 CONTINUE
+ GO TO 900
+
+! CODING FOR FACTOR 8
+! -------------------
+ 800 CONTINUE
+ IBAD=3
+ IF (ILA.NE.M) GO TO 910
+ IA=1
+ IB=IA+IINK
+ IC=IB+IINK
+ ID=IC+IINK
+ IE=ID+IINK
+ IF=IE+IINK
+ IG=IF+IINK
+ IH=IG+IINK
+ JA=1
+ JB=JA+ILA*INC2
+ JC=JB+2*M*INC2
+ JD=JC+2*M*INC2
+ JE=JD+2*M*INC2
+ Z=1.0/FLOAT(N)
+ ZSIN45=Z*SQRT(0.5)
+
+ DO 820 JL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 810 IJK=1,ILOT
+ C(JA+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))+&
+ ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I))))
+ C(JE+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))-&
+ ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I))))
+ C(JC+J)=Z*((A(IA+I)+A(IE+I))-(A(IC+I)+A(IG+I)))
+ D(JC+J)=Z*((A(ID+I)+A(IH+I))-(A(IB+I)+A(IF+I)))
+ C(JB+J)=Z*(A(IA+I)-A(IE+I))&
+ +ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I)))
+ C(JD+J)=Z*(A(IA+I)-A(IE+I))&
+ -ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I)))
+ D(JB+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))&
+ +Z*(A(IG+I)-A(IC+I))
+ D(JD+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))&
+ -Z*(A(IG+I)-A(IC+I))
+ I=I+INC3
+ J=J+INC4
+ 810 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 820 CONTINUE
+
+! RETURN
+! ------
+ 900 CONTINUE
+ IBAD=0
+ 910 CONTINUE
+ IERR=IBAD
+ RETURN
+ END
diff --git a/Source/Fortran/rpassm.f90 b/Source/Fortran/rpassm.f90
new file mode 100644
index 0000000..cdd642f
--- /dev/null
+++ b/Source/Fortran/rpassm.f90
@@ -0,0 +1,760 @@
+!c Copyright 1981-2016 ECMWF.
+!c
+!c This software is licensed under the terms of the Apache Licence
+!c Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+!c
+!c In applying this licence, ECMWF does not waive the privileges and immunities
+!c granted to it by virtue of its status as an intergovernmental organisation
+!c nor does it submit to any jurisdiction.
+!c
+
+ SUBROUTINE RPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,ILOT,N,IFAC,ILA,IERR)
+ REAL A(*),B(*),C(*),D(*),TRIGS(*)
+!C
+!C SUBROUTINE 'RPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART
+!C OF MULTIPLE REAL FFT (FOURIER SYNTHESIS) ROUTINE
+!C
+!C A IS FIRST REAL INPUT VECTOR
+!C EQUIVALENCE B(1) WITH A (ILA*INC1+1)
+!C C IS FIRST REAL OUTPUT VECTOR
+!C EQUIVALENCE D(1) WITH C(IFAC*ILA*INC2+1)
+!C TRIGS IS A PRECALCULATED LIST OF SINES & COSINES
+!C INC1 IS THE ADDRESSING INCREMENT FOR A
+!C INC2 IS THE ADDRESSING INCREMENT FOR C
+!C INC3 IS THE INCREMENT BETWEEN INPUT VECTORS A
+!C INC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS C
+!C ILOT IS THE NUMBER OF VECTORS
+!C N IS THE LENGTH OF THE VECTORS
+!C IFAC IS THE CURRENT FACTOR OF N
+!C ILA IS THE PRODUCT OF PREVIOUS FACTORS
+!C IERR IS AN ERROR INDICATOR:
+!C 0 - PASS COMPLETED WITHOUT ERROR
+!C 1 - ILOT GREATER THAN 64
+!C 2 - IFAC NOT CATERED FOR
+!C 3 - IFAC ONLY CATERED FOR IF ILA=N/IFAC
+!C
+!C-----------------------------------------------------------------------
+!C
+ REAL A10(512),A11(512),&
+ A20(512),A21(512),&
+ B10(512),B11(512),&
+ B20(512),B21(512)
+!C
+ SAVE SIN36, SIN72, QRT5, SIN60
+!C
+ DATA SIN36/0.587785252292473/,SIN72/0.951056516295154/,&
+ QRT5/0.559016994374947/,SIN60/0.866025403784437/
+!C
+ M=N/IFAC
+ IINK=ILA*INC1
+ JINK=ILA*INC2
+ JUMP=(IFAC-1)*JINK
+ KSTOP=(N-IFAC)/(2*IFAC)
+!C
+ IBAD=1
+ IF (ILOT.GT.512) GO TO 910
+ IBASE=0
+ JBASE=0
+ IGO=IFAC-1
+ IF (IGO.EQ.7) IGO=6
+ IBAD=2
+ IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910
+ GO TO (200,300,400,500,600,800),IGO
+!C
+!C CODING FOR FACTOR 2
+!C -------------------
+ 200 CONTINUE
+ IA=1
+ IB=IA+(2*M-ILA)*INC1
+ JA=1
+ JB=JA+JINK
+!C
+ IF (ILA.EQ.M) GO TO 290
+!C
+ DO 220 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 210 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+A(IB+I)
+ C(JB+J)=A(IA+I)-A(IB+I)
+ I=I+INC3
+ J=J+INC4
+ 210 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 220 CONTINUE
+ IA=IA+IINK
+ IINK=2*IINK
+ IB=IB-IINK
+ IBASE=0
+ JBASE=JBASE+JUMP
+ JUMP=2*JUMP+JINK
+ IF (IA.EQ.IB) GO TO 260
+ DO 250 K=ILA,KSTOP,ILA
+ KB=K+K
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ IBASE=0
+ DO 240 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 230 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+A(IB+I)
+ D(JA+J)=B(IA+I)-B(IB+I)
+ C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)+B(IB+I))
+ D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)+B(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 230 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 240 CONTINUE
+ IA=IA+IINK
+ IB=IB-IINK
+ JBASE=JBASE+JUMP
+ 250 CONTINUE
+ IF (IA.GT.IB) GO TO 900
+ 260 CONTINUE
+ IBASE=0
+ DO 280 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 270 IJK=1,ILOT
+ C(JA+J)=A(IA+I)
+ C(JB+J)=-B(IA+I)
+ I=I+INC3
+ J=J+INC4
+ 270 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 280 CONTINUE
+ GO TO 900
+!C
+ 290 CONTINUE
+ DO 294 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 292 IJK=1,ILOT
+ C(JA+J)=2.0*(A(IA+I)+A(IB+I))
+ C(JB+J)=2.0*(A(IA+I)-A(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 292 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 294 CONTINUE
+ GO TO 900
+!C
+!C CODING FOR FACTOR 3
+!C -------------------
+ 300 CONTINUE
+ IA=1
+ IB=IA+(2*M-ILA)*INC1
+ IC=IB
+ JA=1
+ JB=JA+JINK
+ JC=JB+JINK
+!C
+ IF (ILA.EQ.M) GO TO 390
+!C
+ DO 320 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 310 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+A(IB+I)
+ C(JB+J)=(A(IA+I)-0.5*A(IB+I))-(SIN60*(B(IB+I)))
+ C(JC+J)=(A(IA+I)-0.5*A(IB+I))+(SIN60*(B(IB+I)))
+ I=I+INC3
+ J=J+INC4
+ 310 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 320 CONTINUE
+ IA=IA+IINK
+ IINK=2*IINK
+ IB=IB+IINK
+ IC=IC-IINK
+ JBASE=JBASE+JUMP
+ JUMP=2*JUMP+JINK
+ IF (IA.EQ.IC) GO TO 360
+ DO 350 K=ILA,KSTOP,ILA
+ KB=K+K
+ KC=KB+KB
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ C2=TRIGS(KC+1)
+ S2=TRIGS(KC+2)
+ IBASE=0
+ DO 340 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 330 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
+ D(JA+J)=B(IA+I)+(B(IB+I)-B(IC+I))
+ C(JB+J)=C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)+B(IC+I))))&
+ -S1*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
+ D(JB+J)=S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)+B(IC+I))))&
+ +C1*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))))
+ C(JC+J)=C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)+B(IC+I))))&
+ -S2*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
+ D(JC+J)=S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)+B(IC+I))))&
+ +C2*((B(IA+I)-0.5*(B(IB+I)-B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))))
+ I=I+INC3
+ J=J+INC4
+ 330 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 340 CONTINUE
+ IA=IA+IINK
+ IB=IB+IINK
+ IC=IC-IINK
+ JBASE=JBASE+JUMP
+ 350 CONTINUE
+ IF (IA.GT.IC) GO TO 900
+ 360 CONTINUE
+ IBASE=0
+ DO 380 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 370 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+A(IB+I)
+ C(JB+J)=(0.5*A(IA+I)-A(IB+I))-(SIN60*B(IA+I))
+ C(JC+J)=-(0.5*A(IA+I)-A(IB+I))-(SIN60*B(IA+I))
+ I=I+INC3
+ J=J+INC4
+ 370 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 380 CONTINUE
+ GO TO 900
+!C
+ 390 CONTINUE
+ SSIN60=2.0*SIN60
+ DO 394 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 392 IJK=1,ILOT
+ C(JA+J)=2.0*(A(IA+I)+A(IB+I))
+ C(JB+J)=(2.0*A(IA+I)-A(IB+I))-(SSIN60*B(IB+I))
+ C(JC+J)=(2.0*A(IA+I)-A(IB+I))+(SSIN60*B(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 392 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 394 CONTINUE
+ GO TO 900
+!C
+!C CODING FOR FACTOR 4
+!C -------------------
+ 400 CONTINUE
+ IA=1
+ IB=IA+(2*M-ILA)*INC1
+ IC=IB+2*M*INC1
+ ID=IB
+ JA=1
+ JB=JA+JINK
+ JC=JB+JINK
+ JD=JC+JINK
+!C
+ IF (ILA.EQ.M) GO TO 490
+!C
+ DO 420 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 410 IJK=1,ILOT
+ C(JA+J)=(A(IA+I)+A(IC+I))+A(IB+I)
+ C(JB+J)=(A(IA+I)-A(IC+I))-B(IB+I)
+ C(JC+J)=(A(IA+I)+A(IC+I))-A(IB+I)
+ C(JD+J)=(A(IA+I)-A(IC+I))+B(IB+I)
+ I=I+INC3
+ J=J+INC4
+ 410 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 420 CONTINUE
+ IA=IA+IINK
+ IINK=2*IINK
+ IB=IB+IINK
+ IC=IC-IINK
+ ID=ID-IINK
+ JBASE=JBASE+JUMP
+ JUMP=2*JUMP+JINK
+ IF (IB.EQ.IC) GO TO 460
+ DO 450 K=ILA,KSTOP,ILA
+ KB=K+K
+ KC=KB+KB
+ KD=KC+KB
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ C2=TRIGS(KC+1)
+ S2=TRIGS(KC+2)
+ C3=TRIGS(KD+1)
+ S3=TRIGS(KD+2)
+ IBASE=0
+ DO 440 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 430 IJK=1,ILOT
+ C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))
+ D(JA+J)=(B(IA+I)-B(IC+I))+(B(IB+I)-B(ID+I))
+ C(JC+J)=C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))&
+ -S2*((B(IA+I)-B(IC+I))-(B(IB+I)-B(ID+I)))
+ D(JC+J)=S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))&
+ +C2*((B(IA+I)-B(IC+I))-(B(IB+I)-B(ID+I)))
+ C(JB+J)=C1*((A(IA+I)-A(IC+I))-(B(IB+I)+B(ID+I)))&
+ -S1*((B(IA+I)+B(IC+I))+(A(IB+I)-A(ID+I)))
+ D(JB+J)=S1*((A(IA+I)-A(IC+I))-(B(IB+I)+B(ID+I)))&
+ +C1*((B(IA+I)+B(IC+I))+(A(IB+I)-A(ID+I)))
+ C(JD+J)=C3*((A(IA+I)-A(IC+I))+(B(IB+I)+B(ID+I)))&
+ -S3*((B(IA+I)+B(IC+I))-(A(IB+I)-A(ID+I)))
+ D(JD+J)=S3*((A(IA+I)-A(IC+I))+(B(IB+I)+B(ID+I)))&
+ +C3*((B(IA+I)+B(IC+I))-(A(IB+I)-A(ID+I)))
+ I=I+INC3
+ J=J+INC4
+ 430 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 440 CONTINUE
+ IA=IA+IINK
+ IB=IB+IINK
+ IC=IC-IINK
+ ID=ID-IINK
+ JBASE=JBASE+JUMP
+ 450 CONTINUE
+ IF (IB.GT.IC) GO TO 900
+ 460 CONTINUE
+ IBASE=0
+ SIN45=SQRT(0.5)
+ DO 480 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 470 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+A(IB+I)
+ C(JB+J)=SIN45*((A(IA+I)-A(IB+I))-(B(IA+I)+B(IB+I)))
+ C(JC+J)=B(IB+I)-B(IA+I)
+ C(JD+J)=-SIN45*((A(IA+I)-A(IB+I))+(B(IA+I)+B(IB+I)))
+ I=I+INC3
+ J=J+INC4
+ 470 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 480 CONTINUE
+ GO TO 900
+!C
+ 490 CONTINUE
+ DO 494 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 492 IJK=1,ILOT
+ C(JA+J)=2.0*((A(IA+I)+A(IC+I))+A(IB+I))
+ C(JB+J)=2.0*((A(IA+I)-A(IC+I))-B(IB+I))
+ C(JC+J)=2.0*((A(IA+I)+A(IC+I))-A(IB+I))
+ C(JD+J)=2.0*((A(IA+I)-A(IC+I))+B(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 492 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 494 CONTINUE
+ GO TO 900
+!C
+!C CODING FOR FACTOR 5
+!C -------------------
+ 500 CONTINUE
+ IA=1
+ IB=IA+(2*M-ILA)*INC1
+ IC=IB+2*M*INC1
+ ID=IC
+ IE=IB
+ JA=1
+ JB=JA+JINK
+ JC=JB+JINK
+ JD=JC+JINK
+ JE=JD+JINK
+!C
+ IF (ILA.EQ.M) GO TO 590
+!C
+ DO 520 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 510 IJK=1,ILOT
+ C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
+ C(JB+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))+QRT5*(A(IB+I)-A(IC+I)))&
+ -(SIN72*B(IB+I)+SIN36*B(IC+I))
+ C(JC+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))-QRT5*(A(IB+I)-A(IC+I)))&
+ -(SIN36*B(IB+I)-SIN72*B(IC+I))
+ C(JD+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))-QRT5*(A(IB+I)-A(IC+I)))&
+ +(SIN36*B(IB+I)-SIN72*B(IC+I))
+ C(JE+J)=((A(IA+I)-0.25*(A(IB+I)+A(IC+I)))+QRT5*(A(IB+I)-A(IC+I)))&
+ +(SIN72*B(IB+I)+SIN36*B(IC+I))
+ I=I+INC3
+ J=J+INC4
+ 510 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 520 CONTINUE
+ IA=IA+IINK
+ IINK=2*IINK
+ IB=IB+IINK
+ IC=IC+IINK
+ ID=ID-IINK
+ IE=IE-IINK
+ JBASE=JBASE+JUMP
+ JUMP=2*JUMP+JINK
+ IF (IB.EQ.ID) GO TO 560
+ DO 550 K=ILA,KSTOP,ILA
+ KB=K+K
+ KC=KB+KB
+ KD=KC+KB
+ KE=KD+KB
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ C2=TRIGS(KC+1)
+ S2=TRIGS(KC+2)
+ C3=TRIGS(KD+1)
+ S3=TRIGS(KD+2)
+ C4=TRIGS(KE+1)
+ S4=TRIGS(KE+2)
+ IBASE=0
+ DO 540 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 530 IJK=1,ILOT
+!C
+ A10(IJK)=(A(IA+I)-0.25*((A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))))&
+ +QRT5*((A(IB+I)+A(IE+I))-(A(IC+I)+A(ID+I)))
+ A20(IJK)=(A(IA+I)-0.25*((A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I))))&
+ -QRT5*((A(IB+I)+A(IE+I))-(A(IC+I)+A(ID+I)))
+ B10(IJK)=(B(IA+I)-0.25*((B(IB+I)-B(IE+I))+(B(IC+I)-B(ID+I))))&
+ +QRT5*((B(IB+I)-B(IE+I))-(B(IC+I)-B(ID+I)))
+ B20(IJK)=(B(IA+I)-0.25*((B(IB+I)-B(IE+I))+(B(IC+I)-B(ID+I))))&
+ -QRT5*((B(IB+I)-B(IE+I))-(B(IC+I)-B(ID+I)))
+ A11(IJK)=SIN72*(B(IB+I)+B(IE+I))+SIN36*(B(IC+I)+B(ID+I))
+ A21(IJK)=SIN36*(B(IB+I)+B(IE+I))-SIN72*(B(IC+I)+B(ID+I))
+ B11(IJK)=SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))
+ B21(IJK)=SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))
+!C
+ C(JA+J)=A(IA+I)+((A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I)))
+ D(JA+J)=B(IA+I)+((B(IB+I)-B(IE+I))+(B(IC+I)-B(ID+I)))
+ C(JB+J)=C1*(A10(IJK)-A11(IJK))-S1*(B10(IJK)+B11(IJK))
+ D(JB+J)=S1*(A10(IJK)-A11(IJK))+C1*(B10(IJK)+B11(IJK))
+ C(JE+J)=C4*(A10(IJK)+A11(IJK))-S4*(B10(IJK)-B11(IJK))
+ D(JE+J)=S4*(A10(IJK)+A11(IJK))+C4*(B10(IJK)-B11(IJK))
+ C(JC+J)=C2*(A20(IJK)-A21(IJK))-S2*(B20(IJK)+B21(IJK))
+ D(JC+J)=S2*(A20(IJK)-A21(IJK))+C2*(B20(IJK)+B21(IJK))
+ C(JD+J)=C3*(A20(IJK)+A21(IJK))-S3*(B20(IJK)-B21(IJK))
+ D(JD+J)=S3*(A20(IJK)+A21(IJK))+C3*(B20(IJK)-B21(IJK))
+!C
+ I=I+INC3
+ J=J+INC4
+ 530 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 540 CONTINUE
+ IA=IA+IINK
+ IB=IB+IINK
+ IC=IC+IINK
+ ID=ID-IINK
+ IE=IE-IINK
+ JBASE=JBASE+JUMP
+ 550 CONTINUE
+ IF (IB.GT.ID) GO TO 900
+ 560 CONTINUE
+ IBASE=0
+ DO 580 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 570 IJK=1,ILOT
+ C(JA+J)=(A(IA+I)+A(IB+I))+A(IC+I)
+ C(JB+J)=(QRT5*(A(IA+I)-A(IB+I))+(0.25*(A(IA+I)+A(IB+I))-A(IC+I)))&
+ -(SIN36*B(IA+I)+SIN72*B(IB+I))
+ C(JE+J)=-(QRT5*(A(IA+I)-A(IB+I))+(0.25*(A(IA+I)+A(IB+I))-A(IC+I)))&
+ -(SIN36*B(IA+I)+SIN72*B(IB+I))
+ C(JC+J)=(QRT5*(A(IA+I)-A(IB+I))-(0.25*(A(IA+I)+A(IB+I))-A(IC+I)))&
+ -(SIN72*B(IA+I)-SIN36*B(IB+I))
+ C(JD+J)=-(QRT5*(A(IA+I)-A(IB+I))-(0.25*(A(IA+I)+A(IB+I))-A(IC+I)))&
+ -(SIN72*B(IA+I)-SIN36*B(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 570 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 580 CONTINUE
+ GO TO 900
+!C
+ 590 CONTINUE
+ QQRT5=2.0*QRT5
+ SSIN36=2.0*SIN36
+ SSIN72=2.0*SIN72
+ DO 594 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 592 IJK=1,ILOT
+ C(JA+J)=2.0*(A(IA+I)+(A(IB+I)+A(IC+I)))
+ C(JB+J)=(2.0*(A(IA+I)-0.25*(A(IB+I)+A(IC+I)))&
+ +QQRT5*(A(IB+I)-A(IC+I)))-(SSIN72*B(IB+I)+SSIN36*B(IC+I))
+ C(JC+J)=(2.0*(A(IA+I)-0.25*(A(IB+I)+A(IC+I)))&
+ -QQRT5*(A(IB+I)-A(IC+I)))-(SSIN36*B(IB+I)-SSIN72*B(IC+I))
+ C(JD+J)=(2.0*(A(IA+I)-0.25*(A(IB+I)+A(IC+I)))&
+ -QQRT5*(A(IB+I)-A(IC+I)))+(SSIN36*B(IB+I)-SSIN72*B(IC+I))
+ C(JE+J)=(2.0*(A(IA+I)-0.25*(A(IB+I)+A(IC+I)))&
+ +QQRT5*(A(IB+I)-A(IC+I)))+(SSIN72*B(IB+I)+SSIN36*B(IC+I))
+ I=I+INC3
+ J=J+INC4
+ 592 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 594 CONTINUE
+ GO TO 900
+!C
+!C CODING FOR FACTOR 6
+!C -------------------
+ 600 CONTINUE
+ IA=1
+ IB=IA+(2*M-ILA)*INC1
+ IC=IB+2*M*INC1
+ ID=IC+2*M*INC1
+ IE=IC
+ IF=IB
+ JA=1
+ JB=JA+JINK
+ JC=JB+JINK
+ JD=JC+JINK
+ JE=JD+JINK
+ JF=JE+JINK
+!C
+ IF (ILA.EQ.M) GO TO 690
+!C
+ DO 620 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 610 IJK=1,ILOT
+ C(JA+J)=(A(IA+I)+A(ID+I))+(A(IB+I)+A(IC+I))
+ C(JD+J)=(A(IA+I)-A(ID+I))-(A(IB+I)-A(IC+I))
+ C(JB+J)=((A(IA+I)-A(ID+I))+0.5*(A(IB+I)-A(IC+I)))&
+ -(SIN60*(B(IB+I)+B(IC+I)))
+ C(JF+J)=((A(IA+I)-A(ID+I))+0.5*(A(IB+I)-A(IC+I)))&
+ +(SIN60*(B(IB+I)+B(IC+I)))
+ C(JC+J)=((A(IA+I)+A(ID+I))-0.5*(A(IB+I)+A(IC+I)))&
+ -(SIN60*(B(IB+I)-B(IC+I)))
+ C(JE+J)=((A(IA+I)+A(ID+I))-0.5*(A(IB+I)+A(IC+I)))&
+ +(SIN60*(B(IB+I)-B(IC+I)))
+ I=I+INC3
+ J=J+INC4
+ 610 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 620 CONTINUE
+ IA=IA+IINK
+ IINK=2*IINK
+ IB=IB+IINK
+ IC=IC+IINK
+ ID=ID-IINK
+ IE=IE-IINK
+ IF=IF-IINK
+ JBASE=JBASE+JUMP
+ JUMP=2*JUMP+JINK
+ IF (IC.EQ.ID) GO TO 660
+ DO 650 K=ILA,KSTOP,ILA
+ KB=K+K
+ KC=KB+KB
+ KD=KC+KB
+ KE=KD+KB
+ KF=KE+KB
+ C1=TRIGS(KB+1)
+ S1=TRIGS(KB+2)
+ C2=TRIGS(KC+1)
+ S2=TRIGS(KC+2)
+ C3=TRIGS(KD+1)
+ S3=TRIGS(KD+2)
+ C4=TRIGS(KE+1)
+ S4=TRIGS(KE+2)
+ C5=TRIGS(KF+1)
+ S5=TRIGS(KF+2)
+ IBASE=0
+ DO 640 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 630 IJK=1,ILOT
+!C
+ A11(IJK)= (A(IE+I)+A(IB+I))+(A(IC+I)+A(IF+I))
+ A20(IJK)=(A(IA+I)+A(ID+I))-0.5*A11(IJK)
+ A21(IJK)=SIN60*((A(IE+I)+A(IB+I))-(A(IC+I)+A(IF+I)))
+ B11(IJK)= (B(IB+I)-B(IE+I))+(B(IC+I)-B(IF+I))
+ B20(IJK)=(B(IA+I)-B(ID+I))-0.5*B11(IJK)
+ B21(IJK)=SIN60*((B(IB+I)-B(IE+I))-(B(IC+I)-B(IF+I)))
+!C
+ C(JA+J)=(A(IA+I)+A(ID+I))+A11(IJK)
+ D(JA+J)=(B(IA+I)-B(ID+I))+B11(IJK)
+ C(JC+J)=C2*(A20(IJK)-B21(IJK))-S2*(B20(IJK)+A21(IJK))
+ D(JC+J)=S2*(A20(IJK)-B21(IJK))+C2*(B20(IJK)+A21(IJK))
+ C(JE+J)=C4*(A20(IJK)+B21(IJK))-S4*(B20(IJK)-A21(IJK))
+ D(JE+J)=S4*(A20(IJK)+B21(IJK))+C4*(B20(IJK)-A21(IJK))
+!C
+ A11(IJK)=(A(IE+I)-A(IB+I))+(A(IC+I)-A(IF+I))
+ B11(IJK)=(B(IE+I)+B(IB+I))-(B(IC+I)+B(IF+I))
+ A20(IJK)=(A(IA+I)-A(ID+I))-0.5*A11(IJK)
+ A21(IJK)=SIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I)))
+ B20(IJK)=(B(IA+I)+B(ID+I))+0.5*B11(IJK)
+ B21(IJK)=SIN60*((B(IE+I)+B(IB+I))+(B(IC+I)+B(IF+I)))
+!C
+ C(JD+J)=C3*((A(IA+I)-A(ID+I))+A11(IJK))-S3*((B(IA+I)+B(ID+I))-B11(IJK))
+ D(JD+J)=S3*((A(IA+I)-A(ID+I))+A11(IJK))+C3*((B(IA+I)+B(ID+I))-B11(IJK))
+ C(JB+J)=C1*(A20(IJK)-B21(IJK))-S1*(B20(IJK)-A21(IJK))
+ D(JB+J)=S1*(A20(IJK)-B21(IJK))+C1*(B20(IJK)-A21(IJK))
+ C(JF+J)=C5*(A20(IJK)+B21(IJK))-S5*(B20(IJK)+A21(IJK))
+ D(JF+J)=S5*(A20(IJK)+B21(IJK))+C5*(B20(IJK)+A21(IJK))
+!C
+ I=I+INC3
+ J=J+INC4
+ 630 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 640 CONTINUE
+ IA=IA+IINK
+ IB=IB+IINK
+ IC=IC+IINK
+ ID=ID-IINK
+ IE=IE-IINK
+ IF=IF-IINK
+ JBASE=JBASE+JUMP
+ 650 CONTINUE
+ IF (IC.GT.ID) GO TO 900
+ 660 CONTINUE
+ IBASE=0
+ DO 680 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 670 IJK=1,ILOT
+ C(JA+J)=A(IB+I)+(A(IA+I)+A(IC+I))
+ C(JD+J)=B(IB+I)-(B(IA+I)+B(IC+I))
+ C(JB+J)=(SIN60*(A(IA+I)-A(IC+I)))-(0.5*(B(IA+I)+B(IC+I))+B(IB+I))
+ C(JF+J)=-(SIN60*(A(IA+I)-A(IC+I)))-(0.5*(B(IA+I)+B(IC+I))+B(IB+I))
+ C(JC+J)=SIN60*(B(IC+I)-B(IA+I))+(0.5*(A(IA+I)+A(IC+I))-A(IB+I))
+ C(JE+J)=SIN60*(B(IC+I)-B(IA+I))-(0.5*(A(IA+I)+A(IC+I))-A(IB+I))
+ I=I+INC3
+ J=J+INC4
+ 670 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 680 CONTINUE
+ GO TO 900
+!C
+ 690 CONTINUE
+ SSIN60=2.0*SIN60
+ DO 694 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 692 IJK=1,ILOT
+ C(JA+J)=(2.0*(A(IA+I)+A(ID+I)))+(2.0*(A(IB+I)+A(IC+I)))
+ C(JD+J)=(2.0*(A(IA+I)-A(ID+I)))-(2.0*(A(IB+I)-A(IC+I)))
+ C(JB+J)=(2.0*(A(IA+I)-A(ID+I))+(A(IB+I)-A(IC+I)))&
+ -(SSIN60*(B(IB+I)+B(IC+I)))
+ C(JF+J)=(2.0*(A(IA+I)-A(ID+I))+(A(IB+I)-A(IC+I)))&
+ +(SSIN60*(B(IB+I)+B(IC+I)))
+ C(JC+J)=(2.0*(A(IA+I)+A(ID+I))-(A(IB+I)+A(IC+I)))&
+ -(SSIN60*(B(IB+I)-B(IC+I)))
+ C(JE+J)=(2.0*(A(IA+I)+A(ID+I))-(A(IB+I)+A(IC+I)))&
+ +(SSIN60*(B(IB+I)-B(IC+I)))
+ I=I+INC3
+ J=J+INC4
+ 692 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 694 CONTINUE
+ GO TO 900
+!C
+!C CODING FOR FACTOR 8
+!C -------------------
+ 800 CONTINUE
+ IBAD=3
+ IF (ILA.NE.M) GO TO 910
+ IA=1
+ IB=IA+ILA*INC1
+ IC=IB+2*ILA*INC1
+ ID=IC+2*ILA*INC1
+ IE=ID+2*ILA*INC1
+ JA=1
+ JB=JA+JINK
+ JC=JB+JINK
+ JD=JC+JINK
+ JE=JD+JINK
+ JF=JE+JINK
+ JG=JF+JINK
+ JH=JG+JINK
+ SSIN45=SQRT(2.0)
+!C
+ DO 820 IL=1,ILA
+ I=IBASE
+ J=JBASE
+!DIR$ IVDEP
+!VOCL LOOP,NOVREC
+ DO 810 IJK=1,ILOT
+ C(JA+J)=2.0*(((A(IA+I)+A(IE+I))+A(IC+I))+(A(IB+I)+A(ID+I)))
+ C(JE+J)=2.0*(((A(IA+I)+A(IE+I))+A(IC+I))-(A(IB+I)+A(ID+I)))
+ C(JC+J)=2.0*(((A(IA+I)+A(IE+I))-A(IC+I))-(B(IB+I)-B(ID+I)))
+ C(JG+J)=2.0*(((A(IA+I)+A(IE+I))-A(IC+I))+(B(IB+I)-B(ID+I)))
+ C(JB+J)=2.0*((A(IA+I)-A(IE+I))-B(IC+I))+SSIN45&
+ *((A(IB+I)-A(ID+I))-(B(IB+I)+B(ID+I)))
+ C(JF+J)=2.0*((A(IA+I)-A(IE+I))-B(IC+I))&
+ -SSIN45*((A(IB+I)-A(ID+I))-(B(IB+I)+B(ID+I)))
+ C(JD+J)=2.0*((A(IA+I)-A(IE+I))+B(IC+I))&
+ -SSIN45*((A(IB+I)-A(ID+I))+(B(IB+I)+B(ID+I)))
+ C(JH+J)=2.0*((A(IA+I)-A(IE+I))+B(IC+I))&
+ +SSIN45*((A(IB+I)-A(ID+I))+(B(IB+I)+B(ID+I)))
+ I=I+INC3
+ J=J+INC4
+ 810 CONTINUE
+ IBASE=IBASE+INC1
+ JBASE=JBASE+INC2
+ 820 CONTINUE
+!C
+!C RETURN
+!C ------
+ 900 CONTINUE
+ IBAD=0
+ 910 CONTINUE
+ IERR=IBAD
+ RETURN
+ END
diff --git a/Source/Fortran/set99.f90 b/Source/Fortran/set99.f90
new file mode 100644
index 0000000..5d7a03e
--- /dev/null
+++ b/Source/Fortran/set99.f90
@@ -0,0 +1,65 @@
+!c Copyright 1981-2016 ECMWF.
+!c
+!c This software is licensed under the terms of the Apache Licence
+!c Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+!c
+!c In applying this licence, ECMWF does not waive the privileges and immunities
+!c granted to it by virtue of its status as an intergovernmental organisation
+!c nor does it submit to any jurisdiction.
+!c
+
+ SUBROUTINE SET99(TRIGS,IFAX,N)
+ DIMENSION TRIGS(N),IFAX(*),JFAX(10),LFAX(7)
+!C
+!C SUBROUTINE 'SET99' - COMPUTES FACTORS OF N & TRIGONOMETRIC
+!C FUNCTIONS REQUIRED BY FFT99 & FFT991
+!C
+ DATA LFAX/6,8,5,4,3,2,1/
+ IXXX=1
+!C
+ DEL=4.0*ASIN(1.0)/FLOAT(N)
+ NIL=0
+ NHL=(N/2)-1
+ DO 10 K=NIL,NHL
+ ANGLE=FLOAT(K)*DEL
+ TRIGS(2*K+1)=COS(ANGLE)
+ TRIGS(2*K+2)=SIN(ANGLE)
+ 10 CONTINUE
+!C
+!C FIND FACTORS OF N (8,6,5,4,3,2; ONLY ONE 8 ALLOWED)
+!C LOOK FOR SIXES FIRST, STORE FACTORS IN DESCENDING ORDER
+ NU=N
+ IFAC=6
+ K=0
+ L=1
+ 20 CONTINUE
+ IF (MOD(NU,IFAC).NE.0) GO TO 30
+ K=K+1
+ JFAX(K)=IFAC
+ IF (IFAC.NE.8) GO TO 25
+ IF (K.EQ.1) GO TO 25
+ JFAX(1)=8
+ JFAX(K)=6
+ 25 CONTINUE
+ NU=NU/IFAC
+ IF (NU.EQ.1) GO TO 50
+ IF (IFAC.NE.8) GO TO 20
+ 30 CONTINUE
+ L=L+1
+ IFAC=LFAX(L)
+ IF (IFAC.GT.1) GO TO 20
+!C
+ WRITE(6,40) N
+ 40 FORMAT('1N =',I4,' - CONTAINS ILLEGAL FACTORS')
+ RETURN
+!C
+!C NOW REVERSE ORDER OF FACTORS
+ 50 CONTINUE
+ NFAX=K
+ IFAX(1)=NFAX
+ DO 60 I=1,NFAX
+ IFAX(NFAX+2-I)=JFAX(I)
+ 60 CONTINUE
+ IFAX(10)=N
+ RETURN
+ END
--
GitLab