| blob | cebe84c708ebfe34eeb34e9475c12ed2363a36dd |
1 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2 !
3 ! FFTPACK 5.0
4 ! Copyright (C) 1995-2004, Scientific Computing Division,
5 ! University Corporation for Atmospheric Research
6 ! Licensed under the GNU General Public License (GPL)
7 !
8 ! Authors: Paul N. Swarztrauber and Richard A. Valent
9 !
10 ! $Id: mradfg.f,v 1.2 2004/06/15 21:29:20 rodney Exp $
11 !
12 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
14 SUBROUTINE MRADFG (M,IDO,IP,L1,IDL1,CC,C1,C2,IM1,IN1, &
15 & CH,CH2,IM2,IN2,WA)
16 REAL CH(IN2,IDO,L1,IP) ,CC(IN1,IDO,IP,L1), &
17 & C1(IN1,IDO,L1,IP) ,C2(IN1,IDL1,IP), &
18 & CH2(IN2,IDL1,IP) ,WA(IDO)
19 !
20 M1D = (M-1)*IM1+1
21 M2S = 1-IM2
22 TPI=2.*4.*ATAN(1.0)
23 ARG = TPI/FLOAT(IP)
24 DCP = COS(ARG)
25 DSP = SIN(ARG)
26 IPPH = (IP+1)/2
27 IPP2 = IP+2
28 IDP2 = IDO+2
29 NBD = (IDO-1)/2
30 IF (IDO .EQ. 1) GO TO 119
31 DO 101 IK=1,IDL1
32 M2 = M2S
33 DO 1001 M1=1,M1D,IM1
34 M2 = M2+IM2
35 CH2(M2,IK,1) = C2(M1,IK,1)
36 1001 CONTINUE
37 101 END DO
38 DO 103 J=2,IP
39 DO 102 K=1,L1
40 M2 = M2S
41 DO 1002 M1=1,M1D,IM1
42 M2 = M2+IM2
43 CH(M2,1,K,J) = C1(M1,1,K,J)
44 1002 CONTINUE
45 102 CONTINUE
46 103 END DO
47 IF (NBD .GT. L1) GO TO 107
48 IS = -IDO
49 DO 106 J=2,IP
50 IS = IS+IDO
51 IDIJ = IS
52 DO 105 I=3,IDO,2
53 IDIJ = IDIJ+2
54 DO 104 K=1,L1
55 M2 = M2S
56 DO 1004 M1=1,M1D,IM1
57 M2 = M2+IM2
58 CH(M2,I-1,K,J) = WA(IDIJ-1)*C1(M1,I-1,K,J)+WA(IDIJ) &
59 & *C1(M1,I,K,J)
60 CH(M2,I,K,J) = WA(IDIJ-1)*C1(M1,I,K,J)-WA(IDIJ) &
61 & *C1(M1,I-1,K,J)
62 1004 CONTINUE
63 104 CONTINUE
64 105 CONTINUE
65 106 END DO
66 GO TO 111
67 107 IS = -IDO
68 DO 110 J=2,IP
69 IS = IS+IDO
70 DO 109 K=1,L1
71 IDIJ = IS
72 DO 108 I=3,IDO,2
73 IDIJ = IDIJ+2
74 M2 = M2S
75 DO 1008 M1=1,M1D,IM1
76 M2 = M2+IM2
77 CH(M2,I-1,K,J) = WA(IDIJ-1)*C1(M1,I-1,K,J)+WA(IDIJ) &
78 & *C1(M1,I,K,J)
79 CH(M2,I,K,J) = WA(IDIJ-1)*C1(M1,I,K,J)-WA(IDIJ) &
80 & *C1(M1,I-1,K,J)
81 1008 CONTINUE
82 108 CONTINUE
83 109 CONTINUE
84 110 END DO
85 111 IF (NBD .LT. L1) GO TO 115
86 DO 114 J=2,IPPH
87 JC = IPP2-J
88 DO 113 K=1,L1
89 DO 112 I=3,IDO,2
90 M2 = M2S
91 DO 1012 M1=1,M1D,IM1
92 M2 = M2+IM2
93 C1(M1,I-1,K,J) = CH(M2,I-1,K,J)+CH(M2,I-1,K,JC)
94 C1(M1,I-1,K,JC) = CH(M2,I,K,J)-CH(M2,I,K,JC)
95 C1(M1,I,K,J) = CH(M2,I,K,J)+CH(M2,I,K,JC)
96 C1(M1,I,K,JC) = CH(M2,I-1,K,JC)-CH(M2,I-1,K,J)
97 1012 CONTINUE
98 112 CONTINUE
99 113 CONTINUE
100 114 END DO
101 GO TO 121
102 115 DO 118 J=2,IPPH
103 JC = IPP2-J
104 DO 117 I=3,IDO,2
105 DO 116 K=1,L1
106 M2 = M2S
107 DO 1016 M1=1,M1D,IM1
108 M2 = M2+IM2
109 C1(M1,I-1,K,J) = CH(M2,I-1,K,J)+CH(M2,I-1,K,JC)
110 C1(M1,I-1,K,JC) = CH(M2,I,K,J)-CH(M2,I,K,JC)
111 C1(M1,I,K,J) = CH(M2,I,K,J)+CH(M2,I,K,JC)
112 C1(M1,I,K,JC) = CH(M2,I-1,K,JC)-CH(M2,I-1,K,J)
113 1016 CONTINUE
114 116 CONTINUE
115 117 CONTINUE
116 118 END DO
117 GO TO 121
118 119 DO 120 IK=1,IDL1
119 M2 = M2S
120 DO 1020 M1=1,M1D,IM1
121 M2 = M2+IM2
122 C2(M1,IK,1) = CH2(M2,IK,1)
123 1020 CONTINUE
124 120 END DO
125 121 DO 123 J=2,IPPH
126 JC = IPP2-J
127 DO 122 K=1,L1
128 M2 = M2S
129 DO 1022 M1=1,M1D,IM1
130 M2 = M2+IM2
131 C1(M1,1,K,J) = CH(M2,1,K,J)+CH(M2,1,K,JC)
132 C1(M1,1,K,JC) = CH(M2,1,K,JC)-CH(M2,1,K,J)
133 1022 CONTINUE
134 122 CONTINUE
135 123 END DO
136 !
137 AR1 = 1.
138 AI1 = 0.
139 DO 127 L=2,IPPH
140 LC = IPP2-L
141 AR1H = DCP*AR1-DSP*AI1
142 AI1 = DCP*AI1+DSP*AR1
143 AR1 = AR1H
144 DO 124 IK=1,IDL1
145 M2 = M2S
146 DO 1024 M1=1,M1D,IM1
147 M2 = M2+IM2
148 CH2(M2,IK,L) = C2(M1,IK,1)+AR1*C2(M1,IK,2)
149 CH2(M2,IK,LC) = AI1*C2(M1,IK,IP)
150 1024 CONTINUE
151 124 CONTINUE
152 DC2 = AR1
153 DS2 = AI1
154 AR2 = AR1
155 AI2 = AI1
156 DO 126 J=3,IPPH
157 JC = IPP2-J
158 AR2H = DC2*AR2-DS2*AI2
159 AI2 = DC2*AI2+DS2*AR2
160 AR2 = AR2H
161 DO 125 IK=1,IDL1
162 M2 = M2S
163 DO 1025 M1=1,M1D,IM1
164 M2 = M2+IM2
165 CH2(M2,IK,L) = CH2(M2,IK,L)+AR2*C2(M1,IK,J)
166 CH2(M2,IK,LC) = CH2(M2,IK,LC)+AI2*C2(M1,IK,JC)
167 1025 CONTINUE
168 125 CONTINUE
169 126 CONTINUE
170 127 END DO
171 DO 129 J=2,IPPH
172 DO 128 IK=1,IDL1
173 M2 = M2S
174 DO 1028 M1=1,M1D,IM1
175 M2 = M2+IM2
176 CH2(M2,IK,1) = CH2(M2,IK,1)+C2(M1,IK,J)
177 1028 CONTINUE
178 128 CONTINUE
179 129 END DO
180 !
181 IF (IDO .LT. L1) GO TO 132
182 DO 131 K=1,L1
183 DO 130 I=1,IDO
184 M2 = M2S
185 DO 1030 M1=1,M1D,IM1
186 M2 = M2+IM2
187 CC(M1,I,1,K) = CH(M2,I,K,1)
188 1030 CONTINUE
189 130 CONTINUE
190 131 END DO
191 GO TO 135
192 132 DO 134 I=1,IDO
193 DO 133 K=1,L1
194 M2 = M2S
195 DO 1033 M1=1,M1D,IM1
196 M2 = M2+IM2
197 CC(M1,I,1,K) = CH(M2,I,K,1)
198 1033 CONTINUE
199 133 CONTINUE
200 134 END DO
201 135 DO 137 J=2,IPPH
202 JC = IPP2-J
203 J2 = J+J
204 DO 136 K=1,L1
205 M2 = M2S
206 DO 1036 M1=1,M1D,IM1
207 M2 = M2+IM2
208 CC(M1,IDO,J2-2,K) = CH(M2,1,K,J)
209 CC(M1,1,J2-1,K) = CH(M2,1,K,JC)
210 1036 CONTINUE
211 136 CONTINUE
212 137 END DO
213 IF (IDO .EQ. 1) RETURN
214 IF (NBD .LT. L1) GO TO 141
215 DO 140 J=2,IPPH
216 JC = IPP2-J
217 J2 = J+J
218 DO 139 K=1,L1
219 DO 138 I=3,IDO,2
220 IC = IDP2-I
221 M2 = M2S
222 DO 1038 M1=1,M1D,IM1
223 M2 = M2+IM2
224 CC(M1,I-1,J2-1,K) = CH(M2,I-1,K,J)+CH(M2,I-1,K,JC)
225 CC(M1,IC-1,J2-2,K) = CH(M2,I-1,K,J)-CH(M2,I-1,K,JC)
226 CC(M1,I,J2-1,K) = CH(M2,I,K,J)+CH(M2,I,K,JC)
227 CC(M1,IC,J2-2,K) = CH(M2,I,K,JC)-CH(M2,I,K,J)
228 1038 CONTINUE
229 138 CONTINUE
230 139 CONTINUE
231 140 END DO
232 RETURN
233 141 DO 144 J=2,IPPH
234 JC = IPP2-J
235 J2 = J+J
236 DO 143 I=3,IDO,2
237 IC = IDP2-I
238 DO 142 K=1,L1
239 M2 = M2S
240 DO 1042 M1=1,M1D,IM1
241 M2 = M2+IM2
242 CC(M1,I-1,J2-1,K) = CH(M2,I-1,K,J)+CH(M2,I-1,K,JC)
243 CC(M1,IC-1,J2-2,K) = CH(M2,I-1,K,J)-CH(M2,I-1,K,JC)
244 CC(M1,I,J2-1,K) = CH(M2,I,K,J)+CH(M2,I,K,JC)
245 CC(M1,IC,J2-2,K) = CH(M2,I,K,JC)-CH(M2,I,K,J)
246 1042 CONTINUE
247 142 CONTINUE
248 143 CONTINUE
249 144 END DO
250 RETURN
251 END