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