37
38
39
40
41
42#include "implicit_f.inc"
43#include "mvsiz_p.inc"
44
45
46
47 INTEGER JFT,JLT,NPLAT ,IPLAT(*),ISROT
49 . vqn(mvsiz,9,4),vf(mvsiz,3,4),vm(mvsiz,2,4),vq(mvsiz,3,3),
50 . corel(mvsiz,3,4),di(mvsiz,6),vmz(mvsiz,4)
52 . f11(mvsiz), f12(mvsiz), f13(mvsiz), f14(mvsiz),
53 . f21(mvsiz), f22(mvsiz), f23(mvsiz), f24(mvsiz),
54 . f31(mvsiz), f32(mvsiz), f33(mvsiz), f34(mvsiz),
55 . m11(mvsiz), m12(mvsiz), m13(mvsiz), m14(mvsiz),
56 . m21(mvsiz), m22(mvsiz), m23(mvsiz), m24(mvsiz),
57 . m31(mvsiz), m32(mvsiz), m33(mvsiz
58 . off(*)
59
60
61
62 INTEGER I, J, K,EP
64 . mm(3,4),fl(3,4),ml(2,4),c1,z1,
65 . ar(3),ad(4),alr(3),ald(4),dbad(3),mlz(mvsiz,3,4)
66
67#include "vectorize.inc"
68 DO ep=jft,nplat
69 k=iplat(ep)
70
71 fl(1,1)= vf(k,1,1)+vf(k,1,3)
72 fl(1,2)= vf(k,1,2)+vf(k,1,4)
73 fl(1,3)=-vf(k,1,1)+vf(k,1,3)
74 fl(1,4)=-vf(k,1,2)+vf(k,1,4)
75
76 fl(2,1)= vf(k,2,1)+vf(k,2,3)
77 fl(2,2)= vf(k,2,2)+vf(k,2,4)
78 fl(2,3)=-vf(k,2,1)+vf(k,2,3)
79 fl(2,4)=-vf(k,2,2)+vf(k,2,4)
80
81 fl(3,1)= vf(k,3,1)+vf(k,3,3)
82 fl(3,2)= vf(k,3,2)+vf(k,3,4)
83 fl(3,3)=-vf(k,3,1)+vf(k,3,3)
84 fl(3,4)=-vf(k,3,2)+vf(k,3,4)
85
86
87 ml(1,1)= vm(k,1,1)+vm(k,1,3)
88 ml(1,2)= vm(k,1,2)+vm(k,1,4)
89 ml(1,3)=-vm(k,1,1)+vm(k,1,3)
90 ml(1,4)=-vm(k,1,2)+vm(k,1,4)
91
92 ml(2,1)= vm(k,2,1)+vm(k,2,3)
93 ml(2,2)= vm(k,2,2)+vm(k,2,4)
94 ml(2,3)=-vm(k,2,1)+vm(k,2,3)
95 ml(2,4)=-vm(k,2,2)+vm(k,2,4)
96
97
98
99
100
101
102 f11(k)= vq(k,1,1)*fl(1,1)+vq(k,1,2)*fl(2,1)+vq(k,1,3)*fl(3,1)
103 m11(k)= vq(k,1,1)*ml(1,1)+vq(k,1,2)*ml(2,1)
104
105 f21(k)= vq(k,2,1)*fl(1,1)+vq(k,2,2)*fl(2,1)+vq(k,2,3)*fl(3,1)
106 m21(k)= vq(k,2,1)*ml(1,1)+vq(k,2,2)*ml(2,1)
107
108 f31(k)= vq(k,3,1)*fl(1,1)+vq(k,3,2)*fl(2,1)+vq(k,3,3)*fl(3,1)
109 m31(k)= vq(k,3,1)*ml(1,1)+vq(k,3,2)*ml(2,1)
110
111
112
113 f12(k)= vq(k,1,1)*fl(1,2)+vq(k,1,2)*fl(2,2)+vq(k,1,3)*fl(3,2)
114 m12(k)= vq(k,1,1)*ml(1,2)+vq(k,1,2)*ml(2,2)
115
116 f22(k)= vq(k,2,1)*fl(1,2)+vq(k,2,2)*fl(2,2)+vq(k,2,3)*fl(3,2)
117 m22(k)= vq(k,2,1)*ml(1,2)+vq(k,2,2)*ml(2,2)
118
119 f32(k)= vq(k,3,1)*fl(1,2)+vq(k,3,2)*fl(2,2)+vq(k,3,3)*fl(3,2)
120 m32(k)= vq(k,3,1)*ml(1,2)+vq(k,3,2)*ml(2,2)
121
122
123
124 f13(k)= vq(k,1,1)*fl(1,3)+vq(k,1,2)*fl(2,3)+vq(k,1,3)*fl(3,3)
125 m13(k)= vq(k,1,1)*ml(1,3)+vq(k,1,2)*ml(2,3)
126
127 f23(k)= vq(k,2,1)*fl(1,3)+vq(k,2,2)*fl(2,3)+vq(k,2,3)*fl(3,3)
128 m23(k)= vq(k,2,1)*ml(1,3)+vq(k,2,2)*ml(2,3)
129
130 f33(k)= vq(k,3,1)*fl(1,3)+vq(k,3,2)*fl(2,3)+vq(k,3,3)*fl(3,3)
131 m33(k)= vq(k,3,1)*ml(1,3)+vq(k,3,2)*ml(2,3)
132
133
134
135 f14(k)= vq(k,1,1)*fl(1,4)+vq(k,1,2)*fl(2,4)+vq(k,1,3)*fl(3,4)
136 m14(k)= vq(k,1,1)*ml(1,4)+vq(k,1,2)*ml(2,4)
137
138 f24(k)= vq(k,2,1)*fl(1,4)+vq(k,2,2)*fl(2,4)+vq(k,2,3)*fl(3,4)
139 m24(k)= vq(k,2,1)*ml(1,4)+vq(k,2,2)*ml(2,4)
140
141 f34(k)= vq(k,3,1)*fl(1,4)+vq(k,3,2)*fl(2,4)+vq(k,3,3)*fl(3,4)
142 m34(k)= vq(k,3,1)*ml(1,4)+vq(k,3,2)*ml(2,4)
143
144 ENDDO
145
146 IF (isrot>0) THEN
147#include "vectorize.inc"
148 DO ep=jft,nplat
149 k =iplat(ep)
150 m11(k)= m11(k)+ vq(k,1,3)*vmz(ep,1)
151 m21(k)= m21(k)+ vq(k,2,3)*vmz(ep,1)
152 m31(k)= m31(k)+ vq(k,3,3)*vmz(ep,1)
153
154 m12(k)= m12(k)+ vq(k,1,3)*vmz(ep,2)
155 m22(k)= m22(k)+ vq(k,2,3)*vmz(ep,2)
156 m32(k)= m32(k)+ vq(k,3,3)*vmz(ep,2)
157
158 m13(k)= m13(k)+ vq(k,1,3)*vmz(ep,3)
159 m23(k)= m23(k)+ vq(k,2,3)*vmz(ep,3)
160 m33(k)= m33(k)+ vq(k,3,3)*vmz(ep,3)
161
162 m14(k)= m14(k)+ vq(k,1,3)*vmz(ep,4)
163 m24(k)= m24(k)+ vq(k,2,3)*vmz(ep,4)
164 m34(k)= m34(k)+ vq(k,3,3)*vmz(ep,4)
165 ENDDO
166 DO j=1,4
167#include "vectorize.inc"
168 DO ep=nplat+1,jlt
169 k=iplat(ep)
170 mlz(ep,1,j)= vqn(k,7,j)*vmz(ep,j)
171 mlz(ep,2,j)= vqn(k,8,j)*vmz(ep,j)
172 mlz(ep,3,j)= vqn(k,9,j)*vmz(ep,j)
173 ENDDO
174 ENDDO
175 ELSE
176 DO j=1,4
177 DO ep=nplat+1,jlt
178 mlz(ep,1,j)= zero
179 mlz(ep,2,j)= zero
180 mlz(ep,3,j)= zero
181 ENDDO
182 ENDDO
183 END IF
184
185 DO ep=nplat+1,jlt
186 k=iplat(ep)
187 mm(1,1)= vqn(k,1,1)*vm(k,1,1)+vqn(k,4,1)*vm(k,2,1)+mlz(ep,1,1)
188 mm(2,1)= vqn(k,2,1)*vm(k,1,1)+vqn(k,5,1)*vm(k,2,1)+mlz(ep,2,1)
189 mm(3,1)= vqn(k,3,1)*vm(k,1,1)+vqn(k,6,1)*vm(k,2,1)+mlz(ep,3,1)
190
191 mm(1,2)= vqn(k,1,2)*vm(k,1,2)+vqn(k,4,2)*vm(k,2,2)+mlz(ep,1,2)
192 mm(2,2)= vqn(k,2,2)*vm(k,1,2)+vqn(k,5,2)*vm(k,2,2)+mlz(ep,2,2)
193 mm(3,2)= vqn(k,3,2)*vm(k,1,2)+vqn(k,6,2)*vm(k,2,2)+mlz(ep,3,2)
194
195 mm(1,3)= vqn(k,1,3)*vm(k,1,3)+vqn(k,4,3)*vm(k,2,3)+mlz(ep,1,3)
196 mm(2,3)= vqn(k,2,3)*vm(k,1,3)+vqn(k,5,3)*vm(k,2,3)+mlz(ep,2,3)
197 mm(3,3)= vqn(k,3,3)*vm(k,1,3)+vqn(k,6,3)*vm(k,2,3)+mlz(ep,3,3)
198
199 mm(1,4)= vqn(k,1,4)*vm(k,1,4)+vqn(k,4,4)*vm(k,2,4)+mlz(ep,1,4)
200 mm(2,4)= vqn(k,2,4)*vm(k,1,4)+vqn(k,5,4)*vm(k,2,4)+mlz(ep,2,4)
201 mm(3,4)= vqn(k,3,4)*vm(k,1,4)+vqn(k,6,4)*vm(k,2,4)+mlz(ep,3,4)
202
203 z1 = corel(k,3,1)
204 ar(1)= -z1*(vf(k,2,1)-vf(k,2,2)+vf(k,2,3)-vf(k,2,4))
205 1 +corel(k,2,1)*vf(k,3,1)+mm(1,1)
206 2 +corel(k,2,2)*vf(k,3,2)+mm(1,2)
207 3 +corel(k,2,3)*vf(k,3,3)+mm(1,3)
208 4 +corel(k,2,4)*vf(k,3,4)+mm(1,4)
209 ar(2)= z1*(vf(k,1,1)-vf(k,1,2)+vf(k,1,3)-vf(k,1,4))
210 1 -corel(k,1,1)*vf(k,3,1)+mm(2,1)
211 2 -corel(k,1,2)*vf(k,3,2)+mm(2,2)
212 3 -corel(k,1,3)*vf(k,3,3)+mm(2,3)
213 4 -corel(k,1,4)*vf(k,3,4)+mm(2,4)
214 ar(3)=-corel(k,2,1)*vf(k,1,1)+corel(k,1,1)*vf(k,2,1)+mm(3,1)
215 1 -corel(k,2,2)*vf(k,1,2)+corel(k,1,2)*vf(k,2,2)+mm(3,2)
216 2 -corel(k,2,3)*vf(k,1,3)+corel(k,1,3)*vf(k,2,3)+mm(3,3)
217 3 -corel(k,2,4)*vf(k,1,4)+corel(k,1,4)*vf(k,2,4)+mm(3,4)
218
219 alr(1) =di(k,1)*ar(1)+di(k,4)*ar(2)+di(k,5)*ar(3)
220 alr(2) =di(k,4)*ar(1)+di(k,2)*ar(2)+di(k,6)*ar(3)
221 alr(3) =di(k,5)*ar(1)+di(k,6)*ar(2)+di(k,3)*ar(3)
222
223 c1 =z1*alr(2)
224 vf(k,1,1)= vf(k,1,1)-c1+corel(k,2,1)*alr(3)
225 vf(k,1,2)= vf(k,1,2)+c1+corel(k,2,2)*alr(3)
226 vf(k,1,3)= vf(k,1,3)-c1+corel(k,2,3)*alr(3)
227 vf(k,1,4)= vf(k,1,4)+c1+corel(k,2,4)*alr(3)
228
229 c1 =z1*alr(1)
230 vf(k,2,1)= vf(k,2,1)+c1-corel(k,1,1)*alr(3)
231 vf(k,2,2)= vf(k,2,2)-c1-corel(k,1,2)*alr(3)
232 vf(k,2,3)= vf(k,2,3)+c1-corel(k,1,3)*alr(3)
233 vf(k,2,4)= vf(k,2,4)-c1-corel(k,1,4)*alr(3)
234
235 DO j=1,4
236 vf(k,3,j)= vf(k,3,j)-corel(k,2,j)*alr(1)+corel(k,1,j)*alr(2)
237 mm(1,j)= mm(1,j)-alr(1)
238 mm(2,j)= mm(2,j)-alr(2)
239 mm(3,j)= mm(3,j)-alr(3)
240 ENDDO
241
242
243 f11(k)= vq(k,1,1)*vf(k,1,1)+vq(k,1,2)*vf(k,2,1)
244 1 +vq(k,1,3)*vf(k,3,1)
245 m11(k)= vq(k,1,1)*mm(1,1)+vq(k,1,2)*mm(2,1)+vq(k,1,3)*mm(3,1)
246
247 f21(k)= vq(k,2,1)*vf(k,1,1)+vq(k,2,2)*vf(k,2,1)
248 1 +vq(k,2,3)*vf(k,3,1)
249 m21(k)= vq(k,2,1)*mm(1,1)+vq(k,2,2)*mm(2,1)+vq(k,2,3)*mm(3,1)
250
251 f31(k)= vq(k,3,1)*vf(k,1,1)+vq(k,3,2)*vf(k,2,1)
252 1 +vq(k,3,3)*vf(k,3,1)
253 m31(k)= vq(k,3,1)*mm
254
255
256
257 f12(k)= vq(k,1,1)*vf(k,1,2)+vq(k,1,2)*vf(k,2,2)
258 1 +vq(k,1,3)*vf(k,3,2)
259 m12(k)= vq(k,1,1)*mm(1,2)+vq(k,1,2)*mm(2,2)+vq(k,1,3)*mm(3,2)
260
261 f22(k)= vq(k,2,1)*vf(k,1,2)+vq(k,2,2)*vf(k,2,2)
262 1 +vq(k,2,3)*vf(k,3,2)
263 m22(k)= vq(k,2,1)*mm(1,2)+vq(k,2,2)*mm(2,2)+vq(k,2,3)*mm(3,2)
264
265 f32(k)= vq(k,3,1)*vf(k,1,2)+vq(k,3,2)*vf(k,2,2)
266 1 +vq(k,3,3)*vf(k,3,2)
267 m32(k)= vq(k,3,1)*mm(1,2)+vq(k,3,2)*mm(2,2)+vq(k,3,3)*mm(3,2)
268
269
270
271 f13(k)= vq(k,1,1)*vf(k,1,3)+vq(k,1,2)*vf(k,2,3)
272 1 +vq(k,1,3)*vf(k,3,3)
273 m13(k)= vq(k,1,1)*mm(1,3)+vq(k,1,2)*mm(2,3)+vq(k,1,3)*mm(3,3)
274
275 f23(k)= vq(k,2,1)*vf(k,1,3)+vq(k,2,2)*vf(k,2,3)
276 1 +vq(k,2,3)*vf(k,3,3)
277 m23(k)= vq(k,2,1)*mm(1,3)+vq(k,2,2)*mm(2,3)+vq(k,2,3)*mm(3,3)
278
279 f33(k)= vq(k,3,1)*vf(k,1,3)+vq(k,3,2)*vf(k,2,3)
280 1 +vq(k,3,3)*vf(k,3,3)
281 m33(k)= vq(k,3,1)*mm(1,3)+vq(k,3,2)*mm(2,3)+vq(k,3,3)*mm(3,3)
282
283
284
285 f14(k)= vq(k,1,1)*vf(k,1,4)+vq(k,1,2)*vf(k,2,4)
286 1 +vq(k,1,3)*vf(k,3,4)
287 m14(k)= vq(k,1,1)*mm(1,4)+vq(k,1,2)*mm(2,4)+vq(k,1,3)*mm(3,4)
288
289 f24(k)= vq(k,2,1)*vf(k,1,4)+vq(k,2
290 1 +vq(k,2,3)*vf(k,3,4)
291 m24(k)= vq(k,2,1)*mm(1,4)+vq(k,2,2)*mm(2,4)+vq(k,2,3)*mm(3,4)
292
293 f34(k)= vq(k,3,1)*vf(k,1,4)+vq(k,3,2)*vf(k,2,4)
294 1 +vq(k,3,3)*vf(k,3,4)
295 m34(k)= vq(k,3,1)*mm(1,4)+vq(k,3,2)*mm(2,4)+vq(k,3,3)*mm(3,4)
296
297 ENDDO
298
299 DO ep=jft,jlt
300
301 f11(ep)=f11(ep)*off(ep)
302 f21(ep)=f21(ep)*off(ep)
303 f31(ep)=f31(ep)*off(ep)
304
305 f12(ep)=f12(ep)*off(ep)
306 f22(ep)=f22(ep)*off(ep)
307 f32(ep)=f32(ep)*off(ep)
308
309 f13(ep)= f13(ep)*off(ep)
310 f23(ep)= f23(ep)*off(ep)
311 f33(ep)= f33(ep)*off(ep)
312
313 f14(ep)= f14(ep)*off(ep)
314 f24(ep)= f24(ep)*off(ep)
315 f34(ep)= f34(ep)*off(ep)
316
317 m11(ep)=m11(ep)*off(ep)
318 m21(ep)=m21(ep)*off(ep)
319 m31(ep)=m31(ep)*off(ep)
320
321 m12(ep)=m12(ep)*off(ep)
322 m22(ep)=m22(ep)*off(ep)
323 m32(ep)=m32(ep)*off(ep)
324
325 m13(ep)= m13(ep)*off(ep)
326 m23(ep)= m23(ep)*off(ep)
327 m33(ep)= m33(ep)*off(ep)
328
329 m14(ep)= m14(ep)*off(ep)
330 m24(ep)= m24(ep)*off(ep)
331 m34(ep)= m34(ep)*off(ep)
332
333 ENDDO
334
335 RETURN