29 4 RY ,SX ,SY ,RX2 ,RY2 ,
30 5 SX2 ,SY2 ,RHX ,RHY ,SHX ,
31 6 SHY ,GS ,NPLAT ,IPLAT,
32 9 K11,K12,K13,K14,K22,K23,K24,K33,K34,K44,
33 A M11,M12,M13,M14,M22,M23,M24,M33,M34,M44,
34 B MF11,MF12,MF13,MF14,MF22,MF23,MF24,MF33,
35 C MF34,MF44,FM12,FM13,FM14,FM23,FM24,FM34)
40#include "implicit_f.inc"
46 INTEGER JFT,JLT,NPLAT,IPLAT(*) ,IKGEO
48 . VOL(*),HC(MVSIZ,2),RX(*) ,RY(*) ,SX(*) ,SY(*) ,
49 . RX2(*) ,RY2(*) ,SX2(*) ,SY2(*) ,
50 . RHX(MVSIZ,4) ,RHY(MVSIZ,4) ,SHX(MVSIZ,4) ,SHY(MVSIZ,4) ,
51 . K11(3,3,*),K12(3,3,*),K13(3,3,*),K14(3,3,*),
52 . K22(3,3,*),K23(3,3,*),K24(3,3,*),K33(3,3,*),
53 . m11(3,3,*),m12(3,3,*),m13(3,3,*),m14(3,3,*),
54 . m22(3,3,*),m23(3,3,*),m24(3,3,*),m33(3,3,*),
55 . mf11(3,3,*),mf12(3,3,*),mf13(3,3,*),mf14(3,3,*),
56 . mf22(3,3,*),mf23(3,3,*),mf24(3,3,*),mf33(3,3,*),
57 . fm12(3,3,*),fm13(3,3,*),fm14(3,3,*),
58 . fm23(3,3,*),fm24(3,3,*),fm34(3,3,*),
59 . k34(3,3,*),k44(3,3,*),m34(3,3,*),m44(3,3,*),
60 . mf34(3,3,*),mf44(3,3,*),gs(*)
69 . CS1(MVSIZ),CS2(MVSIZ),CS3(MVSIZ),DS1(MVSIZ),QS1(MVSIZ),
70 . C1,C2,DS2(MVSIZ),QS2(),M1C1(2,2,MVSIZ),
71 . M1C2(2,2,MVSIZ),M1C3(2,2,MVSIZ),M2C1(2,2,MVSIZ),
72 . M2C2(2,2,MVSIZ),M2C3(2,2,MVSIZ),
73 . M1C11(2,2,MVSIZ),M1C12(2,2,MVSIZ),M1C13(2,2,MVSIZ),
74 . m1c14(2,2,mvsiz),m1c22(2,2,mvsiz),m1c23(2,2,mvsiz),
75 . m1c24(2,2,mvsiz),m1c33(2,2,mvsiz),m1c34(2,2,mvsiz),
77 . m2c11(2,2,mvsiz),m2c12(2,2,mvsiz),m2c13(2,2,mvsiz),
78 . m2c14(2,2,mvsiz),m2c22(2,2,mvsiz),m2c23(2,2,mvsiz),
79 . m2c24(2,2,mvsiz),m2c33(2,2,mvsiz),m2c34(2,2,mvsiz),
81 . m3c11(2,2,mvsiz),m3c12(2,2,mvsiz),m3c13(2,2,mvsiz),
82 . m3c14(2,2,mvsiz),m3c22(2,2,mvsiz),m3c23(2,2,mvsiz),
83 . m3c24(2,2,mvsiz),m3c33(2,2,mvsiz),m3c34(2,2,mvsiz),
84 . m3c44(2,2,mvsiz),dcx(mvsiz),dcy(mvsiz),c11,
85 . dh(mvsiz),hs1(mvsiz),hs2(mvsiz)
88#include "vectorize.inc"
98 cs1(ep) = dcx(ep)*sy2(ep)+dcy(ep)*sx2(ep)
99 cs2(ep) = dcx(ep)*ry2(ep)+dcy(ep)*rx2(ep)
100 cs3(ep) =-dcx(ep)*sy(ep)*ry(ep)-dcy(ep)*sx(ep)*rx(ep)
110 hs1(ep) = third*dh(ep)*(sy2(ep)+sx2(ep))
111 hs2(ep) = third*dh(ep)*(ry2(ep)+rx2(ep))
116 k11(3,3,ep) = cs1(ep)+cs2(ep)+cs3(ep)+cs3(ep)
118 k22(3,3,ep) = cs1(ep)+cs2(ep)-cs3(ep)-cs3(ep)
119 k33(3,3,ep) = k11(3,3,ep)
120 k44(3,3,ep) = k22(3,3,ep)
122 k12(3,3,ep) =-cs1(ep)+cs2(ep)+cs3(ep)-cs3(ep)
124 k13(3,3,ep) =-cs1(ep)-cs2(ep)-cs3(ep)-cs3(ep)
126 k14(3,3,ep) = - k12(3,3,ep)
128 k23(3,3,ep) = cs1(ep)-cs2(ep)+cs3(ep)-cs3(ep)
130 k24(3,3,ep) =-k22(3,3,ep)
132 k34(3,3,ep) =k12(3,3,ep)
137 k11(3,3,ep) = k11(3,3,ep)+c11
138 k22(3,3,ep) = k22(3,3,ep)+c11
139 k33(3,3,ep) = k33(3,3,ep)+c11
140 k44(3,3,ep) = k44(3,3,ep)+c11
141 k12(3,3,ep) = k12(3,3,ep)-c11
142 k13(3,3,ep) = k13(3,3,ep)+c11
143 k14(3,3,ep) = k14(3,3,ep)-c11
144 k23(3,3,ep) = k23(3,3,ep)-c11
145 k24(3,3,ep) = k24(3,3,ep)+c11
146 k34(3,3,ep) = k34(3,3,ep)-c11
152 m1c1(1,1,ep) = rhy(ep,1)*rhy(ep,1)
153 m1c1(2,2,ep) = rhx(ep,1)*rhx(ep,1)
154 m1c1(1,2,ep) = -rhy(ep,1)*rhx(ep,1)
155 m1c1(2,1,ep) = m1c1(1,2,ep)
157 m1c2(1,1,ep) = rhy(ep,1)*rhy(ep,3)
158 m1c2(2,2,ep) = rhx(ep,1)*rhx(ep,3)
159 m1c2(1,2,ep) = -rhy(ep,1)*rhx(ep,3)
160 m1c2(2,1,ep) = -rhy(ep,3)*rhx(ep,1)
162 m1c3(1,1,ep) = rhy(ep,3)*rhy(ep,3)
163 m1c3(2,2,ep) = rhx(ep,3)*rhx(ep,3)
164 m1c3(1,2,ep) = -rhy(ep,3)*rhx(ep,3)
165 m1c3(2,1,ep) = m1c3(1,2,ep)
166 ds1(ep) = cs1(ep)-hs1(ep)
167 qs1(ep) = cs1(ep)+hs1(ep)
174 m1c11(i,j,ep) = m1c1(i,j,ep)*qs1(ep)
175 m1c12(i,j,ep) = m1c11(i,j,ep)
176 m1c22(i,j,ep) = m1c11(i,j,ep)
178 m1c33(i,j,ep) = m1c3(i,j,ep)*qs1(ep)
179 m1c34(i,j,ep) = m1c33(i,j,ep)
180 m1c44(i,j,ep) = m1c33(i,j,ep)
185 m1c12(2,1,ep)=m1c12(1,2,ep)
186 m1c34(2,1,ep)=m1c34(1,2,ep)
192 m1c13(i,j,ep) = m1c2(i,j,ep)*ds1(ep)
193 m1c14(i,j,ep) = m1c13(i,j,ep)
194 m1c23(i,j,ep) = m1c13(i,j,ep)
195 m1c24(i,j,ep) = m1c13(i,j,ep)
202 m2c1(1,1,ep) = shy(ep,1)*shy(ep,1)
203 m2c1(2,2,ep) = shx(ep,1)*shx(ep,1)
204 m2c1(1,2,ep) = -shy(ep,1)*shx(ep,1)
205 m2c1(2,1,ep) = m2c1(1,2,ep)
207 m2c2(1,1,ep) = shy(ep,1)*shy(ep,2)
208 m2c2(2,2,ep) = shx(ep,1)*shx(ep,2)
210 m2c2(1,2,ep) = -shy(ep,1)*shx(ep,2)
211 m2c2(2,1,ep) = -shy(ep,2)*shx(ep,1)
213 m2c3(1,1,ep) = shy(ep,2)*shy(ep,2)
214 m2c3(2,2,ep) = shx(ep,2)*shx(ep,2)
215 m2c3(1,2,ep) = -shy(ep,2)*shx(ep,2)
216 m2c3(2,1,ep) = m2c3(1,2,ep)
217 ds2(ep) = cs2(ep)-hs2(ep)
218 qs2(ep) = cs2(ep)+hs2(ep)
224 m2c11(i,j,ep) = m2c1(i,j,ep)*qs2(ep)
225 m2c14(i,j,ep) = m2c11(i,j,ep)
226 m2c44(i,j,ep) = m2c11(i,j,ep)
228 m2c22(i,j,ep) = m2c3(i,j,ep)*qs2(ep)
229 m2c23(i,j,ep) = m2c22(i,j,ep)
230 m2c33(i,j,ep) = m2c22(i,j,ep)
235 m2c14(2,1,ep)=m2c14(1,2,ep)
236 m2c23(2,1,ep)=m2c23(1,2,ep)
242 m2c12(i,j,ep) = m2c2(i,j,ep)*ds2(ep)
243 m2c13(i,j,ep) = m2c12(i,j,ep)
251 m2c24(i,j,ep) = m2c12(j,i,ep)
252 m2c34(i,j,ep) = m2c12(j,i,ep)
259 m3c11(1,1,ep) =(rhy(ep,1)*shy(ep,1)+rhy(ep,1)*shy(ep,1))*cs3(ep)
260 m3c12(1,1,ep) =(rhy(ep,1)*shy(ep,2)+rhy(ep,2)*shy(ep,1))*cs3(ep)
261 m3c13(1,1,ep) =(rhy(ep,1)*shy(ep,3)+rhy(ep,3)*shy(ep,1))*cs3(ep)
262 m3c14(1,1,ep) =(rhy(ep,1)*shy(ep,4)+rhy(ep,4)*shy(ep,1))*cs3(ep)
263 m3c22(1,1,ep) =(rhy(ep,2)*shy(ep,2)+rhy(ep,2)*shy(ep,2))*cs3
265 m3c24(1,1,ep) =(rhy(ep,2)*shy(ep,4)+rhy(ep,4)*shy(ep,2))*cs3(ep)
268 m3c44(1,1,ep) =(rhy(ep,4)*shy(ep,4)+rhy(ep,4)*shy(ep,4))*cs3(ep)
271 m3c11(2,2,ep) =(rhx(ep,1)*shx(ep,1)+rhx(ep
272 m3c12(2,2,ep) =(rhx(ep,1)*shx(ep,2)+rhx(ep,2)*shx(ep,1))*cs3(ep)
273 m3c13(2,2,ep) =(rhx(ep,1)*shx(ep,3)+rhx(ep,3)*shx(ep,1))*cs3(ep)
274 m3c14(2,2,ep) =(rhx(ep,1)*shx(ep,4)+rhx(ep,4)*shx(ep,1))*cs3(ep)
275 m3c22(2,2,ep) =(rhx(ep,2)*shx(ep,2)+rhx(ep,2)*shx(ep,2))*cs3(ep)
276 m3c23(2,2,ep) =(rhx(ep,2)*shx(ep,3)+rhx(ep,3)*shx(ep,2))*cs3(ep)
277 m3c24(2,2,ep) =(rhx(ep,2)*shx(ep,4)+rhx(ep,4)*shx(ep,2))*cs3(ep)
278 m3c33(2,2,ep) =(rhx(ep,3)*shx(ep,3)+rhx(ep,3)*shx(ep,3))*cs3(ep)
279 m3c34(2,2,ep) =(rhx(ep,3)*shx(ep,4)+rhx(ep,4)*shx(ep,3))*cs3(ep)
280 m3c44(2,2,ep) =(rhx(ep,4)*shx(ep,4)+rhx(ep,4)*shx(ep,4))*cs3(ep)
283 m3c11(1,2,ep)=(-rhy(ep,1)*shx(ep,1)-rhx(ep,1)*shy(ep,1))*cs3(ep)
284 m3c12(1,2,ep)=(-rhy(ep,1)*shx(ep,2)-rhx(ep,2)*shy(ep,1))*cs3(ep)
285 m3c13(1,2,ep)=(-rhy(ep,1)*shx(ep,3)-rhx(ep,3)*shy(ep,1))*cs3(ep)
286 m3c14(1,2,ep)=(-rhy(ep,1)*shx(ep,4)-rhx(ep,4)*shy(ep,1))*cs3(ep)
287 m3c22(1,2,ep)=(-rhy(ep,2)*shx(ep,2)-rhx(ep,2)*shy(ep,2))*cs3(ep)
288 m3c23(1,2,ep)=(-rhy(ep,2)*shx(ep,3)-rhx(ep,3)*shy(ep,2))*cs3(ep)
289 m3c24(1,2,ep)=(-rhy(ep,2)*shx(ep,4)-rhx(ep,4)*shy(ep,2))*cs3(ep)
290 m3c33(1,2,ep)=(-rhy(ep,3)*shx(ep,3)-rhx(ep,3)*shy(ep,3))*cs3(ep)
291 m3c34(1,2,ep)=(-rhy(ep,3)*shx(ep,4)-rhx(ep,4)*shy(ep,3))*cs3(ep)
292 m3c44(1,2,ep)=(-rhy(ep,4)*shx(ep,4)-rhx(ep,4)*shy(ep,4))*cs3(ep)
296 m3c12(2,1,ep)=(-rhx(ep,1)*shy(ep,2)-rhy(ep,2)*shx(ep,1))*cs3(ep)
297 m3c13(2,1,ep)=(-rhx(ep,1)*shy(ep,3)-rhy(ep,3)*shx(ep,1))*cs3(ep)
298 m3c14(2,1,ep)=(-rhx(ep,1)*shy(ep,4)-rhy(ep,4)*shx(ep,1))*cs3(ep)
299 m3c22(2,1,ep) = m3c22(1,2,ep)
300 m3c23(2,1,ep)=(-rhx(ep,2)*shy(ep,3)-rhy(ep,3)*shx(ep,2))*cs3(ep)
301 m3c24(2,1,ep)=(-rhx(ep,2)*shy(ep,4)-rhy(ep,4)*shx(ep,2))*cs3(ep)
302 m3c33(2,1,ep) = m3c33(1,2,ep)
303 m3c34(2,1,ep)=(-rhx(ep,3)*shy(ep,4)-rhy(ep,4)*shx(ep,3))*cs3(ep)
304 m3c44(2,1,ep) = m3c44(1,2,ep)
310 m11(i,j,ep)=m11(i,j,ep)+
311 1 m1c11(i,j,ep)+m2c11(i,j,ep)+m3c11(i,j,ep)
313 1 m1c22(i,j,ep)+m2c22(i,j,ep)+m3c22(i,j,ep)
314 m33(i,j,ep)=m33(i,j,ep)+
315 1 m1c33(i,j,ep)+m2c33(i,j,ep)+m3c33(i,j,ep)
316 m44(i,j,ep)=m44(i,j,ep)+
317 1 m1c44(i,j,ep)+m2c44(i,j,ep)+m3c44(i,j,ep)
325 m12(i,j,ep)=m12(i,j,ep)+
326 1 m1c12(i,j,ep)+m2c12(i,j,ep)+m3c12(i,j,ep)
327 m13(i,j,ep)=m13(i,j,ep)+
328 1 m1c13(i,j,ep)+m2c13(i,j,ep)+m3c13(i,j,ep)
329 m14(i,j,ep)=m14(i,j,ep)+
330 1 m1c14(i,j,ep)+m2c14(i,j,ep)+m3c14(i,j,ep)
331 m23(i,j,ep)=m23(i,j,ep)+
332 1 m1c23(i,j,ep)+m2c23(i,j,ep)+m3c23(i,j,ep)
333 m24(i,j,ep)=m24(i,j,ep)+
334 1 m1c24(i,j,ep)+m2c24(i,j,ep)+m3c24(i,j,ep)
335 m34(i,j,ep)=m34(i,j,ep)+
336 1 m1c34(i,j,ep)+m2c34(i,j,ep)+m3c34(i,j,ep)
343 mf11(3,1,ep)= qs1(ep)*rhy(ep,1)+qs2(ep)*shy(ep,1)
344 1 +cs3(ep)*(shy(ep,1)+rhy(ep,1))
346 1 +cs3(ep)*(shy(ep,2)+rhy(ep,2))
348 1 +cs3(ep)*(shy(ep,3)+rhy(ep,3))
349 mf14(3,1,ep)= ds1(ep)*rhy(ep,4)+qs2(ep)*shy(ep,4)
350 1 +cs3(ep)*(shy(ep,4)+rhy(ep,4))
351 mf22(3,1,ep)=-qs1(ep)*rhy(ep,2)+qs2(ep)*shy(ep,2)
352 1 -cs3(ep)*(shy(ep,2)-rhy(ep,2))
353 mf23(3,1,ep)=-ds1(ep)*rhy(ep,3)+qs2(ep)*shy(ep,3)
354 1 -cs3(ep)*(shy(ep,3)-rhy(ep,3))
355 mf24(3,1,ep)=-ds1(ep)*rhy(ep,4)+ds2(ep)*shy(ep,4)
356 1 -cs3(ep)*(shy(ep,4)-rhy(ep,4))
357 mf33(3,1,ep)=-qs1(ep)*rhy(ep,3)-qs2(ep)*shy(ep,3)
358 1 -cs3(ep)*(shy(ep,3)+rhy(ep,3))
359 mf34(3,1,ep)=-qs1(ep)*rhy(ep,4)-ds2(ep)*shy(ep,4)
360 1 -cs3(ep)*(shy(ep,4)+rhy(ep,4))
361 mf44(3,1,ep)= qs1(ep)*rhy(ep,4)-qs2(ep)*shy(ep,4)
362 1 +cs3(ep)*(shy(ep,4)-rhy(ep,4))
366 mf11(3,2,ep)=-qs1(ep)*rhx(ep,1)-qs2(ep)*shx(ep,1)
367 1 -cs3(ep)*(shx(ep,1)+rhx(ep,1))
368 mf12(3,2,ep)=-qs1(ep)*rhx(ep,2)-ds2(ep)*shx(ep,2)
369 1 -cs3(ep)*(shx(ep,2)+rhx(ep,2))
370 mf13(3,2,ep)=-ds1(ep)*rhx(ep,3)-ds2(ep)*shx(ep,3)
371 1 -cs3(ep)*(shx(ep,3)+rhx(ep,3))
372 mf14(3,2,ep)=-ds1(ep)*rhx(ep,4)-qs2(ep)*shx(ep,4)
373 1 -cs3(ep)*(shx(ep,4)+rhx(ep,4))
374 mf22(3,2,ep)= qs1(ep)*rhx(ep,2)-qs2(ep)*shx(ep,2)
376 mf23(3,2,ep)= ds1(ep)*rhx(ep,3)-qs2(ep)*shx(ep,3)
377 1 +cs3(ep)*(shx(ep,3)-rhx(ep,3))
378 mf24(3,2,ep)= ds1(ep)*rhx(ep,4)-ds2(ep)*shx(ep,4)
379 1 +cs3(ep)*(shx(ep,4)-rhx(ep,4))
380 mf33(3,2,ep)= qs1(ep)*rhx(ep,3)+qs2(ep)*shx(ep,3)
381 1 +cs3(ep)*(shx(ep,3)+rhx(ep,3))
382 mf34(3,2,ep)= qs1(ep)*rhx(ep,4)+ds2(ep)*shx(ep,4)
383 1 +cs3(ep)*(shx(ep,4)+rhx(ep,4))
384 mf44(3,2,ep)=-qs1(ep)*rhx(ep,4)+qs2(ep)*shx(ep,4)
385 1 -cs3(ep)*(shx(ep,4)-rhx(ep,4))
390 fm12(1,3,ep)=-qs1(ep)*rhy(ep,1)+ds2(ep)*shy(ep,1)
391 1 +cs3(ep)*(rhy(ep,1)-shy(ep,1))
392 fm13(1,3,ep)=-ds1(ep)*rhy(ep,1)-ds2(ep)*shy(ep,1)
393 1 -cs3(ep)*(rhy(ep,1)+shy(ep,1))
394 fm23(1,3,ep)=-ds1(ep)*rhy(ep,2)-qs2(ep)*shy(ep,2)
395 1 -cs3(ep)*(rhy(ep,2)+shy(ep,2))
396 fm14(1,3,ep)= ds1(ep)*rhy(ep,1)-qs2(ep)*shy(ep,1)
397 1 -cs3(ep)*(rhy(ep,1)-shy(ep,1))
398 fm24(1,3,ep)= ds1(ep)*rhy(ep,2)-ds2(ep)*shy(ep,2)
399 1 -cs3(ep)*(rhy(ep,2)-shy(ep,2))
400 fm34(1,3,ep)= qs1(ep)*rhy(ep,3)-ds2(ep)*shy(ep,3)
401 1 -cs3(ep)*(rhy(ep,3)-shy(ep,3))
403 fm12(2,3,ep)= qs1(ep)*rhx(ep,1)-ds2(ep)*shx(ep,1)
404 1 -cs3(ep)*(rhx(ep,1)-shx(ep,1))
405 fm13(2,3,ep)= ds1(ep)*rhx(ep,1)+ds2(ep)*shx(ep,1)
406 1 +cs3(ep)*(rhx(ep,1)+shx(ep,1))
407 fm23(2,3,ep)= ds1(ep)*rhx(ep,2)+qs2(ep)*shx(ep,2)
408 1 +cs3(ep)*(rhx(ep,2)+shx(ep,2))
409 fm14(2,3,ep)=-ds1(ep)*rhx(ep,1)+qs2(ep)*shx(ep,1)
410 1 +cs3(ep)*(rhx(ep,1)-shx(ep,1))
411 fm24(2,3,ep)=-ds1(ep)*rhx(ep,2)+ds2(ep)*shx(ep,2)
412 1 +cs3(ep)*(rhx(ep,2)-shx(ep,2))
413 fm34(2,3,ep)=-qs1(ep)*rhx(ep,3)+ds2(ep)*shx(ep,3)
414 1 +cs3(ep)*(rhx(ep,3)-shx(ep,3))