Merge remote-tracking branch 'origin/release-v4.6.1'
[WRF.git] / external / fftpack / fftpack5 / cmf3kf.F
blob9f991293c706943ed763bbc26cc42b93f19819f4
1 subroutine cmf3kf ( lot, ido, l1, na, cc, im1, in1, ch, im2, in2, wa )
3 !*****************************************************************************80
5 !! CMF3KF is an FFTPACK5 auxiliary routine.
8 !    Copyright (C) 1995-2004, Scientific Computing Division,
9 !    University Corporation for Atmospheric Research
11 !  Modified:
13 !    27 March 2009
15 !  Author:
17 !    Paul Swarztrauber
18 !    Richard Valent
20 !  Reference:
22 !    Paul Swarztrauber,
23 !    Vectorizing the Fast Fourier Transforms,
24 !    in Parallel Computations,
25 !    edited by G. Rodrigue,
26 !    Academic Press, 1982.
28 !    Paul Swarztrauber,
29 !    Fast Fourier Transform Algorithms for Vector Computers,
30 !    Parallel Computing, pages 45-63, 1984.
32 !  Parameters:
34   implicit none
36   integer ( kind = 4 ) ido
37   integer ( kind = 4 ) in1
38   integer ( kind = 4 ) in2
39   integer ( kind = 4 ) l1
41   real ( kind = 4 ) cc(2,in1,l1,ido,3)
42   real ( kind = 4 ) ch(2,in2,l1,3,ido)
43   real ( kind = 4 ) ci2
44   real ( kind = 4 ) ci3
45   real ( kind = 4 ) cr2
46   real ( kind = 4 ) cr3
47   real ( kind = 4 ) di2
48   real ( kind = 4 ) di3
49   real ( kind = 4 ) dr2
50   real ( kind = 4 ) dr3
51   integer ( kind = 4 ) i
52   integer ( kind = 4 ) im1
53   integer ( kind = 4 ) im2
54   integer ( kind = 4 ) k
55   integer ( kind = 4 ) lot
56   integer ( kind = 4 ) m1
57   integer ( kind = 4 ) m1d
58   integer ( kind = 4 ) m2
59   integer ( kind = 4 ) m2s
60   integer ( kind = 4 ) na
61   real ( kind = 4 ) sn
62   real ( kind = 4 ), parameter :: taui = -0.866025403784439E+00
63   real ( kind = 4 ), parameter :: taur = -0.5E+00
64   real ( kind = 4 ) ti2
65   real ( kind = 4 ) tr2
66   real ( kind = 4 ) wa(ido,2,2)
68   m1d = ( lot - 1 ) * im1 + 1
69   m2s = 1 - im2
71   if ( 1 < ido ) then
73     do k = 1, l1
74       m2 = m2s
75       do m1 = 1, m1d, im1
76         m2 = m2 + im2
77         tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,3)
78         cr2 = cc(1,m1,k,1,1)+taur*tr2
79         ch(1,m2,k,1,1) = cc(1,m1,k,1,1)+tr2
80         ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,3)
81         ci2 = cc(2,m1,k,1,1)+taur*ti2
82         ch(2,m2,k,1,1) = cc(2,m1,k,1,1)+ti2
83         cr3 = taui*(cc(1,m1,k,1,2)-cc(1,m1,k,1,3))
84         ci3 = taui*(cc(2,m1,k,1,2)-cc(2,m1,k,1,3))
85         ch(1,m2,k,2,1) = cr2-ci3
86         ch(1,m2,k,3,1) = cr2+ci3
87         ch(2,m2,k,2,1) = ci2+cr3
88         ch(2,m2,k,3,1) = ci2-cr3
89       end do
90     end do
92     do i = 2, ido
93       do k = 1, l1
94         m2 = m2s
95         do m1 = 1, m1d, im1
96           m2 = m2 + im2
97           tr2 = cc(1,m1,k,i,2)+cc(1,m1,k,i,3)
98           cr2 = cc(1,m1,k,i,1)+taur*tr2
99           ch(1,m2,k,1,i) = cc(1,m1,k,i,1)+tr2
100           ti2 = cc(2,m1,k,i,2)+cc(2,m1,k,i,3)
101           ci2 = cc(2,m1,k,i,1)+taur*ti2
102           ch(2,m2,k,1,i) = cc(2,m1,k,i,1)+ti2
103           cr3 = taui*(cc(1,m1,k,i,2)-cc(1,m1,k,i,3))
104           ci3 = taui*(cc(2,m1,k,i,2)-cc(2,m1,k,i,3))
105           dr2 = cr2-ci3
106           dr3 = cr2+ci3
107           di2 = ci2+cr3
108           di3 = ci2-cr3
109           ch(2,m2,k,2,i) = wa(i,1,1)*di2-wa(i,1,2)*dr2
110           ch(1,m2,k,2,i) = wa(i,1,1)*dr2+wa(i,1,2)*di2
111           ch(2,m2,k,3,i) = wa(i,2,1)*di3-wa(i,2,2)*dr3
112           ch(1,m2,k,3,i) = wa(i,2,1)*dr3+wa(i,2,2)*di3
113         end do
114       end do
115     end do
117   else if ( na == 1 ) then
119     sn = 1.0E+00 / real ( 3 * l1, kind = 4 )
121     do k = 1, l1
122       m2 = m2s
123       do m1 = 1, m1d, im1
124         m2 = m2 + im2
125         tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,3)
126         cr2 = cc(1,m1,k,1,1)+taur*tr2
127         ch(1,m2,k,1,1) = sn*(cc(1,m1,k,1,1)+tr2)
128         ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,3)
129         ci2 = cc(2,m1,k,1,1)+taur*ti2
130         ch(2,m2,k,1,1) = sn*(cc(2,m1,k,1,1)+ti2)
131         cr3 = taui*(cc(1,m1,k,1,2)-cc(1,m1,k,1,3))
132         ci3 = taui*(cc(2,m1,k,1,2)-cc(2,m1,k,1,3))
133         ch(1,m2,k,2,1) = sn*(cr2-ci3)
134         ch(1,m2,k,3,1) = sn*(cr2+ci3)
135         ch(2,m2,k,2,1) = sn*(ci2+cr3)
136         ch(2,m2,k,3,1) = sn*(ci2-cr3)
137       end do
138     end do
140   else
142     sn = 1.0E+00 / real ( 3 * l1, kind = 4 )
144     do k = 1, l1
145       do m1 = 1, m1d, im1
146         tr2 = cc(1,m1,k,1,2)+cc(1,m1,k,1,3)
147         cr2 = cc(1,m1,k,1,1)+taur*tr2
148         cc(1,m1,k,1,1) = sn*(cc(1,m1,k,1,1)+tr2)
149         ti2 = cc(2,m1,k,1,2)+cc(2,m1,k,1,3)
150         ci2 = cc(2,m1,k,1,1)+taur*ti2
151         cc(2,m1,k,1,1) = sn*(cc(2,m1,k,1,1)+ti2)
152         cr3 = taui*(cc(1,m1,k,1,2)-cc(1,m1,k,1,3))
153         ci3 = taui*(cc(2,m1,k,1,2)-cc(2,m1,k,1,3))
154         cc(1,m1,k,1,2) = sn*(cr2-ci3)
155         cc(1,m1,k,1,3) = sn*(cr2+ci3)
156         cc(2,m1,k,1,2) = sn*(ci2+cr3)
157         cc(2,m1,k,1,3) = sn*(ci2-cr3)
158       end do
159     end do
161   end if
163   return