Skip to content
Snippets Groups Projects
Commit 4f247989 authored by Anne Tipka's avatar Anne Tipka :headphones:
Browse files

elimination of emoslib in fortran code. See ticket #312

parent 81c180c9
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Showing
with 2017 additions and 7 deletions
Source/Fortran/changelog.txt
+ 10
− 1
View file @ 4f247989
......@@ -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
Source/Fortran/fft99.f90 0 → 100644
+ 218
− 0
View file @ 4f247989
!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
Source/Fortran/jsppole.f90 0 → 100644
+ 182
− 0
View file @ 4f247989
! 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
Source/Fortran/makefile_cray
+ 1
− 1
View file @ 4f247989
......@@ -14,7 +14,7 @@
EXE = calc_etadot
LIB = $(ECCODES_LIB) $(EMOSLIB)
LIB = $(ECCODES_LIB)
INC = -I. -I$(ECCODES_INCLUDE_DIR)
FC = ftn
......
Source/Fortran/makefile_debug
+ 1
− 2
View file @ 4f247989
......@@ -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)
......
Source/Fortran/makefile_ecgate
+ 1
− 1
View file @ 4f247989
......@@ -14,7 +14,7 @@
EXE = calc_etadot_fast.out
LIB = $(ECCODES_LIB) $(EMOSLIB)
LIB = $(ECCODES_LIB)
INC = -I. -I$(ECCODES_INCLUDE_DIR)
FC=gfortran
......
Source/Fortran/makefile_fast
+ 1
− 2
View file @ 4f247989
......@@ -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
......
Source/Fortran/qpassm.f90 0 → 100644
+ 778
− 0
View file @ 4f247989
!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
Source/Fortran/rpassm.f90 0 → 100644
+ 760
− 0
View file @ 4f247989
This diff is collapsed.
Source/Fortran/set99.f90 0 → 100644
+ 65
− 0
View file @ 4f247989
!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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment