OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cbalke3.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com04_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine cbalke3 (jft, jlt, cdet, thk0, thk2, hm, hf, hc, hz, bm, bmf, bf, bc, tc, bzz, nplat, iplat, vol, ikgeo, for, mom, k11, k12, k13, k14, k22, k23, k24, k33, k34, k44, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44, mf11, mf12, mf13, mf14, mf22, mf23, mf24, mf33, mf34, mf44, fm12, fm13, fm14, fm23, fm24, fm34, iorth, hmor, hfor, idril, hmfor, x13, x24, y13, y24, nel)
subroutine cbalkec3 (jft, jlt, vol, x13, x24, y13, y24, hm, k11, k12, k13, k14, k22, k23, k24, k33, k34, k44, nplat, iplat, ikgeo, for, m11, m22, m33, m44, iorth, nel)
subroutine cbalkerz (jft, jlt, vol, thk0, hm, hz, bm, k11, k12, k13, k14, k22, k23, k24, k33, k34, k44, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44, mf11, mf12, mf13, mf14, mf22, mf23, mf24, mf33, mf34, mf44, fm12, fm13, fm14, fm23, fm24, fm34, iorth, hmor, hfor, iplat, nplat, pmrz, brz, frz, ikgeo, ng, hmfor, bf, bmf, nel)
subroutine cbalkeormf (jft, jlt, kmfij, dmf, tij)
subroutine cbalkeorfm (jft, jlt, kfmij, dmf, tij)
subroutine cbalkeormf1 (jft, jlt, kmfij, dmf, tij, t0ij)
subroutine cbalkeorfm1 (jft, jlt, kfmij, dmf, tij, t0ij)

Function/Subroutine Documentation

◆ cbalke3()

subroutine cbalke3 ( integer jft,
integer jlt,
cdet,
thk0,
thk2,
hm,
hf,
hc,
hz,
bm,
bmf,
bf,
bc,
tc,
bzz,
integer nplat,
integer, dimension(*) iplat,
vol,
integer ikgeo,
for,
mom,
k11,
k12,
k13,
k14,
k22,
k23,
k24,
k33,
k34,
k44,
m11,
m12,
m13,
m14,
m22,
m23,
m24,
m33,
m34,
m44,
mf11,
mf12,
mf13,
mf14,
mf22,
mf23,
mf24,
mf33,
mf34,
mf44,
fm12,
fm13,
fm14,
fm23,
fm24,
fm34,
integer iorth,
hmor,
hfor,
integer idril,
hmfor,
x13,
x24,
y13,
y24,
integer nel )

Definition at line 33 of file cbalke3.F.

42C--------------------------------------------------------------------------------------------------
43C-----------------------------------------------
44C I M P L I C I T T Y P E S
45C-----------------------------------------------
46#include "implicit_f.inc"
47#include "mvsiz_p.inc"
48#include "com04_c.inc"
49C-----------------------------------------------
50C D U M M Y A R G U M E N T S
51C-----------------------------------------------
52C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
53 INTEGER JFT,JLT,NPLAT,IPLAT(*),IKGEO,IORTH,IDRIL,NEL
54 my_real
55 . cdet(*),vol(*),
56 . bm(mvsiz,3,3,4),bmf(mvsiz,3,3,4),bf(mvsiz,3,2,4),bc(mvsiz,2,5,4),
57 . thk0(*),thk2(*),bzz(mvsiz,8),tc(mvsiz,2,2),
58 . hm(mvsiz,4),hf(mvsiz,4),hc(mvsiz,2),hz(*),
59 . k11(3,3,*),k12(3,3,*),k13(3,3,*),k14(3,3,*),
60 . k22(3,3,*),k23(3,3,*),k24(3,3,*),k33(3,3,*),
61 . m11(3,3,*),m12(3,3,*),m13(3,3,*),m14(3,3,*),
62 . m22(3,3,*),m23(3,3,*),m24(3,3,*),m33(3,3,*),
63 . mf11(3,3,*),mf12(3,3,*),mf13(3,3,*),mf14(3,3,*),
64 . mf22(3,3,*),mf23(3,3,*),mf24(3,3,*),mf33(3,3,*),
65 . fm12(3,3,*),fm13(3,3,*),fm14(3,3,*),
66 . fm23(3,3,*),fm24(3,3,*),fm34(3,3,*),hmor(mvsiz,2),hfor(mvsiz,2),
67 . k34(3,3,*),k44(3,3,*),m34(3,3,*),m44(3,3,*),
68 . mf34(3,3,*),mf44(3,3,*),for(nel,5),mom(nel,3),hmfor(mvsiz,6),
69 . x13(*),x24(*),y13(*),y24(*)
70C-----------------------------------------------
71C L O C A L V A R I A B L E S
72C-----------------------------------------------
73 INTEGER EP,I,J,L,K,M,I1,J1,IN(2),NF
74 my_real
75 . dm(2,2,mvsiz),df(2,2,mvsiz),dc(2,2,mvsiz),
76 . t11(2,2,mvsiz),t12(2,2,mvsiz),t13(2,2,mvsiz),t14(2,2,mvsiz),
77 . t22(2,2,mvsiz),t23(2,2,mvsiz),t24(2,2,mvsiz),t33(2,2,mvsiz),
78 . t44(2,2,mvsiz),t34(2,2,mvsiz),
79 . s11(2,2,mvsiz),s12(2,2,mvsiz),s13(2,2,mvsiz),s14(2,2,mvsiz),
80 . s22(2,2,mvsiz),s23(2,2,mvsiz),s24(2,2,mvsiz),s33(2,2,mvsiz),
81 . s44(2,2,mvsiz),s34(2,2,mvsiz),c1,c2,dhz(mvsiz),c3,c4,
82 . bb(2,4,mvsiz),bbc(2,3,4,mvsiz),gm(mvsiz),gf(mvsiz),
83 . fxx(mvsiz),fyy(mvsiz),dm5(mvsiz),dm6(mvsiz),df5(mvsiz),
84 . fxy(mvsiz),fxy2(mvsiz),km12,km21,
85 . df6(mvsiz),dm5_2(mvsiz),dm6_2(mvsiz),df5_2(mvsiz),
86 . df6_2(mvsiz),dmf(3,3,mvsiz),fac2z,cbr1(4,mvsiz),
87 . cbr2(4,mvsiz),cbr3(4,mvsiz),bb0(2,2,mvsiz),t0ij(2,2,mvsiz),
88 . t1ij(2,2,mvsiz)
89 DATA in/2,1/
90C-----------Attention Matrice sym Kii ne calcul que la moitie---------72
91#include "vectorize.inc"
92 DO m=jft,jlt
93 ep=iplat(m)
94 c2=vol(ep)
95 c1=thk2(ep)*c2
96 dm(1,1,m)=hm(ep,1)*c2
97 dm(2,2,m)=hm(ep,2)*c2
98 dm(1,2,m)=hm(ep,3)*c2
99 dm(2,1,m)=dm(1,2,m)
100 gm(m) =hm(ep,4)*c2
101 df(1,1,m)=hf(ep,1)*c1
102 df(2,2,m)=hf(ep,2)*c1
103 df(1,2,m)=hf(ep,3)*c1
104 df(2,1,m)=df(1,2,m)
105 gf(m) =hf(ep,4)*c1
106 dhz(m)=hz(ep)*c1
107 ENDDO
108#include "vectorize.inc"
109 DO m=jft,nplat
110 ep=iplat(m)
111 c2=vol(ep)
112 dc(1,1,m)=hc(ep,1)*c2
113 dc(2,2,m)=hc(ep,2)*c2
114 ENDDO
115C
116 nf =nplat+1
117C
118 DO ep=jft,nplat
119 bb(1,1,ep)=bm(ep,1,1,1)
120 bb(1,2,ep)=bm(ep,2,1,1)
121 bb(1,3,ep)=bm(ep,3,1,1)
122 bb(1,4,ep)=bm(ep,1,2,1)
123 bb(2,1,ep)=bm(ep,2,2,1)
124 bb(2,2,ep)=bm(ep,3,2,1)
125 bb(2,3,ep)=bm(ep,1,3,1)
126 bb(2,4,ep)=bm(ep,2,3,1)
127 ENDDO
128C
129 DO i=1,2
130 DO j=i,2
131 DO ep=jft,nplat
132 t11(i,j,ep)=bb(i,1,ep)*bb(j,1,ep)
133 t11(j,i,ep)=t11(i,j,ep)
134 t22(i,j,ep)=bb(i,2,ep)*bb(j,2,ep)
135 t22(j,i,ep)=t22(i,j,ep)
136 t33(i,j,ep)=bb(i,3,ep)*bb(j,3,ep)
137 t33(j,i,ep)=t33(i,j,ep)
138 t44(i,j,ep)=bb(i,4,ep)*bb(j,4,ep)
139 t44(j,i,ep)=t44(i,j,ep)
140 s11(in(i),in(j),ep)=t11(i,j,ep)
141 s22(in(i),in(j),ep)=t22(i,j,ep)
142 s33(in(i),in(j),ep)=t33(i,j,ep)
143 s44(in(i),in(j),ep)=t44(i,j,ep)
144 ENDDO
145 ENDDO
146 ENDDO
147C
148 DO i=1,2
149 DO j=1,2
150 DO ep=jft,nplat
151 t12(i,j,ep)=bb(i,1,ep)*bb(j,2,ep)
152 t13(i,j,ep)=bb(i,1,ep)*bb(j,3,ep)
153 t14(i,j,ep)=bb(i,1,ep)*bb(j,4,ep)
154 t23(i,j,ep)=bb(i,2,ep)*bb(j,3,ep)
155 t24(i,j,ep)=bb(i,2,ep)*bb(j,4,ep)
156 t34(i,j,ep)=bb(i,3,ep)*bb(j,4,ep)
157 s12(in(i),in(j),ep)=t12(i,j,ep)
158 s13(in(i),in(j),ep)=t13(i,j,ep)
159 s14(in(i),in(j),ep)=t14(i,j,ep)
160 s23(in(i),in(j),ep)=t23(i,j,ep)
161 s24(in(i),in(j),ep)=t24(i,j,ep)
162 s34(in(i),in(j),ep)=t34(i,j,ep)
163 ENDDO
164 ENDDO
165 ENDDO
166 DO ep=jft,nplat
167 s11(2,1,ep)=-s11(2,1,ep)
168 s11(1,2,ep)= s11(2,1,ep)
169 s22(2,1,ep)=-s22(2,1,ep)
170 s22(1,2,ep)= s22(2,1,ep)
171 s33(2,1,ep)=-s33(2,1,ep)
172 s33(1,2,ep)= s33(2,1,ep)
173 s44(2,1,ep)=-s44(2,1,ep)
174 s44(1,2,ep)= s44(2,1,ep)
175 s12(1,2,ep)=-s12(1,2,ep)
176 s12(2,1,ep)=-s12(2,1,ep)
177 s13(1,2,ep)=-s13(1,2,ep)
178 s13(2,1,ep)=-s13(2,1,ep)
179 s14(1,2,ep)=-s14(1,2,ep)
180 s14(2,1,ep)=-s14(2,1,ep)
181 s23(1,2,ep)=-s23(1,2,ep)
182 s23(2,1,ep)=-s23(2,1,ep)
183 s24(1,2,ep)=-s24(1,2,ep)
184 s24(2,1,ep)=-s24(2,1,ep)
185 s34(1,2,ep)=-s34(1,2,ep)
186 s34(2,1,ep)=-s34(2,1,ep)
187 ENDDO
188C
189 DO i=1,2
190 DO j=i,2
191 DO ep=jft,nplat
192 k11(i,j,ep)=k11(i,j,ep)+t11(i,j,ep)*dm(i,j,ep)
193 k22(i,j,ep)=k22(i,j,ep)+t22(i,j,ep)*dm(i,j,ep)
194 k33(i,j,ep)=k33(i,j,ep)+t33(i,j,ep)*dm(i,j,ep)
195 k44(i,j,ep)=k44(i,j,ep)+t44(i,j,ep)*dm(i,j,ep)
196 ENDDO
197 ENDDO
198 ENDDO
199 DO i=1,2
200 i1 = in(i)
201 DO j=i,2
202 j1 = in(j)
203#include "vectorize.inc"
204 DO ep=jft,nplat
205 m11(i,j,ep)=m11(i,j,ep)+s11(i,j,ep)*df(i1,j1,ep)+
206 1 s11(i1,j1,ep)*gf(ep)
207 m22(i,j,ep)=m22(i,j,ep)+s22(i,j,ep)*df(i1,j1,ep)+
208 1 s22(i1,j1,ep)*gf(ep)
209 m33(i,j,ep)=m33(i,j,ep)+s33(i,j,ep)*df(i1,j1,ep)+
210 1 s33(i1,j1,ep)*gf(ep)
211 m44(i,j,ep)=m44(i,j,ep)+s44(i,j,ep)*df(i1,j1,ep)+
212 1 s44(i1,j1,ep)*gf(ep)
213 ENDDO
214 ENDDO
215 ENDDO
216C
217 DO i=1,2
218 DO j=1,2
219 DO ep=jft,nplat
220 k12(i,j,ep)=k12(i,j,ep)+t12(i,j,ep)*dm(i,j,ep)
221 k13(i,j,ep)=k13(i,j,ep)+t13(i,j,ep)*dm(i,j,ep)
222 k14(i,j,ep)=k14(i,j,ep)+t14(i,j,ep)*dm(i,j,ep)
223 k23(i,j,ep)=k23(i,j,ep)+t23(i,j,ep)*dm(i,j,ep)
224 k24(i,j,ep)=k24(i,j,ep)+t24(i,j,ep)*dm(i,j,ep)
225 k34(i,j,ep)=k34(i,j,ep)+t34(i,j,ep)*dm(i,j,ep)
226 ENDDO
227 ENDDO
228 ENDDO
229 DO i=1,2
230 i1 = in(i)
231 DO j=1,2
232 j1 = in(j)
233#include "vectorize.inc"
234 DO ep=jft,nplat
235 m12(i,j,ep)=m12(i,j,ep)+s12(i,j,ep)*df(i1,j1,ep)+
236 1 s12(i1,j1,ep)*gf(ep)
237 m13(i,j,ep)=m13(i,j,ep)+s13(i,j,ep)*df(i1,j1,ep)+
238 1 s13(i1,j1,ep)*gf(ep)
239 m14(i,j,ep)=m14(i,j,ep)+s14(i,j,ep)*df(i1,j1,ep)+
240 1 s14(i1,j1,ep)*gf(ep)
241 m23(i,j,ep)=m23(i,j,ep)+s23(i,j,ep)*df(i1,j1,ep)+
242 1 s23(i1,j1,ep)*gf(ep)
243 m24(i,j,ep)=m24(i,j,ep)+s24(i,j,ep)*df(i1,j1,ep)+
244 1 s24(i1,j1,ep)*gf(ep)
245 m34(i,j,ep)=m34(i,j,ep)+s34(i,j,ep)*df(i1,j1,ep)+
246 1 s34(i1,j1,ep)*gf(ep)
247 ENDDO
248 ENDDO
249 ENDDO
250C ----add shear part with drilling dof-----
251 IF (idril>0) THEN
252 DO i=1,2
253 i1 = in(i)
254 DO j=i,2
255 j1 = in(j)
256#include "vectorize.inc"
257 DO ep=jft,nplat
258 k11(i,j,ep)=k11(i,j,ep)+t11(i1,j1,ep)*gm(ep)
259 k22(i,j,ep)=k22(i,j,ep)+t22(i1,j1,ep)*gm(ep)
260 k33(i,j,ep)=k33(i,j,ep)+t33(i1,j1,ep)*gm(ep)
261 k44(i,j,ep)=k44(i,j,ep)+t44(i1,j1,ep)*gm(ep)
262 ENDDO
263 ENDDO
264 ENDDO
265C
266 DO i=1,2
267 i1 = in(i)
268 DO j=1,2
269 j1 = in(j)
270#include "vectorize.inc"
271 DO ep=jft,nplat
272 k12(i,j,ep)=k12(i,j,ep)+t12(i1,j1,ep)*gm(ep)
273 k13(i,j,ep)=k13(i,j,ep)+t13(i1,j1,ep)*gm(ep)
274 k14(i,j,ep)=k14(i,j,ep)+t14(i1,j1,ep)*gm(ep)
275 k23(i,j,ep)=k23(i,j,ep)+t23(i1,j1,ep)*gm(ep)
276 k24(i,j,ep)=k24(i,j,ep)+t24(i1,j1,ep)*gm(ep)
277 k34(i,j,ep)=k34(i,j,ep)+t34(i1,j1,ep)*gm(ep)
278 ENDDO
279 ENDDO
280 ENDDO
281C---+----1----+----2----warped-+----5----+----6----+----7----+----8
282 DO i=1,3
283 DO j=i,3
284 DO ep=nf,jlt
285 k11(i,j,ep)=k11(i,j,ep)+bm(ep,3,i,1)*gm(ep)*bm(ep,3,j,1)
286 k22(i,j,ep)=k22(i,j,ep)+bm(ep,3,i,2)*gm(ep)*bm(ep,3,j,2)
287 k33(i,j,ep)=k33(i,j,ep)+bm(ep,3,i,3)*gm(ep)*bm(ep,3,j,3)
288 k44(i,j,ep)=k44(i,j,ep)+bm(ep,3,i,4)*gm(ep)*bm(ep,3,j,4)
289 ENDDO
290 ENDDO
291 ENDDO
292C
293 DO i=1,3
294 DO j=1,3
295 DO ep=nf,jlt
296 k12(i,j,ep)=k12(i,j,ep)+bm(ep,3,i,1)*gm(ep)*bm(ep,3,j,2)
297 k13(i,j,ep)=k13(i,j,ep)+bm(ep,3,i,1)*gm(ep)*bm(ep,3,j,3)
298 k14(i,j,ep)=k14(i,j,ep)+bm(ep,3,i,1)*gm(ep)*bm(ep,3,j,4)
299 k23(i,j,ep)=k23(i,j,ep)+bm(ep,3,i,2)*gm(ep)*bm(ep,3,j,3)
300 k24(i,j,ep)=k24(i,j,ep)+bm(ep,3,i,2)*gm(ep)*bm(ep,3,j,4)
301 k34(i,j,ep)=k34(i,j,ep)+bm(ep,3,i,3)*gm(ep)*bm(ep,3,j,4)
302 ENDDO
303 ENDDO
304 ENDDO
305 END IF
306C ----ajoute termes sup. due au orthotrope-----
307 IF (iorth>0) THEN
308C ----add mem-bending coupling terms of orthotrope--coplanar first---
309#include "vectorize.inc"
310 DO m=jft,jlt
311 ep=iplat(m)
312 c2=vol(ep)*thk0(ep)
313 dmf(1,1,m)=hmfor(ep,1)*c2
314 dmf(2,2,m)=hmfor(ep,2)*c2
315 dmf(1,2,m)=hmfor(ep,3)*c2
316 dmf(1,3,m)=hmfor(ep,5)*c2
317 dmf(2,3,m)=hmfor(ep,6)*c2
318 dmf(2,1,m)=dmf(1,2,m)
319 dmf(3,1,m)=dmf(1,3,m)
320 dmf(3,2,m)=dmf(2,3,m)
321 dmf(3,3,m)=hmfor(ep,4)*c2
322 ENDDO
323C------add here before TIJ(1,2,*)... have been modified --TIJ(i,j)=B(i,I)*B(j,J)----i,j->1:x;2:y
324 IF (idril >0 ) THEN
325 CALL cbalkeormf(jft ,nplat ,mf11 ,dmf ,t11 )
326 CALL cbalkeormf(jft ,nplat ,mf12 ,dmf ,t12 )
327 CALL cbalkeormf(jft ,nplat ,mf13 ,dmf ,t13 )
328 CALL cbalkeormf(jft ,nplat ,mf14 ,dmf ,t14 )
329 CALL cbalkeormf(jft ,nplat ,mf22 ,dmf ,t22 )
330 CALL cbalkeormf(jft ,nplat ,mf23 ,dmf ,t23 )
331 CALL cbalkeormf(jft ,nplat ,mf24 ,dmf ,t24 )
332 CALL cbalkeormf(jft ,nplat ,mf33 ,dmf ,t33 )
333 CALL cbalkeormf(jft ,nplat ,mf34 ,dmf ,t34 )
334 CALL cbalkeormf(jft ,nplat ,mf44 ,dmf ,t44 )
335C
336 CALL cbalkeorfm(jft ,nplat ,fm12 ,dmf ,t12 )
337 CALL cbalkeorfm(jft ,nplat ,fm13 ,dmf ,t13 )
338 CALL cbalkeorfm(jft ,nplat ,fm14 ,dmf ,t14 )
339 CALL cbalkeorfm(jft ,nplat ,fm23 ,dmf ,t23 )
340 CALL cbalkeorfm(jft ,nplat ,fm24 ,dmf ,t24 )
341 CALL cbalkeorfm(jft ,nplat ,fm34 ,dmf ,t34 )
342 ELSE
343C-----------due to constant in-plane shear --------------
344#include "vectorize.inc"
345 DO m=jft,nplat
346 ep=iplat(m)
347 bb0(1,1,m)=y24(ep)
348 bb0(1,2,m)=-y13(ep)
349 bb0(2,1,m)=-x24(ep)
350 bb0(2,2,m)=x13(ep)
351 ENDDO
352C
353 DO ep=jft,nplat
354 t0ij(1,1,ep)=bb0(1,1,ep)*bb(1,1,ep)
355 t0ij(1,2,ep)=bb0(1,1,ep)*bb(2,1,ep)
356 t0ij(2,1,ep)=bb0(2,1,ep)*bb(1,1,ep)
357 t0ij(2,2,ep)=bb0(2,1,ep)*bb(2,1,ep)
358 t1ij(1,1,ep)=-t0ij(1,1,ep)
359 t1ij(1,2,ep)=-t0ij(2,1,ep)
360 t1ij(2,1,ep)=-t0ij(1,2,ep)
361 t1ij(2,2,ep)=-t0ij(2,2,ep)
362 ENDDO
363 CALL cbalkeormf1(jft ,nplat ,mf11 ,dmf ,t11 ,t0ij )
364 CALL cbalkeorfm1(jft ,nplat ,fm13 ,dmf ,t13 ,t1ij )
365 DO ep=jft,nplat
366 t0ij(1,1,ep)=bb0(1,1,ep)*bb(1,2,ep)
367 t0ij(1,2,ep)=bb0(1,1,ep)*bb(2,2,ep)
368 t0ij(2,1,ep)=bb0(2,1,ep)*bb(1,2,ep)
369 t0ij(2,2,ep)=bb0(2,1,ep)*bb(2,2,ep)
370 t1ij(1,1,ep)=-t0ij(1,1,ep)
371 t1ij(1,2,ep)=-t0ij(2,1,ep)
372 t1ij(2,1,ep)=-t0ij(1,2,ep)
373 t1ij(2,2,ep)=-t0ij(2,2,ep)
374 ENDDO
375 CALL cbalkeormf1(jft ,nplat ,mf12 ,dmf ,t12 ,t0ij )
376 CALL cbalkeorfm1(jft ,nplat ,fm23 ,dmf ,t23 ,t1ij )
377 DO ep=jft,nplat
378 t0ij(1,1,ep)=bb0(1,1,ep)*bb(1,3,ep)
379 t0ij(1,2,ep)=bb0(1,1,ep)*bb(2,3,ep)
380 t0ij(2,1,ep)=bb0(2,1,ep)*bb(1,3,ep)
381 t0ij(2,2,ep)=bb0(2,1,ep)*bb(2,3,ep)
382 t1ij(1,1,ep)=-t0ij(1,1,ep)
383 t1ij(1,2,ep)=-t0ij(1,2,ep)
384 t1ij(2,1,ep)=-t0ij(2,1,ep)
385 t1ij(2,2,ep)=-t0ij(2,2,ep)
386 ENDDO
387 CALL cbalkeormf1(jft ,nplat ,mf13 ,dmf ,t13 ,t0ij )
388 CALL cbalkeormf1(jft ,nplat ,mf33 ,dmf ,t33 ,t1ij )
389 DO ep=jft,nplat
390 t0ij(1,1,ep)=bb0(1,1,ep)*bb(1,4,ep)
391 t0ij(1,2,ep)=bb0(1,1,ep)*bb(2,4,ep)
392 t0ij(2,1,ep)=bb0(2,1,ep)*bb(1,4,ep)
393 t0ij(2,2,ep)=bb0(2,1,ep)*bb(2,4,ep)
394 t1ij(1,1,ep)=-t0ij(1,1,ep)
395 t1ij(1,2,ep)=-t0ij(1,2,ep)
396 t1ij(2,1,ep)=-t0ij(2,1,ep)
397 t1ij(2,2,ep)=-t0ij(2,2,ep)
398 ENDDO
399 CALL cbalkeormf1(jft ,nplat ,mf14 ,dmf ,t14 ,t0ij )
400 CALL cbalkeormf1(jft ,nplat ,mf34 ,dmf ,t34 ,t1ij )
401 DO ep=jft,nplat
402 t0ij(1,1,ep)=bb0(1,2,ep)*bb(1,2,ep)
403 t0ij(1,2,ep)=bb0(1,2,ep)*bb(2,2,ep)
404 t0ij(2,1,ep)=bb0(2,2,ep)*bb(1,2,ep)
405 t0ij(2,2,ep)=bb0(2,2,ep)*bb(2,2,ep)
406 t1ij(1,1,ep)=-t0ij(1,1,ep)
407 t1ij(1,2,ep)=-t0ij(2,1,ep)
408 t1ij(2,1,ep)=-t0ij(1,2,ep)
409 t1ij(2,2,ep)=-t0ij(2,2,ep)
410 ENDDO
411 CALL cbalkeormf1(jft ,nplat ,mf22 ,dmf ,t22 ,t0ij )
412 CALL cbalkeorfm1(jft ,nplat ,fm24 ,dmf ,t24 ,t1ij )
413 DO ep=jft,nplat
414 t0ij(1,1,ep)=bb0(1,2,ep)*bb(1,3,ep)
415 t0ij(1,2,ep)=bb0(1,2,ep)*bb(2,3,ep)
416 t0ij(2,1,ep)=bb0(2,2,ep)*bb(1,3,ep)
417 t0ij(2,2,ep)=bb0(2,2,ep)*bb(2,3,ep)
418 t1ij(1,1,ep)=-t0ij(1,1,ep)
419 t1ij(1,2,ep)=-t0ij(2,1,ep)
420 t1ij(2,1,ep)=-t0ij(1,2,ep)
421 t1ij(2,2,ep)=-t0ij(2,2,ep)
422 ENDDO
423 CALL cbalkeormf1(jft ,nplat ,mf23 ,dmf ,t23 ,t0ij )
424 CALL cbalkeorfm1(jft ,nplat ,fm34 ,dmf ,t34 ,t1ij )
425 DO ep=jft,nplat
426 t0ij(1,1,ep)=bb0(1,2,ep)*bb(1,4,ep)
427 t0ij(1,2,ep)=bb0(1,2,ep)*bb(2,4,ep)
428 t0ij(2,1,ep)=bb0(2,2,ep)*bb(1,4,ep)
429 t0ij(2,2,ep)=bb0(2,2,ep)*bb(2,4,ep)
430 t1ij(1,1,ep)=-t0ij(1,1,ep)
431 t1ij(1,2,ep)=-t0ij(1,2,ep)
432 t1ij(2,1,ep)=-t0ij(2,1,ep)
433 t1ij(2,2,ep)=-t0ij(2,2,ep)
434 ENDDO
435 CALL cbalkeormf1(jft ,nplat ,mf24 ,dmf ,t24 ,t0ij )
436 CALL cbalkeormf1(jft ,nplat ,mf44 ,dmf ,t44 ,t1ij )
437C
438 DO ep=jft,nplat
439 t0ij(1,1,ep)=bb(1,1,ep)*bb0(1,2,ep)
440 t0ij(1,2,ep)=bb(1,1,ep)*bb0(2,2,ep)
441 t0ij(2,1,ep)=bb(2,1,ep)*bb0(1,2,ep)
442 t0ij(2,2,ep)=bb(2,1,ep)*bb0(2,2,ep)
443 t1ij(1,1,ep)=-t0ij(1,1,ep)
444 t1ij(1,2,ep)=-t0ij(1,2,ep)
445 t1ij(2,1,ep)=-t0ij(2,1,ep)
446 t1ij(2,2,ep)=-t0ij(2,2,ep)
447 ENDDO
448 CALL cbalkeorfm1(jft ,nplat ,fm12 ,dmf ,t12 ,t0ij )
449 CALL cbalkeorfm1(jft ,nplat ,fm14 ,dmf ,t14 ,t1ij )
450 END IF !(IDRIL >0 ) THEN
451C
452#include "vectorize.inc"
453 DO m=jft,jlt
454 ep=iplat(m)
455 c2=vol(ep)
456 c1=thk2(ep)*c2
457 dm5(m)=hmor(ep,1)*c2
458 dm6(m)=hmor(ep,2)*c2
459 df5(m)=hfor(ep,1)*c1
460 df6(m)=hfor(ep,2)*c1
461 ENDDO
462 DO m=jft,jlt
463 dm5_2(m)=two*dm5(m)
464 dm6_2(m)=two*dm6(m)
465 df5_2(m)=two*df5(m)
466 df6_2(m)=two*df6(m)
467 ENDDO
468 DO ep=jft,nplat
469 t12(1,2,ep)=half*(t12(1,2,ep)+t12(2,1,ep))
470 t13(1,2,ep)=half*(t13(1,2,ep)+t13(2,1,ep))
471 t14(1,2,ep)=half*(t14(1,2,ep)+t14(2,1,ep))
472 t23(1,2,ep)=half*(t23(1,2,ep)+t23(2,1,ep))
473 t24(1,2,ep)=half*(t24(1,2,ep)+t24(2,1,ep))
474 t34(1,2,ep)=half*(t34(1,2,ep)+t34(2,1,ep))
475 ENDDO
476 IF (idril >0 ) THEN
477 DO ep=jft,nplat
478 k11(1,1,ep)=k11(1,1,ep)+t11(1,2,ep)*dm5_2(ep)
479 km12 =t11(1,1,ep)*dm5(ep)+t11(2,2,ep)*dm6(ep)
480 k11(1,2,ep)=k11(1,2,ep)+km12
481 k11(2,2,ep)=k11(2,2,ep)+t11(1,2,ep)*dm6_2(ep)
482 k12(1,1,ep)=k12(1,1,ep)+t12(1,2,ep)*dm5_2(ep)
483 km12 =t12(1,1,ep)*dm5(ep)+t12(2,2,ep)*dm6(ep)
484 k12(1,2,ep)=k12(1,2,ep)+km12
485 k12(2,1,ep)=k12(2,1,ep)+km12
486 k12(2,2,ep)=k12(2,2,ep)+t12(1,2,ep)*dm6_2(ep)
487 k13(1,1,ep)=k13(1,1,ep)+t13(1,2,ep)*dm5_2(ep)
488 km12 =t13(1,1,ep)*dm5(ep)+t13(2,2,ep)*dm6(ep)
489 k13(1,2,ep)=k13(1,2,ep)+km12
490 k13(2,1,ep)=k13(2,1,ep)+km12
491 k13(2,2,ep)=k13(2,2,ep)+t13(1,2,ep)*dm6_2(ep)
492 k14(1,1,ep)=k14(1,1,ep)+t14(1,2,ep)*dm5_2(ep)
493 km12 =t14(1,1,ep)*dm5(ep)+t14(2,2,ep)*dm6(ep)
494 k14(1,2,ep)=k14(1,2,ep)+km12
495 k14(2,1,ep)=k14(2,1,ep)+km12
496 k14(2,2,ep)=k14(2,2,ep)+t14(1,2,ep)*dm6_2(ep)
497 k22(1,1,ep)=k22(1,1,ep)+t22(1,2,ep)*dm5_2(ep)
498 km12 =t22(1,1,ep)*dm5(ep)+t22(2,2,ep)*dm6(ep)
499 k22(1,2,ep)=k22(1,2,ep)+km12
500 k22(2,2,ep)=k22(2,2,ep)+t22(1,2,ep)*dm6_2(ep)
501 k23(1,1,ep)=k23(1,1,ep)+t23(1,2,ep)*dm5_2(ep)
502 km12 =t23(1,1,ep)*dm5(ep)+t23(2,2,ep)*dm6(ep)
503 k23(1,2,ep)=k23(1,2,ep)+km12
504 k23(2,1,ep)=k23(2,1,ep)+km12
505 k23(2,2,ep)=k23(2,2,ep)+t23(1,2,ep)*dm6_2(ep)
506 k24(1,1,ep)=k24(1,1,ep)+t24(1,2,ep)*dm5_2(ep)
507 km12 =t24(1,1,ep)*dm5(ep)+t24(2,2,ep)*dm6(ep)
508 k24(1,2,ep)=k24(1,2,ep)+km12
509 k24(2,1,ep)=k24(2,1,ep)+km12
510 k24(2,2,ep)=k24(2,2,ep)+t24(1,2,ep)*dm6_2(ep)
511 k33(1,1,ep)=k33(1,1,ep)+t33(1,2,ep)*dm5_2(ep)
512 km12 =t33(1,1,ep)*dm5(ep)+t33(2,2,ep)*dm6(ep)
513 k33(1,2,ep)=k33(1,2,ep)+km12
514 k33(2,2,ep)=k33(2,2,ep)+t33(1,2,ep)*dm6_2(ep)
515 k34(1,1,ep)=k34(1,1,ep)+t34(1,2,ep)*dm5_2(ep)
516 km12 =t34(1,1,ep)*dm5(ep)+t34(2,2,ep)*dm6(ep)
517 k34(1,2,ep)=k34(1,2,ep)+km12
518 k34(2,1,ep)=k34(2,1,ep)+km12
519 k34(2,2,ep)=k34(2,2,ep)+t34(1,2,ep)*dm6_2(ep)
520 k44(1,1,ep)=k44(1,1,ep)+t44(1,2,ep)*dm5_2(ep)
521 km12 =t44(1,1,ep)*dm5(ep)+t44(2,2,ep)*dm6(ep)
522 k44(1,2,ep)=k44(1,2,ep)+km12
523 k44(2,2,ep)=k44(2,2,ep)+t44(1,2,ep)*dm6_2(ep)
524 ENDDO
525 ELSE
526C K11 ,K13 ,K33
527 DO ep=jft,nplat
528 t0ij(1,1,ep)=bb(1,1,ep)*bb0(1,1,ep)
529 t0ij(1,2,ep)=bb(1,1,ep)*bb0(2,1,ep)
530 t0ij(2,1,ep)=bb(2,1,ep)*bb0(1,1,ep)
531 t0ij(2,2,ep)=bb(2,1,ep)*bb0(2,1,ep)
532 t1ij(1,1,ep)=t0ij(1,1,ep)
533 t1ij(1,2,ep)=t0ij(2,1,ep)
534 t1ij(2,1,ep)=t0ij(1,2,ep)
535 t1ij(2,2,ep)=t0ij(2,2,ep)
536 k11(1,1,ep)=k11(1,1,ep)+(t0ij(1,2,ep)+t1ij(2,1,ep))*dm5(ep)
537 km12 =t0ij(1,1,ep)*dm5(ep)+t1ij(2,2,ep)*dm6(ep)
538 k11(1,2,ep)=k11(1,2,ep)+km12
539 k11(2,2,ep)=k11(2,2,ep)+(t0ij(2,1,ep)+t1ij(1,2,ep))*dm6(ep)
540C K13 T0IJ(13)=-T0IJ(11)
541 t1ij(1,1,ep)=bb0(1,1,ep)*bb(1,3,ep)
542 t1ij(1,2,ep)=bb0(1,1,ep)*bb(2,3,ep)
543 t1ij(2,1,ep)=bb0(2,1,ep)*bb(1,3,ep)
544 t1ij(2,2,ep)=bb0(2,1,ep)*bb(2,3,ep)
545 k13(1,1,ep)=k13(1,1,ep)-(t0ij(1,2,ep)-t1ij(2,1,ep))*dm5(ep)
546 km12 =-t0ij(1,1,ep)*dm5(ep)+t1ij(2,2,ep)*dm6(ep)
547 k13(1,2,ep)=k13(1,2,ep)+km12
548 km21 =t1ij(1,1,ep)*dm5(ep)-t0ij(2,2,ep)*dm6(ep)
549 k13(2,1,ep)=k13(2,1,ep)+km21
550 k13(2,2,ep)=k13(2,2,ep)-(t0ij(2,1,ep)-t1ij(1,2,ep))*dm6(ep)
551C-----K33, T1IJ(33)=-T1IJ(13),T0IJ(i,j)=T1IJ(33)(j,i)=-T0IJ(13)
552 t0ij(1,1,ep)=t1ij(1,1,ep)
553 t0ij(1,2,ep)=t1ij(2,1,ep)
554 t0ij(2,1,ep)=t1ij(1,2,ep)
555 t0ij(2,2,ep)=t1ij(2,2,ep)
556 k33(1,1,ep)=k33(1,1,ep)-(t0ij(1,2,ep)+t1ij(2,1,ep))*dm5(ep)
557 km12 =-t0ij(1,1,ep)*dm5(ep)-t1ij(2,2,ep)*dm6(ep)
558 k33(1,2,ep)=k33(1,2,ep)+km12
559 k33(2,2,ep)=k33(2,2,ep)-(t0ij(2,1,ep)+t1ij(1,2,ep))*dm6(ep)
560 ENDDO
561C K12,K14,K34 ,K23
562 DO ep=jft,nplat
563 t0ij(1,1,ep)=bb(1,1,ep)*bb0(1,2,ep)
564 t0ij(1,2,ep)=bb(1,1,ep)*bb0(2,2,ep)
565 t0ij(2,1,ep)=bb(2,1,ep)*bb0(1,2,ep)
566 t0ij(2,2,ep)=bb(2,1,ep)*bb0(2,2,ep)
567 t1ij(1,1,ep)=bb0(1,1,ep)*bb(1,2,ep)
568 t1ij(1,2,ep)=bb0(1,1,ep)*bb(2,2,ep)
569 t1ij(2,1,ep)=bb0(2,1,ep)*bb(1,2,ep)
570 t1ij(2,2,ep)=bb0(2,1,ep)*bb(2,2,ep)
571 k12(1,1,ep)=k12(1,1,ep)+(t0ij(1,2,ep)+t1ij(2,1,ep))*dm5(ep)
572 km12 =t0ij(1,1,ep)*dm5(ep)+t1ij(2,2,ep)*dm6(ep)
573 k12(1,2,ep)=k12(1,2,ep)+km12
574 km21 =t1ij(1,1,ep)*dm5(ep)+t0ij(2,2,ep)*dm6(ep)
575 k12(2,1,ep)=k12(2,1,ep)+km21
576 k12(2,2,ep)=k12(2,2,ep)+(t0ij(2,1,ep)+t1ij(1,2,ep))*dm6(ep)
577C-----K14, T0IJ(14)=-T0IJ(12),
578 t1ij(1,1,ep)=bb0(1,1,ep)*bb(1,4,ep)
579 t1ij(1,2,ep)=bb0(1,1,ep)*bb(2,4,ep)
580 t1ij(2,1,ep)=bb0(2,1,ep)*bb(1,4,ep)
581 t1ij(2,2,ep)=bb0(2,1,ep)*bb(2,4,ep)
582 k14(1,1,ep)=k14(1,1,ep)-(t0ij(1,2,ep)-t1ij(2,1,ep))*dm5(ep)
583 km12 =-t0ij(1,1,ep)*dm5(ep)+t1ij(2,2,ep)*dm6(ep)
584 k14(1,2,ep)=k14(1,2,ep)+km12
585 km21 =t1ij(1,1,ep)*dm5(ep)-t0ij(2,2,ep)*dm6(ep)
586 k14(2,1,ep)=k14(2,1,ep)+km21
587 k14(2,2,ep)=k14(2,2,ep)-(t0ij(2,1,ep)-t1ij(1,2,ep))*dm6(ep)
588C-----K34, T1IJ(34)=-T1IJ(14), T0IJ(34)=-T0IJ(32)
589 t0ij(1,1,ep)=bb(1,3,ep)*bb0(1,2,ep)
590 t0ij(1,2,ep)=bb(1,3,ep)*bb0(2,2,ep)
591 t0ij(2,1,ep)=bb(2,3,ep)*bb0(1,2,ep)
592 t0ij(2,2,ep)=bb(2,3,ep)*bb0(2,2,ep)
593 k34(1,1,ep)=k34(1,1,ep)-(t0ij(1,2,ep)+t1ij(2,1,ep))*dm5(ep)
594 km12 =-t0ij(1,1,ep)*dm5(ep)-t1ij(2,2,ep)*dm6(ep)
595 k34(1,2,ep)=k34(1,2,ep)+km12
596 km21 =-t1ij(1,1,ep)*dm5(ep)-t0ij(2,2,ep)*dm6(ep)
597 k34(2,1,ep)=k34(2,1,ep)+km21
598 k34(2,2,ep)=k34(2,2,ep)-(t0ij(2,1,ep)+t1ij(1,2,ep))*dm6(ep)
599C-----K23, T1IJ(i,j)=-T0IJ(34)(j,i)=T0IJ(32)(j,i),
600 t1ij(1,1,ep)=t0ij(1,1,ep)
601 t1ij(1,2,ep)=t0ij(2,1,ep)
602 t1ij(2,1,ep)=t0ij(1,2,ep)
603 t1ij(2,2,ep)=t0ij(2,2,ep)
604 t0ij(1,1,ep)=-bb(1,2,ep)*bb0(1,1,ep)
605 t0ij(1,2,ep)=-bb(1,2,ep)*bb0(2,1,ep)
606 t0ij(2,1,ep)=-bb(2,2,ep)*bb0(1,1,ep)
607 t0ij(2,2,ep)=-bb(2,2,ep)*bb0(2,1,ep)
608 k23(1,1,ep)=k23(1,1,ep)+(t0ij(1,2,ep)+t1ij(2,1,ep))*dm5(ep)
609 km12 =t0ij(1,1,ep)*dm5(ep)+t1ij(2,2,ep)*dm6(ep)
610 k23(1,2,ep)=k23(1,2,ep)+km12
611 km21 =t1ij(1,1,ep)*dm5(ep)+t0ij(2,2,ep)*dm6(ep)
612 k23(2,1,ep)=k23(2,1,ep)+km21
613 k23(2,2,ep)=k23(2,2,ep)+(t0ij(2,1,ep)+t1ij(1,2,ep))*dm6(ep)
614 ENDDO
615C K22 ,K24 ,K44
616 DO ep=jft,nplat
617 t1ij(1,1,ep)=bb0(1,2,ep)*bb(1,2,ep)
618 t1ij(1,2,ep)=bb0(1,2,ep)*bb(2,2,ep)
619 t1ij(2,1,ep)=bb0(2,2,ep)*bb(1,2,ep)
620 t1ij(2,2,ep)=bb0(2,2,ep)*bb(2,2,ep)
621 t0ij(1,1,ep)=t1ij(1,1,ep)
622 t0ij(1,2,ep)=t1ij(2,1,ep)
623 t0ij(2,1,ep)=t1ij(1,2,ep)
624 t0ij(2,2,ep)=t1ij(2,2,ep)
625 k22(1,1,ep)=k22(1,1,ep)+(t0ij(1,2,ep)+t1ij(2,1,ep))*dm5(ep)
626 km12 =t0ij(1,1,ep)*dm5(ep)+t1ij(2,2,ep)*dm6(ep)
627 k22(1,2,ep)=k22(1,2,ep)+km12
628 k22(2,2,ep)=k22(2,2,ep)+(t0ij(2,1,ep)+t1ij(1,2,ep))*dm6(ep)
629C---K24, T0IJ(24)=-T0IJ(22)
630 t1ij(1,1,ep)=bb0(1,2,ep)*bb(1,4,ep)
631 t1ij(1,2,ep)=bb0(1,2,ep)*bb(2,4,ep)
632 t1ij(2,1,ep)=bb0(2,2,ep)*bb(1,4,ep)
633 t1ij(2,2,ep)=bb0(2,2,ep)*bb(2,4,ep)
634 k24(1,1,ep)=k24(1,1,ep)-(t0ij(1,2,ep)-t1ij(2,1,ep))*dm5(ep)
635 km12 =-t0ij(1,1,ep)*dm5(ep)+t1ij(2,2,ep)*dm6(ep)
636 k24(1,2,ep)=k24(1,2,ep)+km12
637 km21 =t1ij(1,1,ep)*dm5(ep)-t0ij(2,2,ep)*dm6(ep)
638 k24(2,1,ep)=k24(2,1,ep)+km21
639 k24(2,2,ep)=k24(2,2,ep)-(t0ij(2,1,ep)-t1ij(1,2,ep))*dm6(ep)
640C---K44, T1IJ(44)=-T1IJ(24) ; T0IJ(44)=-T0IJ(42)
641 t0ij(1,1,ep)=bb(1,4,ep)*bb0(1,2,ep)
642 t0ij(1,2,ep)=bb(1,4,ep)*bb0(2,2,ep)
643 t0ij(2,1,ep)=bb(2,4,ep)*bb0(1,2,ep)
644 t0ij(2,2,ep)=bb(2,4,ep)*bb0(2,2,ep)
645 k44(1,1,ep)=k44(1,1,ep)-(t0ij(1,2,ep)+t1ij(2,1,ep))*dm5(ep)
646 km12 =-t0ij(1,1,ep)*dm5(ep)-t1ij(2,2,ep)*dm6(ep)
647 k44(1,2,ep)=k44(1,2,ep)+km12
648 k44(2,2,ep)=k44(2,2,ep)-(t0ij(2,1,ep)+t1ij(1,2,ep))*dm6(ep)
649 ENDDO
650 END IF !(IDRIL >0 ) THEN
651C
652 DO ep=jft,nplat
653 m11(1,1,ep)=m11(1,1,ep)+t11(1,2,ep)*df6_2(ep)
654 km12 =t11(1,1,ep)*df5(ep)+t11(2,2,ep)*df6(ep)
655 m11(1,2,ep)=m11(1,2,ep)-km12
656 m11(2,2,ep)=m11(2,2,ep)+t11(1,2,ep)*df5_2(ep)
657
658 m12(1,1,ep)=m12(1,1,ep)+t12(1,2,ep)*df6_2(ep)
659 km12 =t12(1,1,ep)*df5(ep)+t12(2,2,ep)*df6(ep)
660 m12(1,2,ep)=m12(1,2,ep)-km12
661 m12(2,1,ep)=m12(2,1,ep)-km12
662 m12(2,2,ep)=m12(2,2,ep)+t12(1,2,ep)*df5_2(ep)
663
664 m13(1,1,ep)=m13(1,1,ep)+t13(1,2,ep)*df6_2(ep)
665 km12 =t13(1,1,ep)*df5(ep)+t13(2,2,ep)*df6(ep)
666 m13(1,2,ep)=m13(1,2,ep)-km12
667 m13(2,1,ep)=m13(2,1,ep)-km12
668 m13(2,2,ep)=m13(2,2,ep)+t13(1,2,ep)*df5_2(ep)
669
670 m14(1,1,ep)=m14(1,1,ep)+t14(1,2,ep)*df6_2(ep)
671 km12 =t14(1,1,ep)*df5(ep)+t14(2,2,ep)*df6(ep)
672 m14(1,2,ep)=m14(1,2,ep)-km12
673 m14(2,1,ep)=m14(2,1,ep)-km12
674 m14(2,2,ep)=m14(2,2,ep)+t14(1,2,ep)*df5_2(ep)
675
676 m22(1,1,ep)=m22(1,1,ep)+t22(1,2,ep)*df6_2(ep)
677 km12 =t22(1,1,ep)*df5(ep)+t22(2,2,ep)*df6(ep)
678 m22(1,2,ep)=m22(1,2,ep)-km12
679 m22(2,2,ep)=m22(2,2,ep)+t22(1,2,ep)*df5_2(ep)
680
681 m23(1,1,ep)=m23(1,1,ep)+t23(1,2,ep)*df6_2(ep)
682 km12 =t23(1,1,ep)*df5(ep)+t23(2,2,ep)*df6(ep)
683 m23(1,2,ep)=m23(1,2,ep)-km12
684 m23(2,1,ep)=m23(2,1,ep)-km12
685 m23(2,2,ep)=m23(2,2,ep)+t23(1,2,ep)*df5_2(ep)
686
687 m24(1,1,ep)=m24(1,1,ep)+t24(1,2,ep)*df6_2(ep)
688 km12 =t24(1,1,ep)*df5(ep)+t24(2,2,ep)*df6(ep)
689 m24(1,2,ep)=m24(1,2,ep)-km12
690 m24(2,1,ep)=m24(2,1,ep)-km12
691 m24(2,2,ep)=m24(2,2,ep)+t24(1,2,ep)*df5_2(ep)
692
693 m33(1,1,ep)=m33(1,1,ep)+t33(1,2,ep)*df6_2(ep)
694 km12 =t33(1,1,ep)*df5(ep)+t33(2,2,ep)*df6(ep)
695 m33(1,2,ep)=m33(1,2,ep)-km12
696 m33(2,2,ep)=m33(2,2,ep)+t33(1,2,ep)*df5_2(ep)
697
698 m34(1,1,ep)=m34(1,1,ep)+t34(1,2,ep)*df6_2(ep)
699 km12 =t34(1,1,ep)*df5(ep)+t34(2,2,ep)*df6(ep)
700 m34(1,2,ep)=m34(1,2,ep)-km12
701 m34(2,1,ep)=m34(2,1,ep)-km12
702 m34(2,2,ep)=m34(2,2,ep)+t34(1,2,ep)*df5_2(ep)
703
704 m44(1,1,ep)=m44(1,1,ep)+t44(1,2,ep)*df6_2(ep)
705 km12 =t44(1,1,ep)*df5(ep)+t44(2,2,ep)*df6(ep)
706 m44(1,2,ep)=m44(1,2,ep)-km12
707 m44(2,2,ep)=m44(2,2,ep)+t44(1,2,ep)*df5_2(ep)
708 ENDDO
709 ENDIF
710C
711 DO ep=jft,nplat
712 bbc(1,1,1,ep)=bc(ep,1,1,1)
713 bbc(2,1,1,ep)=bc(ep,2,1,1)
714 bbc(1,2,1,ep)=bc(ep,1,2,1)
715 bbc(2,2,1,ep)=bc(ep,2,2,1)
716 bbc(1,3,1,ep)=bc(ep,1,3,1)
717 bbc(2,3,1,ep)=bc(ep,2,3,1)
718 bbc(1,1,2,ep)=bc(ep,1,4,1)
719 bbc(2,1,2,ep)=bc(ep,2,4,1)
720 bbc(1,2,2,ep)=bc(ep,1,5,1)
721 bbc(2,2,2,ep)=bc(ep,2,5,1)
722 bbc(1,3,2,ep)=bc(ep,1,1,2)
723 bbc(2,3,2,ep)=bc(ep,2,1,2)
724 bbc(1,1,3,ep)=bc(ep,1,2,2)
725 bbc(2,1,3,ep)=bc(ep,2,2,2)
726 bbc(1,2,3,ep)=bc(ep,1,3,2)
727 bbc(2,2,3,ep)=bc(ep,2,3,2)
728 bbc(1,3,3,ep)=bc(ep,1,4,2)
729 bbc(2,3,3,ep)=bc(ep,2,4,2)
730 bbc(1,1,4,ep)=bc(ep,1,5,2)
731 bbc(2,1,4,ep)=bc(ep,2,5,2)
732 bbc(1,2,4,ep)=bc(ep,1,1,3)
733 bbc(2,2,4,ep)=bc(ep,2,1,3)
734 bbc(1,3,4,ep)=bc(ep,1,2,3)
735 bbc(2,3,4,ep)=bc(ep,2,2,3)
736 ENDDO
737C
738C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
739 DO i=1,2
740 i1 =i+1
741 DO j=i,2
742 j1 =j+1
743 DO ep=jft,nplat
744 m11(i,j,ep)=m11(i,j,ep)+
745 1 bbc(1,i1,1,ep)*dc(1,1,ep)*bbc(1,j1,1,ep)+
746 2 bbc(2,i1,1,ep)*dc(2,2,ep)*bbc(2,j1,1,ep)
747 m22(i,j,ep)=m22(i,j,ep)+
748 1 bbc(1,i1,2,ep)*dc(1,1,ep)*bbc(1,j1,2,ep)+
749 2 bbc(2,i1,2,ep)*dc(2,2,ep)*bbc(2,j1,2,ep)
750 m33(i,j,ep)=m33(i,j,ep)+
751 1 bbc(1,i1,3,ep)*dc(1,1,ep)*bbc(1,j1,3,ep)+
752 2 bbc(2,i1,3,ep)*dc(2,2,ep)*bbc(2,j1,3,ep)
753 m44(i,j,ep)=m44(i,j,ep)+
754 1 bbc(1,i1,4,ep)*dc(1,1,ep)*bbc(1,j1,4,ep)+
755 2 bbc(2,i1,4,ep)*dc(2,2,ep)*bbc(2,j1,4,ep)
756 ENDDO
757 ENDDO
758 ENDDO
759C
760 DO i=1,2
761 i1 =i+1
762 DO j=1,2
763 j1 =j+1
764 DO ep=jft,nplat
765 m12(i,j,ep)=m12(i,j,ep)+
766 1 bbc(1,i1,1,ep)*dc(1,1,ep)*bbc(1,j1,2,ep)+
767 2 bbc(2,i1,1,ep)*dc(2,2,ep)*bbc(2,j1,2,ep)
768 m13(i,j,ep)=m13(i,j,ep)+
769 1 bbc(1,i1,1,ep)*dc(1,1,ep)*bbc(1,j1,3,ep)+
770 2 bbc(2,i1,1,ep)*dc(2,2,ep)*bbc(2,j1,3,ep)
771 m14(i,j,ep)=m14(i,j,ep)+
772 1 bbc(1,i1,1,ep)*dc(1,1,ep)*bbc(1,j1,4,ep)+
773 2 bbc(2,i1,1,ep)*dc(2,2,ep)*bbc(2,j1,4,ep)
774 m23(i,j,ep)=m23(i,j,ep)+
775 1 bbc(1,i1,2,ep)*dc(1,1,ep)*bbc(1,j1,3,ep)+
776 2 bbc(2,i1,2,ep)*dc(2,2,ep)*bbc(2,j1,3,ep)
777 m24(i,j,ep)=m24(i,j,ep)+
778 1 bbc(1,i1,2,ep)*dc(1,1,ep)*bbc(1,j1,4,ep)+
779 2 bbc(2,i1,2,ep)*dc(2,2,ep)*bbc(2,j1,4,ep)
780 m34(i,j,ep)=m34(i,j,ep)+
781 1 bbc(1,i1,3,ep)*dc(1,1,ep)*bbc(1,j1,4,ep)+
782 2 bbc(2,i1,3,ep)*dc(2,2,ep)*bbc(2,j1,4,ep)
783 ENDDO
784 ENDDO
785 ENDDO
786C
787 DO ep=jft,nplat
788 k11(3,3,ep)=k11(3,3,ep)+
789 1 bbc(1,1,1,ep)*dc(1,1,ep)*bbc(1,1,1,ep)+
790 2 bbc(2,1,1,ep)*dc(2,2,ep)*bbc(2,1,1,ep)
791 k12(3,3,ep)=k12(3,3,ep)+
792 1 bbc(1,1,1,ep)*dc(1,1,ep)*bbc(1,1,2,ep)+
793 2 bbc(2,1,1,ep)*dc(2,2,ep)*bbc(2,1,2,ep)
794 k13(3,3,ep)=k13(3,3,ep)+
795 1 bbc(1,1,1,ep)*dc(1,1,ep)*bbc(1,1,3,ep)+
796 2 bbc(2,1,1,ep)*dc(2,2,ep)*bbc(2,1,3,ep)
797 k14(3,3,ep)=k14(3,3,ep)+
798 1 bbc(1,1,1,ep)*dc(1,1,ep)*bbc(1,1,4,ep)+
799 2 bbc(2,1,1,ep)*dc(2,2,ep)*bbc(2,1,4,ep)
800 k22(3,3,ep)=k22(3,3,ep)+
801 1 bbc(1,1,2,ep)*dc(1,1,ep)*bbc(1,1,2,ep)+
802 2 bbc(2,1,2,ep)*dc(2,2,ep)*bbc(2,1,2,ep)
803 k23(3,3,ep)=k23(3,3,ep)+
804 1 bbc(1,1,2,ep)*dc(1,1,ep)*bbc(1,1,3,ep)+
805 2 bbc(2,1,2,ep)*dc(2,2,ep)*bbc(2,1,3,ep)
806 k24(3,3,ep)=k24(3,3,ep)+
807 1 bbc(1,1,2,ep)*dc(1,1,ep)*bbc(1,1,4,ep)+
808 2 bbc(2,1,2,ep)*dc(2,2,ep)*bbc(2,1,4,ep)
809 k33(3,3,ep)=k33(3,3,ep)+
810 1 bbc(1,1,3,ep)*dc(1,1,ep)*bbc(1,1,3,ep)+
811 2 bbc(2,1,3,ep)*dc(2,2,ep)*bbc(2,1,3,ep)
812 k34(3,3,ep)=k34(3,3,ep)+
813 1 bbc(1,1,3,ep)*dc(1,1,ep)*bbc(1,1,4,ep)+
814 2 bbc(2,1,3,ep)*dc(2,2,ep)*bbc(2,1,4,ep)
815 k44(3,3,ep)=k44(3,3,ep)+
816 1 bbc(1,1,4,ep)*dc(1,1,ep)*bbc(1,1,4,ep)+
817 2 bbc(2,1,4,ep)*dc(2,2,ep)*bbc(2,1,4,ep)
818 ENDDO
819 DO j=1,2
820 j1 =j+1
821 DO ep=jft,nplat
822 mf11(3,j,ep)=mf11(3,j,ep)+
823 1 bbc(1,1,1,ep)*dc(1,1,ep)*bbc(1,j1,1,ep)+
824 2 bbc(2,1,1,ep)*dc(2,2,ep)*bbc(2,j1,1,ep)
825 mf22(3,j,ep)=mf22(3,j,ep)+
826 1 bbc(1,1,2,ep)*dc(1,1,ep)*bbc(1,j1,2,ep)+
827 2 bbc(2,1,2,ep)*dc(2,2,ep)*bbc(2,j1,2,ep)
828 mf33(3,j,ep)=mf33(3,j,ep)+
829 1 bbc(1,1,3,ep)*dc(1,1,ep)*bbc(1,j1,3,ep)+
830 2 bbc(2,1,3,ep)*dc(2,2,ep)*bbc(2,j1,3,ep)
831 mf44(3,j,ep)=mf44(3,j,ep)+
832 1 bbc(1,1,4,ep)*dc(1,1,ep)*bbc(1,j1,4,ep)+
833 2 bbc(2,1,4,ep)*dc(2,2,ep)*bbc(2,j1,4,ep)
834 mf12(3,j,ep)=mf12(3,j,ep)+
835 1 bbc(1,1,1,ep)*dc(1,1,ep)*bbc(1,j1,2,ep)+
836 2 bbc(2,1,1,ep)*dc(2,2,ep)*bbc(2,j1,2,ep)
837 mf13(3,j,ep)=mf13(3,j,ep)+
838 1 bbc(1,1,1,ep)*dc(1,1,ep)*bbc(1,j1,3,ep)+
839 2 bbc(2,1,1,ep)*dc(2,2,ep)*bbc(2,j1,3,ep)
840 mf14(3,j,ep)=mf14(3,j,ep)+
841 1 bbc(1,1,1,ep)*dc(1,1,ep)*bbc(1,j1,4,ep)+
842 2 bbc(2,1,1,ep)*dc(2,2,ep)*bbc(2,j1,4,ep)
843 mf23(3,j,ep)=mf23(3,j,ep)+
844 1 bbc(1,1,2,ep)*dc(1,1,ep)*bbc(1,j1,3,ep)+
845 2 bbc(2,1,2,ep)*dc(2,2,ep)*bbc(2,j1,3,ep)
846 mf24(3,j,ep)=mf24(3,j,ep)+
847 1 bbc(1,1,2,ep)*dc(1,1,ep)*bbc(1,j1,4,ep)+
848 2 bbc(2,1,2,ep)*dc(2,2,ep)*bbc(2,j1,4,ep)
849 mf34(3,j,ep)=mf34(3,j,ep)+
850 1 bbc(1,1,3,ep)*dc(1,1,ep)*bbc(1,j1,4,ep)+
851 2 bbc(2,1,3,ep)*dc(2,2,ep)*bbc(2,j1,4,ep)
852 fm12(j,3,ep)=fm12(j,3,ep)+
853 1 bbc(1,1,2,ep)*dc(1,1,ep)*bbc(1,j1,1,ep)+
854 2 bbc(2,1,2,ep)*dc(2,2,ep)*bbc(2,j1,1,ep)
855 fm13(j,3,ep)=fm13(j,3,ep)+
856 1 bbc(1,1,3,ep)*dc(1,1,ep)*bbc(1,j1,1,ep)+
857 2 bbc(2,1,3,ep)*dc(2,2,ep)*bbc(2,j1,1,ep)
858 fm14(j,3,ep)=fm14(j,3,ep)+
859 1 bbc(1,1,4,ep)*dc(1,1,ep)*bbc(1,j1,1,ep)+
860 2 bbc(2,1,4,ep)*dc(2,2,ep)*bbc(2,j1,1,ep)
861 fm23(j,3,ep)=fm23(j,3,ep)+
862 1 bbc(1,1,3,ep)*dc(1,1,ep)*bbc(1,j1,2,ep)+
863 2 bbc(2,1,3,ep)*dc(2,2,ep)*bbc(2,j1,2,ep)
864 fm24(j,3,ep)=fm24(j,3,ep)+
865 1 bbc(1,1,4,ep)*dc(1,1,ep)*bbc(1,j1,2,ep)+
866 2 bbc(2,1,4,ep)*dc(2,2,ep)*bbc(2,j1,2,ep)
867 fm34(j,3,ep)=fm34(j,3,ep)+
868 1 bbc(1,1,4,ep)*dc(1,1,ep)*bbc(1,j1,3,ep)+
869 2 bbc(2,1,4,ep)*dc(2,2,ep)*bbc(2,j1,3,ep)
870 ENDDO
871 ENDDO
872C------ Mzz pour tous----x:1,3,5,7;y:2,4,6,8---
873 IF (idril==0) THEN
874 DO ep=jft,jlt
875 m11(3,3,ep)=m11(3,3,ep)+dhz(ep)
876 1 *(bzz(ep,1)*bzz(ep,1)+bzz(ep,2)*bzz(ep,2))
877 m12(3,3,ep)=m12(3,3,ep)+dhz(ep)
878 1 *(bzz(ep,1)*bzz(ep,3)+bzz(ep,2)*bzz(ep,4))
879 m13(3,3,ep)=m13(3,3,ep)+dhz(ep)
880 1 *(bzz(ep,1)*bzz(ep,5)+bzz(ep,2)*bzz(ep,6))
881 m14(3,3,ep)=m14(3,3,ep)+dhz(ep)
882 1 *(bzz(ep,1)*bzz(ep,7)+bzz(ep,2)*bzz(ep,8))
883 m22(3,3,ep)=m22(3,3,ep)+dhz(ep)
884 1 *(bzz(ep,3)*bzz(ep,3)+bzz(ep,4)*bzz(ep,4))
885 m23(3,3,ep)=m23(3,3,ep)+dhz(ep)
886 1 *(bzz(ep,3)*bzz(ep,5)+bzz(ep,4)*bzz(ep,6))
887 m24(3,3,ep)=m24(3,3,ep)+dhz(ep)
888 1 *(bzz(ep,3)*bzz(ep,7)+bzz(ep,4)*bzz(ep,8))
889 m33(3,3,ep)=m33(3,3,ep)+dhz(ep)
890 1 *(bzz(ep,5)*bzz(ep,5)+bzz(ep,6)*bzz(ep,6))
891 m34(3,3,ep)=m34(3,3,ep)+dhz(ep)
892 1 *(bzz(ep,5)*bzz(ep,7)+bzz(ep,6)*bzz(ep,8))
893 m44(3,3,ep)=m44(3,3,ep)+dhz(ep)
894 1 *(bzz(ep,7)*bzz(ep,7)+bzz(ep,8)*bzz(ep,8))
895 ENDDO
896C------ renforce positive ----------
897 IF (neig==0) THEN
898 DO ep=jft,jlt
899 c1 = em8*m11(3,3,ep)
900 c2 = em8*m22(3,3,ep)
901 c3 = em8*m33(3,3,ep)
902 c4 = em8*m44(3,3,ep)
903 m11(3,3,ep)=m11(3,3,ep)+c1
904 m22(3,3,ep)=m22(3,3,ep)+c2
905 m33(3,3,ep)=m33(3,3,ep)+c3
906 m44(3,3,ep)=m44(3,3,ep)+c4
907 ENDDO
908 ENDIF
909 END IF !(IDRIL==0) THEN
910C---+----1----+----2----+warped elements----+----5----+----6----+----7----+----8
911#include "vectorize.inc"
912 DO m=nf,jlt
913 ep=iplat(m)
914 c2=vol(ep)
915 c1=thk2(ep)*c2
916 dc(1,1,m)=(hc(ep,1)*tc(ep,1,1)*tc(ep,1,1)+
917 1 hc(ep,2)*tc(ep,1,2)*tc(ep,1,2))*c2
918 dc(2,2,m)=(hc(ep,1)*tc(ep,2,1)*tc(ep,2,1)+
919 1 hc(ep,2)*tc(ep,2,2)*tc(ep,2,2))*c2
920 dc(1,2,m)=(hc(ep,1)*tc(ep,1,1)*tc(ep,2,1)+
921 1 hc(ep,2)*tc(ep,2,2)*tc(ep,1,2))*c2
922 dc(2,1,m)=dc(1,2,m)
923 ENDDO
924 DO i=1,3
925 DO j=i,3
926 DO l=1,2
927 DO k=1,2
928 DO ep=nf,jlt
929 k11(i,j,ep)=k11(i,j,ep)+bm(ep,k,i,1)*dm(k,l,ep)*bm(ep,l,j,1)+
930 1 bmf(ep,k,i,1)*df(k,l,ep)*bmf(ep,l,j,1)+
931 2 bc(ep,k,i,1)*dc(k,l,ep)*bc(ep,l,j,1)
932 k22(i,j,ep)=k22(i,j,ep)+bm(ep,k,i,2)*dm(k,l,ep)*bm(ep,l,j,2)+
933 1 bmf(ep,k,i,2)*df(k,l,ep)*bmf(ep,l,j,2)+
934 2 bc(ep,k,i,2)*dc(k,l,ep)*bc(ep,l,j,2)
935 k33(i,j,ep)=k33(i,j,ep)+bm(ep,k,i,3)*dm(k,l,ep)*bm(ep,l,j,3)+
936 1 bmf(ep,k,i,3)*df(k,l,ep)*bmf(ep,l,j,3)+
937 2 bc(ep,k,i,3)*dc(k,l,ep)*bc(ep,l,j,3)
938 k44(i,j,ep)=k44(i,j,ep)+bm(ep,k,i,4)*dm(k,l,ep)*bm(ep,l,j,4)+
939 1 bmf(ep,k,i,4)*df(k,l,ep)*bmf(ep,l,j,4)+
940 2 bc(ep,k,i,4)*dc(k,l,ep)*bc(ep,l,j,4)
941 ENDDO
942 ENDDO
943 ENDDO
944 ENDDO
945 ENDDO
946C
947 DO i=1,2
948 i1 =i+3
949 DO j=i,2
950 j1 =j+3
951 DO l=1,2
952 DO k=1,2
953 DO ep=nf,jlt
954 m11(i,j,ep)=m11(i,j,ep)+bf(ep,k,i,1)*df(k,l,ep)*bf(ep,l,j,1)+
955 1 bc(ep,k,i1,1)*dc(k,l,ep)*bc(ep,l,j1,1)
956 m22(i,j,ep)=m22(i,j,ep)+bf(ep,k,i,2)*df(k,l,ep)*bf(ep,l,j,2)+
957 1 bc(ep,k,i1,2)*dc(k,l,ep)*bc(ep,l,j1,2)
958 m33(i,j,ep)=m33(i,j,ep)+bf(ep,k,i,3)*df(k,l,ep)*bf(ep,l,j,3)+
959 1 bc(ep,k,i1,3)*dc(k,l,ep)*bc(ep,l,j1,3)
960 m44(i,j,ep)=m44(i,j,ep)+bf(ep,k,i,4)*df(k,l,ep)*bf(ep,l,j,4)+
961 1 bc(ep,k,i1,4)*dc(k,l,ep)*bc(ep,l,j1,4)
962 ENDDO
963 ENDDO
964 ENDDO
965 ENDDO
966 ENDDO
967C---+----1----+----2----+shear contribution-+----5----+----6----+----7----+----8
968 DO i=1,3
969 DO j=i,3
970 DO ep=nf,jlt
971 k11(i,j,ep)=k11(i,j,ep)+bmf(ep,3,i,1)*gf(ep)*bmf(ep,3,j,1)
972 k22(i,j,ep)=k22(i,j,ep)+bmf(ep,3,i,2)*gf(ep)*bmf(ep,3,j,2)
973 k33(i,j,ep)=k33(i,j,ep)+bmf(ep,3,i,3)*gf(ep)*bmf(ep,3,j,3)
974 k44(i,j,ep)=k44(i,j,ep)+bmf(ep,3,i,4)*gf(ep)*bmf(ep,3,j,4)
975 ENDDO
976 ENDDO
977 ENDDO
978 DO i=1,2
979 DO j=i,2
980 DO ep=nf,jlt
981 m11(i,j,ep)=m11(i,j,ep)+bf(ep,3,i,1)*gf(ep)*bf(ep,3,j,1)
982 m22(i,j,ep)=m22(i,j,ep)+bf(ep,3,i,2)*gf(ep)*bf(ep,3,j,2)
983 m33(i,j,ep)=m33(i,j,ep)+bf(ep,3,i,3)*gf(ep)*bf(ep,3,j,3)
984 m44(i,j,ep)=m44(i,j,ep)+bf(ep,3,i,4)*gf(ep)*bf(ep,3,j,4)
985 ENDDO
986 ENDDO
987 ENDDO
988C
989 DO i=1,3
990 DO j=1,3
991 DO l=1,2
992 DO k=1,2
993 DO ep=nf,jlt
994 k12(i,j,ep)=k12(i,j,ep)+bm(ep,k,i,1)*dm(k,l,ep)*bm(ep,l,j,2)+
995 1 bmf(ep,k,i,1)*df(k,l,ep)*bmf(ep,l,j,2)+
996 1 bc(ep,k,i,1)*dc(k,l,ep)*bc(ep,l,j,2)
997 k13(i,j,ep)=k13(i,j,ep)+bm(ep,k,i,1)*dm(k,l,ep)*bm(ep,l,j,3)+
998 1 bmf(ep,k,i,1)*df(k,l,ep)*bmf(ep,l,j,3)+
999 1 bc(ep,k,i,1)*dc(k,l,ep)*bc(ep,l,j,3)
1000 k14(i,j,ep)=k14(i,j,ep)+bm(ep,k,i,1)*dm(k,l,ep)*bm(ep,l,j,4)+
1001 1 bmf(ep,k,i,1)*df(k,l,ep)*bmf(ep,l,j,4)+
1002 1 bc(ep,k,i,1)*dc(k,l,ep)*bc(ep,l,j,4)
1003 k23(i,j,ep)=k23(i,j,ep)+bm(ep,k,i,2)*dm(k,l,ep)*bm(ep,l,j,3)+
1004 1 bmf(ep,k,i,2)*df(k,l,ep)*bmf(ep,l,j,3)+
1005 1 bc(ep,k,i,2)*dc(k,l,ep)*bc(ep,l,j,3)
1006 k24(i,j,ep)=k24(i,j,ep)+bm(ep,k,i,2)*dm(k,l,ep)*bm(ep,l,j,4)+
1007 1 bmf(ep,k,i,2)*df(k,l,ep)*bmf(ep,l,j,4)+
1008 1 bc(ep,k,i,2)*dc(k,l,ep)*bc(ep,l,j,4)
1009 k34(i,j,ep)=k34(i,j,ep)+bm(ep,k,i,3)*dm(k,l,ep)*bm(ep,l,j,4)+
1010 1 bmf(ep,k,i,3)*df(k,l,ep)*bmf(ep,l,j,4)+
1011 1 bc(ep,k,i,3)*dc(k,l,ep)*bc(ep,l,j,4)
1012 ENDDO
1013 ENDDO
1014 ENDDO
1015 ENDDO
1016 ENDDO
1017C
1018 DO i=1,2
1019 i1 =i+3
1020 DO j=1,2
1021 j1 =j+3
1022 DO l=1,2
1023 DO k=1,2
1024 DO ep=nf,jlt
1025 m12(i,j,ep)=m12(i,j,ep)+bf(ep,k,i,1)*df(k,l,ep)*bf(ep,l,j,2)+
1026 1 bc(ep,k,i1,1)*dc(k,l,ep)*bc(ep,l,j1,2)
1027 m13(i,j,ep)=m13(i,j,ep)+bf(ep,k,i,1)*df(k,l,ep)*bf(ep,l,j,3)+
1028 1 bc(ep,k,i1,1)*dc(k,l,ep)*bc(ep,l,j1,3)
1029 m14(i,j,ep)=m14(i,j,ep)+bf(ep,k,i,1)*df(k,l,ep)*bf(ep,l,j,4)+
1030 1 bc(ep,k,i1,1)*dc(k,l,ep)*bc(ep,l,j1,4)
1031 m23(i,j,ep)=m23(i,j,ep)+bf(ep,k,i,2)*df(k,l,ep)*bf(ep,l,j,3)+
1032 1 bc(ep,k,i1,2)*dc(k,l,ep)*bc(ep,l,j1,3)
1033 m24(i,j,ep)=m24(i,j,ep)+bf(ep,k,i,2)*df(k,l,ep)*bf(ep,l,j,4)+
1034 1 bc(ep,k,i1,2)*dc(k,l,ep)*bc(ep,l,j1,4)
1035 m34(i,j,ep)=m34(i,j,ep)+bf(ep,k,i,3)*df(k,l,ep)*bf(ep,l,j,4)+
1036 1 bc(ep,k,i1,3)*dc(k,l,ep)*bc(ep,l,j1,4)
1037 ENDDO
1038 ENDDO
1039 ENDDO
1040 ENDDO
1041 ENDDO
1042C---+----1----+----2----+shear contribution-+----5----+----6----+----7----+----8
1043 DO i=1,3
1044 DO j=1,3
1045 DO ep=nf,jlt
1046 k12(i,j,ep)=k12(i,j,ep)+bmf(ep,3,i,1)*gf(ep)*bmf(ep,3,j,2)
1047 k13(i,j,ep)=k13(i,j,ep)+bmf(ep,3,i,1)*gf(ep)*bmf(ep,3,j,3)
1048 k14(i,j,ep)=k14(i,j,ep)+bmf(ep,3,i,1)*gf(ep)*bmf(ep,3,j,4)
1049 k23(i,j,ep)=k23(i,j,ep)+bmf(ep,3,i,2)*gf(ep)*bmf(ep,3,j,3)
1050 k24(i,j,ep)=k24(i,j,ep)+bmf(ep,3,i,2)*gf(ep)*bmf(ep,3,j,4)
1051 k34(i,j,ep)=k34(i,j,ep)+bmf(ep,3,i,3)*gf(ep)*bmf(ep,3,j,4)
1052 ENDDO
1053 ENDDO
1054 ENDDO
1055 DO i=1,2
1056 DO j=1,2
1057 DO ep=nf,jlt
1058 m12(i,j,ep)=m12(i,j,ep)+bf(ep,3,i,1)*gf(ep)*bf(ep,3,j,2)
1059 m13(i,j,ep)=m13(i,j,ep)+bf(ep,3,i,1)*gf(ep)*bf(ep,3,j,3)
1060 m14(i,j,ep)=m14(i,j,ep)+bf(ep,3,i,1)*gf(ep)*bf(ep,3,j,4)
1061 m23(i,j,ep)=m23(i,j,ep)+bf(ep,3,i,2)*gf(ep)*bf(ep,3,j,3)
1062 m24(i,j,ep)=m24(i,j,ep)+bf(ep,3,i,2)*gf(ep)*bf(ep,3,j,4)
1063 m34(i,j,ep)=m34(i,j,ep)+bf(ep,3,i,3)*gf(ep)*bf(ep,3,j,4)
1064 ENDDO
1065 ENDDO
1066 ENDDO
1067C
1068C---+----1----+----2----+couplage part--[MF](3x2)[FM](2x3)-6----+----7----+----8
1069 DO i=1,3
1070 DO j=1,2
1071 j1=j+3
1072 DO l=1,2
1073 DO k=1,2
1074 DO ep=nf,jlt
1075 mf11(i,j,ep)=mf11(i,j,ep)+bmf(ep,k,i,1)*df(k,l,ep)*bf(ep,l,j,1)
1076 1 +bc(ep,k,i,1)*dc(k,l,ep)*bc(ep,l,j1,1)
1077 mf12(i,j,ep)=mf12(i,j,ep)+bmf(ep,k,i,1)*df(k,l,ep)*bf(ep,l,j,2)
1078 1 +bc(ep,k,i,1)*dc(k,l,ep)*bc(ep,l,j1,2)
1079 mf13(i,j,ep)=mf13(i,j,ep)+bmf(ep,k,i,1)*df(k,l,ep)*bf(ep,l,j,3)
1080 1 +bc(ep,k,i,1)*dc(k,l,ep)*bc(ep,l,j1,3)
1081 mf14(i,j,ep)=mf14(i,j,ep)+bmf(ep,k,i,1)*df(k,l,ep)*bf(ep,l,j,4)
1082 1 +bc(ep,k,i,1)*dc(k,l,ep)*bc(ep,l,j1,4)
1083 mf22(i,j,ep)=mf22(i,j,ep)+bmf(ep,k,i,2)*df(k,l,ep)*bf(ep,l,j,2)
1084 1 +bc(ep,k,i,2)*dc(k,l,ep)*bc(ep,l,j1,2)
1085 mf23(i,j,ep)=mf23(i,j,ep)+bmf(ep,k,i,2)*df(k,l,ep)*bf(ep,l,j,3)
1086 1 +bc(ep,k,i,2)*dc(k,l,ep)*bc(ep,l,j1,3)
1087 mf24(i,j,ep)=mf24(i,j,ep)+bmf(ep,k,i,2)*df(k,l,ep)*bf(ep,l,j,4)
1088 1 +bc(ep,k,i,2)*dc(k,l,ep)*bc(ep,l,j1,4)
1089 mf33(i,j,ep)=mf33(i,j,ep)+bmf(ep,k,i,3)*df(k,l,ep)*bf(ep,l,j,3)
1090 1 +bc(ep,k,i,3)*dc(k,l,ep)*bc(ep,l,j1,3)
1091 mf34(i,j,ep)=mf34(i,j,ep)+bmf(ep,k,i,3)*df(k,l,ep)*bf(ep,l,j,4)
1092 1 +bc(ep,k,i,3)*dc(k,l,ep)*bc(ep,l,j1,4)
1093 mf44(i,j,ep)=mf44(i,j,ep)+bmf(ep,k,i,4)*df(k,l,ep)*bf(ep,l,j,4)
1094 1 +bc(ep,k,i,4)*dc(k,l,ep)*bc(ep,l,j1,4)
1095C
1096 fm12(j,i,ep)=fm12(j,i,ep)+bmf(ep,k,i,2)*df(k,l,ep)*bf(ep,l,j,1)
1097 1 +bc(ep,k,i,2)*dc(k,l,ep)*bc(ep,l,j1,1)
1098 fm13(j,i,ep)=fm13(j,i,ep)+bmf(ep,k,i,3)*df(k,l,ep)*bf(ep,l,j,1)
1099 1 +bc(ep,k,i,3)*dc(k,l,ep)*bc(ep,l,j1,1)
1100 fm14(j,i,ep)=fm14(j,i,ep)+bmf(ep,k,i,4)*df(k,l,ep)*bf(ep,l,j,1)
1101 1 +bc(ep,k,i,4)*dc(k,l,ep)*bc(ep,l,j1,1)
1102 fm23(j,i,ep)=fm23(j,i,ep)+bmf(ep,k,i,3)*df(k,l,ep)*bf(ep,l,j,2)
1103 1 +bc(ep,k,i,3)*dc(k,l,ep)*bc(ep,l,j1,2)
1104 fm24(j,i,ep)=fm24(j,i,ep)+bmf(ep,k,i,4)*df(k,l,ep)*bf(ep,l,j,2)
1105 1 +bc(ep,k,i,4)*dc(k,l,ep)*bc(ep,l,j1,2)
1106 fm34(j,i,ep)=fm34(j,i,ep)+bmf(ep,k,i,4)*df(k,l,ep)*bf(ep,l,j,3)
1107 1 +bc(ep,k,i,4)*dc(k,l,ep)*bc(ep,l,j1,3)
1108 ENDDO
1109 ENDDO
1110 ENDDO
1111 ENDDO
1112 ENDDO
1113C
1114 DO i=1,3
1115 DO j=1,2
1116 DO ep=nf,jlt
1117 mf11(i,j,ep)=mf11(i,j,ep)+bmf(ep,3,i,1)*gf(ep)*bf(ep,3,j,1)
1118 mf12(i,j,ep)=mf12(i,j,ep)+bmf(ep,3,i,1)*gf(ep)*bf(ep,3,j,2)
1119 mf13(i,j,ep)=mf13(i,j,ep)+bmf(ep,3,i,1)*gf(ep)*bf(ep,3,j,3)
1120 mf14(i,j,ep)=mf14(i,j,ep)+bmf(ep,3,i,1)*gf(ep)*bf(ep,3,j,4)
1121 mf22(i,j,ep)=mf22(i,j,ep)+bmf(ep,3,i,2)*gf(ep)*bf(ep,3,j,2)
1122 mf23(i,j,ep)=mf23(i,j,ep)+bmf(ep,3,i,2)*gf(ep)*bf(ep,3,j,3)
1123 mf24(i,j,ep)=mf24(i,j,ep)+bmf(ep,3,i,2)*gf(ep)*bf(ep,3,j,4)
1124 mf33(i,j,ep)=mf33(i,j,ep)+bmf(ep,3,i,3)*gf(ep)*bf(ep,3,j,3)
1125 mf34(i,j,ep)=mf34(i,j,ep)+bmf(ep,3,i,3)*gf(ep)*bf(ep,3,j,4)
1126 mf44(i,j,ep)=mf44(i,j,ep)+bmf(ep,3,i,4)*gf(ep)*bf(ep,3,j,4)
1127C
1128 fm12(j,i,ep)=fm12(j,i,ep)+bmf(ep,3,i,2)*gf(ep)*bf(ep,3,j,1)
1129 fm13(j,i,ep)=fm13(j,i,ep)+bmf(ep,3,i,3)*gf(ep)*bf(ep,3,j,1)
1130 fm14(j,i,ep)=fm14(j,i,ep)+bmf(ep,3,i,4)*gf(ep)*bf(ep,3,j,1)
1131 fm23(j,i,ep)=fm23(j,i,ep)+bmf(ep,3,i,3)*gf(ep)*bf(ep,3,j,2)
1132 fm24(j,i,ep)=fm24(j,i,ep)+bmf(ep,3,i,4)*gf(ep)*bf(ep,3,j,2)
1133 fm34(j,i,ep)=fm34(j,i,ep)+bmf(ep,3,i,4)*gf(ep)*bf(ep,3,j,3)
1134 ENDDO
1135 ENDDO
1136 ENDDO
1137 IF (iorth==1) THEN
1138 l=3
1139 k=1
1140 DO i=1,3
1141 DO j=i,3
1142 DO ep=nf,jlt
1143 k11(i,j,ep)=k11(i,j,ep)+bm(ep,k,i,1)*dm5(ep)*bm(ep,l,j,1)+
1144 1 bm(ep,l,i,1)*dm5(ep)*bm(ep,k,j,1)+
1145 1 bmf(ep,k,i,1)*df5(ep)*bmf(ep,l,j,1)+
1146 1 bmf(ep,l,i,1)*df5(ep)*bmf(ep,k,j,1)
1147 k22(i,j,ep)=k22(i,j,ep)+bm(ep,k,i,2)*dm5(ep)*bm(ep,l,j,2)+
1148 1 bm(ep,l,i,2)*dm5(ep)*bm(ep,k,j,2)+
1149 1 bmf(ep,k,i,2)*df5(ep)*bmf(ep,l,j,2)+
1150 1 bmf(ep,l,i,2)*df5(ep)*bmf(ep,k,j,2)
1151 k33(i,j,ep)=k33(i,j,ep)+bm(ep,k,i,3)*dm5(ep)*bm(ep,l,j,3)+
1152 1 bm(ep,l,i,3)*dm5(ep)*bm(ep,k,j,3)+
1153 1 bmf(ep,k,i,3)*df5(ep)*bmf(ep,l,j,3)+
1154 1 bmf(ep,l,i,3)*df5(ep)*bmf(ep,k,j,3)
1155 k44(i,j,ep)=k44(i,j,ep)+bm(ep,k,i,4)*dm5(ep)*bm(ep,l,j,4)+
1156 1 bm(ep,l,i,4)*dm5(ep)*bm(ep,k,j,4)+
1157 1 bmf(ep,k,i,4)*df5(ep)*bmf(ep,l,j,4)+
1158 1 bmf(ep,l,i,4)*df5(ep)*bmf(ep,k,j,4)
1159 ENDDO
1160 ENDDO
1161 ENDDO
1162 DO i=1,3
1163 DO j=1,3
1164 DO ep=nf,jlt
1165 k12(i,j,ep)=k12(i,j,ep)+bm(ep,k,i,1)*dm5(ep)*bm(ep,l,j,2)+
1166 1 bm(ep,l,i,1)*dm5(ep)*bm(ep,k,j,2)+
1167 1 bmf(ep,k,i,1)*df5(ep)*bmf(ep,l,j,2)+
1168 1 bmf(ep,l,i,1)*df5(ep)*bmf(ep,k,j,2)
1169 k13(i,j,ep)=k13(i,j,ep)+bm(ep,k,i,1)*dm5(ep)*bm(ep,l,j,3)+
1170 1 bm(ep,l,i,1)*dm5(ep)*bm(ep,k,j,3)+
1171 1 bmf(ep,k,i,1)*df5(ep)*bmf(ep,l,j,3)+
1172 1 bmf(ep,l,i,1)*df5(ep)*bmf(ep,k,j,3)
1173 k14(i,j,ep)=k14(i,j,ep)+bm(ep,k,i,1)*dm5(ep)*bm(ep,l,j,4)+
1174 1 bm(ep,l,i,1)*dm5(ep)*bm(ep,k,j,4)+
1175 1 bmf(ep,k,i,1)*df5(ep)*bmf(ep,l,j,4)+
1176 1 bmf(ep,l,i,1)*df5(ep)*bmf(ep,k,j,4)
1177 k23(i,j,ep)=k23(i,j,ep)+bm(ep,k,i,2)*dm5(ep)*bm(ep,l,j,3)+
1178 1 bm(ep,l,i,2)*dm5(ep)*bm(ep,k,j,3)+
1179 1 bmf(ep,k,i,2)*df5(ep)*bmf(ep,l,j,3)+
1180 1 bmf(ep,l,i,2)*df5(ep)*bmf(ep,k,j,3)
1181 k24(i,j,ep)=k24(i,j,ep)+bm(ep,k,i,2)*dm5(ep)*bm(ep,l,j,4)+
1182 1 bm(ep,l,i,2)*dm5(ep)*bm(ep,k,j,4)+
1183 1 bmf(ep,k,i,2)*df5(ep)*bmf(ep,l,j,4)+
1184 1 bmf(ep,l,i,2)*df5(ep)*bmf(ep,k,j,4)
1185 k34(i,j,ep)=k34(i,j,ep)+bm(ep,k,i,3)*dm5(ep)*bm(ep,l,j,4)+
1186 1 bm(ep,l,i,3)*dm5(ep)*bm(ep,k,j,4)+
1187 1 bmf(ep,k,i,3)*df5(ep)*bmf(ep,l,j,4)+
1188 1 bmf(ep,l,i,3)*df5(ep)*bmf(ep,k,j,4)
1189 ENDDO
1190 ENDDO
1191 ENDDO
1192 DO i=1,2
1193 DO j=i,2
1194 DO ep=nf,jlt
1195 m11(i,j,ep)=m11(i,j,ep)+bf(ep,k,i,1)*df5(ep)*bf(ep,l,j,1)
1196 1 +bf(ep,l,i,1)*df5(ep)*bf(ep,k,j,1)
1197 m22(i,j,ep)=m22(i,j,ep)+bf(ep,k,i,2)*df5(ep)*bf(ep,l,j,2)
1198 1 +bf(ep,l,i,2)*df5(ep)*bf(ep,k,j,2)
1199 m33(i,j,ep)=m33(i,j,ep)+bf(ep,k,i,3)*df5(ep)*bf(ep,l,j,3)
1200 1 +bf(ep,l,i,3)*df5(ep)*bf(ep,k,j,3)
1201 m44(i,j,ep)=m44(i,j,ep)+bf(ep,k,i,4)*df5(ep)*bf(ep,l,j,4)
1202 1 +bf(ep,l,i,4)*df5(ep)*bf(ep,k,j,4)
1203 ENDDO
1204 ENDDO
1205 ENDDO
1206 DO i=1,2
1207 DO j=1,2
1208 DO ep=nf,jlt
1209 m12(i,j,ep)=m12(i,j,ep)+bf(ep,k,i,1)*df5(ep)*bf(ep,l,j,2)
1210 1 +bf(ep,l,i,1)*df5(ep)*bf(ep,k,j,2)
1211 m13(i,j,ep)=m13(i,j,ep)+bf(ep,k,i,1)*df5(ep)*bf(ep,l,j,3)
1212 1 +bf(ep,l,i,1)*df5(ep)*bf(ep,k,j,3)
1213 m14(i,j,ep)=m14(i,j,ep)+bf(ep,k,i,1)*df5(ep)*bf(ep,l,j,4)
1214 1 +bf(ep,l,i,1)*df5(ep)*bf(ep,k,j,4)
1215 m23(i,j,ep)=m23(i,j,ep)+bf(ep,k,i,2)*df5(ep)*bf(ep,l,j,3)
1216 1 +bf(ep,l,i,2)*df5(ep)*bf(ep,k,j,3)
1217 m24(i,j,ep)=m24(i,j,ep)+bf(ep,k,i,2)*df5(ep)*bf(ep,l,j,4)
1218 1 +bf(ep,l,i,2)*df5(ep)*bf(ep,k,j,4)
1219 m34(i,j,ep)=m34(i,j,ep)+bf(ep,k,i,3)*df5(ep)*bf(ep,l,j,4)
1220 1 +bf(ep,l,i,3)*df5(ep)*bf(ep,k,j,4)
1221 ENDDO
1222 ENDDO
1223 ENDDO
1224 DO i=1,3
1225 DO j=1,2
1226 DO ep=nf,jlt
1227 mf11(i,j,ep)=mf11(i,j,ep)+bmf(ep,k,i,1)*df5(ep)*bf(ep,l,j,1)
1228 1 +bmf(ep,l,i,1)*df5(ep)*bf(ep,k,j,1)
1229 mf12(i,j,ep)=mf12(i,j,ep)+bmf(ep,k,i,1)*df5(ep)*bf(ep,l,j,2)
1230 1 +bmf(ep,l,i,1)*df5(ep)*bf(ep,k,j,2)
1231 mf13(i,j,ep)=mf13(i,j,ep)+bmf(ep,k,i,1)*df5(ep)*bf(ep,l,j,3)
1232 1 +bmf(ep,l,i,1)*df5(ep)*bf(ep,k,j,3)
1233 mf14(i,j,ep)=mf14(i,j,ep)+bmf(ep,k,i,1)*df5(ep)*bf(ep,l,j,4)
1234 1 +bmf(ep,l,i,1)*df5(ep)*bf(ep,k,j,4)
1235 mf22(i,j,ep)=mf22(i,j,ep)+bmf(ep,k,i,2)*df5(ep)*bf(ep,l,j,2)
1236 1 +bmf(ep,l,i,2)*df5(ep)*bf(ep,k,j,2)
1237 mf23(i,j,ep)=mf23(i,j,ep)+bmf(ep,k,i,2)*df5(ep)*bf(ep,l,j,3)
1238 1 +bmf(ep,l,i,2)*df5(ep)*bf(ep,k,j,3)
1239 mf24(i,j,ep)=mf24(i,j,ep)+bmf(ep,k,i,2)*df5(ep)*bf(ep,l,j,4)
1240 1 +bmf(ep,l,i,2)*df5(ep)*bf(ep,k,j,4)
1241 mf33(i,j,ep)=mf33(i,j,ep)+bmf(ep,k,i,3)*df5(ep)*bf(ep,l,j,3)
1242 1 +bmf(ep,l,i,3)*df5(ep)*bf(ep,k,j,3)
1243 mf34(i,j,ep)=mf34(i,j,ep)+bmf(ep,k,i,3)*df5(ep)*bf(ep,l,j,4)
1244 1 +bmf(ep,l,i,3)*df5(ep)*bf(ep,k,j,4)
1245 mf44(i,j,ep)=mf44(i,j,ep)+bmf(ep,k,i,4)*df5(ep)*bf(ep,l,j,4)
1246 1 +bmf(ep,l,i,4)*df5(ep)*bf(ep,k,j,4)
1247C
1248 fm12(j,i,ep)=fm12(j,i,ep)+bmf(ep,k,i,2)*df5(ep)*bf(ep,l,j,1)
1249 1 +bmf(ep,l,i,2)*df5(ep)*bf(ep,k,j,1)
1250 fm13(j,i,ep)=fm13(j,i,ep)+bmf(ep,k,i,3)*df5(ep)*bf(ep,l,j,1)
1251 1 +bmf(ep,l,i,3)*df5(ep)*bf(ep,k,j,1)
1252 fm14(j,i,ep)=fm14(j,i,ep)+bmf(ep,k,i,4)*df5(ep)*bf(ep,l,j,1)
1253 1 +bmf(ep,l,i,4)*df5(ep)*bf(ep,k,j,1)
1254 fm23(j,i,ep)=fm23(j,i,ep)+bmf(ep,k,i,3)*df5(ep)*bf(ep,l,j,2)
1255 1 +bmf(ep,l,i,3)*df5(ep)*bf(ep,k,j,2)
1256 fm24(j,i,ep)=fm24(j,i,ep)+bmf(ep,k,i,4)*df5(ep)*bf(ep,l,j,2)
1257 1 +bmf(ep,l,i,4)*df5(ep)*bf(ep,k,j,2)
1258 fm34(j,i,ep)=fm34(j,i,ep)+bmf(ep,k,i,4)*df5(ep)*bf(ep,l,j,3)
1259 1 +bmf(ep,l,i,4)*df5(ep)*bf(ep,k,j,3)
1260 ENDDO
1261 ENDDO
1262 ENDDO
1263C
1264 k=2
1265 DO i=1,3
1266 DO j=i,3
1267 DO ep=nf,jlt
1268 k11(i,j,ep)=k11(i,j,ep)+bm(ep,k,i,1)*dm6(ep)*bm(ep,l,j,1)+
1269 1 bm(ep,l,i,1)*dm6(ep)*bm(ep,k,j,1)+
1270 1 bmf(ep,k,i,1)*df6(ep)*bmf(ep,l,j,1)+
1271 1 bmf(ep,l,i,1)*df6(ep)*bmf(ep,k,j,1)
1272 k22(i,j,ep)=k22(i,j,ep)+bm(ep,k,i,2)*dm6(ep)*bm(ep,l,j,2)+
1273 1 bm(ep,l,i,2)*dm6(ep)*bm(ep,k,j,2)+
1274 1 bmf(ep,k,i,2)*df6(ep)*bmf(ep,l,j,2)+
1275 1 bmf(ep,l,i,2)*df6(ep)*bmf(ep,k,j,2)
1276 k33(i,j,ep)=k33(i,j,ep)+bm(ep,k,i,3)*dm6(ep)*bm(ep,l,j,3)+
1277 1 bm(ep,l,i,3)*dm6(ep)*bm(ep,k,j,3)+
1278 1 bmf(ep,k,i,3)*df6(ep)*bmf(ep,l,j,3)+
1279 1 bmf(ep,l,i,3)*df6(ep)*bmf(ep,k,j,3)
1280 k44(i,j,ep)=k44(i,j,ep)+bm(ep,k,i,4)*dm6(ep)*bm(ep,l,j,4)+
1281 1 bm(ep,l,i,4)*dm6(ep)*bm(ep,k,j,4)+
1282 1 bmf(ep,k,i,4)*df6(ep)*bmf(ep,l,j,4)+
1283 1 bmf(ep,l,i,4)*df6(ep)*bmf(ep,k,j,4)
1284 ENDDO
1285 ENDDO
1286 ENDDO
1287 DO i=1,3
1288 DO j=1,3
1289 DO ep=nf,jlt
1290 k12(i,j,ep)=k12(i,j,ep)+bm(ep,k,i,1)*dm6(ep)*bm(ep,l,j,2)+
1291 1 bm(ep,l,i,1)*dm6(ep)*bm(ep,k,j,2)+
1292 1 bmf(ep,k,i,1)*df6(ep)*bmf(ep,l,j,2)+
1293 1 bmf(ep,l,i,1)*df6(ep)*bmf(ep,k,j,2)
1294 k13(i,j,ep)=k13(i,j,ep)+bm(ep,k,i,1)*dm6(ep)*bm(ep,l,j,3)+
1295 1 bm(ep,l,i,1)*dm6(ep)*bm(ep,k,j,3)+
1296 1 bmf(ep,k,i,1)*df6(ep)*bmf(ep,l,j,3)+
1297 1 bmf(ep,l,i,1)*df6(ep)*bmf(ep,k,j,3)
1298 k14(i,j,ep)=k14(i,j,ep)+bm(ep,k,i,1)*dm6(ep)*bm(ep,l,j,4)+
1299 1 bm(ep,l,i,1)*dm6(ep)*bm(ep,k,j,4)+
1300 1 bmf(ep,k,i,1)*df6(ep)*bmf(ep,l,j,4)+
1301 1 bmf(ep,l,i,1)*df6(ep)*bmf(ep,k,j,4)
1302 k23(i,j,ep)=k23(i,j,ep)+bm(ep,k,i,2)*dm6(ep)*bm(ep,l,j,3)+
1303 1 bm(ep,l,i,2)*dm6(ep)*bm(ep,k,j,3)+
1304 1 bmf(ep,k,i,2)*df6(ep)*bmf(ep,l,j,3)+
1305 1 bmf(ep,l,i,2)*df6(ep)*bmf(ep,k,j,3)
1306 k24(i,j,ep)=k24(i,j,ep)+bm(ep,k,i,2)*dm6(ep)*bm(ep,l,j,4)+
1307 1 bm(ep,l,i,2)*dm6(ep)*bm(ep,k,j,4)+
1308 1 bmf(ep,k,i,2)*df6(ep)*bmf(ep,l,j,4)+
1309 1 bmf(ep,l,i,2)*df6(ep)*bmf(ep,k,j,4)
1310 k34(i,j,ep)=k34(i,j,ep)+bm(ep,k,i,3)*dm6(ep)*bm(ep,l,j,4)+
1311 1 bm(ep,l,i,3)*dm6(ep)*bm(ep,k,j,4)+
1312 1 bmf(ep,k,i,3)*df6(ep)*bmf(ep,l,j,4)+
1313 1 bmf(ep,l,i,3)*df6(ep)*bmf(ep,k,j,4)
1314 ENDDO
1315 ENDDO
1316 ENDDO
1317C
1318 DO i=1,2
1319 DO j=i,2
1320 DO ep=nf,jlt
1321 m11(i,j,ep)=m11(i,j,ep)+bf(ep,k,i,1)*df6(ep)*bf(ep,l,j,1)
1322 1 +bf(ep,l,i,1)*df6(ep)*bf(ep,k,j,1)
1323 m22(i,j,ep)=m22(i,j,ep)+bf(ep,k,i,2)*df6(ep)*bf(ep,l,j,2)
1324 1 +bf(ep,l,i,2)*df6(ep)*bf(ep,k,j,2)
1325 m33(i,j,ep)=m33(i,j,ep)+bf(ep,k,i,3)*df6(ep)*bf(ep,l,j,3)
1326 1 +bf(ep,l,i,3)*df6(ep)*bf(ep,k,j,3)
1327 m44(i,j,ep)=m44(i,j,ep)+bf(ep,k,i,4)*df6(ep)*bf(ep,l,j,4)
1328 1 +bf(ep,l,i,4)*df6(ep)*bf(ep,k,j,4)
1329 ENDDO
1330 ENDDO
1331 ENDDO
1332 DO i=1,2
1333 DO j=1,2
1334 DO ep=nf,jlt
1335 m12(i,j,ep)=m12(i,j,ep)+bf(ep,k,i,1)*df6(ep)*bf(ep,l,j,2)
1336 1 +bf(ep,l,i,1)*df6(ep)*bf(ep,k,j,2)
1337 m13(i,j,ep)=m13(i,j,ep)+bf(ep,k,i,1)*df6(ep)*bf(ep,l,j,3)
1338 1 +bf(ep,l,i,1)*df6(ep)*bf(ep,k,j,3)
1339 m14(i,j,ep)=m14(i,j,ep)+bf(ep,k,i,1)*df6(ep)*bf(ep,l,j,4)
1340 1 +bf(ep,l,i,1)*df6(ep)*bf(ep,k,j,4)
1341 m23(i,j,ep)=m23(i,j,ep)+bf(ep,k,i,2)*df6(ep)*bf(ep,l,j,3)
1342 1 +bf(ep,l,i,2)*df6(ep)*bf(ep,k,j,3)
1343 m24(i,j,ep)=m24(i,j,ep)+bf(ep,k,i,2)*df6(ep)*bf(ep,l,j,4)
1344 1 +bf(ep,l,i,2)*df6(ep)*bf(ep,k,j,4)
1345 m34(i,j,ep)=m34(i,j,ep)+bf(ep,k,i,3)*df6(ep)*bf(ep,l,j,4)
1346 1 +bf(ep,l,i,3)*df6(ep)*bf(ep,k,j,4)
1347 ENDDO
1348 ENDDO
1349 ENDDO
1350 DO i=1,3
1351 DO j=1,2
1352 DO ep=nf,jlt
1353 mf11(i,j,ep)=mf11(i,j,ep)+bmf(ep,k,i,1)*df6(ep)*bf(ep,l,j,1)
1354 1 +bmf(ep,l,i,1)*df6(ep)*bf(ep,k,j,1)
1355 mf12(i,j,ep)=mf12(i,j,ep)+bmf(ep,k,i,1)*df6(ep)*bf(ep,l,j,2)
1356 1 +bmf(ep,l,i,1)*df6(ep)*bf(ep,k,j,2)
1357 mf13(i,j,ep)=mf13(i,j,ep)+bmf(ep,k,i,1)*df6(ep)*bf(ep,l,j,3)
1358 1 +bmf(ep,l,i,1)*df6(ep)*bf(ep,k,j,3)
1359 mf14(i,j,ep)=mf14(i,j,ep)+bmf(ep,k,i,1)*df6(ep)*bf(ep,l,j,4)
1360 1 +bmf(ep,l,i,1)*df6(ep)*bf(ep,k,j,4)
1361 mf22(i,j,ep)=mf22(i,j,ep)+bmf(ep,k,i,2)*df6(ep)*bf(ep,l,j,2)
1362 1 +bmf(ep,l,i,2)*df6(ep)*bf(ep,k,j,2)
1363 mf23(i,j,ep)=mf23(i,j,ep)+bmf(ep,k,i,2)*df6(ep)*bf(ep,l,j,3)
1364 1 +bmf(ep,l,i,2)*df6(ep)*bf(ep,k,j,3)
1365 mf24(i,j,ep)=mf24(i,j,ep)+bmf(ep,k,i,2)*df6(ep)*bf(ep,l,j,4)
1366 1 +bmf(ep,l,i,2)*df6(ep)*bf(ep,k,j,4)
1367 mf33(i,j,ep)=mf33(i,j,ep)+bmf(ep,k,i,3)*df6(ep)*bf(ep,l,j,3)
1368 1 +bmf(ep,l,i,3)*df6(ep)*bf(ep,k,j,3)
1369 mf34(i,j,ep)=mf34(i,j,ep)+bmf(ep,k,i,3)*df6(ep)*bf(ep,l,j,4)
1370 1 +bmf(ep,l,i,3)*df6(ep)*bf(ep,k,j,4)
1371 mf44(i,j,ep)=mf44(i,j,ep)+bmf(ep,k,i,4)*df6(ep)*bf(ep,l,j,4)
1372 1 +bmf(ep,l,i,4)*df6(ep)*bf(ep,k,j,4)
1373C
1374 fm12(j,i,ep)=fm12(j,i,ep)+bmf(ep,k,i,2)*df6(ep)*bf(ep,l,j,1)
1375 1 +bmf(ep,l,i,2)*df6(ep)*bf(ep,k,j,1)
1376 fm13(j,i,ep)=fm13(j,i,ep)+bmf(ep,k,i,3)*df6(ep)*bf(ep,l,j,1)
1377 1 +bmf(ep,l,i,3)*df6(ep)*bf(ep,k,j,1)
1378 fm14(j,i,ep)=fm14(j,i,ep)+bmf(ep,k,i,4)*df6(ep)*bf(ep,l,j,1)
1379 1 +bmf(ep,l,i,4)*df6(ep)*bf(ep,k,j,1)
1380 fm23(j,i,ep)=fm23(j,i,ep)+bmf(ep,k,i,3)*df6(ep)*bf(ep,l,j,2)
1381 1 +bmf(ep,l,i,3)*df6(ep)*bf(ep,k,j,2)
1382 fm24(j,i,ep)=fm24(j,i,ep)+bmf(ep,k,i,4)*df6(ep)*bf(ep,l,j,2)
1383 1 +bmf(ep,l,i,4)*df6(ep)*bf(ep,k,j,2)
1384 fm34(j,i,ep)=fm34(j,i,ep)+bmf(ep,k,i,4)*df6(ep)*bf(ep,l,j,3)
1385 1 +bmf(ep,l,i,4)*df6(ep)*bf(ep,k,j,3)
1386 ENDDO
1387 ENDDO
1388 ENDDO
1389C--add mem-bending coupling terms of orthotrope--warped elements
1390 DO i=1,3
1391 DO j=i,3
1392 DO l=1,3
1393 DO k=1,3
1394 DO ep=nf,jlt
1395 k11(i,j,ep)=k11(i,j,ep)+
1396 1 bm(ep,k,i,1)*dmf(k,l,ep)*bmf(ep,l,j,1)+
1397 2 bmf(ep,k,i,1)*dmf(k,l,ep)*bm(ep,l,j,1)
1398 k22(i,j,ep)=k22(i,j,ep)+
1399 1 bm(ep,k,i,2)*dmf(k,l,ep)*bmf(ep,l,j,2)+
1400 2 bmf(ep,k,i,2)*dmf(k,l,ep)*bm(ep,l,j,2)
1401 k33(i,j,ep)=k33(i,j,ep)+
1402 1 bm(ep,k,i,3)*dmf(k,l,ep)*bmf(ep,l,j,3)+
1403 2 bmf(ep,k,i,3)*dmf(k,l,ep)*bm(ep,l,j,3)
1404 k44(i,j,ep)=k44(i,j,ep)+
1405 1 bm(ep,k,i,4)*dmf(k,l,ep)*bmf(ep,l,j,4)+
1406 2 bmf(ep,k,i,4)*dmf(k,l,ep)*bm(ep,l,j,4)
1407 ENDDO
1408 ENDDO
1409 ENDDO
1410 ENDDO
1411 ENDDO
1412C
1413 DO i=1,3
1414 DO j=1,3
1415 DO l=1,3
1416 DO k=1,3
1417 DO ep=nf,jlt
1418 k12(i,j,ep)=k12(i,j,ep)+
1419 1 bm(ep,k,i,1)*dmf(k,l,ep)*bmf(ep,l,j,2)+
1420 2 bmf(ep,k,i,1)*dmf(k,l,ep)*bm(ep,l,j,2)
1421 k13(i,j,ep)=k13(i,j,ep)+
1422 1 bm(ep,k,i,1)*dmf(k,l,ep)*bmf(ep,l,j,3)+
1423 2 bmf(ep,k,i,1)*dmf(k,l,ep)*bm(ep,l,j,3)
1424 k14(i,j,ep)=k14(i,j,ep)+
1425 1 bm(ep,k,i,1)*dmf(k,l,ep)*bmf(ep,l,j,4)+
1426 2 bmf(ep,k,i,1)*dmf(k,l,ep)*bm(ep,l,j,4)
1427 k23(i,j,ep)=k23(i,j,ep)+
1428 1 bm(ep,k,i,2)*dmf(k,l,ep)*bmf(ep,l,j,3)+
1429 2 bmf(ep,k,i,2)*dmf(k,l,ep)*bm(ep,l,j,3)
1430 k24(i,j,ep)=k24(i,j,ep)+
1431 1 bm(ep,k,i,2)*dmf(k,l,ep)*bmf(ep,l,j,4)+
1432 2 bmf(ep,k,i,2)*dmf(k,l,ep)*bm(ep,l,j,4)
1433 k34(i,j,ep)=k34(i,j,ep)+
1434 1 bm(ep,k,i,3)*dmf(k,l,ep)*bmf(ep,l,j,4)+
1435 2 bmf(ep,k,i,3)*dmf(k,l,ep)*bm(ep,l,j,4)
1436 ENDDO
1437 ENDDO
1438 ENDDO
1439 ENDDO
1440 ENDDO
1441C--------------------------
1442 DO i=1,3
1443 DO j=1,2
1444 DO l=1,3
1445 DO k=1,3
1446 DO ep=nf,jlt
1447 mf11(i,j,ep)=mf11(i,j,ep)+
1448 1 bm(ep,k,i,1)*dmf(k,l,ep)*bf(ep,l,j,1)
1449 mf12(i,j,ep)=mf12(i,j,ep)+
1450 1 bm(ep,k,i,1)*dmf(k,l,ep)*bf(ep,l,j,2)
1451 mf13(i,j,ep)=mf13(i,j,ep)+
1452 1 bm(ep,k,i,1)*dmf(k,l,ep)*bf(ep,l,j,3)
1453 mf14(i,j,ep)=mf14(i,j,ep)+
1454 1 bm(ep,k,i,1)*dmf(k,l,ep)*bf(ep,l,j,4)
1455 mf22(i,j,ep)=mf22(i,j,ep)+
1456 1 bm(ep,k,i,2)*dmf(k,l,ep)*bf(ep,l,j,2)
1457 mf23(i,j,ep)=mf23(i,j,ep)+
1458 1 bm(ep,k,i,2)*dmf(k,l,ep)*bf(ep,l,j,3)
1459 mf24(i,j,ep)=mf24(i,j,ep)+
1460 1 bm(ep,k,i,2)*dmf(k,l,ep)*bf(ep,l,j,4)
1461 mf33(i,j,ep)=mf33(i,j,ep)+
1462 1 bm(ep,k,i,3)*dmf(k,l,ep)*bf(ep,l,j,3)
1463 mf34(i,j,ep)=mf34(i,j,ep)+
1464 1 bm(ep,k,i,3)*dmf(k,l,ep)*bf(ep,l,j,4)
1465 mf44(i,j,ep)=mf44(i,j,ep)+
1466 1 bm(ep,k,i,4)*dmf(k,l,ep)*bf(ep,l,j,4)
1467 ENDDO
1468 ENDDO
1469 ENDDO
1470 ENDDO
1471 ENDDO
1472C
1473 DO i=1,2
1474 DO j=1,3
1475 DO l=1,3
1476 DO k=1,3
1477 DO ep=nf,jlt
1478 fm12(i,j,ep)=fm12(i,j,ep)+
1479 1 bf(ep,k,i,1)*dmf(k,l,ep)*bm(ep,l,j,2)
1480 fm13(i,j,ep)=fm13(i,j,ep)+
1481 1 bf(ep,k,i,1)*dmf(k,l,ep)*bm(ep,l,j,3)
1482 fm14(i,j,ep)=fm14(i,j,ep)+
1483 1 bf(ep,k,i,1)*dmf(k,l,ep)*bm(ep,l,j,4)
1484 fm23(i,j,ep)=fm23(i,j,ep)+
1485 1 bf(ep,k,i,2)*dmf(k,l,ep)*bm(ep,l,j,3)
1486 fm24(i,j,ep)=fm24(i,j,ep)+
1487 1 bf(ep,k,i,2)*dmf(k,l,ep)*bm(ep,l,j,4)
1488 fm34(i,j,ep)=fm34(i,j,ep)+
1489 1 bf(ep,k,i,3)*dmf(k,l,ep)*bm(ep,l,j,4)
1490 ENDDO
1491 ENDDO
1492 ENDDO
1493 ENDDO
1494 ENDDO
1495C
1496 ENDIF !(IORTH==1)
1497C
1498 IF (ikgeo==1) THEN
1499#include "vectorize.inc"
1500 DO m=jft,jlt
1501 ep=iplat(m)
1502 c2=vol(ep)
1503 fxx(m)=for(ep,1)*c2
1504 fyy(m)=for(ep,2)*c2
1505 ENDDO
1506C
1507 DO i=1,2
1508 DO ep=jft,nplat
1509 t11(i,i,ep)=bb(i,1,ep)*bb(i,1,ep)
1510 t22(i,i,ep)=bb(i,2,ep)*bb(i,2,ep)
1511 t33(i,i,ep)=bb(i,3,ep)*bb(i,3,ep)
1512 t44(i,i,ep)=bb(i,4,ep)*bb(i,4,ep)
1513 t12(i,i,ep)=bb(i,1,ep)*bb(i,2,ep)
1514 t13(i,i,ep)=bb(i,1,ep)*bb(i,3,ep)
1515 t14(i,i,ep)=bb(i,1,ep)*bb(i,4,ep)
1516 t23(i,i,ep)=bb(i,2,ep)*bb(i,3,ep)
1517 t24(i,i,ep)=bb(i,2,ep)*bb(i,4,ep)
1518 t34(i,i,ep)=bb(i,3,ep)*bb(i,4,ep)
1519 ENDDO
1520 ENDDO
1521C
1522 DO i=1,2
1523 DO ep=nf,jlt
1524 t11(i,i,ep)=bm(ep,i,i,1)*bm(ep,i,i,1)
1525 t22(i,i,ep)=bm(ep,i,i,2)*bm(ep,i,i,2)
1526 t33(i,i,ep)=bm(ep,i,i,3)*bm(ep,i,i,3)
1527 t44(i,i,ep)=bm(ep,i,i,4)*bm(ep,i,i,4)
1528 t12(i,i,ep)=bm(ep,i,i,1)*bm(ep,i,i,2)
1529 t13(i,i,ep)=bm(ep,i,i,1)*bm(ep,i,i,3)
1530 t14(i,i,ep)=bm(ep,i,i,1)*bm(ep,i,i,4)
1531 t23(i,i,ep)=bm(ep,i,i,2)*bm(ep,i,i,3)
1532 t24(i,i,ep)=bm(ep,i,i,2)*bm(ep,i,i,4)
1533 t34(i,i,ep)=bm(ep,i,i,3)*bm(ep,i,i,4)
1534 ENDDO
1535 ENDDO
1536 DO ep=jft,jlt
1537 t11(1,2,ep)=fxx(ep)*t11(1,1,ep)+fyy(ep)*t11(2,2,ep)
1538 t22(1,2,ep)=fxx(ep)*t22(1,1,ep)+fyy(ep)*t22(2,2,ep)
1539 t33(1,2,ep)=fxx(ep)*t33(1,1,ep)+fyy(ep)*t33(2,2,ep)
1540 t44(1,2,ep)=fxx(ep)*t44(1,1,ep)+fyy(ep)*t44(2,2,ep)
1541 t12(1,2,ep)=fxx(ep)*t12(1,1,ep)+fyy(ep)*t12(2,2,ep)
1542 t13(1,2,ep)=fxx(ep)*t13(1,1,ep)+fyy(ep)*t13(2,2,ep)
1543 t14(1,2,ep)=fxx(ep)*t14(1,1,ep)+fyy(ep)*t14(2,2,ep)
1544 t23(1,2,ep)=fxx(ep)*t23(1,1,ep)+fyy(ep)*t23(2,2,ep)
1545 t24(1,2,ep)=fxx(ep)*t24(1,1,ep)+fyy(ep)*t24(2,2,ep)
1546 t34(1,2,ep)=fxx(ep)*t34(1,1,ep)+fyy(ep)*t34(2,2,ep)
1547 ENDDO
1548C
1549 IF (idril>0) THEN
1550#include "vectorize.inc"
1551 DO m=jft,jlt
1552 ep=iplat(m)
1553 fxy(m)=for(ep,3)*vol(ep)
1554 fxy2(m)=two*fxy(m)
1555 ENDDO
1556C---------distinque plat and warped
1557 DO ep=jft,nplat
1558 t11(1,2,ep)=t11(1,2,ep)+fxy2(ep)*bb(1,1,ep)*bb(2,1,ep)
1559 t22(1,2,ep)=t22(1,2,ep)+fxy2(ep)*bb(1,2,ep)*bb(2,2,ep)
1560 t33(1,2,ep)=t33(1,2,ep)+fxy2(ep)*bb(1,3,ep)*bb(2,3,ep)
1561 t44(1,2,ep)=t44(1,2,ep)+fxy2(ep)*bb(1,4,ep)*bb(2,4,ep)
1562 t12(1,2,ep)=t12(1,2,ep)+fxy(ep)*(bb(1,1,ep)*bb(2,2,ep)+
1563 . bb(2,1,ep)*bb(1,2,ep))
1564 t13(1,2,ep)=t13(1,2,ep)+fxy(ep)*(bb(1,1,ep)*bb(2,3,ep)+
1565 . bb(2,1,ep)*bb(1,3,ep))
1566 t14(1,2,ep)=t14(1,2,ep)+fxy(ep)*(bb(1,1,ep)*bb(2,4,ep)+
1567 . bb(2,1,ep)*bb(1,4,ep))
1568 t23(1,2,ep)=t23(1,2,ep)+fxy(ep)*(bb(1,2,ep)*bb(2,3,ep)+
1569 . bb(2,2,ep)*bb(1,3,ep))
1570 t24(1,2,ep)=t24(1,2,ep)+fxy(ep)*(bb(1,2,ep)*bb(2,4,ep)+
1571 . bb(2,2,ep)*bb(1,4,ep))
1572 t34(1,2,ep)=t34(1,2,ep)+fxy(ep)*(bb(1,3,ep)*bb(2,4,ep)+
1573 . bb(2,3,ep)*bb(1,4,ep))
1574 ENDDO
1575 DO ep=nf,jlt
1576 t11(1,2,ep)=t11(1,2,ep)+fxy2(ep)*bm(ep,3,1,1)*bm(ep,3,2,1)
1577 t22(1,2,ep)=t22(1,2,ep)+fxy2(ep)*bm(ep,3,1,2)*bm(ep,3,2,2)
1578 t33(1,2,ep)=t33(1,2,ep)+fxy2(ep)*bm(ep,3,1,3)*bm(ep,3,2,3)
1579 t44(1,2,ep)=t44(1,2,ep)+fxy2(ep)*bm(ep,3,1,4)*bm(ep,3,2,4)
1580 t12(1,2,ep)=t12(1,2,ep)+fxy(ep)*(bm(ep,3,1,1)*bm(ep,3,2,2)+
1581 . bm(ep,3,2,1)*bm(ep,3,1,2))
1582 t13(1,2,ep)=t13(1,2,ep)+fxy(ep)*(bm(ep,3,1,1)*bm(ep,3,2,3)+
1583 . bm(ep,3,2,1)*bm(ep,3,1,3))
1584 t14(1,2,ep)=t14(1,2,ep)+fxy(ep)*(bm(ep,3,1,1)*bm(ep,3,2,4)+
1585 . bm(ep,3,2,1)*bm(ep,3,1,4))
1586 t23(1,2,ep)=t23(1,2,ep)+fxy(ep)*(bm(ep,3,1,2)*bm(ep,3,2,3)+
1587 . bm(ep,3,2,2)*bm(ep,3,1,3))
1588 t24(1,2,ep)=t24(1,2,ep)+fxy(ep)*(bm(ep,3,1,2)*bm(ep,3,2,4)+
1589 . bm(ep,3,2,2)*bm(ep,3,1,4))
1590 t34(1,2,ep)=t34(1,2,ep)+fxy(ep)*(bm(ep,3,1,3)*bm(ep,3,2,4)+
1591 . bm(ep,3,2,3)*bm(ep,3,1,4))
1592 ENDDO
1593C
1594 ENDIF
1595C
1596 DO i=1,3
1597 DO ep=jft,jlt
1598 k11(i,i,ep)=k11(i,i,ep)+t11(1,2,ep)
1599 k22(i,i,ep)=k22(i,i,ep)+t22(1,2,ep)
1600 k33(i,i,ep)=k33(i,i,ep)+t33(1,2,ep)
1601 k44(i,i,ep)=k44(i,i,ep)+t44(1,2,ep)
1602 k12(i,i,ep)=k12(i,i,ep)+t12(1,2,ep)
1603 k13(i,i,ep)=k13(i,i,ep)+t13(1,2,ep)
1604 k14(i,i,ep)=k14(i,i,ep)+t14(1,2,ep)
1605 k23(i,i,ep)=k23(i,i,ep)+t23(1,2,ep)
1606 k24(i,i,ep)=k24(i,i,ep)+t24(1,2,ep)
1607 k34(i,i,ep)=k34(i,i,ep)+t34(1,2,ep)
1608 ENDDO
1609 ENDDO
1610C------ renforce positive ----------
1611 IF (neig==0.AND.idril==0.AND. iorth >0 ) THEN
1612 DO ep=jft,jlt
1613 c1 = min(t11(1,2,ep),t22(1,2,ep),t33(1,2,ep),t44(1,2,ep))
1614 IF (c1 < zero) THEN
1615 c2 =min(m11(3,3,ep),m22(3,3,ep),m33(3,3,ep),m44(3,3,ep))
1616 c1 = min(-c1,five*c2)
1617 m11(3,3,ep)=m11(3,3,ep) + c1
1618 m22(3,3,ep)=m22(3,3,ep) + c1
1619 m33(3,3,ep)=m33(3,3,ep) + c1
1620 m44(3,3,ep)=m44(3,3,ep) + c1
1621 END IF
1622 ENDDO
1623 ENDIF
1624C
1625 END IF !IF (IKGEO==1)
1626C
1627 RETURN
cbalkeormf1
subroutine cbalkeormf1(jft, jlt, kmfij, dmf, tij, t0ij)
Definition cbalke3.F:2611
cbalkeorfm
subroutine cbalkeorfm(jft, jlt, kfmij, dmf, tij)
Definition cbalke3.F:2560
cbalkeormf
subroutine cbalkeormf(jft, jlt, kmfij, dmf, tij)
Definition cbalke3.F:2510
cbalkeorfm1
subroutine cbalkeorfm1(jft, jlt, kfmij, dmf, tij, t0ij)
Definition cbalke3.F:2661
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91

◆ cbalkec3()

subroutine cbalkec3 ( integer jft,
integer jlt,
vol,
x13,
x24,
y13,
y24,
hm,
k11,
k12,
k13,
k14,
k22,
k23,
k24,
k33,
k34,
k44,
integer nplat,
integer, dimension(*) iplat,
integer ikgeo,
for,
m11,
m22,
m33,
m44,
integer iorth,
integer nel )

Definition at line 1634 of file cbalke3.F.

1638C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1639C CALCUL 'membrane shear traitement' use only PARTIE CONSTANTE
1640C-----------------------------------------------
1641C I M P L I C I T T Y P E S
1642C-----------------------------------------------
1643#include "implicit_f.inc"
1644#include "mvsiz_p.inc"
1645#include "com04_c.inc"
1646C-----------------------------------------------
1647C D U M M Y A R G U M E N T S
1648C-----------------------------------------------
1649 INTEGER JFT ,JLT,NPLAT,IPLAT(*),IKGEO,IORTH,NEL
1650 my_real
1651 . vol(*), x13(*) ,x24(*) ,y13(*) ,y24(*) ,
1652 . k11(3,3,*),k12(3,3,*),k13(3,3,*),k14(3,3,*),
1653 . k22(3,3,*),k23(3,3,*),k24(3,3,*),k33(3,3,*),
1654 . k34(3,3,*),k44(3,3,*),hm(mvsiz,4),for(nel,5),
1655 . m11(3,3,*),m22(3,3,*),m33(3,3,*),m44(3,3,*)
1656C-----------------------------------------------
1657C L O C A L V A R I A B L E S
1658C-----------------------------------------------
1659 INTEGER I ,J,EP,M
1660 my_real
1661 . bb(2,2,mvsiz),gvol(mvsiz),fxy(mvsiz),fxy2(mvsiz),
1662 . t11(2,2,mvsiz),t22(2,2,mvsiz),t12(2,2,mvsiz),c1,c2
1663C
1664#include "vectorize.inc"
1665 DO m=jft,jlt
1666 ep=iplat(m)
1667 gvol(m)=hm(ep,4)*vol(ep)
1668 bb(2,1,m)=y24(ep)
1669 bb(2,2,m)=-y13(ep)
1670 bb(1,1,m)=-x24(ep)
1671 bb(1,2,m)=x13(ep)
1672 ENDDO
1673C
1674 DO i=1,2
1675 DO j=1,2
1676 DO ep=jft,jlt
1677 t11(i,j,ep)=gvol(ep)*bb(i,1,ep)*bb(j,1,ep)
1678 t22(i,j,ep)=gvol(ep)*bb(i,2,ep)*bb(j,2,ep)
1679 t12(i,j,ep)=gvol(ep)*bb(i,1,ep)*bb(j,2,ep)
1680 ENDDO
1681 ENDDO
1682 ENDDO
1683C
1684 DO i=1,2
1685 DO j=i,2
1686 DO ep=jft,jlt
1687 k11(i,j,ep)=k11(i,j,ep)+t11(i,j,ep)
1688 k22(i,j,ep)=k22(i,j,ep)+t22(i,j,ep)
1689 k33(i,j,ep)=k33(i,j,ep)+t11(i,j,ep)
1690 k44(i,j,ep)=k44(i,j,ep)+t22(i,j,ep)
1691 ENDDO
1692 ENDDO
1693 ENDDO
1694 DO i=1,2
1695 DO j=1,2
1696 DO ep=jft,jlt
1697 k12(i,j,ep)=k12(i,j,ep)+t12(i,j,ep)
1698 k13(i,j,ep)=k13(i,j,ep)-t11(i,j,ep)
1699 k14(i,j,ep)=k14(i,j,ep)-t12(i,j,ep)
1700 k23(i,j,ep)=k23(i,j,ep)-t12(j,i,ep)
1701 k24(i,j,ep)=k24(i,j,ep)-t22(i,j,ep)
1702 k34(i,j,ep)=k34(i,j,ep)+t12(i,j,ep)
1703 ENDDO
1704 ENDDO
1705 ENDDO
1706 IF (ikgeo==1) THEN
1707#include "vectorize.inc"
1708 DO m=jft,jlt
1709 ep=iplat(m)
1710 fxy(m)=for(ep,3)*vol(ep)
1711 fxy2(m)=two*fxy(m)
1712 ENDDO
1713C
1714 DO ep=jft,jlt
1715 t11(1,2,ep)=fxy2(ep)*bb(1,1,ep)*bb(2,1,ep)
1716 t22(1,2,ep)=fxy2(ep)*bb(1,2,ep)*bb(2,2,ep)
1717 t12(1,2,ep)=fxy(ep)*(bb(1,1,ep)*bb(2,2,ep)+
1718 . bb(2,1,ep)*bb(1,2,ep))
1719 ENDDO
1720 DO i=1,3
1721 DO ep=jft,jlt
1722 k11(i,i,ep)=k11(i,i,ep)+t11(1,2,ep)
1723 k22(i,i,ep)=k22(i,i,ep)+t22(1,2,ep)
1724 k33(i,i,ep)=k33(i,i,ep)+t11(1,2,ep)
1725 k44(i,i,ep)=k44(i,i,ep)+t22(1,2,ep)
1726 k12(i,i,ep)=k12(i,i,ep)+t12(1,2,ep)
1727 k13(i,i,ep)=k13(i,i,ep)-t11(1,2,ep)
1728 k14(i,i,ep)=k14(i,i,ep)-t12(1,2,ep)
1729 k23(i,i,ep)=k23(i,i,ep)-t12(1,2,ep)
1730 k24(i,i,ep)=k24(i,i,ep)-t22(1,2,ep)
1731 k34(i,i,ep)=k34(i,i,ep)+t12(1,2,ep)
1732 ENDDO
1733 ENDDO
1734C------ renforce positive ----------
1735 IF (neig == 0.AND. iorth >0 ) THEN
1736 DO ep=jft,jlt
1737 c1 = min(t11(1,2,ep),t22(1,2,ep))
1738 IF (c1 < zero) THEN
1739 c2 =min(m11(3,3,ep),m22(3,3,ep),m33(3,3,ep),m44(3,3,ep))
1740 c1 = min(-c1,two*c2)
1741 m11(3,3,ep)=m11(3,3,ep)+c1
1742 m22(3,3,ep)=m22(3,3,ep)+c1
1743 m33(3,3,ep)=m33(3,3,ep)+c1
1744 m44(3,3,ep)=m44(3,3,ep)+c1
1745 END IF
1746 ENDDO
1747 ENDIF
1748 ENDIF
1749C
1750 RETURN

◆ cbalkeorfm()

subroutine cbalkeorfm ( integer jft,
integer jlt,
kfmij,
dmf,
tij )

Definition at line 2559 of file cbalke3.F.

2560C-----------------------------------------------
2561C I M P L I C I T T Y P E S
2562C-----------------------------------------------
2563#include "implicit_f.inc"
2564#include "mvsiz_p.inc"
2565C-----------------------------------------------
2566C D U M M Y A R G U M E N T S
2567C-----------------------------------------------
2568C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2569 INTEGER JFT,JLT
2570 my_real
2571 . kfmij(3,3,*),tij(2,2,*),dmf(3,3,*)
2572C-----------------------------------------------
2573C L O C A L V A R I A B L E S
2574C-----------------------------------------------
2575 INTEGER EP,I,J
2576 my_real
2577 . fmij(2,2,mvsiz)
2578C-------FMIJ=[BbI]^t[DMF][BmJ]=MFJI^t;TIJ(i,j)=B(i,I)*B(j,J)--TJI(i,j)=B(i,J)*B(j,I)
2579C-------I<J->TJI(i,j)=TIJ(j,i)
2580 DO ep=jft,jlt
2581 fmij(1,1,ep)=
2582 1 -dmf(1,2,ep)*tij(2,1,ep)-dmf(1,3,ep)*tij(1,1,ep)
2583 2 -dmf(3,2,ep)*tij(2,2,ep)-dmf(3,3,ep)*tij(1,2,ep)
2584 fmij(2,1,ep)=
2585 1 +dmf(1,1,ep)*tij(1,1,ep)+dmf(1,3,ep)*tij(2,1,ep)
2586 2 +dmf(3,1,ep)*tij(1,2,ep)+dmf(3,3,ep)*tij(2,2,ep)
2587 fmij(1,2,ep)=
2588 1 -dmf(2,2,ep)*tij(2,2,ep)-dmf(2,3,ep)*tij(1,2,ep)
2589 2 -dmf(3,2,ep)*tij(2,1,ep)-dmf(3,3,ep)*tij(1,1,ep)
2590 fmij(2,2,ep)=
2591 1 +dmf(2,1,ep)*tij(1,2,ep)+dmf(2,3,ep)*tij(2,2,ep)
2592 2 +dmf(3,1,ep)*tij(1,1,ep)+dmf(3,3,ep)*tij(2,1,ep)
2593 ENDDO
2594C
2595 DO i=1,2
2596 DO j=1,2
2597 DO ep=jft,jlt
2598 kfmij(i,j,ep)=kfmij(i,j,ep)+fmij(i,j,ep)
2599 ENDDO
2600 ENDDO
2601 ENDDO
2602C
2603 RETURN

◆ cbalkeorfm1()

subroutine cbalkeorfm1 ( integer jft,
integer jlt,
kfmij,
dmf,
tij,
t0ij )

Definition at line 2660 of file cbalke3.F.

2661C-----------------------------------------------
2662C I M P L I C I T T Y P E S
2663C-----------------------------------------------
2664#include "implicit_f.inc"
2665#include "mvsiz_p.inc"
2666C-----------------------------------------------
2667C D U M M Y A R G U M E N T S
2668C-----------------------------------------------
2669C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2670 INTEGER JFT,JLT
2671 my_real
2672 . kfmij(3,3,*),tij(2,2,*),dmf(3,3,*),t0ij(2,2,*)
2673C-----------------------------------------------
2674C L O C A L V A R I A B L E S
2675C-----------------------------------------------
2676 INTEGER EP,I,J
2677 my_real
2678 . fmij(2,2,mvsiz)
2679C-------FMIJ=[BbI]^t[DMF][BmJ]=MFJI^t;TIJ(i,j)=B(i,I)*B(j,J)--TJI(i,j)=B(i,J)*B(j,I)
2680C-------I<J->TJI(i,j)=TIJ(j,i)
2681 DO ep=jft,jlt
2682 fmij(1,1,ep)=
2683 1 -dmf(1,2,ep)*tij(2,1,ep)-dmf(1,3,ep)*tij(1,1,ep)
2684 2 -dmf(3,2,ep)*t0ij(2,2,ep)-dmf(3,3,ep)*t0ij(1,2,ep)
2685 fmij(2,1,ep)=
2686 1 +dmf(1,1,ep)*tij(1,1,ep)+dmf(1,3,ep)*tij(2,1,ep)
2687 2 +dmf(3,1,ep)*t0ij(1,2,ep)+dmf(3,3,ep)*t0ij(2,2,ep)
2688 fmij(1,2,ep)=
2689 1 -dmf(2,2,ep)*tij(2,2,ep)-dmf(2,3,ep)*tij(1,2,ep)
2690 2 -dmf(3,2,ep)*t0ij(2,1,ep)-dmf(3,3,ep)*t0ij(1,1,ep)
2691 fmij(2,2,ep)=
2692 1 +dmf(2,1,ep)*tij(1,2,ep)+dmf(2,3,ep)*tij(2,2,ep)
2693 2 +dmf(3,1,ep)*t0ij(1,1,ep)+dmf(3,3,ep)*t0ij(2,1,ep)
2694 ENDDO
2695C
2696 DO i=1,2
2697 DO j=1,2
2698 DO ep=jft,jlt
2699 kfmij(i,j,ep)=kfmij(i,j,ep)+fmij(i,j,ep)
2700 ENDDO
2701 ENDDO
2702 ENDDO
2703C
2704 RETURN

◆ cbalkeormf()

subroutine cbalkeormf ( integer jft,
integer jlt,
kmfij,
dmf,
tij )

Definition at line 2509 of file cbalke3.F.

2510C-----------------------------------------------
2511C I M P L I C I T T Y P E S
2512C-----------------------------------------------
2513#include "implicit_f.inc"
2514#include "mvsiz_p.inc"
2515C-----------------------------------------------
2516C D U M M Y A R G U M E N T S
2517C-----------------------------------------------
2518C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2519 INTEGER JFT,JLT
2520 my_real
2521 . kmfij(3,3,*),tij(2,2,*),dmf(3,3,*)
2522C-----------------------------------------------
2523C L O C A L V A R I A B L E S
2524C-----------------------------------------------
2525 INTEGER EP,I,J
2526 my_real
2527 . mfij(2,2,mvsiz)
2528C-------MFIJ=[BmI]^t[DMF][BbJ] ;TIJ(i,j)=B(i,I)*B(j,J)----i,j->1:x;2:y
2529 DO ep=jft,jlt
2530 mfij(1,1,ep)=
2531 1 -dmf(1,2,ep)*tij(1,2,ep)-dmf(1,3,ep)*tij(1,1,ep)
2532 2 -dmf(3,2,ep)*tij(2,2,ep)-dmf(3,3,ep)*tij(2,1,ep)
2533 mfij(1,2,ep)=
2534 1 +dmf(1,1,ep)*tij(1,1,ep)+dmf(1,3,ep)*tij(1,2,ep)
2535 2 +dmf(3,1,ep)*tij(2,1,ep)+dmf(3,3,ep)*tij(2,2,ep)
2536 mfij(2,1,ep)=
2537 1 -dmf(2,2,ep)*tij(2,2,ep)-dmf(2,3,ep)*tij(2,1,ep)
2538 2 -dmf(3,2,ep)*tij(1,2,ep)-dmf(3,3,ep)*tij(1,1,ep)
2539 mfij(2,2,ep)=
2540 1 +dmf(2,1,ep)*tij(2,1,ep)+dmf(2,3,ep)*tij(2,2,ep)
2541 2 +dmf(3,1,ep)*tij(1,1,ep)+dmf(3,3,ep)*tij(1,2,ep)
2542 ENDDO
2543C
2544 DO i=1,2
2545 DO j=1,2
2546 DO ep=jft,jlt
2547 kmfij(i,j,ep)=kmfij(i,j,ep)+mfij(i,j,ep)
2548 ENDDO
2549 ENDDO
2550 ENDDO
2551C
2552 RETURN

◆ cbalkeormf1()

subroutine cbalkeormf1 ( integer jft,
integer jlt,
kmfij,
dmf,
tij,
t0ij )

Definition at line 2610 of file cbalke3.F.

2611C-----------------------------------------------
2612C I M P L I C I T T Y P E S
2613C-----------------------------------------------
2614#include "implicit_f.inc"
2615#include "mvsiz_p.inc"
2616C-----------------------------------------------
2617C D U M M Y A R G U M E N T S
2618C-----------------------------------------------
2619C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2620 INTEGER JFT,JLT
2621 my_real
2622 . kmfij(3,3,*),tij(2,2,*),dmf(3,3,*),t0ij(2,2,*)
2623C-----------------------------------------------
2624C L O C A L V A R I A B L E S
2625C-----------------------------------------------
2626 INTEGER EP,I,J
2627 my_real
2628 . mfij(2,2,mvsiz)
2629C----cas Idril=0:MFIJ=[BmI]^t[DMF][BbJ] ;TIJ(i,j)=B(i,I)*B(j,J)----i,j->1:x;2:y
2630 DO ep=jft,jlt
2631 mfij(1,1,ep)=
2632 1 -dmf(1,2,ep)*tij(1,2,ep)-dmf(1,3,ep)*tij(1,1,ep)
2633 2 -dmf(3,2,ep)*t0ij(2,2,ep)-dmf(3,3,ep)*t0ij(2,1,ep)
2634 mfij(1,2,ep)=
2635 1 +dmf(1,1,ep)*tij(1,1,ep)+dmf(1,3,ep)*tij(1,2,ep)
2636 2 +dmf(3,1,ep)*t0ij(2,1,ep)+dmf(3,3,ep)*t0ij(2,2,ep)
2637 mfij(2,1,ep)=
2638 1 -dmf(2,2,ep)*tij(2,2,ep)-dmf(2,3,ep)*tij(2,1,ep)
2639 2 -dmf(3,2,ep)*t0ij(1,2,ep)-dmf(3,3,ep)*t0ij(1,1,ep)
2640 mfij(2,2,ep)=
2641 1 +dmf(2,1,ep)*tij(2,1,ep)+dmf(2,3,ep)*tij(2,2,ep)
2642 2 +dmf(3,1,ep)*t0ij(1,1,ep)+dmf(3,3,ep)*t0ij(1,2,ep)
2643 ENDDO
2644C
2645 DO i=1,2
2646 DO j=1,2
2647 DO ep=jft,jlt
2648 kmfij(i,j,ep)=kmfij(i,j,ep)+mfij(i,j,ep)
2649 ENDDO
2650 ENDDO
2651 ENDDO
2652C
2653 RETURN

◆ cbalkerz()

subroutine cbalkerz ( integer jft,
integer jlt,
vol,
thk0,
hm,
hz,
bm,
k11,
k12,
k13,
k14,
k22,
k23,
k24,
k33,
k34,
k44,
m11,
m12,
m13,
m14,
m22,
m23,
m24,
m33,
m34,
m44,
mf11,
mf12,
mf13,
mf14,
mf22,
mf23,
mf24,
mf33,
mf34,
mf44,
fm12,
fm13,
fm14,
fm23,
fm24,
fm34,
integer iorth,
hmor,
hfor,
integer, dimension(*) iplat,
integer nplat,
pmrz,
brz,
frz,
integer ikgeo,
integer ng,
hmfor,
bf,
bmf,
integer nel )

Definition at line 1757 of file cbalke3.F.

1766C--------------------------------------------------------------------------------------------------
1767C-----------------------------------------------
1768C I M P L I C I T T Y P E S
1769C-----------------------------------------------
1770#include "implicit_f.inc"
1771#include "mvsiz_p.inc"
1772C-----------------------------------------------
1773C D U M M Y A R G U M E N T S
1774C-----------------------------------------------
1775C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1776 INTEGER JFT,JLT,IPLAT(*),IORTH,NPLAT,NG,IKGEO,NEL
1777 my_real
1778 . vol(*),thk0(*),hm(mvsiz,4),hz(*),bm(mvsiz,3,3,4),
1779 . k11(3,3,*),k12(3,3,*),k13(3,3,*),k14(3,3,*),
1780 . k22(3,3,*),k23(3,3,*),k24(3,3,*),k33(3,3,*),
1781 . m11(3,3,*),m12(3,3,*),m13(3,3,*),m14(3,3,*),
1782 . m22(3,3,*),m23(3,3,*),m24(3,3,*),m33(3,3,*),
1783 . mf11(3,3,*),mf12(3,3,*),mf13(3,3,*),mf14(3,3,*),
1784 . mf22(3,3,*),mf23(3,3,*),mf24(3,3,*),mf33(3,3,*),
1785 . fm12(3,3,*),fm13(3,3,*),fm14(3,3,*),
1786 . fm23(3,3,*),fm24(3,3,*),fm34(3,3,*),
1787 . k34(3,3,*),k44(3,3,*),m34(3,3,*),m44(3,3,*),
1788 . mf34(3,3,*),mf44(3,3,*),hmor(mvsiz,2),hfor(mvsiz,2),
1789 . pmrz(mvsiz,3,4),brz(mvsiz,4,4),frz(nel,5),
1790 . bf(mvsiz,3,2,4),hmfor(mvsiz,6),bmf(mvsiz,3,3,4)
1791C---------------|[KIJ][MFIJ]|----
1792C-----KE(6x6)= | |
1793C---------------|[FMIJ]{MIJ]|----
1794C-----------------------------------------------
1795C L O C A L V A R I A B L E S
1796C-----------------------------------------------
1797 INTEGER EP,I,J,NF,M
1798 my_real
1799 . t11(3,3,mvsiz),t12(3,3,mvsiz),t13(3,3,mvsiz),t14(3,3,mvsiz),
1800 . t22(3,3,mvsiz),t23(3,3,mvsiz),t24(3,3,mvsiz),t33(3,3,mvsiz),
1801 . t44(3,3,mvsiz),t34(3,3,mvsiz),
1802 . bb(2,4,mvsiz),gm(mvsiz),bmrz1(mvsiz,4,4),bmrz2(mvsiz,4,4),
1803 . bmrz3(mvsiz,4,4),dhz(mvsiz),
1804 . dm(3,3,mvsiz),dprz(4,mvsiz),px(4,mvsiz),py(4,mvsiz),
1805 . c1,c2,cbrx(4,mvsiz),cbry(4,mvsiz),cbrz(4,mvsiz),frzv(mvsiz)
1806 my_real
1807 . prx(4,mvsiz),pry(4,mvsiz),prxy(4,mvsiz),prz(4,mvsiz),
1808 . dmf(3,3,mvsiz),cbr1(4,mvsiz),cbr2(4,mvsiz),cbr3(4,mvsiz)
1809C-----------Attention Matrice sym Kii ne calcul que la moitie---------72
1810 nf =nplat+1
1811#include "vectorize.inc"
1812 DO m=jft,jlt
1813 ep=iplat(m)
1814 c2=vol(ep)
1815 dm(1,1,m)=hm(ep,1)*c2
1816 dm(2,2,m)=hm(ep,2)*c2
1817 dm(1,2,m)=hm(ep,3)*c2
1818 dm(3,3,m)=hm(ep,4)*c2
1819 dhz(m)=hz(ep)*fourth*c2
1820 dm(1,3,m)=zero
1821 dm(2,3,m)=zero
1822 ENDDO
1823 IF (iorth>0) THEN
1824 DO m=jft,jlt
1825 ep=iplat(m)
1826 c2=vol(ep)
1827 dm(1,3,m)=hmor(ep,1)*c2
1828 dm(2,3,m)=hmor(ep,2)*c2
1829 ENDDO
1830 ENDIF
1831C
1832 DO m=jft,jlt
1833 dm(2,1,m)=dm(1,2,m)
1834 dm(3,1,m)=dm(1,3,m)
1835 dm(3,2,m)=dm(2,3,m)
1836 ENDDO
1837C
1838 DO j=1,4
1839 DO m=jft,jlt
1840 prz(j,m)=brz(m,4,j)
1841 dprz(j,m)=prz(j,m)*dhz(m)
1842 ENDDO
1843 ENDDO
1844C-------0.5*[-By Bx 0 0 0 BRZ]^tKG[-By Bx 0 0 0 BRZ]*0.5
1845C---warped----0.5*[Pmrz1 Pmrz2 Pmrz3 0 0 BRZ]^tKG[Pmrz1 Pmrz2 Pmrz3 0 0 BRZ]*0.5
1846C------ BRZ(1-3,J,EP):Pmrz1 Pmrz2 Pmrz3 BRZ(4,J,EP)= (-BRXY+BRYX)-2NI :BB(1,J,I):-ByJ; BB(2,J,I):BxJ;
1847 DO ep=jft,nplat
1848 bb(2,1,ep)=bm(ep,1,1,1)
1849 bb(2,2,ep)=bm(ep,2,1,1)
1850 bb(2,3,ep)=bm(ep,3,1,1)
1851 bb(2,4,ep)=bm(ep,1,2,1)
1852 bb(1,1,ep)=-bm(ep,2,2,1)
1853 bb(1,2,ep)=-bm(ep,3,2,1)
1854 bb(1,3,ep)=-bm(ep,1,3,1)
1855 bb(1,4,ep)=-bm(ep,2,3,1)
1856 ENDDO
1857C
1858 DO j=1,4
1859 DO ep=jft,nplat
1860 px(j,ep) = bb(2,j,ep)
1861 py(j,ep) =-bb(1,j,ep)
1862 ENDDO
1863 ENDDO
1864C
1865 DO i=1,2
1866 DO j=i,2
1867 DO ep=jft,nplat
1868 t11(i,j,ep)=bb(i,1,ep)*bb(j,1,ep)
1869 t11(j,i,ep)=t11(i,j,ep)
1870 t22(i,j,ep)=bb(i,2,ep)*bb(j,2,ep)
1871 t22(j,i,ep)=t22(i,j,ep)
1872 t33(i,j,ep)=bb(i,3,ep)*bb(j,3,ep)
1873 t33(j,i,ep)=t33(i,j,ep)
1874 t44(i,j,ep)=bb(i,4,ep)*bb(j,4,ep)
1875 t44(j,i,ep)=t44(i,j,ep)
1876 ENDDO
1877 ENDDO
1878 ENDDO
1879C
1880 DO i=1,2
1881 DO j=1,2
1882 DO ep=jft,nplat
1883 t12(i,j,ep)=bb(i,1,ep)*bb(j,2,ep)
1884 t13(i,j,ep)=bb(i,1,ep)*bb(j,3,ep)
1885 t14(i,j,ep)=bb(i,1,ep)*bb(j,4,ep)
1886 t23(i,j,ep)=bb(i,2,ep)*bb(j,3,ep)
1887 t24(i,j,ep)=bb(i,2,ep)*bb(j,4,ep)
1888 t34(i,j,ep)=bb(i,3,ep)*bb(j,4,ep)
1889 ENDDO
1890 ENDDO
1891 ENDDO
1892C
1893 DO i=1,2
1894 DO j=i,2
1895 DO ep=jft,nplat
1896 k11(i,j,ep)=k11(i,j,ep)+t11(i,j,ep)*dhz(ep)
1897 k22(i,j,ep)=k22(i,j,ep)+t22(i,j,ep)*dhz(ep)
1898 k33(i,j,ep)=k33(i,j,ep)+t33(i,j,ep)*dhz(ep)
1899 k44(i,j,ep)=k44(i,j,ep)+t44(i,j,ep)*dhz(ep)
1900 ENDDO
1901 ENDDO
1902 ENDDO
1903C
1904 DO i=1,2
1905 DO j=1,2
1906 DO ep=jft,nplat
1907 k12(i,j,ep)=k12(i,j,ep)+t12(i,j,ep)*dhz(ep)
1908 k13(i,j,ep)=k13(i,j,ep)+t13(i,j,ep)*dhz(ep)
1909 k14(i,j,ep)=k14(i,j,ep)+t14(i,j,ep)*dhz(ep)
1910 k23(i,j,ep)=k23(i,j,ep)+t23(i,j,ep)*dhz(ep)
1911 k24(i,j,ep)=k24(i,j,ep)+t24(i,j,ep)*dhz(ep)
1912 k34(i,j,ep)=k34(i,j,ep)+t34(i,j,ep)*dhz(ep)
1913 ENDDO
1914 ENDDO
1915 ENDDO
1916C
1917 DO i=1,4
1918 DO j=1,4
1919 DO ep=jft,nplat
1920 bmrz1(ep,i,j)=bb(1,i,ep)*dprz(j,ep)
1921 bmrz2(ep,i,j)=bb(2,i,ep)*dprz(j,ep)
1922 ENDDO
1923 ENDDO
1924 ENDDO
1925C
1926 DO ep=jft,nplat
1927C------MFIJ(1,3)=BB(1,I)*BRZ(J)=BMRZ1(I,J); MFIJ(2,3)=BB(2,I)*BRZ(J)=BMRZ2(I,J)
1928C------FMIJ(3,1)=BB(1,J)*BRZ(I)=BMRZ1(J,I); FMIJ(3,2)=BB(2,J)*BRZ(I)=BMRZ2(J,I);
1929 mf11(1,3,ep)=mf11(1,3,ep)+ bmrz1(ep,1,1)
1930 mf11(2,3,ep)=mf11(2,3,ep)+ bmrz2(ep,1,1)
1931 mf22(1,3,ep)=mf22(1,3,ep)+ bmrz1(ep,2,2)
1932 mf22(2,3,ep)=mf22(2,3,ep)+ bmrz2(ep,2,2)
1933 mf33(1,3,ep)=mf33(1,3,ep)+ bmrz1(ep,3,3)
1934 mf33(2,3,ep)=mf33(2,3,ep)+ bmrz2(ep,3,3)
1935 mf44(1,3,ep)=mf44(1,3,ep)+ bmrz1(ep,4,4)
1936 mf44(2,3,ep)=mf44(2,3,ep)+ bmrz2(ep,4,4)
1937C
1938 mf12(1,3,ep)=mf12(1,3,ep)+ bmrz1(ep,1,2)
1939 mf12(2,3,ep)=mf12(2,3,ep)+ bmrz2(ep,1,2)
1940 fm12(3,1,ep)=fm12(3,1,ep)+ bmrz1(ep,2,1)
1941 fm12(3,2,ep)=fm12(3,2,ep)+ bmrz2(ep,2,1)
1942 mf13(1,3,ep)=mf13(1,3,ep)+ bmrz1(ep,1,3)
1943 mf13(2,3,ep)=mf13(2,3,ep)+ bmrz2(ep,1,3)
1944 fm13(3,1,ep)=fm13(3,1,ep)+ bmrz1(ep,3,1)
1945 fm13(3,2,ep)=fm13(3,2,ep)+ bmrz2(ep,3,1)
1946 mf14(1,3,ep)=mf14(1,3,ep)+ bmrz1(ep,1,4)
1947 mf14(2,3,ep)=mf14(2,3,ep)+ bmrz2(ep,1,4)
1948 fm14(3,1,ep)=fm14(3,1,ep)+ bmrz1(ep,4,1)
1949 fm14(3,2,ep)=fm14(3,2,ep)+ bmrz2(ep,4,1)
1950 mf23(1,3,ep)=mf23(1,3,ep)+ bmrz1(ep,2,3)
1951 mf23(2,3,ep)=mf23(2,3,ep)+ bmrz2(ep,2,3)
1952 fm23(3,1,ep)=fm23(3,1,ep)+ bmrz1(ep,3,2)
1953 fm23(3,2,ep)=fm23(3,2,ep)+ bmrz2(ep,3,2)
1954 mf24(1,3,ep)=mf24(1,3,ep)+ bmrz1(ep,2,4)
1955 mf24(2,3,ep)=mf24(2,3,ep)+ bmrz2(ep,2,4)
1956 fm24(3,1,ep)=fm24(3,1,ep)+ bmrz1(ep,4,2)
1957 fm24(3,2,ep)=fm24(3,2,ep)+ bmrz2(ep,4,2)
1958 mf34(1,3,ep)=mf34(1,3,ep)+ bmrz1(ep,3,4)
1959 mf34(2,3,ep)=mf34(2,3,ep)+ bmrz2(ep,3,4)
1960 fm34(3,1,ep)=fm34(3,1,ep)+ bmrz1(ep,4,3)
1961 fm34(3,2,ep)=fm34(3,2,ep)+ bmrz2(ep,4,3)
1962 ENDDO
1963 IF (iorth>0) THEN
1964C ----add mem-bending coupling terms of orthotrope--coplanar first---
1965#include "vectorize.inc"
1966 DO m=jft,jlt
1967 ep=iplat(m)
1968 c2=vol(ep)*thk0(ep)
1969 dmf(1,1,m)=hmfor(ep,1)*c2
1970 dmf(2,2,m)=hmfor(ep,2)*c2
1971 dmf(1,2,m)=hmfor(ep,3)*c2
1972 dmf(1,3,m)=hmfor(ep,5)*c2
1973 dmf(2,3,m)=hmfor(ep,6)*c2
1974 dmf(2,1,m)=dmf(1,2,m)
1975 dmf(3,1,m)=dmf(1,3,m)
1976 dmf(3,2,m)=dmf(2,3,m)
1977 dmf(3,3,m)=hmfor(ep,4)*c2
1978 ENDDO
1979C-------[C][BRm];----
1980 DO j=1,4
1981 DO ep=jft,jlt
1982 cbr1(j,ep) =dmf(1,1,ep)*pmrz(ep,1,j)+dmf(1,2,ep)*pmrz(ep,2,j)+
1983 . dmf(1,3,ep)*pmrz(ep,3,j)
1984 cbr2(j,ep) =dmf(2,1,ep)*pmrz(ep,1,j)+dmf(2,2,ep)*pmrz(ep,2,j)+
1985 . dmf(2,3,ep)*pmrz(ep,3,j)
1986 cbr3(j,ep) =dmf(3,1,ep)*pmrz(ep,1,j)+dmf(3,2,ep)*pmrz(ep,2,j)+
1987 . dmf(3,3,ep)*pmrz(ep,3,j)
1988 ENDDO
1989 ENDDO
1990C ----add Rz-bending coupling terms of orthotrope--coplanar first---
1991C MIJ(1,3)=-BY(I)*CBR2(J)-BX(I)*CBR3(J)
1992C MIJ(2,3)=BX(I)*CBR1(J)+BY(I)*CBR3(J)
1993C MIJ(3,1)=-BY(J)*CBR2(I)-BX(J)*CBR3(I)
1994C MIJ(3,2)=BX(J)*CBR1(I)+BY(J)*CBR3(I)
1995C
1996 DO ep=jft,nplat
1997 m11(1,3,ep)=m11(1,3,ep)
1998 . -py(1,ep)*cbr2(1,ep)-px(1,ep)*cbr3(1,ep)
1999 m11(2,3,ep)=m11(2,3,ep)
2000 . +px(1,ep)*cbr1(1,ep)+py(1,ep)*cbr3(1,ep)
2001 m22(1,3,ep)=m22(1,3,ep)
2002 . -py(2,ep)*cbr2(2,ep)-px(2,ep)*cbr3(2,ep)
2003 m22(2,3,ep)=m22(2,3,ep)
2004 . +px(2,ep)*cbr1(2,ep)+py(2,ep)*cbr3(2,ep)
2005 m33(1,3,ep)=m33(1,3,ep)
2006 . -py(3,ep)*cbr2(3,ep)-px(3,ep)*cbr3(3,ep)
2007 m33(2,3,ep)=m33(2,3,ep)
2008 . +px(3,ep)*cbr1(3,ep)+py(3,ep)*cbr3(3,ep)
2009 m44(1,3,ep)=m44(1,3,ep)
2010 . -py(4,ep)*cbr2(4,ep)-px(4,ep)*cbr3(4,ep)
2011 m44(2,3,ep)=m44(2,3,ep)
2012 . +px(4,ep)*cbr1(4,ep)+py(4,ep)*cbr3(4,ep)
2013 ENDDO
2014 DO ep=jft,nplat
2015 m12(1,3,ep)=m12(1,3,ep)
2016 . -py(1,ep)*cbr2(2,ep)-px(1,ep)*cbr3(2,ep)
2017 m12(2,3,ep)=m12(2,3,ep)
2018 . +px(1,ep)*cbr1(2,ep)+py(1,ep)*cbr3(2,ep)
2019 m13(1,3,ep)=m13(1,3,ep)
2020 . -py(1,ep)*cbr2(3,ep)-px(1,ep)*cbr3(3,ep)
2021 m13(2,3,ep)=m13(2,3,ep)
2022 . +px(1,ep)*cbr1(3,ep)+py(1,ep)*cbr3(3,ep)
2023 m14(1,3,ep)=m14(1,3,ep)
2024 . -py(1,ep)*cbr2(4,ep)-px(1,ep)*cbr3(4,ep)
2025 m14(2,3,ep)=m14(2,3,ep)
2026 . +px(1,ep)*cbr1(4,ep)+py(1,ep)*cbr3(4,ep)
2027 m23(1,3,ep)=m23(1,3,ep)
2028 . -py(2,ep)*cbr2(3,ep)-px(2,ep)*cbr3(3,ep)
2029 m23(2,3,ep)=m23(2,3,ep)
2030 . +px(2,ep)*cbr1(3,ep)+py(2,ep)*cbr3(3,ep)
2031 m24(1,3,ep)=m24(1,3,ep)
2032 . -py(2,ep)*cbr2(4,ep)-px(2,ep)*cbr3(4,ep)
2033 m24(2,3,ep)=m24(2,3,ep)
2034 . +px(2,ep)*cbr1(4,ep)+py(2,ep)*cbr3(4,ep)
2035 m34(1,3,ep)=m34(1,3,ep)
2036 . -py(3,ep)*cbr2(4,ep)-px(3,ep)*cbr3(4,ep)
2037 m34(2,3,ep)=m34(2,3,ep)
2038 . +px(3,ep)*cbr1(4,ep)+py(3,ep)*cbr3(4,ep)
2039 ENDDO
2040 DO ep=jft,nplat
2041 m12(3,1,ep)=m12(3,1,ep)
2042 . -py(2,ep)*cbr2(1,ep)-px(2,ep)*cbr3(1,ep)
2043 m12(3,2,ep)=m12(3,2,ep)
2044 . +px(2,ep)*cbr1(1,ep)+py(2,ep)*cbr3(1,ep)
2045 m13(3,1,ep)=m13(3,1,ep)
2046 . -py(3,ep)*cbr2(1,ep)-px(3,ep)*cbr3(1,ep)
2047 m13(3,2,ep)=m13(3,2,ep)
2048 . +px(3,ep)*cbr1(1,ep)+py(3,ep)*cbr3(1,ep)
2049 m14(3,1,ep)=m14(3,1,ep)
2050 . -py(4,ep)*cbr2(1,ep)-px(4,ep)*cbr3(1,ep)
2051 m14(3,2,ep)=m14(3,2,ep)
2052 . +px(4,ep)*cbr1(1,ep)+py(4,ep)*cbr3(1,ep)
2053 m23(3,1,ep)=m23(3,1,ep)
2054 . -py(3,ep)*cbr2(2,ep)-px(3,ep)*cbr3(2,ep)
2055 m23(3,2,ep)=m23(3,2,ep)
2056 . +px(3,ep)*cbr1(2,ep)+py(3,ep)*cbr3(2,ep)
2057 m24(3,1,ep)=m24(3,1,ep)
2058 . -py(4,ep)*cbr2(2,ep)-px(4,ep)*cbr3(2,ep)
2059 m24(3,2,ep)=m24(3,2,ep)
2060 . +px(4,ep)*cbr1(2,ep)+py(4,ep)*cbr3(2,ep)
2061 m34(3,1,ep)=m34(3,1,ep)
2062 . -py(4,ep)*cbr2(3,ep)-px(4,ep)*cbr3(3,ep)
2063 m34(3,2,ep)=m34(3,2,ep)
2064 . +px(4,ep)*cbr1(3,ep)+py(4,ep)*cbr3(3,ep)
2065 ENDDO
2066 ENDIF !IF (IORTH==1)
2067C--------warped elements--------
2068 DO i=1,3
2069 DO j=i,3
2070 DO ep=nf,jlt
2071 t11(i,j,ep)=brz(ep,i,1)*brz(ep,j,1)
2072 t11(j,i,ep)=t11(i,j,ep)
2073 t22(i,j,ep)=brz(ep,i,2)*brz(ep,j,2)
2074 t22(j,i,ep)=t22(i,j,ep)
2075 t33(i,j,ep)=brz(ep,i,3)*brz(ep,j,3)
2076 t33(j,i,ep)=t33(i,j,ep)
2077 t44(i,j,ep)=brz(ep,i,4)*brz(ep,j,4)
2078 t44(j,i,ep)=t44(i,j,ep)
2079 ENDDO
2080 ENDDO
2081 ENDDO
2082C
2083 DO i=1,3
2084 DO j=1,3
2085 DO ep=nf,jlt
2086 t12(i,j,ep)=brz(ep,i,1)*brz(ep,j,2)
2087 t13(i,j,ep)=brz(ep,i,1)*brz(ep,j,3)
2088 t14(i,j,ep)=brz(ep,i,1)*brz(ep,j,4)
2089 t23(i,j,ep)=brz(ep,i,2)*brz(ep,j,3)
2090 t24(i,j,ep)=brz(ep,i,2)*brz(ep,j,4)
2091 t34(i,j,ep)=brz(ep,i,3)*brz(ep,j,4)
2092 ENDDO
2093 ENDDO
2094 ENDDO
2095C
2096 DO i=1,3
2097 DO j=i,3
2098 DO ep=nf,jlt
2099 k11(i,j,ep)=k11(i,j,ep)+t11(i,j,ep)*dhz(ep)
2100 k22(i,j,ep)=k22(i,j,ep)+t22(i,j,ep)*dhz(ep)
2101 k33(i,j,ep)=k33(i,j,ep)+t33(i,j,ep)*dhz(ep)
2102 k44(i,j,ep)=k44(i,j,ep)+t44(i,j,ep)*dhz(ep)
2103 ENDDO
2104 ENDDO
2105 ENDDO
2106C
2107 DO i=1,3
2108 DO j=1,3
2109 DO ep=nf,jlt
2110 k12(i,j,ep)=k12(i,j,ep)+t12(i,j,ep)*dhz(ep)
2111 k13(i,j,ep)=k13(i,j,ep)+t13(i,j,ep)*dhz(ep)
2112 k14(i,j,ep)=k14(i,j,ep)+t14(i,j,ep)*dhz(ep)
2113 k23(i,j,ep)=k23(i,j,ep)+t23(i,j,ep)*dhz(ep)
2114 k24(i,j,ep)=k24(i,j,ep)+t24(i,j,ep)*dhz(ep)
2115 k34(i,j,ep)=k34(i,j,ep)+t34(i,j,ep)*dhz(ep)
2116 ENDDO
2117 ENDDO
2118 ENDDO
2119C
2120 DO i=1,4
2121 DO j=1,4
2122 DO ep=nf,jlt
2123 bmrz1(ep,i,j)=brz(ep,1,i)*dprz(j,ep)
2124 bmrz2(ep,i,j)=brz(ep,2,i)*dprz(j,ep)
2125 bmrz3(ep,i,j)=brz(ep,3,i)*dprz(j,ep)
2126 ENDDO
2127 ENDDO
2128 ENDDO
2129C
2130 DO ep=nf,jlt
2131C------MFIJ(1,3)=BB(1,I)*BRZ(J)=BMRZ1(I,J); MFIJ(2,3)=BB(2,I)*BRZ(J)=BMRZ2(I,J)
2132C------FMIJ(3,1)=BB(1,J)*BRZ(I)=BMRZ1(J,I); FMIJ(3,2)=BB(2,J)*BRZ(I)=BMRZ2(J,I);
2133 mf11(1,3,ep)=mf11(1,3,ep)+ bmrz1(ep,1,1)
2134 mf11(2,3,ep)=mf11(2,3,ep)+ bmrz2(ep,1,1)
2135 mf11(3,3,ep)=mf11(3,3,ep)+ bmrz3(ep,1,1)
2136
2137 mf22(1,3,ep)=mf22(1,3,ep)+ bmrz1(ep,2,2)
2138 mf22(2,3,ep)=mf22(2,3,ep)+ bmrz2(ep,2,2)
2139 mf22(3,3,ep)=mf22(3,3,ep)+ bmrz3(ep,2,2)
2140
2141 mf33(1,3,ep)=mf33(1,3,ep)+ bmrz1(ep,3,3)
2142 mf33(2,3,ep)=mf33(2,3,ep)+ bmrz2(ep,3,3)
2143 mf33(3,3,ep)=mf33(3,3,ep)+ bmrz3(ep,3,3)
2144
2145 mf44(1,3,ep)=mf44(1,3,ep)+ bmrz1(ep,4,4)
2146 mf44(2,3,ep)=mf44(2,3,ep)+ bmrz2(ep,4,4)
2147 mf44(3,3,ep)=mf44(3,3,ep)+ bmrz3(ep,4,4)
2148C
2149 mf12(1,3,ep)=mf12(1,3,ep)+ bmrz1(ep,1,2)
2150 mf12(2,3,ep)=mf12(2,3,ep)+ bmrz2(ep,1,2)
2151 mf12(3,3,ep)=mf12(3,3,ep)+ bmrz3(ep,1,2)
2152
2153 fm12(3,1,ep)=fm12(3,1,ep)+ bmrz1(ep,2,1)
2154 fm12(3,2,ep)=fm12(3,2,ep)+ bmrz2(ep,2,1)
2155 fm12(3,3,ep)=fm12(3,3,ep)+ bmrz3(ep,2,1)
2156
2157 mf13(1,3,ep)=mf13(1,3,ep)+ bmrz1(ep,1,3)
2158 mf13(2,3,ep)=mf13(2,3,ep)+ bmrz2(ep,1,3)
2159 mf13(3,3,ep)=mf13(3,3,ep)+ bmrz3(ep,1,3)
2160
2161 fm13(3,1,ep)=fm13(3,1,ep)+ bmrz1(ep,3,1)
2162 fm13(3,2,ep)=fm13(3,2,ep)+ bmrz2(ep,3,1)
2163 fm13(3,3,ep)=fm13(3,3,ep)+ bmrz3(ep,3,1)
2164
2165 mf14(1,3,ep)=mf14(1,3,ep)+ bmrz1(ep,1,4)
2166 mf14(2,3,ep)=mf14(2,3,ep)+ bmrz2(ep,1,4)
2167 mf14(3,3,ep)=mf14(3,3,ep)+ bmrz3(ep,1,4)
2168
2169 fm14(3,1,ep)=fm14(3,1,ep)+ bmrz1(ep,4,1)
2170 fm14(3,2,ep)=fm14(3,2,ep)+ bmrz2(ep,4,1)
2171 fm14(3,3,ep)=fm14(3,3,ep)+ bmrz3(ep,4,1)
2172
2173 mf23(1,3,ep)=mf23(1,3,ep)+ bmrz1(ep,2,3)
2174 mf23(2,3,ep)=mf23(2,3,ep)+ bmrz2(ep,2,3)
2175 mf23(3,3,ep)=mf23(3,3,ep)+ bmrz3(ep,2,3)
2176
2177 fm23(3,1,ep)=fm23(3,1,ep)+ bmrz1(ep,3,2)
2178 fm23(3,2,ep)=fm23(3,2,ep)+ bmrz2(ep,3,2)
2179 fm23(3,3,ep)=fm23(3,3,ep)+ bmrz3(ep,3,2)
2180
2181 mf24(1,3,ep)=mf24(1,3,ep)+ bmrz1(ep,2,4)
2182 mf24(2,3,ep)=mf24(2,3,ep)+ bmrz2(ep,2,4)
2183 mf24(3,3,ep)=mf24(3,3,ep)+ bmrz3(ep,2,4)
2184
2185 fm24(3,1,ep)=fm24(3,1,ep)+ bmrz1(ep,4,2)
2186 fm24(3,2,ep)=fm24(3,2,ep)+ bmrz2(ep,4,2)
2187 fm24(3,3,ep)=fm24(3,3,ep)+ bmrz3(ep,4,2)
2188
2189 mf34(1,3,ep)=mf34(1,3,ep)+ bmrz1(ep,3,4)
2190 mf34(2,3,ep)=mf34(2,3,ep)+ bmrz2(ep,3,4)
2191 mf34(3,3,ep)=mf34(3,3,ep)+ bmrz3(ep,3,4)
2192
2193 fm34(3,1,ep)=fm34(3,1,ep)+ bmrz1(ep,4,3)
2194 fm34(3,2,ep)=fm34(3,2,ep)+ bmrz2(ep,4,3)
2195 fm34(3,3,ep)=fm34(3,3,ep)+ bmrz3(ep,4,3)
2196 ENDDO
2197C------ Mzz -------
2198 DO m=jft,jlt
2199 m11(3,3,m)=m11(3,3,m)+prz(1,m)*dprz(1,m)
2200 m12(3,3,m)=m12(3,3,m)+prz(1,m)*dprz(2,m)
2201 m13(3,3,m)=m13(3,3,m)+prz(1,m)*dprz(3,m)
2202 m14(3,3,m)=m14(3,3,m)+prz(1,m)*dprz(4,m)
2203 m22(3,3,m)=m22(3,3,m)+prz(2,m)*dprz(2,m)
2204 m23(3,3,m)=m23(3,3,m)+prz(2,m)*dprz(3,m)
2205 m24(3,3,m)=m24(3,3,m)+prz(2,m)*dprz(4,m)
2206 m33(3,3,m)=m33(3,3,m)+prz(3,m)*dprz(3,m)
2207 m34(3,3,m)=m34(3,3,m)+prz(3,m)*dprz(4,m)
2208 m44(3,3,m)=m44(3,3,m)+prz(4,m)*dprz(4,m)
2209 ENDDO
2210C
2211C-------[MFIJ]=[Bm]^t[C][BRm]; [MIJ]=[BRm]^t[C][BRm];----
2212C---------------|0 0 BRX |----
2213C-----BR(3x3)= |0 0 BRY |
2214C---------------|0 0 BRXY+BRYX |----
2215C-- [MFIJ]=|Bxi 0 Byi|{CPRXj} [FMIJ]={CPRXi CPRYi CPRZi}|Bxj 0 0|
2216C-- |0 Byi Bxi|{CPRYj} |0 Byj 0|
2217C-- |0 0 0 |{CPRZj} |ByjBxj 0|
2218C-- warped [MFIJ]=[BMi]^t{CPRXj} [FMIJ]={CPRXi CPRYi CPRZi}[Bmj](3,3)
2219C
2220 DO m=jft,jlt
2221 DO j=1,4
2222 prx(j,m)=pmrz(m,1,j)
2223 pry(j,m)=pmrz(m,2,j)
2224 prxy(j,m)=pmrz(m,3,j)
2225 ENDDO
2226 ENDDO
2227C-------[C][BRm];----
2228 DO j=1,4
2229 DO m=jft,jlt
2230 cbrx(j,m) =dm(1,1,m)*prx(j,m)+dm(1,2,m)*pry(j,m)+
2231 . dm(1,3,m)*prxy(j,m)
2232 cbry(j,m) =dm(2,1,m)*prx(j,m)+dm(2,2,m)*pry(j,m)+
2233 . dm(2,3,m)*prxy(j,m)
2234 cbrz(j,m) =dm(3,1,m)*prx(j,m)+dm(3,2,m)*pry(j,m)+
2235 . dm(3,3,m)*prxy(j,m)
2236 ENDDO
2237 ENDDO
2238 DO m=jft,nplat
2239 mf11(1,3,m)=mf11(1,3,m)+ px(1,m)*cbrx(1,m)+py(1,m)*cbrz(1,m)
2240 mf11(2,3,m)=mf11(2,3,m)+ py(1,m)*cbry(1,m)+px(1,m)*cbrz(1,m)
2241 mf12(1,3,m)=mf12(1,3,m)+ px(1,m)*cbrx(2,m)+py(1,m)*cbrz(2,m)
2242 mf12(2,3,m)=mf12(2,3,m)+ py(1,m)*cbry(2,m)+px(1,m)*cbrz(2,m)
2243 mf13(1,3,m)=mf13(1,3,m)+ px(1,m)*cbrx(3,m)+py(1,m)*cbrz(3,m)
2244 mf13(2,3,m)=mf13(2,3,m)+ py(1,m)*cbry(3,m)+px(1,m)*cbrz(3,m)
2245 mf14(1,3,m)=mf14(1,3,m)+ px(1,m)*cbrx(4,m)+py(1,m)*cbrz(4,m)
2246 mf14(2,3,m)=mf14(2,3,m)+ py(1,m)*cbry(4,m)+px(1,m)*cbrz(4,m)
2247 mf22(1,3,m)=mf22(1,3,m)+ px(2,m)*cbrx(2,m)+py(2,m)*cbrz(2,m)
2248 mf22(2,3,m)=mf22(2,3,m)+ py(2,m)*cbry(2,m)+px(2,m)*cbrz(2,m)
2249 mf23(1,3,m)=mf23(1,3,m)+ px(2,m)*cbrx(3,m)+py(2,m)*cbrz(3,m)
2250 mf23(2,3,m)=mf23(2,3,m)+ py(2,m)*cbry(3,m)+px(2,m)*cbrz(3,m)
2251 mf24(1,3,m)=mf24(1,3,m)+ px(2,m)*cbrx(4,m)+py(2,m)*cbrz(4,m)
2252 mf24(2,3,m)=mf24(2,3,m)+ py(2,m)*cbry(4,m)+px(2,m)*cbrz(4,m)
2253 mf33(1,3,m)=mf33(1,3,m)+ px(3,m)*cbrx(3,m)+py(3,m)*cbrz(3,m)
2254 mf33(2,3,m)=mf33(2,3,m)+ py(3,m)*cbry(3,m)+px(3,m)*cbrz(3,m)
2255 mf34(1,3,m)=mf34(1,3,m)+ px(3,m)*cbrx(4,m)+py(3,m)*cbrz(4,m)
2256 mf34(2,3,m)=mf34(2,3,m)+ py(3,m)*cbry(4,m)+px(3,m)*cbrz(4,m)
2257 mf44(1,3,m)=mf44(1,3,m)+ px(4,m)*cbrx(4,m)+py(4,m)*cbrz(4,m)
2258 mf44(2,3,m)=mf44(2,3,m)+ py(4,m)*cbry(4,m)+px(4,m)*cbrz(4,m)
2259C
2260 fm12(3,1,m)=fm12(3,1,m)+ px(2,m)*cbrx(1,m)+py(2,m)*cbrz(1,m)
2261 fm12(3,2,m)=fm12(3,2,m)+ py(2,m)*cbry(1,m)+px(2,m)*cbrz(1,m)
2262 fm13(3,1,m)=fm13(3,1,m)+ px(3,m)*cbrx(1,m)+py(3,m)*cbrz(1,m)
2263 fm13(3,2,m)=fm13(3,2,m)+ py(3,m)*cbry(1,m)+px(3,m)*cbrz(1,m)
2264 fm14(3,1,m)=fm14(3,1,m)+ px(4,m)*cbrx(1,m)+py(4,m)*cbrz(1,m)
2265 fm14(3,2,m)=fm14(3,2,m)+ py(4,m)*cbry(1,m)+px(4,m)*cbrz(1,m)
2266 fm23(3,1,m)=fm23(3,1,m)+ px(3,m)*cbrx(2,m)+py(3,m)*cbrz(2,m)
2267 fm23(3,2,m)=fm23(3,2,m)+ py(3,m)*cbry(2,m)+px(3,m)*cbrz(2,m)
2268 fm24(3,1,m)=fm24(3,1,m)+ px(4,m)*cbrx(2,m)+py(4,m)*cbrz(2,m)
2269 fm24(3,2,m)=fm24(3,2,m)+ py(4,m)*cbry(2,m)+px(4,m)*cbrz(2,m)
2270 fm34(3,1,m)=fm34(3,1,m)+ px(4,m)*cbrx(3,m)+py(4,m)*cbrz(3,m)
2271 fm34(3,2,m)=fm34(3,2,m)+ py(4,m)*cbry(3,m)+px(4,m)*cbrz(3,m)
2272 ENDDO
2273 DO ep=nf,jlt
2274 DO i=1,3
2275 mf11(i,3,ep)=mf11(i,3,ep)+ bm(ep,1,i,1)*cbrx(1,ep)+
2276 . bm(ep,2,i,1)*cbry(1,ep)+bm(ep,3,i,1)*cbrz(1,ep)
2277 mf12(i,3,ep)=mf12(i,3,ep)+ bm(ep,1,i,1)*cbrx(2,ep)+
2278 . bm(ep,2,i,1)*cbry(2,ep)+bm(ep,3,i,1)*cbrz(2,ep)
2279 mf13(i,3,ep)=mf13(i,3,ep)+ bm(ep,1,i,1)*cbrx(3,ep)+
2280 . bm(ep,2,i,1)*cbry(3,ep)+bm(ep,3,i,1)*cbrz(3,ep)
2281 mf14(i,3,ep)=mf14(i,3,ep)+ bm(ep,1,i,1)*cbrx(4,ep)+
2282 . bm(ep,2,i,1)*cbry(4,ep)+bm(ep,3,i,1)*cbrz(4,ep)
2283 mf22(i,3,ep)=mf22(i,3,ep)+ bm(ep,1,i,2)*cbrx(2,ep)+
2284 . bm(ep,2,i,2)*cbry(2,ep)+bm(ep,3,i,2)*cbrz(2,ep)
2285 mf23(i,3,ep)=mf23(i,3,ep)+ bm(ep,1,i,2)*cbrx(3,ep)+
2286 . bm(ep,2,i,2)*cbry(3,ep)+bm(ep,3,i,2)*cbrz(3,ep)
2287 mf24(i,3,ep)=mf24(i,3,ep)+ bm(ep,1,i,2)*cbrx(4,ep)+
2288 . bm(ep,2,i,2)*cbry(4,ep)+bm(ep,3,i,2)*cbrz(4,ep)
2289 mf33(i,3,ep)=mf33(i,3,ep)+ bm(ep,1,i,3)*cbrx(3,ep)+
2290 . bm(ep,2,i,3)*cbry(3,ep)+bm(ep,3,i,3)*cbrz(3,ep)
2291 mf34(i,3,ep)=mf34(i,3,ep)+ bm(ep,1,i,3)*cbrx(4,ep)+
2292 . bm(ep,2,i,3)*cbry(4,ep)+bm(ep,3,i,3)*cbrz(4,ep)
2293 mf44(i,3,ep)=mf44(i,3,ep)+ bm(ep,1,i,4)*cbrx(4,ep)+
2294 . bm(ep,2,i,4)*cbry(4,ep)+bm(ep,3,i,4)*cbrz(4,ep)
2295C
2296 fm12(3,i,ep)=fm12(3,i,ep)+ bm(ep,1,i,2)*cbrx(1,ep)+
2297 . bm(ep,2,i,2)*cbry(1,ep)+bm(ep,3,i,2)*cbrz(1,ep)
2298 fm13(3,i,ep)=fm13(3,i,ep)+ bm(ep,1,i,3)*cbrx(1,ep)+
2299 . bm(ep,2,i,3)*cbry(1,ep)+bm(ep,3,i,3)*cbrz(1,ep)
2300 fm14(3,i,ep)=fm14(3,i,ep)+ bm(ep,1,i,4)*cbrx(1,ep)+
2301 . bm(ep,2,i,4)*cbry(1,ep)+bm(ep,3,i,4)*cbrz(1,ep)
2302 fm23(3,i,ep)=fm23(3,i,ep)+ bm(ep,1,i,3)*cbrx(2,ep)+
2303 . bm(ep,2,i,3)*cbry(2,ep)+bm(ep,3,i,3)*cbrz(2,ep)
2304 fm24(3,i,ep)=fm24(3,i,ep)+ bm(ep,1,i,4)*cbrx(2,ep)+
2305 . bm(ep,2,i,4)*cbry(2,ep)+bm(ep,3,i,4)*cbrz(2,ep)
2306 fm34(3,i,ep)=fm34(3,i,ep)+ bm(ep,1,i,4)*cbrx(3,ep)+
2307 . bm(ep,2,i,4)*cbry(3,ep)+bm(ep,3,i,4)*cbrz(3,ep)
2308 END DO !I=1,3
2309 ENDDO
2310C------ Mzz ----------
2311 DO m=jft,jlt
2312 m11(3,3,m)=m11(3,3,m)+prx(1,m)*cbrx(1,m)+
2313 . pry(1,m)*cbry(1,m)+prxy(1,m)*cbrz(1,m)
2314 m12(3,3,m)=m12(3,3,m)+prx(1,m)*cbrx(2,m)+
2315 . pry(1,m)*cbry(2,m)+prxy(1,m)*cbrz(2,m)
2316 m13(3,3,m)=m13(3,3,m)+prx(1,m)*cbrx(3,m)+
2317 . pry(1,m)*cbry(3,m)+prxy(1,m)*cbrz(3,m)
2318 m14(3,3,m)=m14(3,3,m)+prx(1,m)*cbrx(4,m)+
2319 . pry(1,m)*cbry(4,m)+prxy(1,m)*cbrz(4,m)
2320 m22(3,3,m)=m22(3,3,m)+prx(2,m)*cbrx(2,m)+
2321 . pry(2,m)*cbry(2,m)+prxy(2,m)*cbrz(2,m)
2322 m23(3,3,m)=m23(3,3,m)+prx(2,m)*cbrx(3,m)+
2323 . pry(2,m)*cbry(3,m)+prxy(2,m)*cbrz(3,m)
2324 m24(3,3,m)=m24(3,3,m)+prx(2,m)*cbrx(4,m)+
2325 . pry(2,m)*cbry(4,m)+prxy(2,m)*cbrz(4,m)
2326 m33(3,3,m)=m33(3,3,m)+prx(3,m)*cbrx(3,m)+
2327 . pry(3,m)*cbry(3,m)+prxy(3,m)*cbrz(3,m)
2328 m34(3,3,m)=m34(3,3,m)+prx(3,m)*cbrx(4,m)+
2329 . pry(3,m)*cbry(4,m)+prxy(3,m)*cbrz(4,m)
2330 m44(3,3,m)=m44(3,3,m)+prx(4,m)*cbrx(4,m)+
2331 . pry(4,m)*cbry(4,m)+prxy(4,m)*cbrz(4,m)
2332 ENDDO
2333C
2334 IF (iorth>0) THEN
2335C ----add Rz-bending coupling terms of orthotrope--warped element---
2336C MIJ(i,3)=BF(I)^t*CBRi(J)
2337C MIJ(3,i)=BF(J)^t*CBRi(I)
2338C FMIJ(i,3)=BMF(I)^t*CBRi(J)
2339C MFIJ(3,i)=BMF(J)^t*CBRi(I)
2340C
2341 DO i=1,2
2342 DO ep=nf,jlt
2343 m11(i,3,ep)=m11(i,3,ep)
2344 . +bf(ep,1,i,1)*cbr1(1,ep)+bf(ep,2,i,1)*cbr2(1,ep)
2345 . +bf(ep,3,i,1)*cbr3(1,ep)
2346 m22(i,3,ep)=m22(i,3,ep)
2347 . +bf(ep,1,i,2)*cbr1(2,ep)+bf(ep,2,i,2)*cbr2(2,ep)
2348 . +bf(ep,3,i,2)*cbr3(2,ep)
2349 m33(i,3,ep)=m33(i,3,ep)
2350 . +bf(ep,1,i,3)*cbr1(3,ep)+bf(ep,2,i,3)*cbr2(3,ep)
2351 . +bf(ep,3,i,3)*cbr3(3,ep)
2352 m44(i,3,ep)=m44(i,3,ep)
2353 . +bf(ep,1,i,4)*cbr1(4,ep)+bf(ep,2,i,4)*cbr2(4,ep)
2354 . +bf(ep,3,i,4)*cbr3(4,ep)
2355 ENDDO
2356 ENDDO
2357C
2358 DO i=1,2
2359 DO ep=nf,jlt
2360 m12(i,3,ep)=m12(i,3,ep)
2361 . +bf(ep,1,i,1)*cbr1(2,ep)+bf(ep,2,i,1)*cbr2(2,ep)
2362 . +bf(ep,3,i,1)*cbr3(2,ep)
2363 m13(i,3,ep)=m13(i,3,ep)
2364 . +bf(ep,1,i,1)*cbr1(3,ep)+bf(ep,2,i,1)*cbr2(3,ep)
2365 . +bf(ep,3,i,1)*cbr3(3,ep)
2366 m14(i,3,ep)=m14(i,3,ep)
2367 . +bf(ep,1,i,1)*cbr1(4,ep)+bf(ep,2,i,1)*cbr2(4,ep)
2368 . +bf(ep,3,i,1)*cbr3(4,ep)
2369 m23(i,3,ep)=m23(i,3,ep)
2370 . +bf(ep,1,i,2)*cbr1(3,ep)+bf(ep,2,i,2)*cbr2(3,ep)
2371 . +bf(ep,3,i,2)*cbr3(3,ep)
2372 m24(i,3,ep)=m24(i,3,ep)
2373 . +bf(ep,1,i,2)*cbr1(4,ep)+bf(ep,2,i,2)*cbr2(4,ep)
2374 . +bf(ep,3,i,2)*cbr3(4,ep)
2375 m34(i,3,ep)=m34(i,3,ep)
2376 . +bf(ep,1,i,3)*cbr1(4,ep)+bf(ep,2,i,3)*cbr2(4,ep)
2377 . +bf(ep,3,i,3)*cbr3(4,ep)
2378 ENDDO
2379 ENDDO
2380 DO i=1,2
2381 DO ep=nf,jlt
2382 m12(3,i,ep)=m12(3,i,ep)
2383 . +bf(ep,1,i,2)*cbr1(1,ep)+bf(ep,2,i,2)*cbr2(1,ep)
2384 . +bf(ep,3,i,2)*cbr3(1,ep)
2385 m13(3,i,ep)=m13(3,i,ep)
2386 . +bf(ep,1,i,3)*cbr1(1,ep)+bf(ep,2,i,3)*cbr2(1,ep)
2387 . +bf(ep,3,i,3)*cbr3(1,ep)
2388 m14(3,i,ep)=m14(3,i,ep)
2389 . +bf(ep,1,i,4)*cbr1(1,ep)+bf(ep,2,i,4)*cbr2(1,ep)
2390 . +bf(ep,3,i,4)*cbr3(1,ep)
2391 m23(3,i,ep)=m23(3,i,ep)
2392 . +bf(ep,1,i,3)*cbr1(2,ep)+bf(ep,2,i,3)*cbr2(2,ep)
2393 . +bf(ep,3,i,3)*cbr3(2,ep)
2394 m24(3,i,ep)=m24(3,i,ep)
2395 . +bf(ep,1,i,4)*cbr1(2,ep)+bf(ep,2,i,4)*cbr2(2,ep)
2396 . +bf(ep,3,i,4)*cbr3(2,ep)
2397 m34(3,i,ep)=m34(3,i,ep)
2398 . +bf(ep,1,i,4)*cbr1(3,ep)+bf(ep,2,i,4)*cbr2(3,ep)
2399 . +bf(ep,3,i,4)*cbr3(3,ep)
2400 ENDDO
2401 ENDDO
2402 DO i=1,3
2403 DO ep=nf,jlt
2404 mf11(i,3,ep)=mf11(i,3,ep)
2405 . +bmf(ep,1,i,1)*cbr1(1,ep)+bmf(ep,2,i,1)*cbr2(1,ep)
2406 . +bmf(ep,3,i,1)*cbr3(1,ep)
2407 mf12(i,3,ep)=mf12(i,3,ep)
2408 . +bmf(ep,1,i,1)*cbr1(2,ep)+bmf(ep,2,i,1)*cbr2(2,ep)
2409 . +bmf(ep,3,i,1)*cbr3(2,ep)
2410 mf13(i,3,ep)=mf13(i,3,ep)
2411 . +bmf(ep,1,i,1)*cbr1(3,ep)+bmf(ep,2,i,1)*cbr2(3,ep)
2412 . +bmf(ep,3,i,1)*cbr3(3,ep)
2413 mf14(i,3,ep)=mf14(i,3,ep)
2414 . +bmf(ep,1,i,1)*cbr1(4,ep)+bmf(ep,2,i,1)*cbr2(4,ep)
2415 . +bmf(ep,3,i,1)*cbr3(4,ep)
2416 mf22(i,3,ep)=mf22(i,3,ep)
2417 . +bmf(ep,1,i,2)*cbr1(2,ep)+bmf(ep,2,i,2)*cbr2(2,ep)
2418 . +bmf(ep,3,i,2)*cbr3(2,ep)
2419 mf23(i,3,ep)=mf23(i,3,ep)
2420 . +bmf(ep,1,i,2)*cbr1(3,ep)+bmf(ep,2,i,2)*cbr2(3,ep)
2421 . +bmf(ep,3,i,2)*cbr3(3,ep)
2422 mf24(i,3,ep)=mf24(i,3,ep)
2423 . +bmf(ep,1,i,2)*cbr1(4,ep)+bmf(ep,2,i,2)*cbr2(4,ep)
2424 . +bmf(ep,3,i,2)*cbr3(4,ep)
2425 mf33(i,3,ep)=mf33(i,3,ep)
2426 . +bmf(ep,1,i,3)*cbr1(3,ep)+bmf(ep,2,i,3)*cbr2(3,ep)
2427 . +bmf(ep,3,i,3)*cbr3(3,ep)
2428 mf34(i,3,ep)=mf34(i,3,ep)
2429 . +bmf(ep,1,i,3)*cbr1(4,ep)+bmf(ep,2,i,3)*cbr2(4,ep)
2430 . +bmf(ep,3,i,3)*cbr3(4,ep)
2431 mf44(i,3,ep)=mf44(i,3,ep)
2432 . +bmf(ep,1,i,4)*cbr1(4,ep)+bmf(ep,2,i,4)*cbr2(4,ep)
2433 . +bmf(ep,3,i,4)*cbr3(4,ep)
2434C
2435 fm12(3,i,ep)=fm12(3,i,ep)
2436 . +bmf(ep,1,i,2)*cbr1(1,ep)+bmf(ep,2,i,2)*cbr2(1,ep)
2437 . +bmf(ep,3,i,2)*cbr3(1,ep)
2438 fm13(3,i,ep)=fm13(3,i,ep)
2439 . +bmf(ep,1,i,3)*cbr1(1,ep)+bmf(ep,2,i,3)*cbr2(1,ep)
2440 . +bmf(ep,3,i,3)*cbr3(1,ep)
2441 fm14(3,i,ep)=fm14(3,i,ep)
2442 . +bmf(ep,1,i,4)*cbr1(1,ep)+bmf(ep,2,i,4)*cbr2(1,ep)
2443 . +bmf(ep,3,i,4)*cbr3(1,ep)
2444 fm23(3,i,ep)=fm23(3,i,ep)
2445 . +bmf(ep,1,i,3)*cbr1(2,ep)+bmf(ep,2,i,3)*cbr2(2,ep)
2446 . +bmf(ep,3,i,3)*cbr3(2,ep)
2447 fm24(3,i,ep)=fm24(3,i,ep)
2448 . +bmf(ep,1,i,4)*cbr1(2,ep)+bmf(ep,2,i,4)*cbr2(2,ep)
2449 . +bmf(ep,3,i,4)*cbr3(2,ep)
2450 fm34(3,i,ep)=fm34(3,i,ep)
2451 . +bmf(ep,1,i,4)*cbr1(3,ep)+bmf(ep,2,i,4)*cbr2(3,ep)
2452 . +bmf(ep,3,i,4)*cbr3(3,ep)
2453 ENDDO
2454 END DO !I=1,3
2455 END IF !(IORTH==1) THEN
2456C------ Geometrical stiffness associated w/ Rz ----------
2457 IF (ikgeo==1) THEN
2458#include "vectorize.inc"
2459 DO m=jft,jlt
2460 ep=iplat(m)
2461 frzv(m)=frz(ep,ng)*vol(ep)
2462 ENDDO
2463 DO ep=jft,jlt
2464 t11(1,2,ep)=frzv(ep)*(prx(1,ep)*prx(1,ep)+pry(1,ep)*pry(1,ep)+
2465 . prxy(1,ep)*prxy(1,ep))
2466 t22(1,2,ep)=frzv(ep)*(prx(2,ep)*prx(2,ep)+pry(2,ep)*pry(2,ep)+
2467 . prxy(2,ep)*prxy(2,ep))
2468 t33(1,2,ep)=frzv(ep)*(prx(3,ep)*prx(3,ep)+pry(3,ep)*pry(3,ep)+
2469 . prxy(3,ep)*prxy(3,ep))
2470 t44(1,2,ep)=frzv(ep)*(prx(4,ep)*prx(4,ep)+pry(4,ep)*pry(4,ep)+
2471 . prxy(4,ep)*prxy(4,ep))
2472 t12(1,2,ep)=frzv(ep)*(prx(1,ep)*prx(2,ep)+pry(1,ep)*pry(2,ep)+
2473 . prxy(1,ep)*prxy(2,ep))
2474 t13(1,2,ep)=frzv(ep)*(prx(1,ep)*prx(3,ep)+pry(1,ep)*pry(3,ep)+
2475 . prxy(1,ep)*prxy(3,ep))
2476 t14(1,2,ep)=frzv(ep)*(prx(1,ep)*prx(4,ep)+pry(1,ep)*pry(4,ep)+
2477 . prxy(1,ep)*prxy(4,ep))
2478 t23(1,2,ep)=frzv(ep)*(prx(2,ep)*prx(3,ep)+pry(2,ep)*pry(3,ep)+
2479 . prxy(2,ep)*prxy(3,ep))
2480 t24(1,2,ep)=frzv(ep)*(prx(2,ep)*prx(4,ep)+pry(2,ep)*pry(4,ep)+
2481 . prxy(2,ep)*prxy(4,ep))
2482 t34(1,2,ep)=frzv(ep)*(prx(3,ep)*prx(4,ep)+pry(3,ep)*pry(4,ep)+
2483 . prxy(3,ep)*prxy(4,ep))
2484 ENDDO
2485C
2486 DO i=1,3
2487 DO ep=jft,jlt
2488 k11(i,i,ep)=k11(i,i,ep)+t11(1,2,ep)
2489 k22(i,i,ep)=k22(i,i,ep)+t22(1,2,ep)
2490 k33(i,i,ep)=k33(i,i,ep)+t33(1,2,ep)
2491 k44(i,i,ep)=k44(i,i,ep)+t44(1,2,ep)
2492 k12(i,i,ep)=k12(i,i,ep)+t12(1,2,ep)
2493 k13(i,i,ep)=k13(i,i,ep)+t13(1,2,ep)
2494 k14(i,i,ep)=k14(i,i,ep)+t14(1,2,ep)
2495 k23(i,i,ep)=k23(i,i,ep)+t23(1,2,ep)
2496 k24(i,i,ep)=k24(i,i,ep)+t24(1,2,ep)
2497 k34(i,i,ep)=k34(i,i,ep)+t34(1,2,ep)
2498 ENDDO
2499 ENDDO
2500 ENDIF
2501C
2502 RETURN