28 SUBROUTINE czlken3(JFT ,JLT ,VOL ,THK0 ,THK2 ,
29 2 HM ,HZ ,A_I ,PX1 ,PX2 ,
30 3 PY1 ,PY2 ,HXX ,HYY ,HXY ,
31 4 PH1 ,PH2 ,Z1 ,NPLAT,IPLAT,DHZ ,
32 5 K11,K12,K13,K14,K22,K23,K24,K33,K34,K44,
33 6 M11,M12,M13,M14,M22,M23,M24,M33,M34,M44,
34 7 MF11,MF12,MF13,MF14,MF22,MF23,MF24,MF33,
35 8 MF34,MF44,FM12,FM13,FM14,FM23,FM24,FM34,
41#include
"implicit_f.inc"
47 INTEGER JFT,JLT,NPLAT,IPLAT(*),IDRIL
49 . VOL(*),THK0(*),THK2(*),A_I(*),Z1(*),HM(MVSIZ,4),HZ(*),
50 . PX1(*) ,PX2(*) ,PY1(*) ,PY2(*) ,PH1(*) ,PH2(*) ,
51 . HXX(*),HYY(*),HXY(*),
52 . K11(3,3,*),K12(3,3,*),K13(3,3,*),K14(3,3,*),
53 . K22(3,3,*),K23(3,3,*),(3,3,*),K33(3,3,*),
54 . M11(3,3,*),M12(3,3,*),M13(3,3,*),M14(3,3,*),
55 . m22(3,3,*),m23(3,3,*),m24(3,3,*),m33(3,3,*),
56 . mf11(3,3,*),mf12(3,3,*),mf13(3,3,*),mf14(3,3,*),
57 . mf22(3,3,*),mf23(3,3,*),mf24(3,3,*),mf33(3,3,*),
58 . fm12(3,3,*),fm13(3,3,*),fm14(3,3,*),
59 . fm23(3,3,*),fm24(3,3,*),fm34(3,3,*),
60 . k34(3,3,*),k44(3,3,*),m34(3,3,*),m44(3,3,*),
61 . mf34(3,3,*),mf44(3,3,*),dhz(*)
70 . DM(2,2,MVSIZ),C1,C2,GM(MVSIZ),
71 . (MVSIZ),G12(MVSIZ),G13(MVSIZ),
72 . G14(MVSIZ),G22(MVSIZ),G23(MVSIZ),
73 . G24(),G33(MVSIZ),G34(MVSIZ),G44(MVSIZ),
74 . CHM(2,2,MVSIZ),CHF(2,2,MVSIZ),FACF(MVSIZ),FAC,ZR(MVSIZ),
75 . GAMA(4,MVSIZ),CXX,CYY,CXY,C11,,C12,CX1,CX2,CY1,CY2
79#include "vectorize.inc"
89 facf(m)=one_over_12*thk2(ep)
95 chm(1,1,ep)=dm(1,1,ep)*hxx(ep)
97 chm(2,2,ep)=dm(2,2,ep)*hyy(ep)
99 chm(1,2,ep)=dm(1,2,ep)*hxy(ep)
101 chm(2,1,ep)=chm(1,2,ep)
103 chf(1,1,ep)=facf(ep)*chm(2,2,ep)
105 chf(2,2,ep)=facf(ep)*chm(1,1,ep)
107 chf(1,2,ep)=-facf(ep)*chm(1,2,ep)
114 chf(2,1,ep)=chf(1,2,ep)
119 chm(1,1,ep)=chm(1,1,ep)+gm(ep)*hyy(ep)
120 chm(2,2,ep)=chm(2,2,ep)+gm(ep)*hxx(ep)
121 chm(1,2,ep)=chm(1,2,ep)+gm(ep)*hxy(ep)
122 chm(2,1,ep)=chm(1,2,ep)
124 fac = facf(ep)*gm(ep)
125 chf(1,1,ep)=chf(1,1,ep)+fac*hxx(ep)
126 chf(2,2,ep)=chf(2,2,ep)+fac*hyy(ep)
127 chf(1,2,ep)=chf(1,2,ep)-fac*hxy(ep)
128 chf(2,1,ep)=chf(1,2,ep)
133 gama(1,ep)=fourth-ph1(ep)
134 gama(2,ep)=-fourth-ph2(ep)
135 gama(3,ep)=fourth+ph1(ep)
136 gama(4,ep)=-fourth+ph2(ep)
137 g11(ep) =gama(1,ep)*gama(1,ep)
138 g12(ep) =gama(1,ep)*gama(2,ep)
139 g13(ep) =gama(1,ep)*gama(3,ep)
140 g14(ep) =gama(1,ep)*gama(4,ep)
141 g22(ep) =gama(2,ep)*gama(2,ep)
142 g23(ep) =gama(2,ep)*gama(3,ep)
143 g24(ep) =gama(2,ep)*gama(4,ep)
144 g33(ep) =gama(3,ep)*gama(3,ep)
145 g34(ep) =gama(3,ep)*gama(4,ep)
146 g44(ep) =gama(4,ep)*gama(4,ep)
152 k11(i,j,ep)=k11(i,j,ep)+chm(i,j,ep)*g11(ep)
153 k22(i,j,ep)=k22(i,j,ep)+chm(i,j,ep)*g22(ep)
154 k33(i,j,ep)=k33(i,j,ep)+chm(i,j,ep)*g33(ep)
155 k44(i,j,ep)=k44(i,j,ep)+chm(i,j,ep)*g44(ep)
156 m11(i,j,ep)=m11(i,j,ep)+chf(i,j,ep)*g11(ep)
157 m22(i,j,ep)=m22(i,j,ep)+chf(i,j,ep)*g22(ep)
158 m33(i,j,ep)=m33(i,j,ep)+chf(i,j,ep)*g33(ep)
159 m44(i,j,ep)=m44(i,j,ep)+chf(i,j,ep)*g44(ep)
166 k12(i,j,ep)=k12(i,j,ep)+chm(i,j,ep)*g12(ep)
167 k13(i,j,ep)=k13(i,j,ep)+chm(i,j,ep)*g13(ep)
168 k14(i,j,ep)=k14(i,j,ep)+chm(i,j,ep)*g14(ep)
169 k23(i,j,ep)=k23(i,j,ep)+chm(i,j,ep)*g23(ep)
170 k24(i,j,ep)=k24(i,j,ep)+chm(i,j,ep)*g24(ep)
171 k34(i,j,ep)=k34(i,j,ep)+chm(i,j,ep)*g34(ep)
172 m12(i,j,ep)=m12(i,j,ep)+chf(i,j,ep)*g12(ep)
173 m13(i,j,ep)=m13(i,j,ep)+chf(i,j,ep)*g13(ep)
174 m14(i,j,ep)=m14(i,j,ep)+chf(i,j,ep)*g14(ep)
175 m23(i,j,ep)=m23(i,j,ep)+chf(i,j,ep)*g23(ep)
176 m24(i,j,ep)=m24(i,j,ep)+chf(i,j,ep)*g24(ep)
177 m34(i,j,ep)=m34(i,j,ep)+chf(i,j,ep)*g34(ep)
183 cxx =zr(ep)*chm(1,1,ep)
184 cyy =zr(ep)*chm(2,2,ep)
185 cxy =zr(ep)*chm(1,2,ep)
186 cx1 = cxx*px1(ep)+cxy*py1(ep)
187 cx2 = cxx*px2(ep)+cxy*py2(ep)
188 cy1 = cxy*px1(ep)+cyy*py1(ep)
189 cy2 = cxy*px2(ep)+cyy*py2(ep)
190 k11(1,3,ep) = cx1*gama(1,ep)
191 k11(2,3,ep) = cy1*gama(1,ep)
192 k11(3,1,ep) = k11(1,3,ep)
193 k11(3,2,ep) = k11(2,3,ep)
194 k22(1,3,ep) = cx2*gama(2,ep)
195 k22(2,3,ep) = cy2*gama(2,ep)
196 k22(3,1,ep) = k22(1,3,ep)
197 k22(3,2,ep) = k22(2,3,ep)
198 k33(1,3,ep) = -cx1*gama(3,ep)
199 k33(2,3,ep) = -cy1*gama(3,ep)
200 k33(3,1,ep) = k33(1,3,ep)
201 k33(3,2,ep) = k33(2,3,ep)
202 k44(1,3,ep) = -cx2*gama(4,ep)
203 k44(2,3,ep) = -cy2*gama(4,ep)
204 k44(3,1,ep) = k44(1,3,ep)
205 k44(3,2,ep) = k44(2,3,ep)
206 k12(1,3,ep) = cx2*gama(1,ep)
207 k12(2,3,ep) = cy2*gama(1,ep)
208 k12(3,1,ep) = cx1*gama(2,ep)
209 k12(3,2,ep) = cy1*gama(2,ep)
210 k13(1,3,ep) = -k11(1,3,ep)
211 k13(2,3,ep) = -k11(2,3,ep)
212 k13(3,1,ep) = -k33(1,3,ep)
213 k13(3,2,ep) = -k33(2,3,ep)
214 k14(1,3,ep) = -k12(1,3,ep)
215 k14(2,3,ep) = -k12(2,3,ep)
216 k14(3,1,ep) = cx1*gama(4,ep)
217 k14(3,2,ep) = cy1*gama(4,ep)
218 k23(1,3,ep) = -k12(3,1,ep)
219 k23(2,3,ep) = -k12(3,2,ep)
220 k23(3,1,ep) = cx2*gama(3,ep)
221 k23(3,2,ep) = cy2*gama(3,ep)
222 k24(1,3,ep) = -k22(1,3,ep)
223 k24(2,3,ep) = -k22(2,3,ep)
224 k24(3,1,ep) = -k44(1,3,ep)
225 k24(3,2,ep) = -k44(2,3,ep)
226 k34(1,3,ep) = -k23(3,1,ep)
227 k34(2,3,ep) = -k23(3,2,ep)
228 k34(3,1,ep) = -k14(3,1,ep)
229 k34(3,2,ep) = -k14(3,2,ep)
236 c11 = cxx*px1(ep)*px1(ep)+cyy*py1(ep)*py1(ep)+
237 1 cxy*(px1(ep)*py1(ep)+py1(ep)*px1(ep))
238 c12 = cxx*px1(ep)*px2(ep)+cyy*py1(ep)*py2(ep)+
239 1 cxy*(px1(ep)*py2(ep)+py1(ep)*px2(ep))
240 c22 = cxx*px2(ep)*px2(ep)+cyy*py2(ep)*py2(ep)+
241 1 cxy*(px2(ep)*py2(ep)+py2(ep)*px2(ep))
242 k11(3,3,ep) = k11(3,3,ep)+c11
243 k12(3,3,ep) = k12(3,3,ep)+c12
244 k13(3,3,ep) = k13(3,3,ep)-c11
245 k14(3,3,ep) = k14(3,3,ep)-c12
246 k22(3,3,ep) = k22(3,3,ep)+c22
247 k23(3,3,ep) = k23(3,3,ep)-c12
248 k24(3,3,ep) = k24(3,3,ep)-c22
249 k33(3,3,ep) = k33(3,3,ep)+c11
250 k34(3,3,ep) = k34(3,3,ep)+c12
251 k44(3,3,ep) = k44(3,3,ep)+c22
256 fac =dhz(ep)*(hxx(ep)+hyy(ep))
257 m11(3,3,ep)=m11(3,3,ep)+fac*g11(ep)
258 m12(3,3,ep)=m12(3,3,ep)+fac*g12(ep)
259 m13(3,3,ep)=m13(3,3,ep)+fac*g13(ep)
260 m14(3,3,ep)=m14(3,3,ep)+fac*g14(ep)
261 m22(3,3,ep)=m22(3,3,ep)+fac*g22(ep)
262 m23(3,3,ep)=m23(3,3,ep)+fac*g23(ep)
263 m24(3,3,ep)=m24(3,3,ep)+fac*g24(ep)
264 m33(3,3,ep)=m33(3,3,ep)+fac*g33(ep)
265 m34(3,3,ep)=m34(3,3,ep)+fac*g34(ep)
266 m44(3,3,ep)=m44(3,3,ep)+fac*g44(ep)
278 2 HM ,HZ ,A_I ,PX1 ,PX2 ,
279 3 PY1 ,PY2 ,HXX ,HYY ,HXY ,
280 4 PH1 ,PH2 ,Z1 ,NPLAT,IPLAT,DHZ ,
281 5 K11,K12,K13,K14,K22,K23,K24,K33,K34,K44,
282 6 M11,M12,M13,M14,M22,M23,M24,M33,M34,M44,
283 7 MF11,MF12,MF13,MF14,MF22,MF23,MF24,MF33,
284 8 MF34,MF44,FM12,FM13,FM14,FM23,FM24,FM34,
285 9 PHKRX,PHKRY,PHKRXY,PHERX,PHERY,PHERXY,
286 A PHKRZ,PHERZ,PHKX ,PHKY ,PHEX ,PHEY )
291#include "implicit_f.inc"
292#include "mvsiz_p.inc"
297 INTEGER JFT,JLT,NPLAT,IPLAT(*)
299 . (*),THK0(*),THK2(*),A_I(*),(*),HM(MVSIZ,4),HZ(*),
300 . PX1(*) ,PX2(*) ,PY1(*) ,PY2(*) ,PH1(*) ,PH2(*) ,
302 . K11(3,3,*),K12(3,3,*),K13(3,3,*),K14(3,3,*),
303 . K22(3,3,*),K23(3,3,*),K24(3,3,*),K33(3,3,*),
304 . M11(3,3,*),M12(3,3,*),(3,3,*),M14(3,3,*),
305 . M22(3,3,*),M23(3,3,*),M24(3,3,*),M33(3,3,*),
306 . MF11(3,3,*),MF12(3,3,*),MF13(3,3,*),MF14(3,3,*),
307 . MF22(3,3,*),MF23(3,3,*),MF24(3,3,*),MF33(3,3,*),
308 . FM12(3,3,*),FM13(3,3,*),FM14(3,3,*),
309 . FM23(3,3,*),FM24(3,3,*),FM34(3,3,*),
310 . K34(3,3,*),K44(3,3,*),M34(3,3,*),M44(3,3,*),
311 . MF34(3,3,*),MF44(3,3,*),DHZ(*)
313 . PHKRX(4,*),PHKRY(4,*),
314 . phkrxy(4,*),pherx(4,*),phery(4,*),pherxy(4,*),
315 . phkrz(4,*),pherz(4,*),phkx(*) ,phky(*) ,phex(*) ,phey(*)
324 . DM(3,3,MVSIZ),C1,C2,GM(MVSIZ),
325 . G11(MVSIZ),G12(MVSIZ),G13(MVSIZ),
326 . G14(MVSIZ),G22(MVSIZ),G23(MVSIZ),
327 . G24(MVSIZ),G33(MVSIZ),G34(MVSIZ),G44(MVSIZ),
328 . CHM(2,2,MVSIZ),CHF(2,2,MVSIZ),FACF(MVSIZ),FAC,
329 . gama(4,mvsiz),cxx,cyy,cxy,c11,c22,c12,rz_3,
330 . cbkrx(4,mvsiz),cbkry(4,mvsiz),cbkrz(4,mvsiz),
331 . cberx(4,mvsiz),cbery(4,mvsiz),cberz(4,mvsiz),
332 . grz3(4,mvsiz),zr(mvsiz),c3,cx1,cx2,cy1,cy2,zr2
335#include "vectorize.inc"
340 dm(1,1,m)=hm(ep,1)*c1
341 dm(2,2,m)=hm(ep,2)*c1
342 dm(1,2,m)=hm(ep,3)*c1
346 dm(3,3,m)=hm(ep,4)*c1
359 chm(1,1,ep)=dhz(ep)*hyy(ep)
360 chm(2,2,ep)=dhz(ep)*hxx(ep)
361 chm(1,2,ep)=-dhz(ep)*hxy(ep)
362 chm(2,1,ep)=chm(1,2,ep)
366 gama(1,ep)=fourth-ph1(ep)
367 gama(2,ep)=-fourth-ph2(ep)
368 gama(3,ep)=fourth+ph1(ep)
369 gama(4,ep)=-fourth+ph2(ep)
370 g11(ep) =gama(1,ep)*gama(1,ep)
371 g12(ep) =gama(1,ep)*gama(2,ep)
372 g13(ep) =gama(1,ep)*gama(3,ep)
373 g14(ep) =gama(1,ep)*gama(4,ep)
374 g22(ep) =gama(2,ep)*gama(2,ep)
375 g23(ep) =gama(2,ep)*gama(3,ep)
376 g24(ep) =gama(2,ep)*gama(4,ep)
377 g33(ep) =gama(3,ep)*gama(3,ep)
378 g34(ep) =gama(3,ep)*gama(4,ep)
379 g44(ep) =gama(4,ep)*gama(4,ep)
385 k11(i,j,ep)=k11(i,j,ep)+chm(i,j,ep)*g11(ep)
386 k22(i,j,ep)=k22(i,j,ep)+chm(i,j,ep)*g22(ep)
387 k33(i,j,ep)=k33(i,j,ep)+chm(i,j,ep)*g33(ep)
388 k44(i,j,ep)=k44(i,j,ep)+chm(i,j,ep)*g44(ep)
395 k12(i,j,ep)=k12(i,j,ep)+chm(i,j,ep)*g12(ep)
396 k13(i,j,ep)=k13(i,j,ep)+chm(i,j,ep)*g13(ep)
397 k14(i,j,ep)=k14(i,j,ep)+chm(i,j,ep)*g14(ep)
398 k23(i,j,ep)=k23(i,j,ep)+chm(i,j,ep)*g23(ep)
399 k24(i,j,ep)=k24(i,j,ep)+chm(i,j,ep)*g24(ep)
400 k34(i,j,ep)=k34(i,j,ep)+chm(i,j,ep)*g34(ep)
408 grz3(j,m)=gama(j,m)*rz_3
415 c1 = -(phky(m)*phkrz(j,m)+phey(m)*pherz(j,m))
416 c2 = (phkx(m)*phkrz(j,m)+phex
417 mf11(1,3,m)=mf11(1,3,m)+grz3(i,m)*c1
418 mf11(2,3,m)=mf11(2,3,m)+grz3(i,m)*c2
420 fm12(3,1,m)=fm12(3,1,m)+ grz3(i,m)*c1
421 fm12(3,2,m)=fm12(3,2,m)+ grz3(i,m)*c2
423 fm13(3,1,m)=fm13(3,1,m)+grz3(i,m)*c1
424 fm13(3,2,m)=fm13(3,2,m)+grz3(i,m)*c2
426 fm14(3,1,m)=fm14(3,1,m)+ grz3(i,m)*c1
427 fm14(3,2,m)=fm14(3,2,m)+ grz3(i,m)*c2
429 c1 = -(phky(m)*phkrz(j,m)+phey(m)*pherz(j,m))
430 c2 = (phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m))
432 mf12(1,3,m)=mf12(1,3,m)+grz3(i,m)*c1
433 mf12(2,3,m)=mf12(2,3,m)+grz3(i,m)*c2
435 mf22(1,3,m)=mf22(1,3,m)+grz3(i,m)*c1
436 mf22(2,3,m)=mf22(2,3,m)+grz3(i,m)*c2
438 fm23(3,1,m)=fm23(3,1,m
439 fm23(3,2,m)=fm23(3,2,m)+grz3(i,m)*c2
441 fm24(3,1,m)=fm24(3,1,m)+ grz3(i,m)*c1
442 fm24(3,2,m)=fm24(3,2,m)+ grz3(i,m)*c2
444 c1 = -(phky(m)*phkrz(j,m)+phey(m)*pherz(j,m))
445 c2 = (phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m))
447 mf13(1,3,m)=mf13(1,3,m)+ grz3(i,m)*c1
448 mf13(2,3,m)=mf13(2,3,m)+ grz3(i,m)*c2
450 mf23(1,3,m)=mf23(1,3,m)+grz3(i,m)*c1
451 mf23(2,3,m)=mf23(2,3,m)+grz3(i,m)*c2
453 mf33(1,3,m)=mf33(1,3,m)+ grz3(i,m)*c1
454 mf33(2,3,m)=mf33(2,3,m)+ grz3(i,m)*c2
456 fm34(3,1,m)=fm34(3,1,m)+ grz3(i,m)*c1
457 fm34(3,2,m)=fm34(3,2,m)+ grz3(i,m)*c2
459 c1 = -(phky(m)*phkrz(j,m)+phey(m)*pherz(j,m))
460 c2 = (phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m))
462 mf14(1,3,m)=mf14(1,3,m)+ grz3(i,m)*c1
463 mf14(2,3,m)=mf14(2,3,m)+ grz3(i,m)*c2
465 mf24(1,3,m)=mf24(1,3,m)+ grz3(i,m)*c1
466 mf24(2,3,m)=mf24(2,3,m)+ grz3(i,m)*c2
468 mf34(1,3,m)=mf34(1,3,m)+ grz3(i,m)*c1
469 mf34(2,3,m)=mf34(2,3,m)+ grz3(i,m)*c2
471 mf44(1,3,m)=mf44(1,3,m)+ grz3(i,m)*c1
472 mf44(2,3,m)=mf44(2,3,m)+ grz3(i,m)*c2
479 m11(3,3,ep)=m11(3,3,ep)+
480 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
482 m12(3,3,ep)=m12(3,3,ep)+
483 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
485 m13(3,3,ep)=m13(3,3,ep)+
486 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
488 m14(3,3,ep)=m14(3,3,ep)+
489 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
492 m22(3,3,ep)=m22(3,3,ep)+
493 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
495 m23(3,3,ep)=m23(3,3,ep)+
496 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
498 m24(3,3,ep)=m24(3,3,ep)+
499 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
502 m33(3,3,ep)=m33(3,3,ep)+
503 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
505 m34(3,3,ep)=m34(3,3,ep)+
506 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
508 m44(3,3,ep)=m44(3,3,ep)+
509 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
518 cbkrx(j,m) =dm(1,1,m)*phkrx(j,m)+dm(1,2,m)*phkry(j,m)+
519 . dm(1,3,m)*phkrxy(j,m)
520 cbkry(j,m) =dm(2,1,m)*phkrx(j,m)+dm(2,2,m)*phkry(j,m)+
521 . dm(2,3,m)*phkrxy(j,m)
522 cbkrz(j,m) =dm(3,1,m)*phkrx(j,m)+dm(3,2,m)*phkry(j,m)+
523 . dm(3,3,m)*phkrxy(j,m)
524 cberx(j,m) =dm(1,1,m)*pherx(j,m)+dm(1,2,m)*phery(j,m)+
525 . dm(1,3,m)*pherxy(j,m)
526 cbery(j,m) =dm(2,1,m)*pherx(j,m)+dm(2,2,m)*phery(j,m)+
527 . dm(2,3,m)*pherxy(j,m)
528 cberz(j,m) =dm(3,1,m)*pherx(j,m)+dm(3,2,m)*phery(j,m)+
529 . dm(3,3,m)*pherxy(j,m)
536 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
537 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
538 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
539 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
540 mf11(1,3,m)=mf11(1,3,m)+gama(i,m)*c1
541 mf11(2,3,m)=mf11(2,3,m)+gama(i,m)*c2
543 fm12(3,1,m)=fm12(3,1,m)+ gama(i,m)*c1
544 fm12(3,2,m)=fm12(3,2,m)+ gama(i,m)*c2
546 fm13(3,1,m)=fm13(3,1,m)+gama(i,m)*c1
547 fm13(3,2,m)=fm13(3,2,m)+gama(i,m)*c2
549 fm14(3,1,m)=fm14(3,1,m)+ gama(i,m)*c1
550 fm14(3,2,m)=fm14(3,2,m)+ gama(i,m)*c2
552 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
553 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
554 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
555 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
557 mf12(1,3,m)=mf12(1,3,m)+gama(i,m)*c1
558 mf12(2,3,m)=mf12(2,3,m)+gama(i,m)*c2
560 mf22(1,3,m)=mf22(1,3,m)+gama(i,m)*c1
561 mf22(2,3,m)=mf22(2,3,m)+gama(i,m)*c2
563 fm23(3,1,m)=fm23(3,1,m)+ gama(i,m)*c1
564 fm23(3,2,m)=fm23(3,2,m)+ gama(i,m)*c2
566 fm24(3,1,m)=fm24(3,1,m)+ gama(i,m)*c1
567 fm24(3,2,m)=fm24(3,2,m)+ gama(i,m)*c2
569 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
570 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
571 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
572 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
574 mf13(1,3,m)=mf13(1,3,m)+ gama(i,m)*c1
575 mf13(2,3,m)=mf13(2,3,m)+ gama(i,m)*c2
577 mf23(1,3,m)=mf23(1,3,m)+gama(i,m)*c1
578 mf23(2,3,m)=mf23(2,3,m)+gama(i,m)*c2
580 mf33(1,3,m)=mf33(1,3,m)+ gama(i,m)*c1
581 mf33(2,3,m)=mf33(2,3,m)+ gama(i,m)*c2
583 fm34(3,1,m)=fm34(3,1,m)+ gama(i,m)*c1
584 fm34(3,2,m)=fm34(3,2,m)+ gama(i,m)*c2
587 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
588 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
589 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
591 mf14(1,3,m)=mf14(1,3,m)+ gama(i,m)*c1
592 mf14(2,3,m)=mf14(2,3,m)+ gama(i,m)*c2
594 mf24(1,3,m)=mf24(1,3,m)+ gama(i,m)*c1
595 mf24(2,3,m)=mf24(2,3,m)+ gama(i,m)*c2
597 mf34(1,3,m)=mf34(1,3,m)+ gama(i,m)*c1
598 mf34(2,3,m)=mf34(2,3,m)+ gama(i,m)*c2
600 mf44(1,3,m)=mf44(1,3,m)+ gama(i,m)*c1
601 mf44(2,3,m)=mf44(2,3,m)+ gama(i,m)*c2
607 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
608 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
609 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
610 m11(3,3,m)=m11(3,3,m)+c1
612 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
613 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
614 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
615 m12(3,3,m)=m12(3,3,m)+c1
617 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
618 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
619 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
620 m13(3,3,m)=m13(3,3,m)+c1
622 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
623 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
624 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
625 m14(3,3,m)=m14(3,3,m)+c1
628 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
629 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
630 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
631 m22(3,3,m)=m22(3,3,m)+c1
633 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
634 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
635 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
636 m23(3,3,m)=m23(3,3,m)+c1
638 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
639 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
640 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
641 m24(3,3,m)=m24(3,3,m)+c1
644 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
646 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
647 m33(3,3,m)=m33(3,3,m)+c1
649 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
650 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
651 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
652 m34(3,3,m)=m34(3,3,m)+c1
654 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
655 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
656 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
657 m44(3,3,m)=m44(3,3,m)+c1
663 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
664 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
665 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
666 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
668 c3 = zr(m)*(px1(m)*c1+py1(m)*c2)
669 mf11(3,3,m)=mf11(3,3,m)+ c3
670 fm13(3,3,m)=fm13(3,3,m)- c3
672 c3 = zr(m)*(px2(m)*c1+py2(m)*c2)
673 fm12(3,3,m)=fm12(3,3,m)+ c3
674 fm14(3,3,m)=fm14(3,3,m)- c3
676 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
677 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
678 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
679 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
681 c3 = zr(m)*(px1(m)*c1+py1(m)*c2)
682 mf12(3,3,m)=mf12(3,3,m)+ c3
683 fm23(3,3,m)=fm23(3,3,m)- c3
685 c3 = zr(m)*(px2(m)*c1+py2(m)*c2)
686 mf22(3,3,m)=mf22(3,3,m)+ c3
687 fm24(3,3,m)=fm24(3,3,m)- c3
689 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
690 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
691 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
692 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
694 c3 = zr(m)*(px1(m)*c1+py1(m)*c2)
695 mf13(3,3,m)=mf13(3,3,m)+ c3
696 mf33(3,3,m)=mf33(3,3,m)- c3
698 c3 = zr(m)*(px2(m)*c1+py2(m)*c2)
699 mf23(3,3,m)=mf23(3,3,m)+ c3
700 fm34(3,3,m)=fm34(3,3,m)- c3
702 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
703 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
704 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
705 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
707 c3 = zr(m)*(px1(m)*c1+py1(m)*c2)
708 mf14(3,3,m)=mf14(3,3,m)+ c3
709 mf34(3,3,m)=mf34(3,3,m)- c3
711 c3 = zr(m)*(px2(m)*c1+py2(m)*c2)
712 mf24(3,3,m)=mf24(3,3,m)+ c3
713 mf44(3,3,m)=mf44(3,3,m)- c3
724 c11 = zr(m)*gama(j,m)
726 c1 = (chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))*c11
727 c2 = (chm(1,2,m)*px1(m)+chm(2,2,m)*py1(m))*c11
728 k11(1,3,m)=k11(1,3,m)+ c1
729 k13(1,3,m)=k13(1,3,m)- c1
730 k11(2,3,m)=k11(2,3,m)+ c2
731 k13(2,3,m)=k13(2,3,m)- c2
733 c1 = (chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))*c11
734 c2 = (chm(1,2,m)*px2(m)+chm(2,2,m)*py2(m))*c11
735 k12(1,3,m)=k12(1,3,m)+ c1
736 k14(1,3,m)=k14(1,3,m)- c1
737 k12(2,3,m)=k12(2,3,m)+ c2
738 k14(2,3,m)=k14(2,3,m)- c2
742 c11 = zr(m)*gama(j,m)
744 c1 = (chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))*c11
745 c2 = (chm(1,2,m)*px1(m)+chm(2,2,m)*py1(m))*c11
746 k12(3,1,m)=k12(3,1,m)+ c1
747 k23(1,3,m)=k23(1,3,m)- c1
748 k12(3,2,m)=k12(3,2,m)+ c2
749 k23(2,3,m)=k23(2,3,m)- c2
751 c1 = (chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))*c11
752 c2 = (chm(1,2,m)*px2(m)+chm(2,2,m)*py2(m))*c11
753 k22(1,3,m)=k22(1,3,m)+ c1
754 k24(1,3,m)=k24(1,3,m)- c1
755 k22(2,3,m)=k22(2,3,m)+ c2
756 k24(2,3,m)=k24(2,3,m)- c2
760 c11 = zr(m)*gama(j,m)
762 c1 = (chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))*c11
763 c2 = (chm(1,2,m)*px1(m)+chm(2,2,m)*py1(m))*c11
764 k13(3,1,m)=k13(3,1,m)+ c1
766 k13(3,2,m)=k13(3,2,m)+ c2
767 k33(2,3,m)=k33(2,3,m)- c2
769 c1 = (chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))*c11
770 c2 = (chm(1,2,m)*px2(m)+chm(2,2,m
771 k23(3,1,m)=k23(3,1,m)+ c1
772 k34(1,3,m)=k34(1,3,m)- c1
773 k23(3,2,m)=k23(3,2,m)+ c2
774 k34(2,3,m)=k34(2,3,m)- c2
778 c11 = zr(m)*gama(j,m)
780 c1 = (chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))*c11
781 c2 = (chm(1,2,m)*px1(m)+chm(2,2,m)*py1(m))*c11
782 k14(3,1,m)=k14(3,1,m)+ c1
783 k34(3,1,m)=k34(3,1,m)- c1
784 k14(3,2,m)=k14(3,2,m)+ c2
785 k34(3,2,m)=k34(3,2,m)- c2
787 c1 = (chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))*c11
789 k24(3,1,m)=k24(3,1,m)+ c1
790 k44(1,3,m)=k44(1,3,m)- c1
791 k24(3,2,m)=k24(3,2,m)+ c2
792 k44(2,3,m)=k44(2,3,m)- c2
800 cx1=zr2*(chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))
804 k11(3,3,m)=k11(3,3,m)+ c3
806 k13(3,3,m)=k13(3,3,m)- c3
807 k33(3,3,m)=k33(3,3,m)+ c3
809 c3=px2(m)*cx1+py2(m)*cy1
811 k23(3,3,m)=k23(3,3,m)- c3
815 cx1=zr2*(chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))
816 cy1=zr2*(chm(1,2,m)*px2(m)+chm(2,2,m)*py2(m))
818 c3=px1(m)*cx1+py1(m)*cy1
819 k12(3,3,m)=k12(3,3,m)+ c3
821 k14(3,3,m)=k14(3,3,m)- c3
822 k34(3,3,m)=k34(3,3,m)+ c3
824 c3=px2(m)*cx1+py2(m)*cy1
825 k22(3,3,m)=k22(3,3,m)+ c3
827 k24(3,3,m)=k24(3,3,m)- c3
828 k44(3,3,m)=k44(3,3,m)+ c3
832 rz_3= zr(m)*dhz(m)*third
836 cx2 = phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m)
837 cy2 = phky(m)*phkrz(j,m)+phey(m)*pherz(j,m)
839 c3 = rz_3*(-px1(m)*cy2+py1(m)*cx2)
840 mf11(3,3,m)=mf11(3,3,m)+ c3
841 fm13(3,3,m)=fm13(3,3,m)- c3
843 c3 = rz_3*(-px2(m)*cy2+py2(m)*cx2)
844 fm12(3,3,m)=fm12(3,3,m)+ c3
845 fm14(3,3,m)=fm14(3,3,m)- c3
849 cx2 = phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m)
850 cy2 = phky(m)*phkrz(j,m)+phey(m)*pherz(j,m)
852 c3 = rz_3*(-px1(m)*cy2+py1(m)*cx2)
853 mf12(3,3,m)=mf12(3,3,m)+ c3
854 fm23(3,3,m)=fm23(3,3,m)- c3
856 c3 = rz_3*(-px2(m)*cy2+py2(m)*cx2)
857 mf22(3,3,m)=mf22(3,3,m)+ c3
858 fm24(3,3,m)=fm24(3,3,m)- c3
862 cx2 = phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m)
863 cy2 = phky(m)*phkrz(j,m)+phey(m)*pherz(j,m)
865 c3 = rz_3*(-px1(m)*cy2+py1(m)*cx2)
866 mf13(3,3,m)=mf13(3,3,m)+ c3
867 mf33(3,3,m)=mf33(3,3,m)- c3
869 c3 = rz_3*(-px2(m)*cy2+py2(m)*cx2)
870 mf23(3,3,m)=mf23(3,3,m)+ c3
871 fm34(3,3,m)=fm34(3,3,m)- c3
875 cx2 = phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m)
876 cy2 = phky(m)*phkrz(j,m)+phey(m)*pherz(j,m)
878 c3 = rz_3*(-px1(m)*cy2+py1(m)*cx2)
879 mf14(3,3,m)=mf14(3,3,m)+ c3
880 mf34(3,3,m)=mf34(3,3,m)- c3
882 c3 = rz_3*(-px2(m)*cy2+py2(m)*cx2)
883 mf24(3,3,m)=mf24(3,3,m)+ c3
884 mf44(3,3,m)=mf44(3,3,m)- c3