C $Header: /u/gcmpack/MITgcm/pkg/zonal_filt/fftpack.F,v 1.7 2003/12/24 00:27:53 jmc Exp $
C $Name: $
#include "ZONAL_FILT_OPTIONS.h"
C SUBROUTINE RFFTB (N,R,WSAVE)
C DIMENSION R(1) ,WSAVE(1)
C IF (N .EQ. 1) RETURN
C CALL RFFTB1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
C RETURN
C END
SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
#ifdef ALLOW_ZONAL_FILT
DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) ,
1 WA1(*)
DO 101 K=1,L1
CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
101 CONTINUE
IF (IDO-2) 107,105,102
102 IDP2 = IDO+2
DO 104 K=1,L1
DO 103 I=3,IDO,2
IC = IDP2-I
CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
TI2 = CC(I,1,K)+CC(IC,2,K)
CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
103 CONTINUE
104 CONTINUE
IF (MOD(IDO,2) .EQ. 1) RETURN
105 DO 106 K=1,L1
CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
106 CONTINUE
107 RETURN
END
SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) ,
1 WA1(*) ,WA2(*)
DATA TAUR,TAUI /-.5D0,.866025403784439D0/
DO 101 K=1,L1
TR2 = CC(IDO,2,K)+CC(IDO,2,K)
CR2 = CC(1,1,K)+TAUR*TR2
CH(1,K,1) = CC(1,1,K)+TR2
CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
CH(1,K,2) = CR2-CI3
CH(1,K,3) = CR2+CI3
101 CONTINUE
IF (IDO .EQ. 1) RETURN
IDP2 = IDO+2
DO 103 K=1,L1
DO 102 I=3,IDO,2
IC = IDP2-I
TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
CR2 = CC(I-1,1,K)+TAUR*TR2
CH(I-1,K,1) = CC(I-1,1,K)+TR2
TI2 = CC(I,3,K)-CC(IC,2,K)
CI2 = CC(I,1,K)+TAUR*TI2
CH(I,K,1) = CC(I,1,K)+TI2
CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
DR2 = CR2-CI3
DR3 = CR2+CI3
DI2 = CI2+CR3
DI3 = CI2-CR3
CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
102 CONTINUE
103 CONTINUE
RETURN
END
SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) ,
1 WA1(*) ,WA2(*) ,WA3(*)
DATA SQRT2 /1.414213562373095D0/
DO 101 K=1,L1
TR1 = CC(1,1,K)-CC(IDO,4,K)
TR2 = CC(1,1,K)+CC(IDO,4,K)
TR3 = CC(IDO,2,K)+CC(IDO,2,K)
TR4 = CC(1,3,K)+CC(1,3,K)
CH(1,K,1) = TR2+TR3
CH(1,K,2) = TR1-TR4
CH(1,K,3) = TR2-TR3
CH(1,K,4) = TR1+TR4
101 CONTINUE
IF (IDO-2) 107,105,102
102 IDP2 = IDO+2
DO 104 K=1,L1
DO 103 I=3,IDO,2
IC = IDP2-I
TI1 = CC(I,1,K)+CC(IC,4,K)
TI2 = CC(I,1,K)-CC(IC,4,K)
TI3 = CC(I,3,K)-CC(IC,2,K)
TR4 = CC(I,3,K)+CC(IC,2,K)
TR1 = CC(I-1,1,K)-CC(IC-1,4,K)
TR2 = CC(I-1,1,K)+CC(IC-1,4,K)
TI4 = CC(I-1,3,K)-CC(IC-1,2,K)
TR3 = CC(I-1,3,K)+CC(IC-1,2,K)
CH(I-1,K,1) = TR2+TR3
CR3 = TR2-TR3
CH(I,K,1) = TI2+TI3
CI3 = TI2-TI3
CR2 = TR1-TR4
CR4 = TR1+TR4
CI2 = TI1+TI4
CI4 = TI1-TI4
CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2
CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2
CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3
CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3
CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4
CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4
103 CONTINUE
104 CONTINUE
IF (MOD(IDO,2) .EQ. 1) RETURN
105 CONTINUE
DO 106 K=1,L1
TI1 = CC(1,2,K)+CC(1,4,K)
TI2 = CC(1,4,K)-CC(1,2,K)
TR1 = CC(IDO,1,K)-CC(IDO,3,K)
TR2 = CC(IDO,1,K)+CC(IDO,3,K)
CH(IDO,K,1) = TR2+TR2
CH(IDO,K,2) = SQRT2*(TR1-TI1)
CH(IDO,K,3) = TI2+TI2
CH(IDO,K,4) = -SQRT2*(TR1+TI1)
106 CONTINUE
107 RETURN
END
SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) ,
1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*)
DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
1-.809016994374947D0,.587785252292473D0/
DO 101 K=1,L1
TI5 = CC(1,3,K)+CC(1,3,K)
TI4 = CC(1,5,K)+CC(1,5,K)
TR2 = CC(IDO,2,K)+CC(IDO,2,K)
TR3 = CC(IDO,4,K)+CC(IDO,4,K)
CH(1,K,1) = CC(1,1,K)+TR2+TR3
CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
CI5 = TI11*TI5+TI12*TI4
CI4 = TI12*TI5-TI11*TI4
CH(1,K,2) = CR2-CI5
CH(1,K,3) = CR3-CI4
CH(1,K,4) = CR3+CI4
CH(1,K,5) = CR2+CI5
101 CONTINUE
IF (IDO .EQ. 1) RETURN
IDP2 = IDO+2
DO 103 K=1,L1
DO 102 I=3,IDO,2
IC = IDP2-I
TI5 = CC(I,3,K)+CC(IC,2,K)
TI2 = CC(I,3,K)-CC(IC,2,K)
TI4 = CC(I,5,K)+CC(IC,4,K)
TI3 = CC(I,5,K)-CC(IC,4,K)
TR5 = CC(I-1,3,K)-CC(IC-1,2,K)
TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
TR4 = CC(I-1,5,K)-CC(IC-1,4,K)
TR3 = CC(I-1,5,K)+CC(IC-1,4,K)
CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
CH(I,K,1) = CC(I,1,K)+TI2+TI3
CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
CR5 = TI11*TR5+TI12*TR4
CI5 = TI11*TI5+TI12*TI4
CR4 = TI12*TR5-TI11*TR4
CI4 = TI12*TI5-TI11*TI4
DR3 = CR3-CI4
DR4 = CR3+CI4
DI3 = CI3+CR4
DI4 = CI3-CR4
DR5 = CR2+CI5
DR2 = CR2-CI5
DI5 = CI2-CR5
DI2 = CI2+CR5
CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4
CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4
CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5
CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5
102 CONTINUE
103 CONTINUE
RETURN
END
SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
1 C1(IDO,L1,IP) ,C2(IDL1,IP),
2 CH2(IDL1,IP) ,WA(*)
DATA TPI/6.28318530717959D0/
ARG = TPI/FLOAT(IP)
DCP = COS(ARG)
DSP = SIN(ARG)
IDP2 = IDO+2
NBD = (IDO-1)/2
IPP2 = IP+2
IPPH = (IP+1)/2
IF (IDO .LT. L1) GO TO 103
DO 102 K=1,L1
DO 101 I=1,IDO
CH(I,K,1) = CC(I,1,K)
101 CONTINUE
102 CONTINUE
GO TO 106
103 DO 105 I=1,IDO
DO 104 K=1,L1
CH(I,K,1) = CC(I,1,K)
104 CONTINUE
105 CONTINUE
106 DO 108 J=2,IPPH
JC = IPP2-J
J2 = J+J
DO 107 K=1,L1
CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
107 CONTINUE
108 CONTINUE
IF (IDO .EQ. 1) GO TO 116
IF (NBD .LT. L1) GO TO 112
DO 111 J=2,IPPH
JC = IPP2-J
DO 110 K=1,L1
DO 109 I=3,IDO,2
IC = IDP2-I
CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
109 CONTINUE
110 CONTINUE
111 CONTINUE
GO TO 116
112 DO 115 J=2,IPPH
JC = IPP2-J
DO 114 I=3,IDO,2
IC = IDP2-I
DO 113 K=1,L1
CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
113 CONTINUE
114 CONTINUE
115 CONTINUE
116 AR1 = 1.D0
AI1 = 0.D0
DO 120 L=2,IPPH
LC = IPP2-L
AR1H = DCP*AR1-DSP*AI1
AI1 = DCP*AI1+DSP*AR1
AR1 = AR1H
DO 117 IK=1,IDL1
C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
C2(IK,LC) = AI1*CH2(IK,IP)
117 CONTINUE
DC2 = AR1
DS2 = AI1
AR2 = AR1
AI2 = AI1
DO 119 J=3,IPPH
JC = IPP2-J
AR2H = DC2*AR2-DS2*AI2
AI2 = DC2*AI2+DS2*AR2
AR2 = AR2H
DO 118 IK=1,IDL1
C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
118 CONTINUE
119 CONTINUE
120 CONTINUE
DO 122 J=2,IPPH
DO 121 IK=1,IDL1
CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
121 CONTINUE
122 CONTINUE
DO 124 J=2,IPPH
JC = IPP2-J
DO 123 K=1,L1
CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
123 CONTINUE
124 CONTINUE
IF (IDO .EQ. 1) GO TO 132
IF (NBD .LT. L1) GO TO 128
DO 127 J=2,IPPH
JC = IPP2-J
DO 126 K=1,L1
DO 125 I=3,IDO,2
CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
125 CONTINUE
126 CONTINUE
127 CONTINUE
GO TO 132
128 DO 131 J=2,IPPH
JC = IPP2-J
DO 130 I=3,IDO,2
DO 129 K=1,L1
CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
129 CONTINUE
130 CONTINUE
131 CONTINUE
132 CONTINUE
IF (IDO .EQ. 1) RETURN
DO 133 IK=1,IDL1
C2(IK,1) = CH2(IK,1)
133 CONTINUE
DO 135 J=2,IP
DO 134 K=1,L1
C1(1,K,J) = CH(1,K,J)
134 CONTINUE
135 CONTINUE
IF (NBD .GT. L1) GO TO 139
IS = -IDO
DO 138 J=2,IP
IS = IS+IDO
IDIJ = IS
DO 137 I=3,IDO,2
IDIJ = IDIJ+2
DO 136 K=1,L1
C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
136 CONTINUE
137 CONTINUE
138 CONTINUE
GO TO 143
139 IS = -IDO
DO 142 J=2,IP
IS = IS+IDO
DO 141 K=1,L1
IDIJ = IS
DO 140 I=3,IDO,2
IDIJ = IDIJ+2
C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
140 CONTINUE
141 CONTINUE
142 CONTINUE
143 RETURN
END
SUBROUTINE RFFTB1 (N,C,CH,WA,IFAC)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*)
NF = IFAC(2)
NA = 0
L1 = 1
IW = 1
DO 116 K1=1,NF
IP = IFAC(K1+2)
L2 = IP*L1
IDO = N/L2
IDL1 = IDO*L1
IF (IP .NE. 4) GO TO 103
IX2 = IW+IDO
IX3 = IX2+IDO
IF (NA .NE. 0) GO TO 101
CALL RADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
GO TO 102
101 CALL RADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
102 NA = 1-NA
GO TO 115
103 IF (IP .NE. 2) GO TO 106
IF (NA .NE. 0) GO TO 104
CALL RADB2 (IDO,L1,C,CH,WA(IW))
GO TO 105
104 CALL RADB2 (IDO,L1,CH,C,WA(IW))
105 NA = 1-NA
GO TO 115
106 IF (IP .NE. 3) GO TO 109
IX2 = IW+IDO
IF (NA .NE. 0) GO TO 107
CALL RADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
GO TO 108
107 CALL RADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
108 NA = 1-NA
GO TO 115
109 IF (IP .NE. 5) GO TO 112
IX2 = IW+IDO
IX3 = IX2+IDO
IX4 = IX3+IDO
IF (NA .NE. 0) GO TO 110
CALL RADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
GO TO 111
110 CALL RADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
111 NA = 1-NA
GO TO 115
112 IF (NA .NE. 0) GO TO 113
CALL RADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
GO TO 114
113 CALL RADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
114 IF (IDO .EQ. 1) NA = 1-NA
115 L1 = L2
IW = IW+(IP-1)*IDO
116 CONTINUE
IF (NA .EQ. 0) RETURN
DO 117 I=1,N
C(I) = CH(I)
117 CONTINUE
RETURN
END
C SUBROUTINE RFFTF (N,R,WSAVE)
C IMPLICIT REAL*8 (A-H,O-Z)
C DIMENSION R(1) ,WSAVE(1)
C IF (N .EQ. 1) RETURN
C CALL RFFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
C RETURN
C END
SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) ,
1 WA1(*)
DO 101 K=1,L1
CH(1,1,K) = CC(1,K,1)+CC(1,K,2)
CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2)
101 CONTINUE
IF (IDO-2) 107,105,102
102 IDP2 = IDO+2
DO 104 K=1,L1
DO 103 I=3,IDO,2
IC = IDP2-I
TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
CH(I,1,K) = CC(I,K,1)+TI2
CH(IC,2,K) = TI2-CC(I,K,1)
CH(I-1,1,K) = CC(I-1,K,1)+TR2
CH(IC-1,2,K) = CC(I-1,K,1)-TR2
103 CONTINUE
104 CONTINUE
IF (MOD(IDO,2) .EQ. 1) RETURN
105 DO 106 K=1,L1
CH(1,2,K) = -CC(IDO,K,2)
CH(IDO,1,K) = CC(IDO,K,1)
106 CONTINUE
107 RETURN
END
SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) ,
1 WA1(*) ,WA2(*)
DATA TAUR,TAUI /-.5D0,.866025403784439D0/
DO 101 K=1,L1
CR2 = CC(1,K,2)+CC(1,K,3)
CH(1,1,K) = CC(1,K,1)+CR2
CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))
CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2
101 CONTINUE
IF (IDO .EQ. 1) RETURN
IDP2 = IDO+2
DO 103 K=1,L1
DO 102 I=3,IDO,2
IC = IDP2-I
DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
CR2 = DR2+DR3
CI2 = DI2+DI3
CH(I-1,1,K) = CC(I-1,K,1)+CR2
CH(I,1,K) = CC(I,K,1)+CI2
TR2 = CC(I-1,K,1)+TAUR*CR2
TI2 = CC(I,K,1)+TAUR*CI2
TR3 = TAUI*(DI2-DI3)
TI3 = TAUI*(DR3-DR2)
CH(I-1,3,K) = TR2+TR3
CH(IC-1,2,K) = TR2-TR3
CH(I,3,K) = TI2+TI3
CH(IC,2,K) = TI3-TI2
102 CONTINUE
103 CONTINUE
RETURN
END
SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) ,
1 WA1(*) ,WA2(*) ,WA3(*)
DATA HSQT2 /.7071067811865475D0/
DO 101 K=1,L1
TR1 = CC(1,K,2)+CC(1,K,4)
TR2 = CC(1,K,1)+CC(1,K,3)
CH(1,1,K) = TR1+TR2
CH(IDO,4,K) = TR2-TR1
CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3)
CH(1,3,K) = CC(1,K,4)-CC(1,K,2)
101 CONTINUE
IF (IDO-2) 107,105,102
102 IDP2 = IDO+2
DO 104 K=1,L1
DO 103 I=3,IDO,2
IC = IDP2-I
CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
TR1 = CR2+CR4
TR4 = CR4-CR2
TI1 = CI2+CI4
TI4 = CI2-CI4
TI2 = CC(I,K,1)+CI3
TI3 = CC(I,K,1)-CI3
TR2 = CC(I-1,K,1)+CR3
TR3 = CC(I-1,K,1)-CR3
CH(I-1,1,K) = TR1+TR2
CH(IC-1,4,K) = TR2-TR1
CH(I,1,K) = TI1+TI2
CH(IC,4,K) = TI1-TI2
CH(I-1,3,K) = TI4+TR3
CH(IC-1,2,K) = TR3-TI4
CH(I,3,K) = TR4+TI3
CH(IC,2,K) = TR4-TI3
103 CONTINUE
104 CONTINUE
IF (MOD(IDO,2) .EQ. 1) RETURN
105 CONTINUE
DO 106 K=1,L1
TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4))
TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4))
CH(IDO,1,K) = TR1+CC(IDO,K,1)
CH(IDO,3,K) = CC(IDO,K,1)-TR1
CH(1,2,K) = TI1-CC(IDO,K,3)
CH(1,4,K) = TI1+CC(IDO,K,3)
106 CONTINUE
107 RETURN
END
SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT INTEGER (I-N)
DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) ,
1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*)
DATA TR11,TI11,TR12,TI12 /.309016994374947D0,.951056516295154D0,
1-.809016994374947D0,.587785252292473D0/
DO 101 K=1,L1
CR2 = CC(1,K,5)+CC(1,K,2)
CI5 = CC(1,K,5)-CC(1,K,2)
CR3 = CC(1,K,4)+CC(1,K,3)
CI4 = CC(1,K,4)-CC(1,<