36
37
38
39
40#include "implicit_f.inc"
41#include "mvsiz_p.inc"
42#include "com04_c.inc"
43
44
45
46
47 INTEGER JFT,JLT,IKGEO,IORTH,IDRIL,NEL
49 .
area(*),vol(*), px1(*),py1(*),py2(*),thk0(*),thk2(*),
50 . hm(mvsiz,4),hf(mvsiz,4),hc(mvsiz,2),hz(*),hmor(mvsiz,2),hfor(mvsiz,2),
51 . k11(3,3,*),k12(3,3,*),k13(3,3,*),
52 . k22(3,3,*),k23(3,3,*),k33(3,3,*),
53 . m11(3,3,*),m12(3,3,*),m13(3,3,*),
54 . m22(3,3,*),m23(3,3,*),m33(3,3,*),
55 . mf11(3,3,*),mf12(3,3,*),mf13(3,3,*),
56 . mf22(3,3,*),mf23(3,3,*),mf33(3,3,*),
57 . fm12(3,3,*),fm13(3,3,*),fm23(3,3,*),
for(nel,5),mom(nel,3),
58 . hmfor(mvsiz,6)
59
60
61
62 INTEGER EP,I,J,L,K,I1,J1,IN(2),NF
64 . dm(2,2,mvsiz),df(2,2,mvsiz),dcx(mvsiz),dcy(mvsiz),
65 . px1px1(mvsiz),px1py1(mvsiz),px1py2(mvsiz),px1py3(mvsiz),
66 . py1py1(mvsiz),py1py2(mvsiz),py1py3(mvsiz),py3(mvsiz),
67 . py2py2(mvsiz),py2py3(mvsiz),py3py3(mvsiz),
68 . c1,c2,dhz(mvsiz),gm(mvsiz),gf(mvsiz),c3,
69 . gpx1px1,gpx1py1,gpx1py2,gpx1py3,bcxx1,bcxx2,bcxy3,bcxy1,
70 . bcxy2,bcyy1,bcyy2,bcyy3,bcyx1,bcyx2,bcyx3,
71 . px1dx,py1dy,py2dy,py3dy,fac,fxx,fyy
72 . h11(mvsiz),h12(mvsiz),h22(mvsiz),h13(mvsiz),h23(mvsiz),
73 . dm5(mvsiz),dm6(mvsiz),df5(mvsiz),df6(mvsiz),dm5_2(mvsiz),
74 . dm6_2(mvsiz),df5_2(mvsiz),df6_2(mvsiz),
75 . dmf(3,3,mvsiz),mfij(2,2,mvsiz)
76
77 DO ep=jft,jlt
78 c2=vol(ep)
79 c1=thk2(ep)*c2
80 dm(1,1,ep)=hm(ep,1)*c2
81 dm(2,2,ep)=hm(ep,2)*c2
82 dm(1,2,ep)=hm(ep,3)*c2
83 dm(2,1,ep)=dm(1,2,ep)
84 gm(ep) =hm(ep,4)*c2
85 df(1,1,ep)=hf(ep,1)*c1
86 df(2,2,ep)=hf(ep,2)*c1
87 df(1,2,ep)=hf(ep,3)*c1
88 df(2,1,ep)=df(1,2,ep)
89 gf(ep) =hf(ep,4)*c1
90 dhz(ep)=hz(ep)*c1
91 dcx(ep)=hc(ep,1)*c2
92 dcy(ep)=hc(ep,2)*c2
93 ENDDO
94 DO ep=jft,jlt
95
96 py3(ep)= -py1(ep)-py2(ep)
97 px1px1(ep)=px1(ep)*px1(ep)
98 px1py1(ep)=px1(ep)*py1(ep)
99 px1py2(ep)=px1(ep)*py2(ep)
100 px1py3(ep)=-px1py1(ep)-px1py2(ep)
101 py1py1(ep)=py1(ep)*py1(ep)
102 py1py2(ep)=py1(ep)*py2(ep)
103 py1py3(ep)=-py1py1(ep)-py1py2(ep)
104 py2py2(ep)=py2(ep)*py2(ep)
105 py2py3(ep)=py2(ep)*py3(ep)
106 py3py3(ep)=py3(ep)*py3(ep)
107 ENDDO
108
109 DO ep=jft,jlt
110 k11(1,1,ep)=px1px1(ep)*dm(1,1,ep)
111 k11(1,2,ep)=px1py1(ep)*dm(1,2,ep)
112 k11(2,2,ep)=py1py1(ep)*dm(2,2,ep)
113 k22(1,1,ep)=k11(1,1,ep)
114 k22(1,2,ep)=-px1py2(ep)*dm(1,2,ep)
115 k22(2,2,ep)=py2py2(ep)*dm(2,2,ep)
116 k33(2,2,ep)=py3py3(ep)*dm(2,2,ep)
117
118 k12(1,1,ep)=-k11(1,1,ep)
119 k12(1,2,ep)=-k22(1,2,ep)
120 k12(2,1,ep)=-k11(1,2,ep)
121 k12(2,2,ep)=py1py2(ep)*dm(2,2,ep)
122 k13(1,2,ep)=px1py3(ep)*dm(1,2,ep)
123 k13(2,2,ep)=py1py3(ep)*dm(2,2,ep)
124 k23(1,2,ep)=-k13(1,2,ep)
125 k23(2,2,ep)=py2py3(ep)*dm(2,2,ep)
126 ENDDO
127
128 DO ep=jft,jlt
129 m11(2,2,ep)=px1px1(ep)*df(1,1,ep)
130 m11(1,2,ep)=-px1py1(ep)*df(1,2,ep)
131 m11(1,1,ep)=py1py1(ep)*df(2,2,ep)
132 m22(2,2,ep)=m11(2,2,ep)
133 m22(1,2,ep)=px1py2(ep)*df(1,2,ep)
134 m22(1,1,ep)=py2py2(ep)*df(2,2,ep)
135 m33(1,1,ep)=py3py3(ep)*df(2,2,ep)
136
137 m12(2,2,ep)=-m11(2,2,ep)
138 m12(1,2,ep)=-m11(1,2,ep)
139 m12(2,1,ep)=-m22(1,2,ep)
140 m12(1,1,ep)=py1py2(ep)*df(2,2,ep)
141 m13(2,1,ep)=-px1py3(ep)*df(1,2,ep)
142 m13(1,1,ep)=py1py3(ep)*df(2,2,ep)
143 m23(2,1,ep)=-m13(2,1,ep)
144 m23(1,1,ep)=py2py3(ep)*df(2,2,ep)
145 ENDDO
146
147 DO ep=jft,jlt
148 gpx1px1=px1px1(ep)*gm(ep)
149 gpx1py1=px1py1(ep)*gm(ep)
150 gpx1py2=px1py2(ep)*gm(ep)
151 gpx1py3=-gpx1py1-gpx1py2
152 k11(1,1,ep)=k11(1,1,ep)+py1py1(ep)*gm(ep)
153 k11(1,2,ep)=k11(1,2,ep)+gpx1py1
154 k11(2,2,ep)=k11(2,2,ep)+gpx1px1
155 k22(1,1,ep)=k22(1,1,ep)+py2py2(ep)*gm(ep)
156 k22(1,2,ep)=k22(1,2,ep)-gpx1py2
157 k22(2,2,ep)=k22(2,2,ep)+gpx1px1
158 k33(1,1,ep)=k33(1,1,ep)+py3py3(ep)*gm(ep)
159
160 k12(1,1,ep)=k12(1,1,ep)+py1py2(ep)*gm(ep)
161 k12(1,2,ep)=k12(1,2,ep)-gpx1py1
162 k12(2,1,ep)=k12(2,1,ep)+gpx1py2
163 k12(2,2,ep)=k12(2,2,ep)-gpx1px1
164 k13(1,1,ep)=k13(1,1,ep)+py1py3(ep)*gm(ep)
165 k13(2,1,ep)=k13(2,1,ep)+gpx1py3
166 k23(1,1,ep)=k23(1,1,ep)+py2py3(ep)*gm(ep)
167 k23(2,1,ep)=k23(2,1,ep)-gpx1py3
168 ENDDO
169
170 DO ep=jft,jlt
171 gpx1px1=px1px1(ep)*gf(ep)
172 gpx1py1=px1py1(ep)*gf(ep)
173 gpx1py2=px1py2(ep)*gf(ep)
174 gpx1py3=-gpx1py1-gpx1py2
175 m11(2,2,ep)=m11(2,2,ep)+py1py1(ep)*gf(ep)
176 m11(1,2,ep)=m11(1,2,ep)-gpx1py1
177 m11(1,1,ep)=m11(1,1,ep)+gpx1px1
178 m22(2,2,ep)=m22(2,2,ep)+py2py2(ep)*gf(ep)
179 m22(1,2,ep)=m22(1,2,ep)+gpx1py2
180 m22(1,1,ep)=m22(1,1,ep)+gpx1px1
181 m33(2,2,ep)=m33(2,2,ep)+py3py3(ep)*gf(ep)
182
183 m12(2,2,ep)=m12(2,2,ep)+py1py2(ep)*gf(ep)
184 m12(1,2,ep)=m12(1,2,ep)-gpx1py2
185 m12(2,1,ep)=m12(2,1,ep)+gpx1py1
186 m12(1,1,ep)=m12(1,1,ep)-gpx1px1
187 m13(1,2,ep)=m13(1,2,ep)-gpx1py3
188 m13(2,2,ep)=m13(2,2,ep)+py1py3(ep)*gf(ep)
189 m23(1,2,ep)=m23(1,2,ep)+gpx1py3
190 m23(2,2,ep)=m23(2,2,ep)+py2py3(ep)*gf(ep)
191 ENDDO
192 IF (iorth==1) THEN
193#include "vectorize.inc"
194 DO ep=jft,jlt
195 c2=vol(ep)
196 c1=thk2(ep)*c2
197 dm5(ep)=hmor(ep,1)*c2
198 dm6(ep)=hmor(ep,2)*c2
199 df5(ep)=hfor(ep,1)*c1
200 df6(ep)=hfor(ep,2)*c1
201 ENDDO
202 DO ep=jft,jlt
203 dm5_2(ep)=two*dm5(ep)
204 dm6_2(ep)=two*dm6(ep)
205 df5_2(ep)=two*df5(ep)
206 df6_2(ep)=two*df6(ep)
207 ENDDO
208
209 DO ep=jft,jlt
210 k11(1,1,ep)=k11(1,1,ep)+px1py1(ep)*dm5_2(ep)
211 k11(1,2,ep)=k11(1,2,ep)+
212 1 px1px1(ep)*dm5(ep)+py1py1(ep)*dm6(ep)
213 k11(2,2,ep)=k11(2,2,ep)+px1py1(ep)*dm6_2(ep)
214 k12(1,1,ep)=k12(1,1,ep)+(px1py2(ep)-px1py1(ep))*dm5(ep)
215 c1 = -px1px1(ep)*dm5(ep)+py1py2(ep)*dm6(ep)
216 k12(1,2,ep)=k12(1,2,ep)+c1
217 k12(2,1,ep)=k12(2,1,ep)+c1
218 k12(2,2,ep)=k12(2,2,ep)+(px1py2(ep)-px1py1(ep))*dm6(ep)
219 k13(1,1,ep)=k13(1,1,ep)+px1py3(ep)*dm5(ep)
220 c1 = py1py3(ep)*dm6(ep)
221 k13(1,2,ep)=k13(1,2,ep)+c1
222 k13(2,1,ep)=k13(2,1,ep)+c1
223 k13(2,2,ep)=k13(2,2,ep)+px1py3(ep)*dm6(ep)
224
225 k22(1,1,ep)=k22(1,1,ep)-px1py2(ep)*dm5_2(ep)
226 k22(1,2,ep)=k22(1,2,ep)+
227 1 px1px1(ep)*dm5(ep)+py2py2(ep)*dm6(ep)
228 k22(2,2,ep)=k22(2,2,ep)-px1py2(ep)*dm6_2(ep)
229 k23(1,1,ep)=k23(1,1,ep)-px1py3(ep)*dm5(ep)
230
231 c1 = py2py3(ep)*dm6(ep)
232 k23(1,2,ep)=k23(1,2,ep)+c1
233 k23(2,1,ep)=k23(2,1,ep)+c1
234 k23(2,2,ep)=k23(2,2,ep)-px1py3(ep)*dm6(ep)
235
236 k33(1,2,ep)=k33(1,2,ep)+py3py3(ep)*dm6(ep)
237
238 m11(1,1,ep)=m11(1,1,ep)+px1py1(ep)*df6_2(ep)
239 m11(1,2,ep)=m11(1,2,ep)-
240 1 px1px1(ep)*df5(ep)-py1py1(ep)*df6(ep)
241 m11(2,2,ep)=m11(2,2,ep)+px1py1(ep)*df5_2(ep)
242 m12(1,1,ep)=m12(1,1,ep)+(px1py2(ep)-px1py1(ep))*df6(ep)
243 c2 = -px1px1(ep)*df5(ep)+py1py2(ep)*df6(ep)
244 m12(1,2,ep)=m12(1,2,ep)-c2
245 m12(2,1,ep)=m12(2,1,ep)-c2
246 m12(2,2,ep)=m12(2,2,ep)+(px1py2(ep)-px1py1(ep))*df5(ep)
247 m13(1,1,ep)=m13(1,1,ep)+px1py3(ep)*df6(ep)
248 m13(1,2,ep)=m13(1,2,ep)-py1py3(ep)*df6(ep)
249 m13(2,1,ep)=m13(2,1,ep)-py1py3(ep)*df6(ep)
250 m13(2,2,ep)=m13(2,2,ep)+px1py3(ep)*df5(ep)
251 m22(1,1,ep)=m22(1,1,ep)-px1py2(ep)*df6_2(ep)
252 m22(1,2,ep)=m22(1,2,ep)-
253 1 px1px1(ep)*df5(ep)-py2py2(ep)*df6(ep)
254 m22(2,2,ep)=m22(2,2,ep)-px1py2(ep)*df5_2(ep)
255 m23(1,1,ep)=m23(1,1,ep)-px1py3(ep)*df6(ep)
256 m23(1,2,ep)=m23(1,2,ep)-py2py3(ep)*df6(ep)
257 m23(2,1,ep)=m23(2,1,ep)-py2py3(ep)*df6(ep)
258 m23(2,2,ep)=m23(2,2,ep)-px1py3(ep)*df5(ep)
259 m33(1,2,ep)=m33(1,2,ep)-py3py3(ep)*df6(ep)
260 ENDDO
261 DO ep=jft,jlt
262 c2=vol(ep)*thk0(ep)
263 dmf(1,1,ep)=hmfor(ep,1)*c2
264 dmf(2,2,ep)=hmfor(ep,2)*c2
265 dmf(1,2,ep)=hmfor(ep,3)*c2
266 dmf(1,3,ep)=hmfor(ep,5)*c2
267 dmf(2,3,ep)=hmfor(ep,6)*c2
268 dmf(2,1,ep)=dmf(1,2,ep)
269 dmf(3,1,ep)=dmf(1,3,ep)
270 dmf(3,2,ep)=dmf(2,3,ep)
271 dmf(3,3,ep)=hmfor(ep,4)*c2
272 ENDDO
273
274 DO ep=jft,jlt
275 mf11(1,1,ep)=
276 1 -dmf(1,2,ep)*px1py1(ep)-dmf(1,3,ep)*px1px1(ep)
277 2 -dmf(2,3,ep)*py1py1(ep)-dmf(3,3,ep)*px1py1(ep)
278 mf11(1,2,ep)=
279 1 dmf(1,1,ep)*px1px1(ep)+dmf(1,3,ep)*px1py1(ep)
280 2 +dmf(1,3,ep)*px1py1(ep)+dmf(3,3,ep)*py1py1(ep)
281 mf11(2,1,ep)=
282 1 -dmf(2,2,ep)*py1py1(ep)-dmf(2,3,ep)*px1py1(ep)
283 2 -dmf(2,3,ep)*px1py1(ep)-dmf(3,3,ep)*px1px1(ep)
284 mf11(2,2,ep)=
285 1 dmf(1,2,ep)*px1py1(ep)+dmf(2,3,ep)*py1py1(ep)
286 2 +dmf(1,3,ep)*px1px1(ep)+dmf(3,3,ep)*px1py1(ep)
287 ENDDO
288 DO ep=jft,jlt
289 mf12(1,1,ep)=
290 1 -dmf(1,2,ep)*px1py2(ep)+dmf(1,3,ep)*px1px1(ep)
291
292 2 -dmf(2,3,ep)*py1py2(ep)+dmf(3,3,ep)*px1py1(ep)
293 mf12(1,2,ep)=
294 1 -dmf(1,1,ep)*px1px1(ep)+dmf(1,3,ep)*px1py2(ep)
295
296 2 -dmf(1,3,ep)*px1py1(ep)+dmf(3,3,ep)*py1py2(ep)
297 mf12(2,1,ep)=
298
299 1 -dmf(2,2,ep)*py1py2(ep)+dmf(2,3,ep)*px1py1(ep)
300 2 -dmf(2,3,ep)*px1py2(ep)+dmf(3,3,ep)*px1px1(ep)
301 mf12(2,2,ep)=
302
303 1 -dmf(1,2,ep)*px1py1(ep)+dmf(2,3,ep)*py1py2(ep)
304 2 -dmf(1,3,ep)*px1px1(ep)+dmf(3,3,ep)*px1py2(ep)
305 ENDDO
306 DO ep=jft,jlt
307 mf22(1,1,ep)=
308 1 dmf(1,2,ep)*px1py2(ep)-dmf(1,3,ep)*px1px1(ep)
309
310 2 -dmf(2,3,ep)*py2py2(ep)+dmf(3,3,ep)*px1py2(ep)
311 mf22(1,2,ep)=
312 1 dmf(1,1,ep)*px1px1(ep)-dmf(1,3,ep)*px1py2(ep)
313
314 2 -dmf(1,3,ep)*px1py2(ep)+dmf(3,3,ep)*py2py2(ep)
315 mf22(2,1,ep)=
316
317 1 -dmf(2,2,ep)*py2py2(ep)+dmf(2,3,ep)*px1py2(ep)
318 2 +dmf(2,3,ep)*px1py2(ep)-dmf(3,3,ep)*px1px1(ep)
319 mf22(2,2,ep)=
320
321 1 -dmf(1,2,ep)*px1py2(ep)+dmf(2,3,ep)*py2py2(ep)
322 2 +dmf(1,3,ep)*px1px1(ep)-dmf(3,3,ep)*px1py2(ep)
323 ENDDO
324 DO ep=jft,jlt
325 mf23(1,1,ep)=
326 1 dmf(1,2,ep)*px1py3(ep)
327
328 2 -dmf(2,3,ep)*py2py3(ep)
329 mf23(1,2,ep)=
330 1 -dmf(1,3,ep)*px1py3(ep)
331
332 2 +dmf(3,3,ep)*py2py3(ep)
333 mf23(2,1,ep)=
334
335 1 -dmf(2,2,ep)*py2py3(ep)
336 2 +dmf(2,3,ep)*px1py3(ep)
337 mf23(2,2,ep)=
338
339 1 dmf(2,3,ep)*py2py3(ep)
340 2 -dmf(3,3,ep)*px1py3(ep)
341 ENDDO
342 DO ep=jft,jlt
343 mf13(1,1,ep)=
344 1 -dmf(1,2,ep)*px1py3(ep)
345
346 2 -dmf(2,3,ep)*py1py3(ep)
347 mf13(1,2,ep)=
348 1 dmf(1,3,ep)*px1py3(ep)
349
350 2 +dmf(3,3,ep)*py1py3(ep)
351 mf13(2,1,ep)=
352
353 1 -dmf(2,2,ep)*py1py3(ep)
354 2 -dmf(2,3,ep)*px1py3(ep)
355 mf13(2,2,ep)=
356
357 1 dmf(2,3,ep)*py1py3(ep)
358 2 +dmf(3,3,ep)*px1py3(ep)
359 ENDDO
360 DO ep=jft,jlt
361 mf33(1,1,ep)=
362 2 -dmf(2,3,ep)*py3py3(ep)
363 mf33(1,2,ep)=
364 2 dmf(3,3,ep)*py3py3(ep)
365 mf33(2,1,ep)=
366 1 -dmf(2,2,ep)*py3py3(ep)
367 mf33(2,2,ep)=
368 1 dmf(2,3,ep)*py3py3(ep)
369 ENDDO
370
371 DO ep=jft,jlt
372 fm12(1,1,ep)=
373 1 dmf(1,2,ep)*px1py1(ep)+dmf(1,3,ep)*px1px1(ep)
374
375 2 -dmf(2,3,ep)*py1py2(ep)-dmf(3,3,ep)*px1py2(ep)
376 fm12(2,1,ep)=
377 1 -dmf(1,1,ep)*px1px1(ep)-dmf(1,3,ep)*px1py1(ep)
378
379
380 2 +dmf(1,3,ep)*px1py2(ep)+dmf(3,3,ep)*py1py2(ep)
381 fm12(1,2,ep)=
382
383 1 -dmf(2,2,ep)*py1py2(ep)-dmf(2,3,ep)*px1py2(ep)
384 2 +dmf(2,3,ep)*px1py1(ep)+dmf(3,3,ep)*px1px1(ep)
385 fm12(2,2,ep)=
386
387 1 dmf(1,2,ep)*px1py2(ep)+dmf(2,3,ep)*py1py2(ep)
388 2 -dmf(1,3,ep)*px1px1(ep)-dmf(3,3,ep)*px1py1(ep)
389 ENDDO
390 DO ep=jft,jlt
391 fm13(1,1,ep)=
392
393
394 2 -dmf(2,3,ep)*py1py3(ep)-dmf(3,3,ep)*px1py3(ep)
395 fm13(2,1,ep)=
396
397
398 2 dmf(1,3,ep)*px1py3(ep)+dmf(3,3,ep)*py1py3(ep)
399 fm13(1,2,ep)=
400
401 1 -dmf(2,2,ep)*py1py3(ep)-dmf(2,3,ep)*px1py3(ep)
402
403 fm13(2,2,ep)=
404
405 1 dmf(1,2,ep)*px1py3(ep)+dmf(2,3,ep)*py1py3(ep)
406
407 ENDDO
408 DO ep=jft,jlt
409 fm23(1,1,ep)=
410
411
412 2 -dmf(2,3,ep)*py2py3(ep)+dmf(3,3,ep)*px1py3(ep)
413 fm23(2,1,ep)=
414
415
416 2 -dmf(1,3,ep)*px1py3(ep)+dmf(3,3,ep)*py2py3(ep)
417 fm23(1,2,ep)=
418
419 1 -dmf(2,2,ep)*py2py3(ep)+dmf(2,3,ep)*px1py3(ep)
420
421 fm23(2,2,ep)=
422
423 1 -dmf(1,2,ep)*px1py3(ep)+dmf(2,3,ep)*py2py3(ep)
424
425 ENDDO
426 ENDIF
427
428 DO ep=jft,jlt
429 gpx1px1=px1px1(ep)*dcx(ep)
430 k11(3,3,ep)=gpx1px1+py1py1(ep)*dcy(ep)
431 k22(3,3,ep)=gpx1px1+py2py2(ep)*dcy(ep)
432 k33(3,3,ep)=py3py3(ep)*dcy(ep)
433 k12(3,3,ep)=-gpx1px1+py1py2(ep)*dcy(ep)
434 k13(3,3,ep)=py1py3(ep)*dcy(ep)
435 k23(3,3,ep)=py2py3(ep)*dcy(ep)
436 ENDDO
437 DO ep=jft,jlt
439 bcxx1=-fac*px1px1(ep)
440 bcxx2= -bcxx1
441 bcxy3= fac*(px1py1(ep)+px1py2(ep))
442 bcxy1= bcxy3+bcxy3+fac*px1py2(ep)
443 bcxy2= bcxy3+bcxy3+fac*px1py1(ep)
444 bcyy1= fac*(py1py1(ep)+py1py2(ep)+py1py2(ep))
445 bcyy2= -fac*(py2py2(ep)+py1py2(ep)+py1py2(ep))
446 bcyy3= -bcyy1-bcyy2
447 bcyx1= -bcxy3-fac*px1py1(ep)
448 bcyx2= -bcxy3-fac*px1py2(ep)
449 bcyx3= bcyx1+bcyx2
450 m11(1,1,ep)=m11(1,1,ep)+bcxx1*bcxx1*dcx(ep)+bcyx1*bcyx1*dcy(ep)
451 m11(1,2,ep)=m11(1,2,ep)+bcxx1*bcxy1*dcx(ep)+bcyx1*bcyy1*dcy(ep)
452 m11(2,2,ep)=m11(2,2,ep)+bcxy1*bcxy1*dcx(ep)+bcyy1*bcyy1*dcy(ep)
453 m22(1,1,ep)=m22(1,1,ep)+bcxx2*bcxx2*dcx(ep)+bcyx2*bcyx2*dcy(ep)
454 m22(1,2,ep)=m22(1,2,ep)+bcxx2*bcxy2*dcx(ep)+bcyx2*bcyy2*dcy(ep)
455 m22(2,2,ep)=m22(2,2,ep)+bcxy2*bcxy2*dcx(ep)+bcyy2*bcyy2*dcy(ep)
456 m33(1,1,ep)=m33(1,1,ep)+bcyx3*bcyx3*dcy(ep)
457 m33(1,2,ep)=m33(1,2,ep)+bcyx3*bcyy3*dcy(ep)
458 m33(2,2,ep)=m33(2,2,ep)+bcxy3*bcxy3*dcx(ep)+bcyy3*bcyy3*dcy(ep)
459 m12(1,1,ep)=m12(1,1,ep)+bcxx1*bcxx2*dcx(ep)+bcyx1*bcyx2*dcy(ep)
460 m12(1,2,ep)=m12(1,2,ep)+bcxx1*bcxy2*dcx(ep)+bcyx1*bcyy2*dcy(ep)
461 m12(2,1,ep)=m12(2,1,ep)+bcxx2*bcxy1*dcx(ep)+bcyx2*bcyy1*dcy(ep)
462 m12(2,2,ep)=m12(2,2,ep)+bcxy1*bcxy2*dcx(ep)+bcyy1*bcyy2*dcy(ep)
463 m13(1,1,ep)=m13(1,1,ep)+bcyx1*bcyx3*dcy(ep)
464 m13(1,2,ep)=m13(1,2,ep)+bcxx1*bcxy3*dcx(ep)+bcyx1*bcyy3*dcy(ep)
465 m13(2,1,ep)=m13(2,1,ep)+bcyx3*bcyy1*dcy(ep)
466 m13(2,2,ep)=m13(2,2,ep)+bcxy1*bcxy3*dcx(ep)+bcyy1*bcyy3*dcy(ep)
467 m23(1,1,ep)=m23(1,1,ep)+bcyx2*bcyx3*dcy(ep)
468 m23(1,2,ep)=m23(1,2,ep)+bcxx2*bcxy3*dcx(ep)+bcyx2*bcyy3*dcy(ep)
469 m23(2,1,ep)=m23(2,1,ep)+bcyx3*bcyy2*dcy(ep)
470 m23(2,2,ep)=m23(2,2,ep)+bcxy2*bcxy3*dcx(ep)+bcyy2*bcyy3*dcy(ep)
471
472 px1dx=px1(ep)*dcx(ep)
473 py1dy=py1(ep)*dcy(ep)
474 py2dy=py2(ep)*dcy(ep)
475 py3dy=py3(ep)*dcy(ep)
476 mf11(3,1,ep)=bcxx1*px1dx+bcyx1*py1dy
477 mf11(3,2,ep)=bcxy1*px1dx+bcyy1*py1dy
478 mf22(3,1,ep)=-bcxx2*px1dx+bcyx2*py2dy
479 mf22(3,2,ep)=-bcxy2*px1dx+bcyy2*py2dy
480 mf33(3,1,ep)=bcyx3*py3dy
481 mf33(3,2,ep)=bcyy3*py3dy
482 mf12(3,1,ep)=bcxx2*px1dx+bcyx2*py1dy
483 mf12(3,2,ep)=bcxy2*px1dx+bcyy2*py1dy
484 mf13(3,1,ep)=bcyx3*py1dy
485 mf13(3,2,ep)=bcxy3*px1dx+bcyy3*py1dy
486 mf23(3,1,ep)=bcyx3*py2dy
487 mf23(3,2,ep)=-bcxy3*px1dx+bcyy3*py2dy
488 fm12(1,3,ep)=-bcxx1*px1dx+bcyx1*py2dy
489 fm12(2,3,ep)=-bcxy1*px1dx+bcyy1*py2dy
490 fm13(1,3,ep)=bcyx1*py3dy
491 fm13(2,3,ep)=bcyy1*py3dy
492 fm23(1,3,ep)=bcyx2*py3dy
493 fm23(2,3,ep)=bcyy2*py3dy
494 ENDDO
495
496
497 DO ep=jft,jlt
498 m11(3,3,ep)=m11(3,3,ep)+dhz(ep)*(px1px1(ep)+py1py1(ep))
499 m12(3,3,ep)=m12(3,3,ep)+dhz(ep)*(py1py2(ep)-px1px1(ep))
500 m13(3,3,ep)=m13(3,3,ep)+dhz(ep)*py1py3(ep)
501 m22(3,3,ep)=m22(3,3,ep)+dhz(ep)*(px1px1(ep)+py2py2(ep))
502 m23(3,3,ep)=m23(3,3,ep)+dhz(ep)*py2py3(ep)
503 m33(3,3,ep)=m33(3,3,ep)+dhz(ep)*py3py3(ep)
504 ENDDO
505
506 IF (neig==0) THEN
507 DO ep=jft,jlt
508 c1 = em8*m11(3,3,ep)
509 c2 = em8*m22(3,3,ep)
510 c3 = em8*m33(3,3,ep)
511 m11(3,3,ep)=m11(3,3,ep)+c1
512 m22(3,3,ep)=m22(3,3,ep)+c2
513 m33(3,3,ep)=m33(3,3,ep)+c3
514 ENDDO
515 ENDIF
516
517 IF (ikgeo == 1 ) THEN
518
519 DO ep=jft,jlt
520 c2=vol(ep)
524 fxy2=two*fxy
525 h11(ep)=fxx*px1px1(ep)+fyy*py1py1(ep)+fxy2*px1py1(ep)
526 h12(ep)=-fxx*px1px1(ep)+fyy*py1py2(ep)
527 . +fxy*(px1py2(ep)-px1py1(ep))
528 h22(ep)=fxx*px1px1(ep)+fyy*py2py2(ep)-fxy2*px1py2(ep)
529 h13(ep)=fyy*py1py3(ep)+fxy*px1py3(ep)
530 h23(ep)=fyy*py2py3(ep)-fxy*px1py3(ep)
531 h33(ep)=fyy*py3py3(ep)
532 ENDDO
533 DO i=1,3
534 DO ep=jft,jlt
535 k11(i,i,ep) = k11(i,i,ep)+h11(ep)
536 k12(i,i,ep) = k12(i,i,ep)+h12(ep)
537 k22(i,i,ep) = k22(i,i,ep)+h22(ep)
538 k13(i,i,ep) = k13(i,i,ep)+h13(ep)
539 k23(i,i,ep) = k23(i,i,ep)+h23(ep)
540 k33(i,i,ep) = k33(i,i,ep)+h33(ep)
541 ENDDO
542 ENDDO
543
544 IF (neig==0.AND.idril==0.AND.iorth>0) THEN
545 DO ep=jft,jlt
546 c1 =
min(h11(ep),h22(ep),h33(ep))
547 IF (c1 < zero) THEN
548 c2 =
min(m11(3,3,ep),m22(3,3,ep),m33(3,3,ep))
550 m11(3,3,ep)=m11(3,3,ep) + c1
551 m22(3,3,ep)=m22(3,3,ep) + c1
552 m33(3,3,ep)=m33(3,3,ep) + c1
553 END IF
554 ENDDO
555 ENDIF
556 END IF
557
558 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
for(i8=*sizetab-1;i8 >=0;i8--)