38
39
40
41
42#include "implicit_f.inc"
43#include "mvsiz_p.inc"
44#include "com04_c.inc"
45
46
47
48
49 INTEGER JFT,JLT,NPLAT,IPLAT(*),IDRIL,IORTH,NEL
51 . vol(*),thk0(*),thk2(*),
52 . rx(*) ,ry(*) ,sx(*) ,sy(*) ,
53 . rx2(*) ,ry2(*) ,sx2(*) ,sy2(*) ,
for(nel,5),mom(nel,3),
54 . rhx(mvsiz,4) ,rhy(mvsiz,4) ,shx(mvsiz,4) ,shy(mvsiz,4) ,
55 . px1(*) ,px2(*) ,py1(*) ,py2(*) ,
56 . k11(3,3,*),k12(3,3,*),k13(3,3,*),k14(3,3,*),
57 . k22(3,3,*),k23(3,3,*),k24(3,3,*),k33(3,3,*),
58 . m11(3,3,*),m12(3,3,*),m13(3,3,*),m14(3,3,*),
59 . m22(3,3,*),m23(3,3,*),m24(3,3,*),m33(3,3,*),
60 . mf11(3,3,*),mf12(3,3,*),mf13(3,3,*),mf14(3,3,*),
61 . mf22(3,3,*),mf23(3,3,*),mf24(3,3,*),mf33(3,3,*),
62 . fm12(3,3,*),fm13(3,3,*),fm14(3,3,*),
63 . fm23(3,3,*),fm24(3,3,*),fm34(3,3,*),
64 . k34(3,3,*),k44(3,3,*),m34(3,3,*),m44(3,3,*),
65 . mf34(3,3,*),mf44(3,3,*)
66
67
68
69
70
71
72 INTEGER EP,I,J,NF,M
74 . c1,c2,
75 . px1px1(mvsiz),px1px2(mvsiz) ,px2px2(mvsiz),
76 . px1py1(mvsiz),px1py2(mvsiz) ,px2py1(mvsiz),
77 . px2py2(mvsiz),py1py1(mvsiz),py1py2(mvsiz),py2py2(mvsiz)
78 . h11(mvsiz),h12(mvsiz),h22(mvsiz),h13(mvsiz),h23(mvsiz),
79 . fxx(mvsiz),fyy(mvsiz),fxy(mvsiz),mxx(mvsiz),myy(mvsiz),
80 . mxy(mvsiz),fxy2(mvsiz),mxy2(mvsiz),scx(mvsiz),scy(mvsiz),
81 . cs1(mvsiz),cs2(mvsiz),cs3(mvsiz)
82
83 nf=nplat+1
84#include "vectorize.inc"
85
86 DO m=jft,jlt
87 ep=iplat(m)
88 c2=vol(ep)
89 c1=thk2(ep)*c2
93 fxy2(m)=two*fxy(m)
94 mxx(m)=mom(ep,1)*c1
95 myy(m)=mom(ep,2)*c1
96 mxy(m)=mom(ep,3)*c1
97 mxy2(m)=two*mxy(m)
100 ENDDO
101
102 DO ep=jft,jlt
103 px1px1(ep) = px1(ep)*px1(ep)
104 px1px2(ep) = px1(ep)*px2(ep)
105 px2px2(ep) = px2(ep)*px2(ep)
106 px1py1(ep) = px1(ep)*py1(ep)
107 px1py2(ep) = px1(ep)*py2(ep)
108 px2py1(ep) = px2(ep)*py1(ep)
109 px2py2(ep) = px2(ep)*py2(ep)
110 py1py1(ep) = py1(ep)*py1(ep)
111 py1py2(ep) = py1(ep)*py2(ep)
112 py2py2(ep) = py2(ep)*py2(ep)
113 cs1(ep) = scx(ep)*sy2(ep)+scy(ep)*sx2(ep)
114 cs2(ep) = scx(ep)*ry2(ep)+scy(ep)*rx2(ep)
115 cs3(ep) =-scx(ep)*sy(ep)*ry(ep)-scy(ep)*sx(ep)*rx(ep)
116 ENDDO
117
118 DO ep=jft,jlt
119 h11(ep)=fxx(ep)*px1px1(ep)+fyy(ep)*py1py1(ep)
120 . +fxy2(ep)*px1py1(ep)
121 h12(ep)=fxx(ep)*px1px2(ep)+fyy(ep)*py1py2(ep)
122 . +fxy(ep)*(px1py2(ep)+px2py1(ep))
123 h22(ep)=fxx(ep)*px2px2(ep)+fyy(ep)*py2py2(ep)
124 . +fxy2(ep)*px2py2(ep)
125 ENDDO
126 DO i=1,3
127 DO ep=jft,jlt
128 k11(i,i,ep) = k11(i,i,ep)+h11(ep)
129 k12(i,i,ep) = k12(i,i,ep)+h12(ep)
130 k22(i,i,ep) = k22(i,i,ep)+h22(ep)
131 k13(i,i,ep) = k13(i,i,ep)-h11(ep)
132 k14(i,i,ep) = k14(i,i,ep)-h12(ep)
133 k23(i,i,ep) = k23(i,i,ep)-h12(ep)
134 k24(i,i,ep) = k24(i,i,ep)-h22(ep)
135 k33(i,i,ep) = k33(i,i,ep)+h11(ep)
136 k34(i,i,ep) = k34(i,i,ep)+h12(ep)
137 k44(i,i,ep) = k44(i,i,ep)+h22(ep)
138 ENDDO
139 ENDDO
140
141 IF (neig==0.AND.idril==0.AND.iorth >0) THEN
142 DO ep=jft,jlt
143 c1 =
min(h11(ep),h22(ep),-h11(ep),-h22(ep))
144 IF (c1 < zero) THEN
145 c2 =
min(m11(3,3,ep),m22(3,3,ep),m33(3,3,ep),m44(3,3,ep))
147 m11(3,3,ep)=m11(3,3,ep) + c1
148 m22(3,3,ep)=m22(3,3,ep) + c1
149 m33(3,3,ep)=m33(3,3,ep) + c1
150 m44(3,3,ep)=m44(3,3,ep) + c1
151 END IF
152 ENDDO
153 ENDIF
154
155 RETURN
156
157 DO ep=jft,jlt
158 h11(ep)=mxx(ep)*py1py1(ep)+myy(ep)*px1px1(ep)
159 . -mxy2(ep)*px1py1(ep)
160 h12(ep)=mxx(ep)*py1py2(ep)+myy(ep)*px1px2(ep)
161 . -mxy(ep)*(px2py1(ep)+px1py2(ep))
162 h22(ep)=mxx(ep)*py2py2(ep)+myy(ep)*px2px2(ep)
163 . -mxy2(ep)*px2py2(ep)
164 ENDDO
165 DO i=1,2
166 DO ep=jft,jlt
167 m11(i,i,ep) = m11(i,i,ep)+h11(ep)
168 m12(i,i,ep) = m12(i,i,ep)+h12(ep)
169 m22(i,i,ep) = m22(i,i,ep)+h22(ep)
170 m13(i,i,ep) = m13(i,i,ep)-h11(ep)
171 m14(i,i,ep) = m14(i,i,ep)-h12(ep)
172 m23(i,i,ep) = m23(i,i,ep)-h12(ep)
173 m24(i,i,ep) = m24(i,i,ep)-h22(ep)
174 m33(i,i,ep) = m33(i,i,ep)+h11(ep)
175 m34(i,i,ep) = m34(i,i,ep)+h12(ep)
176 m44(i,i,ep) = m44(i,i,ep)+h22(ep)
177 ENDDO
178 ENDDO
179
180 DO ep=jft,jlt
181 h11(ep)=cs1(ep)+cs2(ep)+cs3(ep)+cs3(ep)
182 h22(ep)=cs1(ep)+cs2(ep)-cs3(ep)-cs3(ep)
183 h12(ep)=-cs1(ep)+cs2(ep)+cs3(ep)-cs3(ep)
184 h13(ep)=-cs1(ep)-cs2(ep)-cs3(ep)-cs3(ep)
185 h23(ep)= cs1(ep)-cs2(ep)+cs3(ep)-cs3(ep)
186 ENDDO
187 DO i=1,3
188 DO ep=jft,jlt
189 k11(i,i,ep) = k11(i,i,ep)+h11(ep)
190 k12(i,i,ep) = k12(i,i,ep)+h12(ep)
191 k22(i,i,ep) = k22(i,i,ep)+h22(ep)
192 k13(i,i,ep) = k13(i,i,ep)+h13(ep)
193 k14(i,i,ep) = k14(i,i,ep)-h12(ep)
194 k23(i,i,ep) = k23(i,i,ep)+h23(ep)
195 k24(i,i,ep) = k24(i,i,ep)-h22(ep)
196 k33(i,i,ep) = k33(i,i,ep)+h11(ep)
197 k34(i,i,ep) = k34(i,i,ep)+h12(ep)
198 k44(i,i,ep) = k44(i,i,ep)+h22(ep)
199 ENDDO
200 ENDDO
201
202
203 DO ep=jft,jlt
204
205 h11(ep)=(rhy(ep,1)*rhy(ep,1)+rhx(ep,1)*rhx(ep,1))*cs1(ep)
206
207 h12(ep)=(rhy(ep,1)*rhy(ep,3)+rhx(ep,1)*rhx(ep,3))*cs1(ep)
208
209 h13(ep)=(rhy(ep,3)*rhy(ep,3)+rhx(ep,3)*rhx(ep,3))*cs1(ep)
210 ENDDO
211 DO i=1,2
212 DO ep=jft,jlt
213 m11(i,i,ep) = m11(i,i,ep)+h11(ep)
214 m12(i,i,ep) = m12(i,i,ep)+h11(ep)
215 m22(i,i,ep) = m22(i,i,ep)+h11(ep)
216 m13(i,i,ep) = m13(i,i,ep)+h12(ep)
217 m14(i,i,ep) = m14(i,i,ep)+h12(ep)
218 m23(i,i,ep) = m23(i,i,ep)+h12(ep)
219 m24(i,i,ep) = m24(i,i,ep)+h12(ep)
220 m33(i,i,ep) = m33(i,i,ep)+h13(ep)
221 m34(i,i,ep) = m34(i,i,ep)+h13(ep)
222 m44(i,i,ep) = m44(i,i,ep)+h13(ep)
223 ENDDO
224 ENDDO
225 DO ep=jft,jlt
226
227
228 h11(ep)=(shy(ep,1)*shy(ep,1)+shx(ep,1)*shx(ep,1))*cs2(ep)
229
230 h12(ep)=(shy(ep,1)*shy(ep,2)+shx(ep,1)*shx(ep,2))*cs2(ep)
231
232 h22(ep)=(shy(ep,2)*shy(ep,2)+shx(ep,2)*shx(ep,2))*cs2(ep)
233 ENDDO
234 DO i=1,2
235 DO ep=jft,jlt
236 m11(i,i,ep) = m11(i,i,ep)+h11(ep)
237 m12(i,i,ep) = m12(i,i,ep)+h12(ep)
238 m22(i,i,ep) = m22(i,i,ep)+h22(ep)
239 m13(i,i,ep) = m13(i,i,ep)+h12(ep)
240 m14(i,i,ep) = m14(i,i,ep)+h11(ep)
241 m23(i,i,ep) = m23(i,i,ep)+h22(ep)
242 m24(i,i,ep) = m24(i,i,ep)+h12(ep)
243 m33(i,i,ep) = m33(i,i,ep)+h22(ep)
244 m34(i,i,ep) = m34(i,i,ep)+h12(ep)
245 m44(i,i,ep) = m44(i,i,ep)+h11(ep)
246 ENDDO
247 ENDDO
248
249 DO ep=jft,jlt
250 h11(ep)=(rhy(ep,1)*shy(ep,1)+rhy(ep,1)*shy(ep,1)+
251 . rhx(ep,1)*shx(ep,1)+rhx(ep,1)*shx(ep,1))*cs3(ep)
252 h12(ep)=(rhy(ep,1)*shy(ep,2)+rhy(ep,2)*shy(ep,1)+
253 . rhx(ep,1)*shx(ep,2)+rhx(ep,2)*shx(ep,1))*cs3(ep)
254 h22(ep)=(rhy(ep,2)*shy(ep,2)+rhy(ep,2)*shy(ep,2)+
255 . rhx(ep,2)*shx(ep,2)+rhx(ep,2)*shx(ep,2))*cs3(ep)
256 h13(ep)=(rhy(ep,1)*shy(ep,3)+rhy(ep,3)*shy(ep,1)+
257 . rhx(ep,1)*shx(ep,3)+rhx(ep,3)*shx(ep,1))*cs3(ep)
258 h23(ep)=(rhy(ep,2)*shy(ep,3)+rhy(ep,3)*shy(ep,2)+
259 . rhx(ep,2)*shx(ep,3)+rhx(ep,3)*shx(ep,2))*cs3(ep)
260 ENDDO
261 DO i=1,2
262 DO ep=jft,jlt
263 m11(i,i,ep) = m11(i,i,ep)+h11(ep)
264 m12(i,i,ep) = m12(i,i,ep)+h12(ep)
265 m22(i,i,ep) = m22(i,i,ep)+h22(ep)
266 m13(i,i,ep) = m13(i,i,ep)+h12(ep)
267 m23(i,i,ep) = m23(i,i,ep)+h22(ep)
268 ENDDO
269 ENDDO
270
271 DO ep=jft,jlt
272 h11(ep)=(rhy(ep,1)*shy(ep,4)+rhy(ep,4)*shy(ep,1)+
273 . rhx(ep,1)*shx(ep,4)+rhx(ep,4)*shx(ep,1))*cs3(ep)
274 h12(ep)=(rhy(ep,2)*shy(ep,4)+rhy(ep,4)*shy(ep,2)+
275 . rhx(ep,2)*shx(ep,4)+rhx(ep,4)*shx(ep,2))*cs3(ep)
276 h22(ep)=(rhy(ep,3)*shy(ep,3)+rhy(ep,3)*shy(ep,3)+
277 . rhx(ep,3)*shx(ep,3)+rhx(ep,3)*shx(ep,3))*cs3(ep)
278 h13(ep)=(rhy(ep,3)*shy(ep,4)+rhy(ep,4)*shy(ep,3)+
279 . rhx(ep,3)*shx(ep,4)+rhx(ep,4)*shx(ep,3))*cs3(ep)
280 h23(ep)=(rhy(ep,4)*shy(ep,4)+rhy(ep,4)*shy(ep,4)+
281 . rhx(ep,4)*shx(ep,4)+rhx(ep,4)*shx(ep,4))*cs3(ep)
282 ENDDO
283 DO i=1,2
284 DO ep=jft,jlt
285 m14(i,i,ep) = m14
286 m24(i,i,ep) = m24(i,i,ep)+h12(ep)
287 m33(i,i,ep) = m33(i,i,ep)+h22(ep)
288 m34(i,i,ep) = m34(i,i,ep)+h12(ep)
289 m44(i,i,ep) = m44(i,i,ep)+h22(ep)
290 ENDDO
291 ENDDO
292
293 RETURN
for(i8=*sizetab-1;i8 >=0;i8--)