28 SUBROUTINE czlkec3(JFT ,JLT ,VOL ,THK0 ,THK2 ,
29 2 HM ,HF ,HZ ,A_I ,Z1 ,
30 3 PX1 ,PX2 ,PY1 ,PY2 ,NPLAT,
32 6 K11,K12,K13,K14,K22,K23,K24,K33,K34,K44,
33 7 M11,M12,M13,M14,M22,M23,M24,M33,M34,M44,
34 8 MF11,MF12,MF13,MF14,MF22,MF23,MF24,MF33,
35 9 MF34,MF44,FM12,FM13,FM14,FM23,FM24,FM34,
36 A IORTH,HMOR,HFOR,HMFOR)
41#include "implicit_f.inc"
48 INTEGER JFT,JLT,NPLAT,IPLAT(*),IKGEO,IORTH
50 . VOL(*),THK0(*),THK2(*),A_I(*),Z1(*),
51 . HM(MVSIZ,4),HF(MVSIZ,4),HZ(*),
52 . PX1(*) ,PX2(*) ,PY1(*) ,PY2(*) ,
53 . K11(3,3,*),K12(3,3,*),K13(3,3,*),K14(3,3,*),
54 . K22(3,3,*),K23(3,3,*),K24(3,3,*),K33(3,3,*),
55 . M11(3,3,*),M12(3,3,*),M13(3,3,*),M14(3,3,*),
56 . m22(3,3,*),m23(3,3,*),m24
57 . mf11(3,3,*),mf12(3,3,*),mf13(3,3,*),mf14(3,3,*),
58 . mf22(3,3,*),mf23(3,3,*),mf24(3,3,*),mf33(3,3,*),
59 . fm12(3,3,*),fm13(3,3,*),fm14(3,3,*),
60 . fm23(3,3,*),fm24(3,3,*),fm34(3,3,*),
61 . k34(3,3,*),k44(3,3,*),m34(3,3,*),m44(3,3,*),
62 . mf34(3,3,*),mf44(3,3,*),dhz(*),hmor(mvsiz,2),hfor(mvsiz,2),
72 . DM(2,2,MVSIZ),DF(2,2,MVSIZ),
73 . C1,C2,GM(MVSIZ),GF(MVSIZ),DG(MVSIZ),
74 . PX1PX1(MVSIZ),PX1PX2(MVSIZ) ,PX2PX2(MVSIZ),
79 . df6_2(mvsiz),dmf(3,3,mvsiz),
80 . bxibxj(2,2,mvsiz),bxibyj(2,2,mvsiz),byibxj(2,2,mvsiz),
81 . byibyj(2,2,mvsiz),fac2z
84#include "vectorize.inc"
104 px1px1(ep) = px1(ep)*px1(ep)
105 px1px2(ep) = px1(ep)*px2(ep)
106 px2px2(ep) = px2(ep)*px2(ep)
107 px1py1(ep) = px1(ep)*py1(ep)
108 px1py2(ep) = px1(ep)*py2(ep)
109 px2py1(ep) = px2(ep)*py1(ep)
110 px2py2(ep) = px2(ep)*py2(ep)
111 py1py1(ep) = py1(ep)*py1(ep)
112 py1py2(ep) = py1(ep)*py2(ep)
113 py2py2(ep) = py2(ep)*py2(ep)
117 k11(1,1,ep)=dm(1,1,ep)*px1px1(ep)+gm(ep)*py1py1(ep)
118 k12(1,1,ep)=dm(1,1,ep)*px1px2(ep)+gm(ep)*py1py2(ep)
119 k22(1,1,ep)=dm(1,1,ep)*px2px2(ep)+gm(ep)*py2py2(ep)
121 m11(2,2,ep)=df(1,1,ep)*px1px1(ep)+gf(ep)*py1py1(ep)
122 m12(2,2,ep)=df(1,1,ep)*px1px2(ep)+gf(ep)*py1py2(ep)
123 m22(2,2,ep)=df(1,1,ep)*px2px2(ep)+gf(ep)*py2py2(ep)
127 k11(1,2,ep)=(dm(1,2,ep)+gm(ep))*px1py1(ep)
128 k12(1,2,ep)=dm(1,2,ep)*px1py2(ep)+gm(ep)*px2py1(ep)
129 k22(1,2,ep)=(dm(1,2,ep)+gm(ep))*px2py2(ep)
130 k12(2,1,ep)=dm(1,2,ep)*px2py1(ep)+gm(ep)*px1py2(ep)
132 m11(1,2,ep)=-(df(1,2,ep)+gf(ep))*px1py1(ep)
133 m12(1,2,ep)=-df(1,2,ep)*px2py1(ep)-gf(ep)*px1py2(ep)
134 m22(1,2,ep)=-df(1,2,ep)*px2py2(ep)-gf(ep)*px2py2(ep)
135 m12(2,1,ep)=-df(1,2,ep)*px1py2(ep)-gf(ep)*px2py1
139 k11(2,2,ep)=dm(2,2,ep)*py1py1(ep)+gm(ep)*px1px1(ep)
140 k12(2,2,ep)=dm(2,2,ep)*py1py2(ep)+gm(ep)*px1px2(ep)
141 k22(2,2,ep)=dm(2,2,ep)*py2py2(ep)+gm(ep)*px2px2(ep)
143 m11(1,1,ep)=df(2,2,ep)*py1py1(ep)+gf(ep)*px1px1(ep)
144 m12(1,1,ep)=df(2,2,ep)*py1py2(ep)+gf(ep)*px1px2(ep)
145 m22(1,1,ep)=df(2,2,ep)*py2py2(ep)+gf(ep)*px2px2(ep)
148#include
"vectorize.inc"
165 k11(1,1,ep)=k11(1,1,ep)+px1py1(ep)*dm5_2(ep)
166 k11(1,2,ep)=k11(1,2,ep)+
167 1 px1px1(ep)*dm5(ep)+py1py1(ep)*dm6(ep)
168 k11(2,2,ep)=k11(2,2,ep)+px1py1(ep)*dm6_2(ep)
169 k12(1,1,ep)=k12(1,1,ep)+(px1py2(ep)+px2py1(ep))*dm5(ep)
170 c1 = px1px2(ep)*dm5(ep)+py1py2(ep)*dm6(ep)
171 k12(1,2,ep)=k12(1,2,ep)+c1
172 k12(2,1,ep)=k12(2,1,ep)+c1
173 k12(2,2,ep)=k12(2,2,ep)+(px1py2(ep)+px2py1(ep))*dm6(ep)
174 k22(1,1,ep)=k22(1,1,ep)+px2py2(ep)*dm5_2(ep)
175 k22(1,2,ep)=k22(1,2,ep)+
176 1 px2px2(ep)*dm5(ep)+py2py2(ep)*dm6(ep)
177 k22(2,2,ep)=k22(2,2,ep)+px2py2(ep)*dm6_2(ep)
179 m11(1,1,ep)=m11(1,1,ep)+px1py1(ep)*df6_2(ep)
180 m11(1,2,ep)=m11(1,2,ep)-
181 1 px1px1(ep)*df5(ep)-py1py1(ep)*df6(ep)
182 m11(2,2,ep)=m11(2,2,ep)+px1py1(ep)*df5_2(ep)
183 m12(1,1,ep)=m12(1,1,ep)+(px1py2(ep)+px2py1(ep))*df6(ep)
184 c2 = px1px2(ep)*df5(ep)+py1py2(ep)*df6(ep)
185 m12(1,2,ep)=m12(1,2,ep)-c2
186 m12(2,1,ep)=m12(2,1,ep)-c2
187 m12(2,2,ep)=m12(2,2,ep)+(px1py2(ep)+px2py1(ep))*df5(ep)
188 m22(1,1,ep)=m22(1,1,ep)+px2py2(ep)*df6_2(ep)
189 m22(1,2,ep)=m22(1,2,ep)-
190 1 px2px2(ep)*df5(ep)-py2py2(ep)*df6(ep)
191 m22(2,2,ep)=m22(2,2,ep)+px2py2(ep)*df5_2(ep)
196 dmf(1,1,m)=hmfor(ep,1)*c2
197 dmf(2,2,m)=hmfor(ep,2)*c2
198 dmf(1,2,m)=hmfor(ep,3)*c2
199 dmf(1,3,m)=hmfor(ep,5)*c2
200 dmf(2,3,m)=hmfor(ep,6)*c2
201 dmf(2,1,m)=dmf(1,2,m)
202 dmf(3,1,m)=dmf(1,3,m)
203 dmf(3,2,m)=dmf(2,3,m)
204 dmf(3,3,m)=hmfor(ep,4)*c2
209 1 -dmf(1,2,ep)*px1py1(ep)-dmf(1,3,ep)*px1px1(ep)
210 2 -dmf(2,3,ep)*py1py1(ep)-dmf(3,3,ep)*px1py1(ep)
212 1 dmf(1,1,ep)*px1px1(ep)+dmf(1,3,ep)*px1py1(ep)
213 2 +dmf(1,3,ep)*px1py1(ep)+dmf(3,3,ep)*py1py1(ep
215 1 -dmf(2,2,ep)*py1py1(ep)-dmf(2,3,ep)*px1py1(ep)
216 2 -dmf(2,3,ep)*px1py1(ep)-dmf(3,3,ep)*px1px1(ep)
218 1 dmf(1,2,ep)*px1py1(ep)+dmf(2,3,ep)*py1py1(ep)
219 2 +dmf(1,3,ep)*px1px1(ep)+dmf(3,3,ep)*px1py1(ep)
223 1 -dmf(1,2,ep)*px1py2(ep)-dmf(1,3,ep)*px1px2(ep)
227 1 dmf(1,1,ep)*px1px2(ep)+dmf(1,3,ep)*px1py2(ep)
229 2 +dmf(1,3,ep)*px2py1(ep)+dmf(3,3,ep)*py1py2(ep)
232 1 -dmf(2,2,ep)*py1py2(ep)-dmf(2,3,ep)*px2py1(ep)
233 2 -dmf(2,3,ep)*px1py2(ep)-dmf(3,3,ep)*px1px2(ep)
236 1 dmf(1,2,ep)*px2py1(ep)+dmf(2,3,ep)*py1py2(ep)
237 2 +dmf(1,3,ep)*px1px2(ep)+dmf(3,3,ep)*px1py2(ep)
241 1 -dmf(1,2,ep)*px2py2(ep)-dmf(1,3,ep)*px2px2(ep)
243 2 -dmf(2,3,ep)*py2py2(ep)-dmf
245 1 dmf(1,1,ep)*px2px2(ep)+dmf(1,3,ep)*px2py2(ep)
247 2 +dmf(1,3,ep)*px2py2(ep)+dmf(3,3,ep)*py2py2
250 1 -dmf(2,2,ep)*py2py2(ep)-dmf(2,3,ep)*px2py2(ep)
251 2 -dmf(2,3,ep)*px2py2(ep)-dmf(3,3,ep)*px2px2(ep)
254 1 dmf(1,2,ep)*px2py2(ep)+dmf(2,3,ep)*py2py2(ep)
255 2 +dmf(1,3,ep)*px2px2(ep)+dmf(3,3,ep)*px2py2(ep)
259 1 dmf(1,2,ep)*px2py1(ep)+dmf(1,3,ep)*px1px2(ep)
261 2 +dmf(2,3,ep)*py1py2(ep)+dmf(3,3,ep)*px1py2(ep)
263 1 -dmf(1,1,ep)*px1px2(ep)-dmf(1,3,ep)*px2py1(ep)
265 2 -dmf(1,3,ep)*px1py2(ep)-dmf(3,3,ep)*py1py2(ep)
268 1 dmf(2,2,ep)*py1py2(ep)+dmf(2,3,ep)*px1py2(ep)
269 2 +dmf(2,3,ep)*px2py1(ep)+dmf(3,3,ep)*px1px2(ep)
272 1 -dmf(1,2,ep)*px1py2(ep)-dmf(2,3,ep)*py1py2(ep)
273 2 -dmf(1,3,ep)*px1px2(ep)-dmf(3,3,ep)*px2py1(ep)
278 mf13(i,j,ep)=-mf11(i,j,ep)
279 mf14(i,j,ep)=-mf12(i,j,ep)
280 mf24(i,j,ep)=-mf22(i,j,ep)
281 mf33(i,j,ep)=mf11(i,j,ep)
282 mf34(i,j,ep)=mf12(i,j,ep)
283 mf44(i,j,ep)=mf22(i,j,ep)
291 fm12(i,j,ep)=-mf23(j,i,ep)
292 fm13(i,j,ep)=-mf11(j,i,ep)
293 fm14(i,j,ep)= mf23(j,i,ep)
294 fm23(i,j,ep)=-mf12(j,i,ep)
295 fm24(i,j,ep)=-mf22(j,i,ep)
296 fm34(i,j,ep)=-mf23(j,i,ep)
302#include "vectorize.inc"
305 facz(m) = four*z1(ep)*a_i(ep)
306 facz2(m) = facz(m)*facz(m)
312 dfz =df(1,1,ep)*facz2(ep)
313 k11(1,1,ep)=k11(1,1,ep)+dfz*py2py2(ep)+dg(ep)*px2px2(ep)
314 k12(1,1,ep)=k12(1,1,ep)+dfz
315 k22(1,1,ep)=k22(1,1,ep)+dfz*py1py1(ep)+dg(ep)*px1px1(ep)
319 dfz =df(1,2,ep)*facz2(ep)
320 k11(1,2,ep)=k11(1,2,ep)-(dfz+dg(ep))*px2py2(ep)
321 k12(1,2,ep)=k12(1,2,ep)-dfz*px1py2(ep)-dg(ep)*px2py1(ep)
322 k22(1,2,ep)=k22(1,2,ep)-(dfz+dg(ep))*px1py1(ep)
323 k12(2,1,ep)=k12(2,1,ep)-dfz*px2py1(ep)-dg(ep)*px1py2(ep)
327 dfz =df(2,2,ep)*facz2(ep)
328 k11(2,2,ep)=k11(2,2,ep)+dfz*px2px2(ep)+dg(ep)*py2py2(ep)
329 k12(2,2,ep)=k12(2,2,ep)+dfz*px1px2(ep)+dg(ep)*py1py2(ep)
330 k22(2,2,ep)=k22(2,2,ep)+dfz*px1px1(ep)+dg(ep)*py1py1(ep)
337 k11(1,1,ep)=k11(1,1,ep)-px2py2(ep)*df5_2(ep)*facz2(ep)
338 k11(1,2,ep)=k11(1,2,ep)+facz2(ep)*(
339 1 py2py2(ep)*df5(ep)+px2px2(ep)*df6(ep))
340 k11(2,2,ep)=k11(2,2,ep)-px2py2(ep)*df6_2(ep)*facz2(ep)
341 k12(1,1,ep)=k12(1,1,ep)-facz2(ep)*(
342 1 px1py2(ep)+px2py1(ep))*df5(ep)
343 c1 = (py1py2(ep)*df5(ep)+px1px2(ep)*df6(ep))*facz2(ep)
344 k12(1,2,ep)=k12(1,2,ep)+c1
345 k12(2,1,ep)=k12(2,1,ep)+c1
346 k12(2,2,ep)=k12(2,2,ep)-facz2(ep)*(
347 1 px1py2(ep)+px2py1(ep))*df6(ep)
348 k22(1,1,ep)=k22(1,1,ep)-px1py1(ep)*df5_2(ep)*facz2(ep)
349 k22(1,2,ep)=k22(1,2,ep)+facz2(ep)*(
350 1 py1py1(ep)*df5(ep)+px1px1(ep)*df6(ep))
351 k22(2,2,ep)=k22(2,2,ep)-px1py1(ep)*df6_2(ep)*facz2(ep)
359 bxibxj(1,1,ep)=fac2z*px1py2(ep)
360 bxibxj(1,2,ep)=facz(ep)*(px1py1(ep)+px2py2(ep))
361 bxibxj(2,1,ep)=bxibxj(1,2,ep)
362 bxibxj(2,2,ep)=fac2z*px2py1(ep)
364 bxibyj(1,1,ep)=facz(ep)*(-px1px2(ep)+py1py2(ep))
365 bxibyj(1,2,ep)=facz(ep)*(-px1px1(ep)+py2py2(ep))
366 bxibyj(2,1,ep)=facz(ep)*(-px2px2(ep)+py1py1(ep))
367 bxibyj(2,2,ep)=bxibyj(1,1,ep)
369 byibyj(1,1,ep)=-fac2z*px2py1(ep)
370 byibyj(1,2,ep)=facz(ep)*(-px1py1(ep)-px2py2(ep))
372 byibyj(2,2,ep)=-fac2z*px1py2(ep)
374 byibxj(1,1,ep)=facz(ep)*(py1py2(ep)-px1px2(ep))
375 byibxj(1,2,ep)=facz(ep)*(py1py1(ep)-px2px2(ep))
376 byibxj(2,1,ep)=facz(ep)*(py2py2(ep)-px1px1(ep))
377 byibxj(2,2,ep)=byibxj(1,1,ep)
381 k11(1,1,ep)=k11(1,1,ep)+dmf(1,1,ep)*bxibxj(1,1,ep)+
382 1 dmf(1,3,ep)*(bxibyj(1,1,ep)+byibxj(1,1,ep))+
383 2 dmf(3,3,ep)*byibyj(1,1,ep)
384 k11(1,2,ep)=k11(1,2,ep)+dmf(1,2,ep)*bxibyj(1,1,ep)+
385 1 dmf(1,3,ep)*bxibxj(1,1,ep)+
386 2 dmf(2,3,ep)*byibyj(1,1,ep)+
387 3 dmf(3,3,ep)*byibxj(1,1,ep)
388 k11(2,2,ep)=k11(2,2,ep)+dmf(2,2,ep)*byibyj(1,1,ep)+
389 1 dmf(2,3,ep)*(bxibyj(1,1,ep
390 2 dmf(3,3,ep)*bxibxj(1,1,ep)
392 k12(1,1,ep)=k12(1,1,ep)+dmf(1,1,ep)*bxibxj(1,2,ep)+
393 1 dmf(1,3,ep)*(bxibyj(1,2,ep)+byibxj(1,2,ep))+
394 2 dmf(3,3,ep)*byibyj(1,2,ep)
395 k12(1,2,ep)=k12(1,2,ep)+dmf(1,2,ep)*bxibyj(1,2,ep)+
396 1 dmf(1,3,ep)*bxibxj(1,2,ep)+
397 2 dmf(2,3,ep)*byibyj(1,2,ep)+
398 3 dmf(3,3,ep)*byibxj(1,2,ep)
399 k12(2,1,ep)=k12(2,1,ep)+dmf(1,2,ep)*byibxj(1,2,ep)+
400 1 dmf(1,3,ep)*bxibxj(1,2,ep)+
401 2 dmf(2,3,ep)*byibyj(1,2,ep)+
402 3 dmf(3,3,ep)*bxibyj(1,2,ep)
403 k12(2,2,ep)=k12(2,2,ep)+dmf(2,2,ep)*byibyj(1,2,ep)+
404 1 dmf(2,3,ep)*(bxibyj(1,2,ep)+byibxj(1,2,ep))+
405 2 dmf(3,3,ep)*bxibxj(1,2,ep)
407 k22(1,1,ep)=k22(1,1,ep)+dmf(1,1,ep)*bxibxj(2,2,ep)+
408 1 dmf(1,3,ep)*(bxibyj(2,2,ep)+byibxj(2,2,ep))+
409 2 dmf(3,3,ep)*byibyj(2,2,ep)
410 k22(1,2,ep)=k22(1,2,ep)+dmf(1,2,ep)*bxibyj(2,2,ep)+
411 1 dmf(1,3,ep)*bxibxj(2,2,ep)+
412 2 dmf(2,3,ep)*byibyj(2,2,ep)+
413 3 dmf(3,3,ep)*byibxj(2,2,ep)
414 k22(2,2,ep)=k22(2,2,ep)+dmf(2,2,ep)*byibyj(2,2,ep)+
415 1 dmf(2,3,ep)*(bxibyj(2,2,ep)+byibxj(2,2,ep))+
416 2 dmf(3,3,ep)*bxibxj(2,2,ep)
425 k13(i,j,ep)=-k11(i,j,ep)
426 k14(i,j,ep)=-k12(i,j,ep)
427 k23(i,j,ep)=-k12(j,i,ep)
428 k24(i,j,ep)=-k22(i,j,ep)
429 k33(i,j,ep)=k11(i,j,ep)
430 k34(i,j,ep)=k12(i,j,ep)
431 k44(i,j,ep)=k22(i,j,ep)
432 m13(i,j,ep)=-m11(i,j,ep)
433 m14(i,j,ep)=-m12(i,j,ep)
434 m23(i,j,ep)=-m12(j,i,ep)
435 m24(i,j,ep)=-m22(i,j,ep)
436 m33(i,j,ep)=m11(i,j,ep)
437 m34(i,j,ep)=m12(i,j,ep)
438 m44(i,j,ep)=m22(i,j,ep)
444 k13(2,1,ep)=-k11(1,2,ep)
445 k14(2,1,ep)=-k12(2,1,ep)
446 k23(2,1,ep)=-k12(1,2,ep)
447 k24(2,1,ep)=-k22(1,2,ep)
448 k34(2,1,ep)=k12(2,1,ep)
449 m13(2,1,ep)=-m11(1,2,ep)
450 m14(2,1,ep)=-m12(2,1,ep)
451 m23(2,1,ep)=-m12(1,2,ep)
452 m24(2,1,ep)=-m22(1,2,ep)
453 m34(2,1,ep)=m12(2,1,ep)
461 dg(ep)=gf(ep)*facz(ep)
462 dfz =df(1,1,ep)*facz(ep)
463 mf11(1,2,ep)=dfz*px1py2(ep)-dg(ep)*px2py1(ep)
464 mf12(1,2,ep)=(dfz-dg(ep))*px2py2(ep)
465 mf22(1,2,ep)=dfz*px2py1(ep)-dg(ep)*px1py2(ep)
466 mf23(1,2,ep)=-(dfz-dg(ep))*px1py1(ep)
469 dfz =df(1,2,ep)*facz(ep)
470 mf11(1,1,ep)=-dfz*py1py2(ep)+dg(ep)*px1px2(ep)
471 mf12(1,1,ep)=-dfz*py2py2(ep)+dg(ep)*px2px2(ep)
472 mf22(1,1,ep)=mf11(1,1,ep)
473 mf23(1,1,ep)=dfz*py1py1(ep)-dg(ep)*px1px1(ep)
474 mf11(2,2,ep)=-dfz*px1px2(ep)+dg(ep)*py1py2(ep)
475 mf12(2,2,ep)=-dfz*px2px2(ep)+dg(ep)*py2py2(ep)
476 mf22(2,2,ep)= mf11(2,2,ep)
477 mf23(2,2,ep)=dfz*px1px1(ep)-dg(ep)*py1py1(ep)
480 dfz =df(2,2,ep)*facz(ep)
481 mf11(2,1,ep)=dfz*px2py1(ep)-dg(ep)*px1py2(ep)
482 mf12(2,1,ep)=dfz*px2py2(ep)-dg(ep)*px2py2(ep)
483 mf22(2,1,ep)=dfz*px1py2(ep)-dg(ep)*px2py1(ep)
484 mf23(2,1,ep)=-(dfz-dg(ep))*px1py1(ep)
488 df5(m)=facz(m)*df5(m)
489 df6(m)=facz(m)*df6(m)
493 dfz=df5(ep)*px1py2(ep)-df6(ep)*px2py1(ep)
494 mf11(1,1,ep)=mf11(1,1,ep)-dfz
495 mf11(1,2,ep)=mf11(1,2,ep)+df5(ep)*(py1py2(ep)-px1px2(ep))
496 mf11(2,1,ep)=mf11(2,1,ep)+df6(ep)*(px1px2(ep)-py1py2(ep))
497 mf11(2,2,ep)=mf11(2,2,ep)+dfz
500 dfz=(df5(ep)-df6(ep))*px2py2(ep)
501 mf12(1,1,ep)=mf12(1,1,ep)-dfz
502 mf12(1,2,ep)=mf12(1,2,ep)+df5(ep)*(py2py2(ep)-px2px2(ep))
503 mf12(2,1,ep)=mf12(2,1,ep)+df6(ep)*(px2px2(ep)-py2py2(ep))
504 mf12(2,2,ep)=mf12(2,2,ep)+dfz
507 dfz=df5(ep)*px2py1(ep)-df6(ep)*px1py2(ep)
508 mf22(1,1,ep)=mf22(1,1,ep)-dfz
509 mf22(1,2,ep)=mf22(1,2,ep)+df5(ep)*(py1py2(ep)-px1px2(ep))
510 mf22(2,1,ep)=mf22(2,1,ep)+df6(ep)*(px1px2(ep)-py1py2(ep))
511 mf22(2,2,ep)=mf22(2,2,ep)+dfz
514 dfz=(df6(ep)-df5(ep))*px1py1(ep)
515 mf23(1,1,ep)=mf23(1,1,ep)-dfz
516 mf23(1,2,ep)=mf23(1,2,ep)-df5(ep)*(py1py1(ep)-px1px1(ep))
517 mf23(2,1,ep)=mf23(2,1,ep)-df6(ep)*(px1px1(ep)-py1py1(ep))
518 mf23(2,2,ep)=mf23(2,2,ep)+dfz
522 mf11(1,1,ep)=mf11(1,1,ep)
523 1 -dmf(1,2,ep)*px1py1(ep)-dmf(1,3,ep)*px1px1(ep)
524 2 -dmf(2,3,ep)*py1py1(ep)-dmf(3,3,ep)*px1py1(ep)
525 mf11(1,2,ep)=mf11(1,2,ep)
526 1 +dmf(1,1,ep)*px1px1(ep)+dmf(1,3,ep)*px1py1(ep)
527 2 +dmf(1,3,ep)*px1py1(ep)+dmf(3,3,ep)*py1py1(ep)
528 mf11(2,1,ep)=mf11(2,1,ep)
529 1 -dmf(2,2,ep)*py1py1(ep)-dmf(2,3,ep)*px1py1(ep)
530 2 -dmf(2,3,ep)*px1py1(ep)-dmf(3,3,ep)*px1px1(ep)
531 mf11(2,2,ep)=mf11(2,2,ep)
532 1 +dmf(1,2,ep)*px1py1(ep)+dmf(2,3,ep)*py1py1(ep)
533 2 +dmf(1,3,ep)*px1px1(ep)+dmf(3,3,ep)*px1py1(ep)
536 mf12(1,1,ep)=mf12(1,1,ep)
537 1 -dmf(1,2,ep)*px1py2(ep)-dmf(1,3,ep)*px1px2(ep)
538 2 -dmf(2,3,ep)*py1py2(ep)-dmf(3,3,ep)*px2py1(ep)
539 mf12(1,2,ep)=mf12(1,2,ep)
540 1 +dmf(1,1,ep)*px1px2(ep)+dmf(1,3,ep)*px1py2(ep)
541 2 +dmf(1,3,ep)*px2py1(ep)+dmf(3,3,ep)*py1py2(ep)
542 mf12(2,1,ep)=mf12(2,1,ep)
543 1 -dmf(2,2,ep)*py1py2(ep)-dmf(2,3,ep)*px2py1(ep)
544 2 -dmf(2,3,ep)*px1py2(ep)-dmf(3,3,ep)*px1px2(ep)
545 mf12(2,2,ep)=mf12(2,2,ep)
546 1 +dmf(1,2,ep)*px2py1(ep)+dmf(2,3,ep)*py1py2(ep)
547 2 +dmf(1,3,ep)*px1px2(ep)+dmf(3,3,ep)*px1py2(ep)
550 mf22(1,1,ep)=mf22(1,1,ep)
551 1 -dmf(1,2,ep)*px2py2(ep)-dmf(1,3,ep)*px2px2(ep)
552 2 -dmf(2,3,ep)*py2py2(ep)-dmf(3,3,ep)*px2py2(ep)
553 mf22(1,2,ep)=mf22(1,2,ep)
554 1 +dmf(1,1,ep)*px2px2(ep)+dmf(1,3,ep)*px2py2(ep)
555 2 +dmf(1,3,ep)*px2py2(ep)+dmf(3,3,ep)*py2py2(ep)
556 mf22(2,1,ep)=mf22(2,1,ep)
557 1 -dmf(2,2,ep)*py2py2(ep)-dmf(2,3,ep)*px2py2(ep)
558 2 -dmf(2,3,ep)*px2py2(ep)-dmf(3,3,ep)*px2px2(ep)
559 mf22(2,2,ep)=mf22(2,2,ep)
560 1 +dmf(1,2,ep)*px2py2(ep)+dmf(2,3,ep)*py2py2(ep)
561 2 +dmf(1,3,ep)*px2px2(ep)+dmf(3,3,ep)*px2py2(ep)
564 mf23(1,1,ep)=mf23(1,1,ep)
565 1 +dmf(1,2,ep)*px2py1(ep)+dmf(1,3,ep)*px1px2(ep)
566 2 +dmf(2,3,ep)*py1py2(ep)+dmf(3,3,ep)*px1py2(ep)
567 mf23(1,2,ep)=mf23(1,2,ep)
568 1 -dmf(1,1,ep)*px1px2(ep)-dmf(1,3,ep)*px2py1(ep)
569 2 -dmf(1,3,ep)*px1py2(ep)-dmf(3,3,ep)*py1py2(ep)
570 mf23(2,1,ep)=mf23(2,1,ep)
571 1 +dmf(2,2,ep)*py1py2
572 2 +dmf(2,3,ep)*px2py1(ep)+dmf(3,3,ep)*px1px2(ep)
573 mf23(2,2,ep)=mf23(2,2,ep)
574 1 -dmf(1,2,ep)*px1py2(ep)-dmf(2,3,ep)*py1py2(ep)
575 2 -dmf(1,3,ep)*px1px2(ep)-dmf(3,3,ep)*px2py1(ep)
582 mf13(i,j,ep)=-mf11(i,j,ep)
583 mf14(i,j,ep)=-mf12(i,j,ep)
584 mf24(i,j,ep)=-mf22(i,j,ep)
585 mf33(i,j,ep)=mf11(i,j,ep)
586 mf34(i,j,ep)=mf12(i,j,ep)
587 mf44(i,j,ep)=mf22(i,j,ep)
595 fm12(i,j,ep)=-mf23(j,i,ep)
596 fm13(i,j,ep)=-mf11(j,i,ep)
597 fm14(i,j,ep)= mf23(j,i,ep)
598 fm23(i,j,ep)=-mf12(j,i,ep)
599 fm24(i,j,ep)=-mf22(j,i,ep)
600 fm34(i,j,ep)=-mf23(j,i,ep)
606 m11(3,3,ep)=dhz(ep)*(px1px1(ep)+py1py1(ep))
607 m12(3,3,ep)=dhz(ep)*(px1px2(ep)+py1py2(ep))
608 m13(3,3,ep)=-m11(3,3,ep)
609 m14(3,3,ep)=-m12(3,3,ep)
610 m22(3,3,ep)=dhz(ep)*(px2px2(ep)+py2py2(ep))
611 m23(3,3,ep)=-m12(3,3,ep)
612 m24(3,3,ep)=-m22(3,3,ep)
613 m33(3,3,ep)=m11(3,3,ep)
615 m44(3,3,ep)=m22(3,3,ep)
622 m11(3,3,ep)=m11(3,3,ep)+c1
623 m22(3,3,ep)=m22(3,3,ep)+c2
624 m33(3,3,ep)=m11(3,3,ep)
625 m44(3,3,ep)=m22(3,3,ep)
637 2 HM ,HF ,HZ ,A_I ,Z1 ,
638 3 PX1 ,PX2 ,PY1 ,PY2 ,
639 6 K11,K12,K13,K14,K22,K23,K24,K33,K34,K44,
640 7 M11,M12,M13,M14,M22,M23,M24,M33,M34,M44,
641 8 MF11,MF12,MF13,MF14,MF22,MF23,MF24,MF33,
642 9 MF34,MF44,FM12,FM13,FM14,FM23,FM24,FM34,
643 A IORTH,HMOR,HFOR ,IPLAT,DHZ ,
644 4 PRX ,PRY ,PRXY ,PRZ ,HMFOR,NPLAT)
649#include "implicit_f.inc"
650#include "mvsiz_p.inc"
655 INTEGER JFT,JLT,IPLAT(*),IORTH,NPLAT
657 . VOL(*),THK0(*),THK2(*),A_I(*),Z1(*),
658 . HM(MVSIZ,4),HF(MVSIZ,4),HZ(*),
659 . PX1(*) ,PX2(*) ,PY1(*) ,PY2(*) ,
660 . K11(3,3,*),K12(3,3,*),K13(3,3,*),K14(3,3,*),
661 . K22(3,3,*),K23(3,3,*),K24(3,3,*),K33(3,3,*),
662 . M11(3,3,*),M12(3,3,*),M13(3,3,*),M14(3,3,*),
666(3,3,*),FM13(3,3,*),FM14(3,3,*),
667 . FM23(3,3,*),FM24(3,3,*),FM34(3,3,*),
668 . (3,3,*),K44(3,3,*),M34(3,3,*),M44(3,3,*),
669 . MF34(3,3,*),MF44(3,3,*),DHZ(*),HMOR(MVSIZ,2),HFOR(MVSIZ,2)
671 . prx(4,*),pry(4,*),prxy(4,*),prz(4,*) ,hmfor(mvsiz,6)
680 . DM(3,3,MVSIZ),DPRZ(4,),
681 . C1,C2,(4,MVSIZ),CBRY(4,MVSIZ),CBRZ(4,),
682 . KZ11(2,2,MVSIZ),KZ12(2,2,MVSIZ),KZ22(2,2,MVSIZ),
683 . DMF(3,3,MVSIZ),PCX1(MVSIZ),PCX2(MVSIZ),PCY1(MVSIZ),
684 . PCY2(MVSIZ),FACZ(MVSIZ)
686#include "vectorize.inc"
691 dm(1,1,m)=hm(ep,1)*c2
692 dm(2,2,m)=hm(ep,2)*c2
693 dm(1,2,m)=hm(ep,3)*c2
694 dm(3,3,m)=hm(ep,4)*c2
695 dhz(m)=hz(ep)*fourth*c2
703 dm(1,3,m)=hmor(ep,1)*c2
704 dm(2,3,m)=hmor(ep,2)*c2
716 kz11(1,1,ep)=py1(ep)*dhz(ep)*py1(ep)
717 kz11(1,2,ep)=-py1(ep)*dhz(ep)*px1(ep)
718 kz11(2,2,ep)=px1(ep)*dhz(ep)*px1(ep)
719 kz12(1,1,ep)=py1(ep)*dhz(ep)*py2(ep)
720 kz12(1,2,ep)=-py1(ep)*dhz(ep)*px2(ep)
721 kz12(2,1,ep)=-px1(ep)*dhz(ep)*py2(ep)
722 kz12(2,2,ep)=px1(ep)*dhz(ep)*px2(ep)
723 kz22(1,1,ep)=py2(ep)*dhz(ep)*py2(ep)
724 kz22(1,2,ep)=-py2(ep)*dhz(ep)*px2(ep)
725 kz22(2,2,ep)=px2(ep)*dhz(ep)*px2(ep)
731 k11(i,j,ep)=k11(i,j,ep)+kz11(i,j,ep)
732 k12(i,j,ep)=k12(i,j,ep)+kz12(i,j,ep)
733 k22(i,j,ep)=k22(i,j,ep)+kz22(i,j,ep)
734 k13(i,j,ep)=k13(i,j,ep)-kz11(i,j,ep)
736 k23(i,j,ep)=k23(i,j,ep)-kz12(j,i,ep)
737 k24(i,j,ep)=k24(i,j,ep)-kz22(i,j,ep)
738 k33(i,j,ep)=k33(i,j,ep)+kz11(i,j,ep)
739 k34(i,j,ep)=k34(i,j,ep)+kz12(i,j,ep)
740 k44(i,j,ep)=k44(i,j,ep)+kz22(i,j,ep)
746 k12(2,1,ep)=k12(2,1,ep)+kz12(2,1,ep)
747 k13(2,1,ep)=k13(2,1,ep)-kz11(1,2,ep)
748 k14(2,1,ep)=k14(2,1,ep)-kz12(2,1,ep)
749 k23(2,1,ep)=k23(2,1,ep)-kz12(1,2,ep)
750 k24(2,1,ep)=k24(2,1,ep)-kz22(1,2,ep)
751 k34(2,1,ep)=k34(2,1,ep)+kz12(2,1,ep)
756 dprz(j,m)=prz(j,m)*dhz(m)
761 c1 = -py2(m)*dprz(1,m)
762 c2 = px2(m)*dprz(1,m)
763 fm12(3,1,m)=fm12(3,1,m)+ c1
764 fm14(3,1,m)=fm14(3,1,m)- c1
765 fm12(3,2,m)=fm12(3,2,m)+ c2
766 fm14(3,2,m)=fm14(3,2,m)- c2
767 c1 = -py1(m)*dprz(1,m)
768 c2 = px1(m)*dprz(1,m)
769 mf11(1,3,m)=mf11(1,3,m)+ c1
770 fm13(3,1,m)=fm13(3,1,m)- c1
771 mf11(2,3,m)=mf11(2,3,m)+ c2
772 fm13(3,2,m)=fm13(3,2,m)- c2
773 c1 = -py1(m)*dprz(2,m)
774 c2 = px1(m)*dprz(2,m)
775 mf12(1,3,m)=mf12(1,3,m)+ c1
776 fm23(3,1,m)=fm23(3,1,m)- c1
777 mf12(2,3,m)=mf12(2,3,m)+ c2
778 fm23(3,2,m)=fm23(3,2,m)- c2
779 c1 = -py1(m)*dprz(3,m)
780 c2 = px1(m)*dprz(3,m)
781 mf13(1,3,m)=mf13(1,3,m)+ c1
782 mf33(1,3,m)=mf33(1,3,m)- c1
783 mf13(2,3,m)=mf13(2,3,m)+ c2
784 mf33(2,3,m)=mf33(2,3,m)- c2
785 c1 = -py1(m)*dprz(4,m)
786 c2 = px1(m)*dprz(4,m)
787 mf14(1,3,m)=mf14(1,3,m)+ c1
788 mf34(1,3,m)=mf34(1,3,m)- c1
789 mf14(2,3,m)=mf14(2,3,m)+ c2
790 mf34(2,3,m)=mf34(2,3,m)- c2
791 c1 = -py2(m)*dprz(2,m)
792 c2 = px2(m)*dprz(2,m)
793 mf22(1,3,m)=mf22(1,3,m)+ c1
794 fm24(3,1,m)=fm24(3,1,m)- c1
795 mf22(2,3,m)=mf22(2,3,m)+ c2
796 fm24(3,2,m)=fm24(3,2,m)- c2
797 c1 = -py2(m)*dprz(3,m)
798 c2 = px2(m)*dprz(3,m)
799 mf23(1,3,m)=mf23(1,3,m)+ c1
800 fm34(3,1,m)=fm34(3,1,m)- c1
801 mf23(2,3,m)=mf23(2,3,m)+ c2
802 fm34(3,2,m)=fm34(3,2,m)- c2
803 c1 = -py2(m)*dprz(4,m)
804 c2 = px2(m)*dprz(4,m)
805 mf24(1,3,m)=mf24(1,3,m)+ c1
807 mf24(2,3,m)=mf24(2,3,m)+ c2
808 mf44(2,3,m)=mf44(2,3,m)- c2
812 m11(3,3,ep)=prz(1,ep)*dprz(1,ep)
813 m12(3,3,ep)=prz(1,ep)*dprz(2,ep)
814 m13(3,3,ep)=prz(1,ep)*dprz(3,ep)
815 m14(3,3,ep)=prz(1,ep)*dprz(4,ep)
816 m22(3,3,ep)=prz(2,ep)*dprz(2,ep)
817 m23(3,3,ep)=prz(2,ep)*dprz(3,ep)
818 m24(3,3,ep)=prz(2,ep)*dprz(4,ep)
819 m33(3,3,ep)=prz(3,ep)*dprz(3,ep)
820 m34(3,3,ep)=prz(3,ep)*dprz(4,ep)
821 m44(3,3,ep)=prz(4,ep)*dprz(4,ep)
829 cbrx(j,m) =dm(1,1,m)*prx(j,m)+dm(1,2,m)*pry(j,m)+
830 . dm(1,3,m)*prxy(j,m)
831 cbry(j,m) =dm(2,1,m)*prx(j,m)+dm(2,2,m)*pry(j,m)+
832 . dm(2,3,m)*prxy(j,m)
833 cbrz(j,m) =dm(3,1,m)*prx(j,m)+dm(3,2,m)*pry(j,m)+
834 . dm(3,3,m)*prxy(j,m)
839 c1 = px2(m)*cbrx(1,m)+py2(m)*cbrz(1,m)
840 c2 = py2(m)*cbry(1,m)+px2(m)*cbrz(1,m)
841 fm12(3,1,m)=fm12(3,1,m)+ c1
842 fm14(3,1,m)=fm14(3,1,m)- c1
843 fm12(3,2,m)=fm12(3,2,m)+ c2
844 fm14(3,2,m)=fm14(3,2,m)- c2
846 c1 = px1(m)*cbrx(1,m)+py1(m)*cbrz(1,m)
847 c2 = py1(m)*cbry(1,m)+px1(m)*cbrz(1,m)
848 mf11(1,3,m)=mf11(1,3,m)+ c1
849 fm13(3,1,m)=fm13(3,1,m)- c1
850 mf11(2,3,m)=mf11(2,3,m)+ c2
851 fm13(3,2,m)=fm13(3,2,m)- c2
852 c1 = px1(m)*cbrx(2,m)+py1(m)*cbrz(2,m)
853 c2 = py1(m)*cbry(2,m)+px1(m)*cbrz(2,m)
854 mf12(1,3,m)=mf12(1,3,m)+ c1
855 fm23(3,1,m)=fm23(3,1,m)- c1
856 mf12(2,3,m)=mf12(2,3,m)+ c2
857 fm23(3,2,m)=fm23(3,2,m)- c2
858 c1 = px1(m)*cbrx(3,m)+py1(m)*cbrz(3,m)
859 c2 = py1(m)*cbry(3,m)+px1(m)*cbrz(3,m)
860 mf13(1,3,m)=mf13(1,3,m)+ c1
861 mf33(1,3,m)=mf33(1,3,m)- c1
862 mf13(2,3,m)=mf13(2,3,m)+ c2
863 mf33(2,3,m)=mf33(2,3,m)- c2
864 c1 = px1(m)*cbrx(4,m)+py1(m)*cbrz(4,m)
865 c2 = py1(m)*cbry(4,m)+px1(m)*cbrz(4,m)
866 mf14(1,3,m)=mf14(1,3,m)+ c1
867 mf34(1,3,m)=mf34(1,3,m)- c1
868 mf14(2,3,m)=mf14(2,3,m)+ c2
869 mf34(2,3,m)=mf34(2,3,m)- c2
870 c1 = px2(m)*cbrx(2,m)+py2(m)*cbrz(2,m)
871 c2 = py2(m)*cbry(2,m)+px2(m)*cbrz(2,m)
872 mf22(1,3,m)=mf22(1,3,m)+ c1
873 fm24(3,1,m)=fm24(3,1,m)- c1
874 mf22(2,3,m)=mf22(2,3,m)+ c2
875 fm24(3,2,m)=fm24(3,2,m)- c2
876 c1 = px2(m)*cbrx(3,m)+py2(m)*cbrz(3,m)
877 c2 = py2(m)*cbry(3,m)+px2(m)*cbrz(3,m)
878 mf23(1,3,m)=mf23(1,3,m)+ c1
879 fm34(3,1,m)=fm34(3,1,m)- c1
880 mf23(2,3,m)=mf23(2,3,m)+ c2
881 fm34(3,2,m)=fm34(3,2,m)- c2
882 c1 = px2(m)*cbrx(4,m)+py2(m)*cbrz(4,m)
883 c2 = py2(m)*cbry(4,m)+px2(m)*cbrz(4,m)
884 mf24(1,3,m)=mf24(1,3,m)+ c1
885 mf44(1,3,m)=mf44(1,3,m)- c1
886 mf24(2,3,m)=mf24(2,3,m)+ c2
887 mf44(2,3,m)=mf44(2,3,m)- c2
891 m11(3,3,m)=m11(3,3,m)+prx(1,m)*cbrx(1,m)+
892 . pry(1,m)*cbry(1,m)+prxy(1,m)*cbrz(1,m)
893 m12(3,3,m)=m12(3,3,m)+prx(1,m)*cbrx(2,m)+
894 . pry(1,m)*cbry(2,m)+prxy(1,m)*cbrz(2,m)
895 m13(3,3,m)=m13(3,3,m)+prx(1,m)*cbrx(3,m)+
896 . pry(1,m)*cbry(3,m)+prxy(1,m)*cbrz(3,m)
897 m14(3,3,m)=m14(3,3,m)+prx(1,m)*cbrx(4,m)+
898 . pry(1,m)*cbry(4,m)+prxy(1,m)*cbrz(4,m)
899 m22(3,3,m)=m22(3,3,m)+prx(2,m)*cbrx(2,m)+
900 . pry(2,m)*cbry(2,m)+prxy(2,m)*cbrz(2,m)
901 m23(3,3,m)=m23(3,3,m)+prx(2,m)*cbrx(3,m)+
902 . pry(2,m)*cbry(3,m)+prxy(2,m)*cbrz(3,m)
903 m24(3,3,m)=m24(3,3,m)+prx(2,m)*cbrx(4,m)+
904 . pry(2,m)*cbry(4,m)+prxy(2,m)*cbrz(4,m)
905 m33(3,3,m)=m33(3,3,m)+prx(3,m)*cbrx(3,m)+
906 . pry(3,m)*cbry(3,m)+prxy(3,m)*cbrz(3,m)
907 m34(3,3,m)=m34(3,3,m)+prx(3,m)*cbrx(4,m)+
908 . pry(3,m)*cbry(4,m)+prxy(3,m)*cbrz(4,m)
909 m44(3,3,m)=m44(3,3,m)+prx(4,m)*cbrx(4,m)+
910 . pry(4,m)*cbry(4,m)+prxy(4,m)*cbrz(4,m)
915#include "vectorize.inc"
919 dmf(1,1,m)=hmfor(ep,1)*c2
920 dmf(2,2,m)=hmfor(ep,2)*c2
921 dmf(1,2,m)=hmfor(ep,3)*c2
922 dmf(1,3,m)=hmfor(ep,5)*c2
923 dmf(2,3,m)=hmfor(ep,6)*c2
924 dmf(2,1,m)=dmf(1,2,m)
925 dmf(3,1,m)=dmf(1,3,m)
926 dmf(3,2,m)=dmf(2,3,m)
927 dmf(3,3,m)=hmfor(ep,4)*c2
932 cbrx(j,ep) =dmf(1,1,ep)*prx(j,ep)+dmf(1,2,ep)*pry(j,ep)+
933 . dmf(1,3,ep)*prxy(j,ep)
934 cbry(j,ep) =dmf(2,1,ep)*prx(j,ep)+dmf(2,2,ep)*pry(j,ep)+
935 . dmf(2,3,ep)*prxy(j,ep)
936 cbrz(j,ep) =dmf(3,1,ep)*prx(j,ep)+dmf(3,2,ep)*pry(j,ep)+
937 . dmf(3,3,ep)*prxy(j,ep)
943 facz(m) = four*z1(ep)*a_i(ep)
947 pcx1(ep)=facz(ep)*py2(ep)
948 pcx2(ep)=facz(ep)*py1(ep)
949 pcy1(ep)=-facz(ep)*px2(ep)
950 pcy2(ep)=-facz(ep)*px1(ep)
959 c1=-py1(ep)*cbry(1,ep)-px1(ep)*cbrz(1,ep)
960 c2= px1(ep)*cbrx(1,ep)+py1(ep)*cbrz(1,ep)
961 m11(1,3,ep)=m11(1,3,ep)+c1
962 m11(2,3,ep)=m11(2,3,ep)+c2
963 m13(3,1,ep)=m13(3,1,ep)-c1
964 m13(3,2,ep)=m13(3,2,ep)-c2
965 c1=-py1(ep)*cbry(2,ep)-px1(ep)*cbrz(2,ep)
966 c2= px1(ep)*cbrx(2,ep)+py1(ep)*cbrz(2,ep)
967 m12(1,3,ep)=m12(1,3,ep)+c1
968 m12(2,3,ep)=m12(2,3,ep)+c2
969 m23(3,1,ep)=m23(3,1,ep)-c1
970 m23(3,2,ep)=m23(3,2,ep)-c2
971 c1=-py1(ep)*cbry(4,ep)-px1(ep)*cbrz(4,ep)
972 c2= px1(ep)*cbrx(4,ep)+py1(ep)*cbrz(4,ep)
973 m14(1,3,ep)=m14(1,3,ep)+c1
974 m14(2,3,ep)=m14(2,3,ep)+c2
975 m34(1,3,ep)=m34(1,3,ep)-c1
976 m34(2,3,ep)=m34(2,3,ep)-c2
977 c1=-py2(ep)*cbry(2,ep)-px2(ep)*cbrz(2,ep)
978 c2= px2(ep)*cbrx(2,ep)+py2(ep)*cbrz(2,ep)
979 m22(1,3,ep)=m22(1,3,ep)+c1
980 m22(2,3,ep)=m22(2,3,ep)+c2
981 m24(3,1,ep)=m24(3,1,ep)-c1
982 m24(3,2,ep)=m24(3,2,ep)-c2
983 c1=-py2(ep)*cbry(3,ep)-px2(ep)*cbrz(3,ep)
985 m23(1,3,ep)=m23(1,3,ep)+c1
986 m23(2,3,ep)=m23(2,3,ep)+c2
987 m34(3,1,ep)=m34(3,1,ep)-c1
988 m34(3,2,ep)=m34(3,2,ep)-c2
989 c1= py1(ep)*cbry(3,ep)+px1(ep)*cbrz(3,ep)
990 c2=-px1(ep)*cbrx(3,ep)-py1(ep)*cbrz(3,ep)
991 m33(1,3,ep)=m33(1,3,ep)+c1
992 m33(2,3,ep)=m33(2,3,ep)+c2
993 m13(1,3,ep)=m13(1,3,ep)-c1
994 m13(2,3,ep)=m13(2,3,ep)-c2
995 c1= py2(ep)*cbry(4,ep)+px2(ep)*cbrz(4,ep)
996 c2=-px2(ep)*cbrx(4,ep)-py2(ep)*cbrz(4,ep)
997 m44(1,3,ep)=m44(1,3,ep)+c1
998 m44(2,3,ep)=m44(2,3,ep)+c2
999 m24(1,3,ep)=m24(1,3,ep)-c1
1000 m24(2,3,ep)=m24(2,3,ep)-c2
1001 c1=-py2(ep)*cbry(1,ep)-px2(ep)*cbrz(1,ep)
1002 c2= px2(ep)*cbrx(1,ep)+py2(ep)*cbrz(1,ep)
1003 m12(3,1,ep)=m12(3,1,ep)+c1
1004 m12(3,2,ep)=m12(3,2,ep)+c2
1005 m14(3,1,ep)=m14(3,1,ep)-c1
1006 m14(3,2,ep)=m14(3,2,ep)-c2
1015 c1= pcx1(ep)*cbrx(1,ep)+pcy1(ep)*cbrz(1,ep)
1016 c2= pcy1(ep)*cbry(1,ep)+pcx1(ep)*cbrz(1,ep)
1017 mf11(1,3,ep)=mf11(1,3,ep)+c1
1018 mf11(2,3,ep)=mf11(2,3,ep)+c2
1019 fm13(3,1,ep)=fm13(3,1,ep)-c1
1020 fm13(3,2,ep)=fm13(3,2,ep)-c2
1021 c1= pcx1(ep)*cbrx(2,ep)+pcy1(ep)*cbrz(2,ep)
1022 c2= pcy1(ep)*cbry(2,ep)+pcx1(ep)*cbrz(2,ep)
1023 mf12(1,3,ep)=mf12(1,3,ep)+c1
1024 mf12(2,3,ep)=mf12(2,3,ep)+c2
1025 fm23(3,1,ep)=fm23(3,1,ep)-c1
1026 fm23(3,2,ep)=fm23(3,2,ep)-c2
1027 c1= pcx1(ep)*cbrx(4,ep)+pcy1(ep)*cbrz(4,ep)
1028 c2= pcy1(ep)*cbry(4,ep)+pcx1(ep)*cbrz(4,ep)
1029 mf14(1,3,ep)=mf14(1,3,ep)+c1
1030 mf14(2,3,ep)=mf14(2,3,ep)+c2
1031 mf34(1,3,ep)=mf34(1,3,ep)-c1
1032 mf34(2,3,ep)=mf34(2,3,ep)-c2
1033 c1= pcx2(ep)*cbrx(2,ep)+pcy2(ep)*cbrz(2,ep)
1034 c2= pcy2(ep)*cbry(2,ep)+pcx2(ep)*cbrz(2,ep)
1035 mf22(1,3,ep)=mf22(1,3,ep)+c1
1036 mf22(2,3,ep)=mf22(2,3,ep)+c2
1037 fm24(3,1,ep)=fm24(3,1,ep)-c1
1039 c1= pcx2(ep)*cbrx(3,ep)+pcy2(ep)*cbrz
1040 c2= pcy2(ep)*cbry(3,ep)+pcx2(ep)*cbrz(3,ep)
1041 mf23(1,3,ep)=mf23(1,3,ep)+c1
1042 mf23(2,3,ep)=mf23(2,3,ep)+c2
1043 fm34(3,1,ep)=fm34(3,1,ep)-c1
1044 fm34(3,2,ep)=fm34(3,2,ep)-c2
1045 c1=-pcx1(ep)*cbrx(3,ep)-pcy1(ep)*cbrz(3,ep)
1046 c2=-pcy1(ep)*cbry(3,ep)-pcx1(ep)*cbrz(3,ep)
1047 mf33(1,3,ep)=mf33(1,3,ep)+c1
1048 mf33(2,3,ep)=mf33(2,3,ep)+c2
1049 mf13(1,3,ep)=mf13(1,3,ep)-c1
1050 mf13(2,3,ep)=mf13(2,3,ep)-c2
1051 c1=-pcx2(ep)*cbrx(4,ep)-pcy2(ep)*cbrz(4,ep)
1052 c2=-pcy2(ep)*cbry(4,ep)-pcx2(ep)*cbrz(4,ep)
1053 mf44(1,3,ep)=mf44(1,3,ep)+c1
1054 mf44(2,3,ep)=mf44(2,3,ep)+c2
1055 mf24(1,3,ep)=mf24(1,3,ep)-c1
1056 mf24(2,3,ep)=mf24(2,3,ep)-c2
1057 c1= pcx2(ep)*cbrx(1,ep)+pcy2(ep)*cbrz(1,ep)
1058 c2= pcy2(ep)*cbry(1,ep)+pcx2(ep)*cbrz(1,ep)
1059 fm12(3,1,ep)=fm12(3,1,ep)+c1
1060 fm12(3,2,ep)=fm12(3,2,ep)+c2
1061 fm14(3,1,ep)=fm14(3,1,ep)-c1
1062 fm14(3,2,ep)=fm14(3,2,ep)-c2