58
59
60
64
65
66
67#include "implicit_f.inc"
68#include "comlock.inc"
69
70
71
72#include "mvsiz_p.inc"
73
74
75
76#include "com01_c.inc"
77#include "com08_c.inc"
78#include "com04_c.inc"
79#include "com06_c.inc"
80#include "scr07_c.inc"
81#include "scr14_c.inc"
82#include "scr16_c.inc"
83#include "scr18_c.inc"
84#include "sms_c.inc"
85#include "parit_c.inc"
86#include "param_c.inc"
87#include "impl1_c.inc"
88
89
90
91 INTEGER JLT,ITIED, NISKYFI,NIN,NELTST,ITYPTST,JTASK
92 INTEGER ISKY(*),ISECIN, NSTRF(*)
93 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
94 . NSVG(MVSIZ), NSMS(MVSIZ), INDEX(*),
95 . NISUB, LISUB(*), ADDSUBS(*), ADDSUBM(*), LISUBS(*),
96 . LISUBM(*),CN_LOC(*),CE_LOC(*),ISKYI_SMS(*),ICONTACT(*),
97 . ISENSINT(*),NFT
98 INTEGER ,INTENT(IN) :: NODADT_THERM
100 . a(3,*), ms(*), fsav(*),x1(*),x2(*),x3(*),x4(*),
101 . y1(*),y2(*),y3(*),y4(*),z1(*),z2(*),z3(*),z4(*),
102 . visc,stifn(*),cand_f(6,*), v(3,*),fskyi(lskyi,nfskyi),
103 . fcont(3,*),
104 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz)
106 . nx1(mvsiz), nx2(mvsiz), nx3(mvsiz), nx4(mvsiz),
107 . ny1(mvsiz), ny2(mvsiz), ny3(mvsiz), ny4(mvsiz),
108 . nz1(mvsiz), nz2(mvsiz), nz3(mvsiz), nz4(mvsiz),
109 . lb1(mvsiz), lb2(mvsiz), lb3(mvsiz), lb4(mvsiz),
110 . lc1(mvsiz), lc2(mvsiz), lc3(mvsiz), lc4(mvsiz),
111 . p1(mvsiz), p2(mvsiz), p3(mvsiz), p4(mvsiz), stif(mvsiz),
112 . gapv(mvsiz),
113 . secfcum(7,numnod,nsect), viscn(*),fsavsub(nthvki,*),
114 . fncont(3,*), ftcont(3,*), mskyi_sms(*),
115 . xi(mvsiz),yi(mvsiz),zi(mvsiz),dt2t,
116 . fsavparit(nisub+1,11,*)
117 TYPE(H3D_DATABASE) :: H3D_DATA
118
119
120
121 INTEGER I, J1, IG, II , K0,NBINTER,K1S,K,J,NN,JG
122 INTEGER NISKYL,NISKYL1,NOINT
123 INTEGER JSUB,KSUB,JJ,KK,IN,IE,NSUB,IBID
125 . fsavsub1(24,nisub),impx,impy,impz
127 . fxi(mvsiz), fyi(mvsiz), fzi(mvsiz), fni(mvsiz),
128 . fx1(mvsiz), fx2(mvsiz), fx3(mvsiz), fx4(mvsiz),
129 . fy1(mvsiz), fy2(mvsiz), fy3(mvsiz), fy4(mvsiz),
130 . fz1(mvsiz), fz2(mvsiz), fz3(mvsiz), fz4(mvsiz),
131 . ft1(mvsiz), ft2(mvsiz),
132 . n1(mvsiz), n2(mvsiz), n3(mvsiz), pene(mvsiz),
133 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
134 . vt1(mvsiz), vt2(mvsiz),
135 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
136 . t1x(mvsiz),t1y(mvsiz),t1z(mvsiz),
137 . t2x(mvsiz),t2y(mvsiz),t2z(mvsiz),
138 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
139 . xp(mvsiz), yp(mvsiz), zp(mvsiz),
140 . s2,d1,d2,d3,d4,a1,a2,a3,a4,la1,la2,la3,la4,h0,
141 . econtt, dt05, norminv, vis, dt1inv,econtdt,
142 . fsav1, fsav2, fsav3,fsav4 , fsav5, fsav6, fsav7, fsav8,
143 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15,
144 . fsav22, fsav23, fsav24
146 . c(mvsiz),c1(mvsiz),c2(mvsiz),c3(mvsiz),c4(mvsiz),
147 . k1(mvsiz),k2(mvsiz),k3(mvsiz),k4(mvsiz),
148 . fxn(mvsiz), fyn(mvsiz), fzn(mvsiz),
149 . fxt(mvsiz), fyt(mvsiz), fzt(mvsiz),bid ,dti
150
151 IF(dt1>zero)THEN
152 dt1inv = one/dt1
153 ELSE
154 dt1inv =zero
155 ENDIF
156 econtt = zero
157 econtdt = zero
158
159 dt05 = half * dt1
160 ibid = 0
161 bid = zero
162
163
164
165 DO i=1,jlt
166
167 d1 = sqrt(p1(i))
168 p1(i) =
max(zero, gapv(i) - d1)
169
170 d2 = sqrt(p2(i))
171 p2(i) =
max(zero, gapv(i) - d2)
172
173 d3 = sqrt(p3(i))
174 p3(i) =
max(zero, gapv(i) - d3)
175
176 d4 = sqrt(p4(i))
177 p4(i) =
max(zero, gapv(i) - d4)
178
179 a1 = p1(i)/
max(em20,d1)
180 a2 = p2(i)/
max(em20,d2)
181 a3 = p3(i)/
max(em20,d3)
182 a4 = p4(i)/
max(em20,d4)
183 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
184 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
185 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
186 ENDDO
187
188 DO i=1,jlt
189 IF(ix3(i)/=ix4(i))THEN
190 pene(i) =
max(p1(i),p2(i),p3(i),p4(i))
191
192 la1 = one - lb1(i) - lc1(i)
193 la2 = one - lb2(i) - lc2(i)
194 la3 = one - lb3(i) - lc3(i)
195 la4 = one - lb4(i) - lc4(i)
196
197 h0 = fourth *
198 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
199 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
200 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
201 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
202 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
203
204 h0 = one/
max(em20,h1(i) + h2(i) + h3(i) + h4(i))
205 h1(i) = h1(i) * h0
206 h2(i) = h2(i) * h0
207 h3(i) = h3(i) * h0
208 h4(i) = h4(i) * h0
209
210 ELSE
211 pene(i) = p1(i
212 n1(i) = nx1(i)
213 n2(i) = ny1(i)
214 n3(i) = nz1(i)
215 h1(i) = lb1(i)
216 h2(i) = lc1(i)
217 h3(i) = one - lb1(i) - lc1(i)
218 h4(i) = zero
219 ENDIF
220 ENDDO
221
222 DO i=1,jlt
223 s2 = one/
max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
224 n1(i) = n1(i)*s2
225 n2(i) = n2(i)*s2
226 n3(i) = n3(i)*s2
227 ENDDO
228
229
230 DO i=1,jlt
231
232 IF(ix3(i)/=ix4(i))THEN
233 h0 = -fourth*(h1(i) - h2(i) + h3(i) - h4(i))
234 h0 =
min(h0,h2(i),h4(i))
235 h0 =
max(h0,-h1(i),-h3(i))
236 h1(i) = h1(i) + h0
237 h2(i) = h2(i) - h0
238 h3(i) = h3(i) + h0
239 h4(i) = h4(i) - h0
240 ENDIF
241 ENDDO
242
243 DO 5 i=1,jlt
244 ii = index(i)
245 IF(cand_f(1,ii)==zero)THEN
246
247
248
249 cand_f(4,ii) = h1(i)
250 cand_f(5,ii) = h2(i)
251 cand_f(6,ii) = h3(i)
252 ELSE
253
254
255
256 h1(i) = cand_f(4,ii)
257 h2(i) = cand_f(5,ii)
258 h3(i) = cand_f(6,ii)
259 h4(i) = one - h1(i) - h2(i) - h3(i)
260 ENDIF
261 5 CONTINUE
262
263 DO 10 i=1,jlt
264 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
265 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
266 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
267 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
268 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
269 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
270 10 CONTINUE
271
272
273 DO 20 i=1,jlt
274 t1x(i) = x3(i) - x1(i)
275 t1y(i) = y3(i) - y1(i)
276 t1z(i) = z3(i) - z1(i)
277 norminv = one/sqrt(t1x(i)**2+t1y(i)**2+t1z(i)**2)
278 t1x(i) = t1x(i)*norminv
279 t1y(i) = t1y(i)*norminv
280 t1z(i) = t1z(i)*norminv
281
282 t2x(i) = x4(i) - x2(i)
283 t2y(i) = y4(i) - y2(i)
284 t2z(i) = z4(i) - z2(i)
285
286 nx(i) = t1y(i)*t2z(i) - t1z(i)*t2y(i)
287 ny(i) = t1z(i)*t2x(i) - t1x(i)*t2z(i)
288 nz(i) = t1x(i)*t2y(i) - t1y(i)*t2x(i)
289 norminv = one/sqrt(nx(i)**2+ny(i)**2+nz(i)**2)
290 nx(i) = nx(i)*norminv
291 ny(i) = ny(i)*norminv
292 nz(i) = nz(i)*norminv
293
294 t2x(i) = ny(i)*t1z(i) - nz(i)*t1y(i)
295 t2y(i) = nz(i)*t1x(i) - nx(i)*t1z(i)
296 t2z(i) = nx(i)*t1y(i) - ny(i)*t1x(i)
297
298 vn(i) = vx(i)*nx(i) + vy(i)*ny(i) + vz(i)*nz(i)
299 vt1(i) = vx(i)*t1x(i) + vy(i)*t1y(i) + vz(i)*t1z(i)
300 vt2(i) = vx(i)*t2x(i) + vy(i)*t2y(i) + vz(i)*t2z(i)
301 20 CONTINUE
302
303 DO 25 i=1,jlt
304 IF(pene(i)==zero.AND.cand_f(1,index(i))==zero)THEN
305
306
307
308 vn(i) = zero
309 vt1(i) = zero
310 vt2(i) = zero
311 ENDIF
312 25 CONTINUE
313
314 DO 40 i=1,jlt
315 ii = index(i)
316 econtt = econtt + cand_f(1,ii) * vn(i) * dt05
317 econtt = econtt + cand_f(2,ii) * vt1(i) * dt05
318 econtt = econtt + cand_f(3,ii) * vt2(i) * dt05
319 fni(i) = cand_f(1,ii) + vn(i) * dt1 * stif(i)
320 ft1(i) = cand_f(2,ii) + vt1(i) * dt1 * stif(i)
321 ft2(i) = cand_f(3,ii) + vt2(i) * dt1 * stif(i)
322 40 CONTINUE
323
324 DO 100 i=1,jlt
325 IF(itied==0)THEN
326 IF(cand_f(1,index(i))*fni(i)<zero)THEN
327
328
329
330 IF(pene(i)==zero)THEN
331
332 cand_f(1,index(i)) =zero
333 ELSE
334 cand_f(1,index(i)) = sign(em30,cand_f(1,index(i)))
335 ENDIF
336 fni(i) = zero
337 ft1(i) = zero
338 ft2(i) = zero
339 vn(i) = zero
340 vt1(i) = zero
341 vt2(i) = zero
342 ELSE
343 IF (inconv==1) cand_f(1,index(i)) = fni(i)
344 ENDIF
345 ELSEIF(fni(i)==zero.AND.pene(i)/=zero)THEN
346 cand_f(1,index(i)) = em30
347 ELSE
348 IF (inconv==1) cand_f(1,index(i)) = fni(i)
349 ENDIF
350 IF (inconv==1) THEN
351 cand_f(2,index(i)) = ft1(i)
352 cand_f(3,index(i)) = ft2(i)
353 ENDIF
354
355 100 CONTINUE
356
357 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
358 DO 120 i=1,jlt
359 vis = visc * sqrt(two * stif(i) * msi(i))
360 fni(i) = fni(i) + vn(i) * vis
361 ft1(i) = ft1(i) + vt1(i) * vis
362 ft2(i) = ft2(i) + vt2(i) * vis
363 stif(i) = stif(i) + two * vis * dt1inv
364
365 econtdt = econtdt
366 . + vis * (vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)) * dt1
367 120 CONTINUE
368 ELSE
369 DO i=1,jlt
370 vis = visc * sqrt(two * stif(i) * msi(i))
371 fni(i) = fni(i) + vn(i) * vis
372 ft1(i) = ft1(i) + vt1(i) * vis
373 ft2(i) = ft2(i) + vt2(i) * vis
374
375 c(i) = vis
376
377 econtdt = econtdt
378 . + vis * (vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)) * dt1
379 ENDDO
380 ENDIF
381
382
383
384 fsav1 = zero
385 fsav2 = zero
386 fsav3 = zero
387 fsav4 = zero
388 fsav5 = zero
389 fsav6 = zero
390 fsav8 = zero
391 fsav9 = zero
392 fsav10= zero
393 fsav11= zero
394 fsav12= zero
395 fsav13= zero
396 fsav14= zero
397 fsav15= zero
398 fsav22= zero
399 fsav23= zero
400 fsav24= zero
401 DO i=1,jlt
402 ii = index(i)
403 econtt = econtt + cand_f(1,ii) * vn(i) * dt05
404 econtt = econtt + cand_f(2,ii) * vt1(i) * dt05
405 econtt = econtt + cand_f(3,ii) * vt2(i) * dt05
406 fxn(i)= nx(i)*fni(i)
407 fyn(i)= ny(i)*fni(i)
408 fzn(i)= nz(i)*fni(i)
409 fxt(i)= t1x(i)*ft1(i) + t2x(i)*ft2(i)
410 fyt(i)= t1y(i)*ft1(i) + t2y(i)*ft2(i)
411 fzt(i)= t1z(i)*ft1(i) + t2z(i)*ft2(i)
412 impx=fxn(i)*dt12
413 impy=fyn(i)*dt12
414 impz=fzn(i)*dt12
415 fsav1=fsav1+impx
416 fsav2=fsav2+impy
417 fsav3=fsav3+impz
418 fsav8 =fsav8 +abs(impx)
419 fsav9 =fsav9 +abs(impy)
420 fsav10=fsav10+abs(impz)
421 fsav11=fsav11+fni(i)*dt12
422 impx=fxt(i)*dt12
423 impy=fyt(i)*dt12
424 impz=fzt(i)*dt12
425 fsav4=fsav4+impx
426 fsav5=fsav5+impy
427 fsav6=fsav6+impz
428 fsav12=fsav12+abs(impx)
429 fsav13=fsav13+abs(impy)
430 fsav14=fsav14+abs(impz)
431 fsav15=fsav15+sqrt(impx*impx+impy*impy+impz*impz)
432 fxi(i) = fxn(i) + fxt(i)
433 fyi(i) = fyn(i) + fyt(i)
434 fzi(i) = fzn(i) + fzt(i)
435 impx=fxi(i)*dt12
436 impy=fyi(i)*dt12
437 impz=fzi(i)*dt12
438 xp(i)=xi(i)+pene(i)*n1(i)
439 yp(i)=yi(i)+pene(i)*n2(i)
440 zp(i)=zi(i)+pene(i)*n3(i)
441 fsav22=fsav22+yp(i)*impz-zp(i)*impy
442 fsav23=fsav23+zp(i)*impx-xp(i)*impz
443 fsav24=fsav24+xp(i)*impy-yp(i)*impx
444 ENDDO
445 IF(isensint(1)/=0) THEN
446 DO i=1,jlt
447 fsavparit(1,1,i+nft) = fxn(i)
448 fsavparit(1,2,i+nft) = fyn(i)
449 fsavparit(1,3,i+nft) = fzn(i)
450 fsavparit(1,4,i+nft) = fxt(i)
451 fsavparit(1,5,i+nft) = fyt(i)
452 fsavparit(1,6,i+nft) = fzt(i)
453 ENDDO
454 ENDIF
455
456 IF (inconv==1) THEN
457#include "lockon.inc"
458 fsav(1)=fsav(1)+fsav1
459 fsav(2)=fsav(2)+fsav2
460 fsav(3)=fsav(3)+fsav3
461 fsav(4)=fsav(4)+fsav4
462 fsav(5)=fsav(5)+fsav5
463 fsav(6)=fsav(6)+fsav6
464 fsav(8) = fsav(8) +fsav8
465
466 fsav(10) = fsav(10) +fsav10
467 fsav(11) = fsav(11) +fsav11
468 fsav(12) = fsav(12) + fsav12
469 fsav(13) = fsav(13) + fsav13
470 fsav(14) = fsav(14) + fsav14
471 fsav(15) = fsav(15) + fsav15
472 fsav(22) = fsav(22) + fsav22
473 fsav(23) = fsav(23) + fsav23
474 fsav(24) = fsav(24) + fsav24
475 econt_cumu = econt_cumu + econtt
476 econtd = econtd + econtdt
477 fsav(26) = fsav(26) + econtt
478 fsav(28) = fsav(28) + econtdt
479#include "lockoff.inc"
480 ENDIF
481
482
483
484 IF(nisub/=0)THEN
485 DO i=1,jlt
486 nn = nsvg(i)
487 IF(nn>0)THEN
488 in=cn_loc(i)
489 ie=ce_loc(i)
490 jj =addsubs(in)
491 kk =addsubm(ie)
492 DO WHILE(jj<addsubs(in+1))
493 jsub=lisubs(jj)
494 DO WHILE(kk<addsubm(ie+1))
495 ksub=lisubm(kk)
496 IF(ksub==jsub)THEN
497 impx=fxn(i)*dt12
498 impy=fyn(i)*dt12
499 impz=fzn(i)*dt12
500
501 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
502 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
503 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
504
505 IF(isensint(jsub+1)/=0) THEN
506 fsavparit(jsub+1,1,i+nft) = fxn(i)
507 fsavparit(jsub+1,2,i+nft) = fyn(i)
508 fsavparit(jsub+1,3,i+nft) = fzn(i)
509 ENDIF
510
511 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
512 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
513 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
514 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
515
516 impx=fxt(i)*dt12
517 impy=fyt(i)*dt12
518 impz=fzt(i)*dt12
519
520 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
521 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
522 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
523
524 IF(isensint(jsub+1)/=0) THEN
525 fsavparit(jsub+1,4,i+nft) = fxt(i)
526 fsavparit(jsub+1,5,i+nft) = fyt(i)
527 fsavparit(jsub+1,6,i+nft) = fzt(i)
528 ENDIF
529
530 impx=fxi(i)*dt12
531 impy=fyi(i)*dt12
532 impz=fzi(i)*dt12
533 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
534 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
535 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
536 fsavsub1(15,jsub)= fsavsub1(15,jsub)
537 . +sqrt(impx*impx+impy*impy+impz*impz)
538 fsavsub1(22,jsub)=fsavsub1(22,jsub)
539 . +yp(i)*impz-zp(i)*impy
540 fsavsub1(23,jsub)=fsavsub1(23,jsub)
541 . +zp(i)*impx-xp(i)*impz
542 fsavsub1(24,jsub)=fsavsub1(24,jsub)
543 . +xp(i)*impy-yp(i)*impx
544 kk=kk+1
545 GO TO 200
546 ELSE IF(ksub<jsub)THEN
547 kk=kk+1
548 ELSE
549 GO TO 200
550 END IF
551 END DO
552 200 CONTINUE
553 jj=jj+1
554 END DO
555 END IF
556 END DO
557 IF(nspmd>1) THEN
558
559 DO i=1,jlt
560 nn = nsvg(i)
561 IF(nn<0)THEN
562 nn = -nn
563 ie=ce_loc(i)
565 kk =addsubm(ie)
568 DO WHILE(kk<addsubm(ie+1))
569 ksub=lisubm(kk)
570 IF(ksub==jsub)THEN
571 impx=fxn(i)*dt12
572 impy=fyn(i)*dt12
573 impz=fzn(i)*dt12
574
575 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
576 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
577 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
578
579 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
580 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
581 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
582 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
583
584 impx=fxt(i)*dt12
585 impy=fyt(i)*dt12
586 impz=fzt(i)*dt12
587
588 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
589 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
590 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
591
592 impx=fxi(i)*dt12
593 impy=fyi(i)*dt12
594 impz=fzi(i)*dt12
595 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
596 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
597 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
598 fsavsub1(15,jsub)= fsavsub1(15,jsub)
599 . +sqrt(impx*impx+impy*impy+impz*impz)
600 fsavsub1(22,jsub)=fsavsub1(22,jsub)
601 . +yp(i)*impz-zp(i)*impy
602 fsavsub1(23,jsub)=fsavsub1(23,jsub)
603 . +zp(i)*impx-xp(i)*impz
604 fsavsub1(24,jsub)=fsavsub1(24,jsub)
605 . +xp(i)*impy-yp(i)*impx
606 kk=kk+1
607 GO TO 250
608 ELSE IF(ksub<jsub)THEN
609 kk=kk+1
610 ELSE
611 GO TO 250
612 END IF
613 END DO
614 250 CONTINUE
615 jj=jj+1
616 END DO
617 END IF
618 END DO
619 END IF
620#include "lockon.inc"
621 DO jsub=1,nisub
622 nsub=lisub(jsub)
623 DO j=1,24
624 fsavsub(j,nsub)=fsavsub(j,nsub)+fsavsub1(j,jsub)
625 END DO
626 END DO
627#include "lockoff.inc"
628 END IF
629
630 DO 160 i=1,jlt
631 fx1(i)=fxi(i)*h1(i)
632 fy1(i)=fyi(i)*h1(i)
633 fz1(i)=fzi(i)*h1(i)
634
635 fx2(i)=fxi(i)*h2(i)
636 fy2(i)=fyi(i)*h2(i)
637 fz2(i)=fzi(i)*h2(i)
638
639 fx3(i)=fxi(i)*h3(i)
640 fy3(i)=fyi(i)*h3(i)
641 fz3(i)=fzi(i)*h3(i)
642
643 fx4(i)=fxi(i)*h4(i)
644 fy4(i)=fyi(i)*h4(i)
645 fz4(i)=fzi(i)*h4(i)
646
647 160 CONTINUE
648
649
650 IF (nspmd>1) THEN
651
652#include "mic_lockon.inc"
653 DO i = 1,jlt
654 nn = nsvg(i)
655 IF(nn<0)THEN
656
658 ENDIF
659 ENDDO
660
661#include "mic_lockoff.inc"
662 ENDIF
663
664 IF(kdtint/=0)THEN
665 DO i=1,jlt
666 k1(i) =stif(i)*abs(h1(i))
667 c1(i) =c(i)*abs(h1(i))
668 k2(i) =stif(i)*abs(h2(i))
669 c2(i) =c(i)*abs(h2(i))
670 k3(i) =stif(i)*abs(h3(i))
671 c3(i) =c(i)*abs(h3(i))
672 k4(i) =stif(i)*abs(h4(i))
673 c4(i) =c(i)*abs(h4(i))
674 ENDDO
675 END IF
676
677 IF(idtmins==2.OR.idtmins_int/=0)THEN
678 dti=dt2t
679 CALL i10sms2(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
680 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
681 3 nin ,noint ,mskyi_sms, iskyi_sms,nsms ,
682 4 stif ,c ,dti )
683 IF(dti<dt2t)THEN
684 dt2t = dti
685 neltst = noint
686 ityptst = 10
687 ENDIF
688 ENDIF
689
690 IF(idtmins_int/=0)THEN
691 stif(1:jlt)=zero
692 END IF
693
694 IF(iparit==0)THEN
695 IF(kdtint==0)THEN
696 CALL i7ass0(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
697 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
698 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
699 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
700 5 fxi ,fyi ,fzi ,a ,stifn ,nin ,
701 6 ibid ,bid ,bid ,bid ,bid ,bid ,
702 7 bid ,bid ,bid ,jtask,ibid ,ibid)
703 ELSE
704 CALL i7ass05(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
705 2 nsvg ,h1 ,h2 ,h3 ,h4 ,
706 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
707 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
708 5 fxi ,fyi ,fzi ,a ,stifn ,viscn,
709 6 stif ,k1 ,k2 ,k3 ,k4 ,c ,
710 7 c1 ,c2 ,c3 ,c4 ,nin ,ibid ,
711 8 bid ,bid ,bid ,bid ,bid ,bid,
712 9 jtask ,bid ,bid ,ibid ,ibid )
713 END IF
714 ELSE
715 IF(kdtint==0)THEN
716 CALL i7ass2(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
717 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
718 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
719 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
720 5 fxi ,fyi ,fzi ,fskyi,isky ,niskyfi,
721 6 nin ,noint ,ibid ,bid ,bid ,bid ,
722 7 bid ,bid ,bid ,bid ,bid ,
723 a ibid ,ibid )
724 ELSE
725 CALL i7ass25(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
726 2 nsvg ,h1 ,h2 ,h3 ,h4 ,
727 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
728 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
729 5 fxi ,fyi ,fzi ,fskyi,niskyfi,nin ,
730 6 stif ,k1 ,k2 ,k3 ,k4 ,c ,
731 7 c1 ,c2 ,c3 ,c4 ,isky ,noint,
732 8 ibid ,bid ,bid ,bid ,bid ,bid ,
733 9 bid ,bid ,bid ,ibid ,ibid )
734 END IF
735 END IF
736
737 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT >0)THEN
738 IF (inconv==1) THEN
739#include "lockon.inc"
740 DO i=1,jlt
741 fcont(1,ix1(i)) =fcont(1,ix1(i)) + fx1(i)
742 fcont(2,ix1(i)) =fcont(2,ix1(i)) + fy1(i)
743 fcont(3,ix1(i)) =fcont(3,ix1(i)) + fz1(i)
744 fcont(1,ix2(i)) =fcont(1,ix2(i)) + fx2(i)
745 fcont(2,ix2(i)) =fcont(2,ix2(i)) + fy2(i)
746 fcont(3,ix2(i)) =fcont(3,ix2(i)) + fz2(i)
747 fcont(1,ix3(i)) =fcont(1,ix3(i)) + fx3(i)
748 fcont(2,ix3(i)) =fcont(2,ix3(i)) + fy3(i)
749 fcont(3,ix3(i)) =fcont(3,ix3(i)) + fz3(i)
750 fcont(1,ix4(i)) =fcont(1,ix4(i)) + fx4(i)
751 fcont(2,ix4(i)) =fcont(2,ix4(i)) + fy4(i)
752 fcont(3,ix4(i)) =fcont(3,ix4(i)) + fz4(i)
753 jg = nsvg(i)
754 IF(jg>0) THEN
755
756 fcont(1,jg)=fcont(1,jg)- fxi(i)
757 fcont(2,jg)=fcont(2,jg)- fyi(i)
758 fcont(3,jg)=fcont(3,jg)- fzi(i)
759 ENDIF
760 ENDDO
761#include "lockoff.inc"
762 END IF
763 ENDIF
764
765 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
766 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
767 . (manim>=4.AND.manim<=15).OR. h3d_data%MH3D/=0))
768 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
769 IF (inconv==1) THEN
770#include "lockon.inc"
771 DO i=1,jlt
772 fncont(1,ix1(i)) =fncont(1,ix1(i)) + fx1(i)
773 fncont(2,ix1(i)) =fncont(2,ix1(i)) + fy1(i)
774 fncont(3,ix1(i)) =fncont(3,ix1(i)) + fz1(i)
775 fncont(1,ix2(i)
776 fncont(2,ix2(i)) =fncont(2,ix2
777 fncont(3,ix2(i)) =fncont(3,ix2(i)) + fz2(i)
778 fncont(1,ix3(i)) =fncont(1,ix3(i)) + fx3(i)
779 fncont(2,ix3(i)) =fncont(2,ix3(i)) + fy3(i)
780 fncont(3,ix3(i)) =fncont(3,ix3(i)) + fz3(i)
781 fncont(1,ix4(i)) =fncont(1,ix4(i)) + fx4(i)
782 fncont(2,ix4(i)) =fncont(2,ix4(i)) + fy4(i)
783 fncont(3,ix4(i)) =fncont(3,ix4(i)) + fz4(i)
784 jg = nsvg(i)
785 IF(jg>0) THEN
786
787 fncont(1,jg)=fncont(1,jg)- fxi(i)
788 fncont(2,jg)=fncont(2,jg)- fyi(i)
789 fncont(3,jg)=fncont(3,jg)- fzi(i)
790 ELSE
791 jg = -jg
795 ENDIF
796 ENDDO
797#include "lockoff.inc"
798 END IF
799 ENDIF
800
801
802 IF(isecin>0.AND.inconv==1)THEN
803 k0=nstrf(25)
804 IF(nstrf(1)+nstrf(2)/=0)THEN
805 DO i=1,nsect
806 nbinter=nstrf(k0+14)
807 k1s=k0+30
808 DO j=1,nbinter
809 IF(nstrf(k1s)==noint)THEN
810 IF(isecut/=0)THEN
811#include "lockon.inc"
812 DO k=1,jlt
813
814
815 IF(secfcum(4,ix1(k),i)==1.)THEN
816 secfcum(1,ix1(k),i)=secfcum(1,ix1(k),i)-fx1(k)
817 secfcum(2,ix1(k),i)
818 secfcum(3,ix1(k),i)=secfcum(3,ix1(k),i)-fz1(k)
819 ENDIF
820 IF(secfcum(4,ix2(k),i)==1.)THEN
821 secfcum(1,ix2(k),i)=secfcum(1,ix2(k),i)-fx2(k)
822 secfcum(2,ix2(k),i
823 secfcum(3,ix2(k),i)=secfcum(3,ix2(k),i)-fz2(k)
824 ENDIF
825 IF(secfcum(4,ix3(k),i)==1.)THEN
826 secfcum(1,ix3(k),i)=secfcum(1,ix3(k),i)-fx3(k)
827 secfcum(2,ix3(k),i)=secfcum(2,ix3(k),i)-fy3(k)
828 secfcum(3,ix3(k),i)=secfcum(3,ix3
829 ENDIF
830 IF(secfcum(4,ix4(k),i)==1.)THEN
831 secfcum(1,ix4(k),i)=secfcum(1,ix4(k),i)-fx4(k)
832 secfcum(2,ix4(k),i)=secfcum(2,ix4(k),i)-fy4(k)
833 secfcum(3,ix4(k),i)=secfcum(3,ix4(k),i)-fz4(k)
834 ENDIF
835 jg = nsvg(k)
836 IF(jg>0) THEN
837
838 IF(secfcum(4,jg,i)==1.)THEN
839 secfcum(1,jg,i)=secfcum(1,jg,i)+fxi(k)
840 secfcum(2,jg,i)=secfcum(2,jg,i)+fyi(k)
841 secfcum(3,jg,i)=secfcum(3,jg,i)+fzi(k)
842 ENDIF
843 ENDIF
844 ENDDO
845#include "lockoff.inc"
846 ENDIF
847
848 ENDIF
849 k1s=k1s+1
850 ENDDO
851 k0=nstrf(k0+24)
852 ENDDO
853 ENDIF
854 ENDIF
855
856 IF (idamp_rdof/=0) THEN
857 DO i=1,jlt
858
859
860 IF(fxi(i)/=zero.OR.fyi(i)/=zero.OR.fzi(i)/=zero)THEN
861 jg = nsvg(i)
862 IF(jg>0) THEN
863
864 icontact(jg)=1
865 ENDIF
866 icontact(ix1(i))=1
867 icontact(ix2(i))=1
868 icontact(ix3(i))=1
869 icontact(ix4(i))=1
870 ENDIF
871 ENDDO
872 ENDIF
873
874 RETURN
subroutine i10sms2(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, nin, noint, mskyi_sms, iskyi_sms, nsms, kt, c, dti)
subroutine i7ass0(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, a, stifn, nin, intth, phi, fthe, phi1, phi2, phi3, phi4, condn, condint, jtask, iform, nodadt_therm)
subroutine i7ass2(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, fskyi, isky, niskyfi, nin, noint, intth, phi, ftheskyi, phi1, phi2, phi3, phi4, condnskyi, condint, iform, nodadt_therm)
subroutine i7ass05(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, a, stifn, viscn, ks, k1, k2, k3, k4, cs, c1, c2, c3, c4, nin, intth, phi, fthe, phi1, phi2, phi3, phi4, jtask, condn, condint, iform, nodadt_therm)
subroutine i7ass25(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, fskyi, niskyfi, nin, ks, k1, k2, k3, k4, cs, c1, c2, c3, c4, isky, noint, intth, phi, ftheskyi, phi1, phi2, phi3, phi4, condnskyi, condint, iform, nodadt_therm)
type(real_pointer2), dimension(:), allocatable fnconti
type(int_pointer), dimension(:), allocatable lisubsfi
type(int_pointer), dimension(:), allocatable nsvfi
type(int_pointer), dimension(:), allocatable addsubsfi