61
62
63
67
68
69
70#include "implicit_f.inc"
71#include "comlock.inc"
72#include "assert.inc"
73
74
75
76#include "mvsiz_p.inc"
77
78
79
80#include "com01_c.inc"
81#include "com04_c.inc"
82#include "com06_c.inc"
83#include "com08_c.inc"
84#include "scr05_c.inc"
85#include "scr11_c.inc"
86#include "scr14_c.inc"
87#include "scr16_c.inc"
88#include "scr18_c.inc"
89#include "param_c.inc"
90#include "parit_c.inc"
91#include "impl1_c.inc"
92#include "sms_c.inc"
93
94
95
96 INTEGER NELTST,ITYPTST,JLT,IBC,IVIS2,INACTI,NEDGE,NISKYFIE,NIN,ILEV,
97 . IFQ
98 INTEGER ICODT(*), ITAB(*), ISKY(*),
99 . NOINT,NEWFRONT,ISECIN, NSTRF(*), ISKYI_SMS(*),
100 . NISUB, LISUB(*), ADDSUBE(*), LISUBE(*), INFLG_SUBE(*),
101 . TYPSUB(*)
102 INTEGER N1(MVSIZ), N2(MVSIZ), M1(MVSIZ), M2(MVSIZ),
103 . CS_LOC(MVSIZ), CM_LOC(MVSIZ), NSMS(MVSIZ),JTASK,
104 . ISENSINT(*), NFT, INDEX(*), EBINFLG(*), IBM(*),IFPEN(*),
105 . TAGNCONT(NLOADP_HYD_INTER,NUMNOD)
106 INTEGER :: EDGE_ID(2,MVSIZ)
107 INTEGER , INTENT(IN) :: NINLOADP,S_LOADPINTER
108 INTEGER , INTENT(IN) :: KLOADPINTER(NINTER+1),LOADPINTER(S_LOADPINTER),
109 . LOADP_HYD_INTER(NLOADP_HYD)
111 . stiglo,
112 . a(3,*), ms(*), v(3,*), fsav(*),fcont(3,*),
113 . stifn(*),fskyi(lskyi,nfskyi),fsavsub(nthvki,*),
114 . mskyi_sms(*), gapve(*), cand_p(*),
115 . gap,fric,visc,viscf,vis,dt2t, startt
117 . hs1(mvsiz), hs2(mvsiz), hm1(mvsiz), hm2(mvsiz),
118 . nx(mvsiz), ny(mvsiz), nz(mvsiz), stif(mvsiz),
119 . secfcum(7,numnod,nsect), viscn(*),
120 . ms1(mvsiz),ms2(mvsiz),mm1(mvsiz),mm2(mvsiz),
121 . vxs1(mvsiz),vys1(mvsiz),vzs1(mvsiz),vxs2(mvsiz),vys2(mvsiz),
122 . vzs2(mvsiz),vxm1(mvsiz),vym1(mvsiz),vzm1(mvsiz),vxm2(mvsiz),
123 . vym2(mvsiz),vzm2(mvsiz),fsavparit(nisub+1,11,*),fricc(mvsiz),
124 . cand_fx(*),cand_fy(*),cand_fz(*)
125 my_real ,
INTENT(IN) :: dgaploadint(s_loadpinter)
126 TYPE(H3D_DATABASE) ::
127
128
129
130 INTEGER I, J, K0, NBINTER, K1S, K, NI, IL, IE, PP, PPL
131
132 INTEGER JSUB,KSUB,NSUB,JJ,KK,ISS1,ISS2,IMS1,,ITYPSUB
133 INTEGER TAGIP(MVSIZ)
135 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
136 . fxi(mvsiz), fyi(mvsiz), fzi(mvsiz), fni(mvsiz),
137 . fx1(mvsiz), fx2(mvsiz), fx3(mvsiz), fx4(mvsiz),
138 . fy1(mvsiz), fy2(mvsiz), fy3(mvsiz), fy4(mvsiz),
139 . fz1(mvsiz), fz2(mvsiz), fz3(mvsiz), fz4(mvsiz),
140 . fxt(mvsiz), fyt(mvsiz), fzt(mvsiz),
141 . pene(mvsiz),dist(mvsiz),
142 . vis2(mvsiz),
143 . vmax,s2,
144 . dt1inv, fac, ff,
145 . fx, fy, fz, mas2,dti,
146 . econtt, econvt,masm,
147 . fsav1, fsav2, fsav3, fsav4, fsav5, fsav6, fsav8,
148 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15,
149 . dgapload,
150 . fsavsub1(24,nisub), impx, impy, impz,ftn ,fn , ft,beta
152 . prec
154 . st1(mvsiz),st2(mvsiz),st3(mvsiz),st4(mvsiz),stif0(mvsiz),
155 . kt(mvsiz),c(mvsiz),cf(mvsiz),
156 . k1(mvsiz),k2(mvsiz),k3(mvsiz),k4(mvsiz),
157 . c1(mvsiz),c2(mvsiz),c3(mvsiz),c4(mvsiz),
158 . cx,cy,cfi,aux,gapp
159
160 INTEGER BITGET
162
163 IF (iresp == 1) THEN
164 prec = fiveem4
165 ELSE
166 prec = em10
167 ENDIF
168 IF(dt1>zero)THEN
169 dt1inv = one/dt1
170 ELSE
171 dt1inv =zero
172 ENDIF
173 econtt = zero
174 econvt = zero
175 DO i=1,jlt
176 stif0(i) = stif(i)
177 ENDDO
178
179 DO i=1,jlt
180 s2 = sqrt(nx(i)**2 + ny(i)**2 + nz(i)**2)
181 dist(i)=s2
182 IF(gapve(i)/=zero)THEN
183 pene(i) =
max(zero,gapve(i) - s2)
184 ELSE
185 pene(i) = s2
186 END IF
187 s2 = one/
max(em30,s2)
188 nx(i) = nx(i)*s2
189 ny(i) = ny(i)*s2
190 nz(i) = nz(i)*s2
191 ENDDO
192
193 DO i=1,jlt
194
195 debug_e2e(edge_id(1,i)==d_em.AND.edge_id(2,i)==d_es,cand_p(index(i)))
196 IF(cand_p(index(i))==zero)cand_p(index(i))=pene(i)
197 ENDDO
198
199 IF(inacti/=-1)THEN
200 DO i=1,jlt
201
202 IF(cand_p(index(i))<zero) THEN
203 IF(startt>zero) THEN
204 cand_p(index(i))=pene(i) ! 1st impact
205 ELSE
206 cand_p(index(i))=-cand_p(index(i))
207 ENDIF
208 ENDIF
209
210 IF(pene(i)/=cand_p(index(i)))
211 . cand_p(index(i))=
min(cand_p(index(i)),
212 . ((one-fiveem2)*cand_p(index(i))+fiveem2*pene(i)) )
213 pene(i)=
max(zero,pene(i)-cand_p(index(i)))
214 IF( pene(i)==zero ) stif(i) = zero
215 ENDDO
216 ELSE
217 DO i=1,jlt
218 IF(cand_p(index(i)) < zero)THEN
219
220
221 cand_p(index(i)) = -cand_p(index(i))
222 IF(pene(i)/=cand_p(index(i)))
223 . cand_p(index(i)) =
min(cand_p(index(i)),
224 . ((one-fiveem2)*cand_p(index(i))+fiveem2*pene(i)) )
225 cand_p(index(i)) = -cand_p(index(i))
226 IF( pene(i)==zero ) stif(i) = zero
227
228 ELSE
229
230
231 IF(pene(i)/=cand_p(index(i)))
232 . cand_p(index(i))=
min(cand_p(index(i)),
233 . ((one-fiveem2)*cand_p(index(i))+fiveem2*pene(i)) )
234
235 pene(i)=
max(zero,pene(i)-cand_p(index(i)))
236 IF( pene(i)==zero ) stif(i) = zero
237 END IF
238 ENDDO
239 ENDIF
240
241 vmax = zero
242 DO i=1,jlt
243 vx(i) = hs1(i)*vxs1(i) + hs2(i)*vxs2(i)
244 . - hm1(i)*vxm1(i) - hm2(i)*vxm2(i)
245 vy(i) = hs1(i)*vys1(i) + hs2(i)*vys2(i)
246 . - hm1(i)*vym1(i) - hm2(i)*vym2(i)
247 vz(i) = hs1(i)*vzs1(i) + hs2(i)*vzs2(i)
248 . - hm1(i)*vzm1(i) - hm2(i)*vzm2(i)
249 vn(i) = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
250 ENDDO
251
252 DO i=1,jlt
253 stif(i)= half*stif(i)
254 fni(i) = -stif(i) * pene(i)
255 econvt = econvt+fni(i)*vn(i)*dt1
256 ENDDO
257
258
259
260 IF(visc/=zero)THEN
261 IF(ivis2==-1)THEN
262 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
263 DO i=1,jlt
264 fac = stif(i) /
max(em30,stif(i))
265 mas2 = ms1(i)*hs1(i)
266 . + ms2(i)*hs2(i)
267 masm = mm1(i)*hm1(i)
268 . + mm2(i)*hm2(i)
269 vis2(i) = two * stif(i) *
min(mas2,masm)
270 vis = sqrt(vis2(i))
271 ff = fac * visc * vis
272 stif(i) = stif0(i) + ff * dt1inv
273 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
274 ff = ff * vn(i)
275 econvt = econvt + ff * vn(i) * dt1
276 fni(i) = fni(i) + ff
277 ENDDO
278
279 ELSE
280 DO i=1,jlt
281 fac = stif(i) /
max(em30,stif(i))
282 mas2 = ms1(i)*hs1(i)
283 . + ms2(i)*hs2(i)
284 masm = mm1(i)*hm1(i)
285 . + mm2(i)*hm2(i)
286 vis2(i) = two * stif(i) *
min(mas2,masm)
287 vis = sqrt(vis2(i))
288 c(i)= fac * visc * vis
289 kt(i)= stif0(i)
290 stif(i) = stif(i) + c(i) * dt1inv
291 ff = c(i) * vn(i)
292 econvt = econvt + ff * vn(i) * dt1
293 fni(i) = fni(i) + ff
294 cf(i) = fac*sqrt(viscf)*vis
295 stif(i) =
max(stif(i) ,cf(i)*dt1inv)
296 ENDDO
297 ENDIF
298 ELSEIF(ivis2==1)THEN
299
300 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
301 DO i=1,jlt
302 fac = stif(i) /
max(em30,stif(i))
303 mas2 = ms1(i)*hs1(i)
304 . + ms2(i)*hs2(i)
305 masm = mm1(i)*hm1(i)
306 . + mm2(i)*hm2(i)
307 vis2(i) = two* stif(i) * masm * mas2 /
308 .
max(em30,masm+mas2)
309 vis = sqrt(vis2(i))
310 ff = fac * visc * vis
311 stif(i) = stif0(i) + ff * dt1inv
312 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
313 ff = ff * vn(i)
314 econvt = econvt + ff * vn(i) * dt1
315 fni(i) = fni(i) + ff
316 ENDDO
317
318 ELSE
319 DO i=1,jlt
320 fac = stif(i) /
max(em30,stif(i))
321 mas2 = ms1(i)*hs1(i)
322 . + ms2(i)*hs2(i)
323 masm = mm1(i)*hm1(i)
324 . + mm2(i)*hm2(i)
325 vis2(i) = two* stif(i) * masm * mas2 /
326 .
max(em30,masm+mas2)
327 vis = sqrt(vis2(i))
328 c(i)= fac * visc * vis
329 kt(i)= stif0(i)
330 stif(i) = stif(i) + c(i) * dt1inv
331 ff = c(i) * vn(i)
332 econvt = econvt + ff * vn(i) * dt1
333 fni(i) = fni(i) + ff
334 cf(i) = fac*sqrt(viscf)*vis
335 stif(i) =
max(stif(i) ,cf(i
336 ENDDO
337 ENDIF
338
339 ELSEIF(ivis2==2)THEN
340
341
342
343 DO i=1,jlt
344 fac = stif(i) /
max(em30,stif(i))
345 mas2 = ms1(i)*hs1(i)
346 . + ms2(i)*hs2(i)
347 masm = mm1(i)*hm1(i)
348 . + mm2(i)*hm2(i)
349 vis2(i) = two * stif(i) *
min(mas2,masm)
350 vis = sqrt(vis2(i))
351 ff = fac * visc * vis
352 stif(i) = stif0(i) + two * ff * dt1inv
353 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv)
354 ff = ff * vn(i)
355 econvt = econvt + ff * vn(i) * dt1
356 fni(i) = fni(i) + ff
357 ENDDO
358 ELSEIF(ivis2==3)THEN
359
360
361
362 DO i=1,jlt
363 fac = stif(i) /
max(em30,stif(i))
364 mas2 = ms1(i)*hs1(i)
365 . + ms2(i)*hs2(i)
366 masm = mm1(i)*hm1(i)
367 . + mm2(i)*hm2(i)
368 vis2(i) = two * stif(i) *
min(mas2,masm)
369 vis = sqrt(vis2(i))
370 ff = fac * visc * vis
371 stif(i) = stif0(i) + two* ff * dt1inv
372 stif(i) =
max(stif(i) ,fac*sqrt(viscf)*vis*dt1inv
373 ff = ff * vn(i)
374 econvt = econvt + ff * vn(i) * dt1
375 fni(i) = fni(i) + ff
376 ENDDO
377 ELSEIF(ivis2==4)THEN
378
379
380
381 DO i=1,jlt
382 fac = stif(i) /
max(em30,stif(i))
383 mas2 = ms1(i)*hs1(i)
384 . + ms2(i)*hs2(i)
385 masm = mm1(i)*hm1(i)
386 . + mm2(i)*hm2(i)
387 vis2(i) = two * stif(i) *
min(mas2,masm)
388 vis = sqrt(vis2(i))
390 ENDDO
391 ELSEIF(ivis2==5)THEN
392
393
394
395
396
397 DO i=1,jlt
398 fac = stif(i) /
max(em30,stif(i))
399 mas2 = ms1(i)*hs1(i)
400 . + ms2(i)*hs2
401 masm = mm1(i)*hm1(i)
402 . + mm2(i)*hm2(i)
403 vis2(i) = two* stif(i) * masm * mas2 /
404 .
max(em30,masm+mas2)
405 vis = 2. * visc * dt1inv * masm * mas2 /
406 .
max(em30,masm+mas2)
407 stif(i) =
max(stif0(i) ,fac*sqrt(viscf*vis2(i))*dt1inv)
408 ff = vis * vn(i)
409 econvt = econvt + ff * vn(i) * dt1
410 fni(i) =
min(fni(i),ff)
411 ENDDO
412 ELSE
413 ENDIF
414 ELSE
415 ENDIF
416
417
418
419 fsav1 = zero
420 fsav2 = zero
421 fsav3 = zero
422 fsav8 = zero
423 fsav9 = zero
424 fsav10= zero
425 fsav11= zero
426 IF(ilev==2)THEN
427 DO i=1,jlt
428 IF(pene(i)==zero) cycle
429 ie=cm_loc(i)
430 ims2 =
bitget(ebinflg(ie),1)
431 fxi(i)=nx(i)*fni(i)
432 fyi(i)=ny(i)*fni(i)
433 fzi(i)=nz(i)*fni(i)
434 impx=fxi(i)*dt12
435 impy=fyi(i)*dt12
436 impz=fzi(i)*dt12
437 IF (ims2 > 0 ) THEN
438 fsav1 =fsav1 -impx
439 fsav2 =fsav2 -impy
440 fsav3 =fsav3 -impz
441 fsav11=fsav11-fni(i)*dt12
442 ELSE
443 fsav1 =fsav1 +impx
444 fsav2 =fsav2 +impy
445 fsav3 =fsav3 +impz
446 fsav11=fsav11+fni(i)*dt12
447 END IF
448 fsav8 =fsav8 +abs(impx)
449 fsav9 =fsav9 +abs(impy)
450 fsav10=fsav10+abs(impz)
451 IF(isensint(1)/=0) THEN
452 IF (ims2 >0 ) THEN
453 fsavparit(1,1,i) = -fxi(i)
454 fsavparit(1,2,i) = -fyi(i)
455 fsavparit(1,3,i) = -fzi(i)
456 ELSE
457 fsavparit(1,1,i) = fxi(i)
458 fsavparit(1,2,i) = fyi(i)
459
460 END IF
461 ENDIF
462 ENDDO
463 ELSE
464 DO i=1,jlt
465 IF
466 fxi(i)=nx(i)*fni(i)
467 fyi(i)=ny(i)*fni(i)
468 fzi(i)=nz(i)*fni(i)
469 impx=fxi(i)*dt12
470 impy=fyi(i)*dt12
471 impz=fzi(i)*dt12
472 fsav1 =fsav1 +impx
473 fsav2 =fsav2 +impy
474 fsav3 =fsav3 +impz
475 fsav8 =fsav8 +abs(impx)
476 fsav9 =fsav9 +abs(impy)
477 fsav10=fsav10+abs(impz)
478 fsav11=fsav11+fni(i)*dt12
479 IF(isensint(1)/=0) THEN
480 fsavparit(1,1,i) = fxi(i)
481 fsavparit(1,2,i) = fyi(i)
482 fsavparit(1,3,i) = fzi(i)
483 ENDIF
484 ENDDO
485 END IF
486 IF (imconv==1) THEN
487#include "lockon.inc"
488 fsav(1)=fsav(1)+fsav1
489 fsav(2)=fsav(2)+fsav2
490 fsav(3)=fsav(3)+fsav3
491 fsav(8)=fsav(8)+fsav8
492 fsav(9)=fsav(9)+fsav9
493 fsav(10)=fsav(10)+fsav10
494 fsav(11)=fsav(11)+fsav11
495#include "lockoff.inc"
496 ENDIF
497
498
499
500 IF(nisub/=0)THEN
501 DO jsub=1,nisub
502 DO j=1,24
503 fsavsub1(j,jsub)=zero
504 END DO
505 ENDDO
506 DO i=1,jlt
507
508 IF(pene(i) == zero)cycle
509
510 il = cs_loc(i)
511 IF(il<=nedge)THEN
512
513 ie = cm_loc(i)
514
515 jj =addsube(il)
516 kk =addsube(ie)
517 DO WHILE(jj<addsube(il+1))
518 jsub=lisube(jj)
519 itypsub = typsub(jsub)
520 IF(itypsub == 1 ) THEN
521 iss1 =
bitget(inflg_sube(jj),0)
522 iss2 =
bitget(inflg_sube(jj),1)
523 IF(kk<addsube(ie+1)) THEN
524 ksub=lisube(kk)
525 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
526 ims1 =
bitget(inflg_sube(kk),0)
527 ims2 =
bitget(inflg_sube(kk),1)
528 IF(ksub==jsub)THEN
529 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
530 . (ims2 == 1 .AND. iss1 == 1))) THEN
531 kk=kk+1
532 IF(kk<addsube(ie+1)) ksub=lisube(kk)
533 cycle
534 END IF
535 impx=fxi(i)*dt12
536 impy=fyi(i)*dt12
537 impz=fzi(i)*dt12
538
539 IF(ims2 > 0)THEN
540 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
541 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
542 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
543 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
544 ELSE
545 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
546 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
547 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
548 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
549 END IF
550
551 IF(isensint(jsub+1)/=0) THEN
552 IF(ims2 > 0)THEN
553 fsavparit(jsub+1,1,i) = -fxi(i)
554 fsavparit(jsub+1,2,i) = -fyi(i)
555 fsavparit(jsub+1,3,i) = -fzi(i)
556 ELSE
557 fsavparit(jsub+1,1,i) = fxi(i)
558 fsavparit(jsub+1,2,i) = fyi(i)
559 fsavparit(jsub+1,3,i) = fzi
560 END IF
561 ENDIF
562
563 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs
564 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
565 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
566
567 ENDIF
568 kk=kk+1
569 ksub=lisube(kk)
570 ENDDO
571 ENDIF
572 jj=jj+1
573
574 ELSEIF(itypsub == 2 ) THEN
575
576 impx=fxi(i)*dt12
577 impy=fyi(i)*dt12
578 impz=fzi(i)*dt12
579
580
581 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
582 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
583 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
584
585 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
586 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
587 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
588
589 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
590
591 IF(isensint(jsub+1)/=0) THEN
592 fsavparit(jsub+1,1,i) = fxi(i)
593 fsavparit(jsub+1,2,i) = fyi(i)
594 fsavparit(jsub+1,3,i) = fzi(i)
595 ENDIF
596
597 jj=jj+1
598
599 ELSEIF(itypsub == 3 ) THEN
600
601 iss2 =
bitget(inflg_sube(jj),0)
602 iss1 =
bitget(inflg_sube(jj),1)
603 IF(kk<addsube(ie+1)) THEN
604 ksub=lisube(kk)
605 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
606 ims2 =
bitget(inflg_sube(kk),0)
607 ims1 =
bitget(inflg_sube(kk),1)
608 IF(ksub==jsub)THEN
609 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
610 . (ims2 == 1 .AND. iss1 == 1))) THEN
611 kk=kk+1
612 IF(kk<addsube(ie+1)) ksub=lisube(kk)
613 cycle
614 END IF
615
616 IF(ims2 > 0)THEN
617 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
618 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
619 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
620 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
621 ELSE
622 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
623 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
624 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
625 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
626 ENDIF
627
628 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
629 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
630 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
631
632 IF(isensint(jsub+1)/=0) THEN
633 IF(ims2 > 0)THEN
634 fsavparit(jsub+1,1,i) = -fxi(i)
635 fsavparit(jsub+1,2,i) = -fyi(i)
636 fsavparit(jsub+1,3,i) = -fzi(i)
637 ELSE
638 fsavparit(jsub+1,1,i) = fxi(i)
639 fsavparit(jsub+1,2,i) = fyi(i)
640 fsavparit(jsub+1,3,i) = fzi(i)
641 END IF
642 ENDIF
643
644
645 ENDIF
646 kk=kk+1
647 ksub=lisube(kk)
648 ENDDO
649 ENDIF
650 jj=jj+1
651 ENDIF
652 END DO
653 END IF
654
655 ie = cm_loc(i)
656
657 kk =addsube(ie)
658 DO WHILE(kk<addsube(ie+1))
659 ksub=lisube(kk)
660 itypsub = typsub(ksub)
661 IF(itypsub == 2 ) THEN
662
663 impx=-fxi(i)*dt12
664 impy=-fyi(i)*dt12
665 impz=-fzi(i)*dt12
666
667 fsavsub1(1,ksub)=fsavsub1(1,ksub)+impx
668 fsavsub1(2,ksub)=fsavsub1(2,ksub)+impy
669 fsavsub1(3,ksub)=fsavsub1(3,ksub)+impz
670
671 fsavsub1(8,ksub) =fsavsub1(8,ksub) +abs(impx)
672 fsavsub1(9,ksub) =fsavsub1(9,ksub) +abs(impy)
673 fsavsub1(10,ksub)=fsavsub1(10,ksub)+abs(impz)
674
675 fsavsub1(11,ksub)=fsavsub1(11,ksub)-fni(i)*dt12
676
677 IF(isensint(ksub+1)/=0) THEN
678 fsavparit(ksub+1,1,i) = -fxi(i)
679 fsavparit(ksub+1,2,i) = -fyi(i)
680 fsavparit(ksub+1,3,i) = -fzi(i)
681 ENDIF
682
683 ENDIF
684 kk=kk+1
685 ENDDO
686
687 END DO
688 IF(nspmd>1) THEN
689#ifdef WITH_ASSERT
690 s_addsubfie = 1
691 s_lisubsfie = 0
692 DO p = 1,nspmd
693 s_addsubfie = s_addsubfie +
nsnfie(nin)%P(p)
694 s_lisubsfie = s_lisubsfie +
nisubsfie(nin)%P(p)
695 END DO
696#endif
697 DO i=1,jlt
698 IF(pene(i) == zero)cycle
699
700
701
702
703
704 il = cs_loc(i)
705 IF(il > nedge)THEN
706 il = il - nedge
707 ie = cm_loc(i)
708
709
710 assert(il <= s_addsubfie)
712 kk =addsube(ie)
713
714
715
716
718
719 assert(jj <= s_lisubsfie)
721 itypsub = typsub(jsub)
722 IF(itypsub == 1 ) THEN
723 assert(jsub <= nisub)
726 IF(kk<addsube(ie+1)) THEN
727 ksub=lisube(kk)
728 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
729 assert(ksub <= nisub)
730 ims1 =
bitget(inflg_sube(kk),0)
731 ims2 =
bitget(inflg_sube(kk),1)
732 IF(ksub==jsub)THEN
733 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
734 . (ims2 == 1 .AND. iss1 == 1))) THEN
735 kk=kk+1
736 IF(kk<addsube(ie+1)) ksub=lisube(kk)
737 cycle
738 END IF
739 impx=fxi(i)*dt12
740 impy=fyi(i)*dt12
741 impz=fzi(i)*dt12
742
743 IF(ims2 > 0)THEN
744 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
745 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
746 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
747 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
748 ELSE
749 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
750 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
751 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
752 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
753 END IF
754
755 IF(isensint(jsub+1)/=0) THEN
756 IF(ims2 > 0)THEN
757 fsavparit(jsub+1,1,i) = -fxi(i)
758 fsavparit(jsub+1,2,i) = -fyi(i)
759 fsavparit(jsub+1,3,i) = -fzi(i)
760 ELSE
761 fsavparit(jsub+1,1,i) = fxi(i)
762 fsavparit(jsub+1,2,i) = fyi(i)
763 fsavparit(jsub+1,3,i) = fzi(i)
764 END IF
765 ENDIF
766
767 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
768 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
769 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
770
771 ENDIF
772 kk=kk+1
773 IF(kk<addsube(ie+1)) ksub=lisube(kk)
774 ENDDO
775 ENDIF
776 jj=jj+1
777
778 ELSEIF(itypsub == 2 ) THEN
779
780 impx=fxi(i)*dt12
781 impy=fyi(i)*dt12
782 impz=fzi(i)*dt12
783
784
785 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
786 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
787 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
788
789 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
790 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
791 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
792
793 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
794
795 IF(isensint(jsub+1)/=0) THEN
796 fsavparit(jsub+1,1,i) = fxi(i)
797 fsavparit(jsub+1,2,i) = fyi(i)
798 fsavparit(jsub+1,3,i) = fzi(i)
799 ENDIF
800
801 jj=jj+1
802
803 ELSEIF(itypsub == 3 ) THEN
804
807 IF(kk<addsube(ie+1)) THEN
808 ksub=lisube(kk
809 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
810 ims2 =
bitget(inflg_sube(kk),0)
811 ims1 =
bitget(inflg_sube(kk),1)
812 IF(ksubTHEN
813 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
814 . (ims2 == 1 .AND. iss1 == 1))) THEN
815 kk=kk+1
816 IF(kk<addsube(ie+1)) ksub=lisube(kk)
817 cycle
818 END IF
819
820 impx=fxi(i)*dt12
821 impy=fyi(i)*dt12
822 impz=fzi(i)*dt12
823
824 IF(ims2 > 0)THEN
825 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
826 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
827 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
828 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i
829 ELSE
830 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
831 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
832 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
833 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
834 ENDIF
835
836 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
837 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
838 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
839
840
841 ENDIF
842 kk=kk+1
843 IF(kk<addsube(ie+1)) ksub=lisube(kk)
844 ENDDO
845 ENDIF
846 jj=jj+1
847
848 ENDIF
849 END DO
850 END IF
851 END DO
852 END IF
853 END IF
854
855
856
857 fxt(1:jlt)=zero
858 fyt(1:jlt)=zero
859 fzt(1:jlt)=zero
860
861 fsav4
862 fsav5 = zero
863 fsav6 = zero
864 fsav12= zero
865 fsav13= zero
866 fsav14= zero
867 fsav15= zero
868
869
870 IF (ifq /= 0) THEN
871 DO i=1,jlt
872
873 IF(pene(i) == zero)cycle
874
875 fx = stif0(i)*vx(i)*dt12
876 fy = stif0(i)*vy(i)*dt12
877 fz = stif0(i)*vz(i)*dt12
878
879 fx = cand_fx(index(i)) + fx
880 fy = cand_fy(index(i)) + fy
881 fz = cand_fz(index(i)) + fz
882
883 ftn = fx*nx(i) + fy
884 fx = fx - ftn*nx(i)
885 fy = fy - ftn*ny(i)
886 fz = fz - ftn*nz(i)
887 ft = fx*fx + fy*fy + fz*fz
889
890 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
891 beta =
min(one,fricc(i)*sqrt(fn/ft))
892 fxt(i) = fx * beta
893 fyt(i) = fy * beta
894 fzt(i) = fz * beta
895
896 cand_fx(index(i)) = fxt(i)
897 cand_fy(index(i)) = fyt(i)
898 cand_fz(index(i)) = fzt(i)
899
900 ifpen(index(i)) = 1
901
902 fxi(i)=fxi(i) + fxt(i)
903 fyi(i)=fyi(i) + fyt(i)
904 fzi(i)=fzi(i) + fzt(i)
905
906 fsav4 = fsav4 + fxt(i)*dt12
907 fsav5 = fsav5 + fyt(i)*dt12
908 fsav6 = fsav6 + fzt(i)*dt12
909
910 fsav12 = fsav12 + abs(fxi(i)*dt12)
911 fsav13 = fsav13 + abs(fyi(i)*dt12)
912 fsav14 = fsav14 + abs(fzi(i)*dt12)
913 fsav15 = fsav15 + sqrt(fxi(i)*fxi(i)+fyi(i)*fyi(i)+fzi(i)*fzi(i))*dt12
914 econvt = econvt
915 . + dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
916
917 ENDDO
918 ENDIF
919
920 IF (inconv==1) THEN
921#include "lockon.inc"
922 fsav(4) = fsav(4) + fsav4
923 fsav(5) = fsav(5) + fsav5
924 fsav(6) = fsav(6) + fsav6
925 fsav(12) = fsav(12) + fsav12
926 fsav(13) = fsav(13) + fsav13
927 fsav(14) = fsav(14) + fsav14
928 fsav(15) = fsav(15) + fsav15
929#include "lockoff.inc"
930 ENDIF
931
932
933
934
935
936
937
938 IF(nisub/=0)THEN
939 DO i=1,jlt
940
941 IF(pene(i) == zero)cycle
942
943 il = cs_loc(i)
944 IF(il<=nedge)THEN
945
946 ie = cm_loc(i)
947
948 jj =addsube(il)
949 kk =addsube(ie)
950 DO WHILE(jj<addsube(il+1))
951 jsub=lisube(jj)
952 itypsub = typsub(jsub)
953
954 IF(itypsub == 1 ) THEN
955 iss1 =
bitget(inflg_sube(jj),0)
956 iss2 =
bitget(inflg_sube(jj),1)
957 IF(kk<addsube(ie+1)) THEN
958 ksub=lisube(kk)
959 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
960 ims1 =
bitget(inflg_sube(kk),0)
961 ims2 =
bitget(inflg_sube(kk),1)
962 IF(ksub==jsub)THEN
963 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
964 . (ims2 == 1 .AND. iss1 == 1))) THEN
965 kk=kk+1
966 IF(kk<addsube(ie+1)) ksub=lisube(kk)
967 cycle
968 END IF
969 impx=fxt(i)*dt12
970 impy=fyt(i)*dt12
971 impz=fzt(i)*dt12
972
973 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
974 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
975 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
976
977 impx=fxi(i)*dt12
978 impy=fyi(i)*dt12
979 impz=fzi(i)*dt12
980 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
981
982 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
983
984 IF(isensint(jsub+1)/=0) THEN
985 fsavparit(jsub+1,4,i) = fxt(i)
986 fsavparit(jsub+1,5,i) = fyt(i)
987 fsavparit(jsub+1,6,i) = fzt(i)
988 ENDIF
989
990 fsavsub1(15,jsub)= fsavsub1(15,jsub)
991 . +sqrt(impx*impx+impy*impy+impz*impz)
992
993
994 ENDIF
995 kk=kk+1
996 IF(kk<addsube(ie+1)) ksub=lisube(kk)
997 ENDDO
998 ENDIF
999 jj=jj+1
1000
1001
1002 ELSEIF(itypsub == 2 ) THEN
1003
1004 impx=fxt(i)*dt12
1005 impy=fyt(i)*dt12
1006 impz=fzt(i)*dt12
1007
1008 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1009 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1010 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1011
1012 impx=fxi(i)*dt12
1013 impy=fyi(i)*dt12
1014 impz=fzi(i)*dt12
1015 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1016 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1017 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1018
1019 IF(isensint(jsub+1)/=0) THEN
1020 fsavparit(jsub+1,4,i) = fxt(i)
1021 fsavparit(jsub+1,5,i) = fyt(i)
1022 fsavparit(jsub+1,6,i) = fzt(i)
1023 ENDIF
1024
1025 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1026 . +sqrt(impx*impx+impy*impy+impz*impz)
1027
1028
1029 jj=jj+1
1030
1031 ELSEIF(itypsub == 3 ) THEN
1032
1033 iss2 =
bitget(inflg_sube(jj),0)
1034 iss1 =
bitget(inflg_sube(jj),1)
1035 IF(kk<addsube(ie+1)) THEN
1036 ksub=lisube(kk)
1037 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
1038 ims2 =
bitget(inflg_sube(kk),0)
1039 ims1 =
bitget(inflg_sube(kk),1)
1040 IF(ksub==jsub)THEN
1041 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1042 . (ims2 == 1 .AND. iss1 == 1))) THEN
1043 kk=kk+1
1044 IF(kk<addsube(ie+1)) ksub=lisube(kk)
1045 cycle
1046 END IF
1047
1048 impx=fxt(i)*dt12
1049 impy=fyt(i)*dt12
1050 impz=fzt(i)*dt12
1051
1052 IF(ims2 > 0) THEN
1053 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
1054 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
1055 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
1056 ELSE
1057
1058 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1059 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1060 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1061 ENDIF
1062
1063 impx=fxi(i)*dt12
1064 impy=fyi(i)*dt12
1065 impz=fzi(i)*dt12
1066 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1067 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1068 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1069
1070 IF(isensint(jsub+1)/=0) THEN
1071 IF(ims2 > 0) THEN
1072 fsavparit(jsub+1,4,i) = -fxt(i)
1073 fsavparit(jsub+1,5,i) = -fyt(i)
1074 fsavparit(jsub+1,6,i) = -fzt(i)
1075 ELSE
1076 fsavparit(jsub+1,4,i) = fxt(i)
1077 fsavparit(jsub+1,5,i) = fyt(i)
1078 fsavparit(jsub+1,6,i) = fzt(i)
1079 ENDIF
1080 ENDIF
1081
1082 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1083 . +sqrt(impx*impx+impy*impy+impz*impz)
1084
1085 ENDIF
1086 kk=kk+1
1087 IF(kk<addsube(ie+1))ksub=lisube(kk)
1088 ENDDO
1089 ENDIF
1090 jj=jj+1
1091
1092 ENDIF
1093
1094 END DO
1095 END IF
1096
1097 ie = cm_loc(i)
1098 kk =addsube(ie)
1099 DO WHILE(kk<addsube(ie+1))
1100 ksub=lisube(kk)
1101 itypsub = typsub(ksub)
1102 IF(itypsub == 2 ) THEN
1103
1104 impx=-fxt(i)*dt12
1105 impy=-fyt(i)*dt12
1106 impz=-fzt(i)*dt12
1107
1108 fsavsub1(4,ksub)=fsavsub1(4,ksub)+impx
1109 fsavsub1(5,ksub)=fsavsub1(5,ksub)+impy
1110 fsavsub1(6,ksub)=fsavsub1(6,ksub)+impz
1111
1112 impx=fxi(i)*dt12
1113 impy=fyi(i)*dt12
1114 impz=fzi(i)*dt12
1115 fsavsub1(12,ksub)=fsavsub1(12,jsub)+abs(impx)
1116 fsavsub1(13,ksub)=fsavsub1(13,jsub)+abs(impy)
1117 fsavsub1(14,ksub)=fsavsub1(14,jsub)+abs(impz)
1118
1119 IF(isensint(ksub+1)/=0) THEN
1120 fsavparit(ksub+1,4,i) = -fxt(i)
1121 fsavparit(ksub+1,5,i) = -fyt(i)
1122 fsavparit(ksub+1,6,i) = -fzt(i)
1123 ENDIF
1124
1125 fsavsub1(15,ksub)= fsavsub1(15,ksub)
1126 . +sqrt(impx*impx+impy*impy+impz*impz)
1127
1128 ENDIF
1129 kk=kk+1
1130 ENDDO
1131 END DO
1132 IF(nspmdTHEN
1133 DO i=1,jlt
1134
1135 IF(pene(i) == zero)cycle
1136
1137 il = cs_loc(i)
1138 IF(il>nedge)THEN
1139 il = il - nedge
1140
1141 ie = cm_loc(i)
1142
1144 kk =addsube(ie)
1147 itypsub = typsub
1148
1149 IF(itypsub == 1 ) THEN
1152 IF(kk<addsube(ie+1)) THEN
1153 ksub=lisube(kk)
1154 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie
1155 ims1 =
bitget(inflg_sube(kk),0)
1156 ims2 =
bitget(inflg_sube(kk),1)
1157 IF(ksub==jsub)THEN
1158 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1159 . (ims2 == 1 .AND. iss1 == 1))) THEN
1160 kk=kk+1
1161 IF(kk<addsube(ie+1)) ksub=lisube(kk)
1162 cycle
1163 END IF
1164 impx=fxt(i)*dt12
1165 impy=fyt(i)*dt12
1166 impz=fzt(i)*dt12
1167
1168 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1169 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1170 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1171
1172 impx=fxi(i)*dt12
1173 impy=fyi(i)*dt12
1174 impz=fzi(i)*dt12
1175 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1176 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1177 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1178
1179 IF(isensint(jsub+1)/=0) THEN
1180 fsavparit(jsub+1,4,i) = fxt(i)
1181 fsavparit(jsub+1,5,i) = fyt(i)
1182 fsavparit(jsub+1,6,i) = fzt(i)
1183 ENDIF
1184
1185 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1186 . +sqrt(impx*impx+impy*impy+impz*impz)
1187
1188
1189 ENDIF
1190 kk=kk+1
1191 IF(kk<addsube(ie+1)) ksub=lisube(kk)
1192 ENDDO
1193 ENDIF
1194 jj=jj+1
1195
1196
1197
1198 ELSEIF(itypsub == 2 ) THEN
1199
1200 impx=fxt(i)*dt12
1201 impy=fyt(i)*dt12
1202 impz=fzt(i)*dt12
1203
1204 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1205 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1206 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1207
1208 impx=fxi(i)*dt12
1209 impy=fyi(i)*dt12
1210 impz=fzi(i)*dt12
1211 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1212 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1213 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1214
1215 IF(isensint(jsub+1)/=0) THEN
1216 fsavparit(jsub+1,4,i) = fxt(i)
1217 fsavparit
1218 fsavparit(jsub+1,6,i) = fzt(i)
1219 ENDIF
1220
1221 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1222 . +sqrt(impx*impx+impy*impy+impz*impz)
1223
1224
1225 jj=jj+1
1226
1227 ELSEIF(itypsub == 3 ) THEN
1228
1231 IF(kk<addsube(ie+1)) THEN
1232 ksub=lisube(kk)
1233 DO WHILE((ksub<=jsub).AND.(kk<addsube(ie+1)))
1234 ims2 =
bitget(inflg_sube(kk),0)
1235 ims1 =
bitget(inflg_sube(kk),1)
1236 IF(ksub==jsub)THEN
1237 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1238 . (ims2 == 1 .AND. iss1 == 1))) THEN
1239 kk=kk+1
1240 IF(kk<addsube(ie+1)) ksub=lisube(kk)
1241 cycle
1242 END IF
1243
1244 impx=fxt(i)*dt12
1245 impy=fyt(i)*dt12
1246 impz=fzt(i)*dt12
1247 IF(ims2 > 0 ) THEN
1248
1249 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
1250 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
1251 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
1252 ELSE
1253
1254 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
1255 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
1256 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
1257 ENDIF
1258
1259 impx=fxi(i)*dt12
1260 impy=fyi(i)*dt12
1261 impz=fzi(i)*dt12
1262 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
1263 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
1264 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
1265
1266 IF(isensint(jsub+1)/=0) THEN
1267 IF(ims2 > 0 ) THEN
1268 fsavparit(jsub+1,4,i) = -fxt(i)
1269 fsavparit(jsub+1,5,i) = -fyt(i)
1270 fsavparit(jsub+1,6,i) = -fzt(i)
1271 ELSE
1272 fsavparit(jsub+1,4,i) = fxt(i)
1273 fsavparit(jsub+1,5,i) = fyt(i)
1274 fsavparit(jsub+1,6,i) = fzt(i)
1275 ENDIF
1276 ENDIF
1277
1278 fsavsub1(15,jsub)= fsavsub1(15,jsub)
1279 . +sqrt(impx*impx+impy*impy+impz*impz)
1280
1281
1282 ENDIF
1283 kk=kk+1
1284 IF(kk<addsube(ie+1)) ksub=lisube(kk)
1285 ENDDO
1286 ENDIF
1287 jj=jj+1
1288
1289 ENDIF
1290
1291
1292
1293 END DO
1294 END IF
1295 END DO
1296 ENDIF
1297#include "lockon.inc"
1298 DO jsub=1,nisub
1299 nsub=lisub(jsub)
1300 DO j=1,15
1301 fsavsub(j,nsub)=fsavsub(j,nsub)+fsavsub1(j,jsub)
1302 END DO
1303 fsavsub(22,nsub)=fsavsub(22,nsub)+fsavsub1(22,jsub)
1304 fsavsub(23,nsub)=fsavsub(23,nsub)+fsavsub1(23,jsub)
1305 fsavsub(24,nsub)=fsavsub(24,nsub)+fsavsub1(24,jsub)
1306 END DO
1307#include "lockoff.inc"
1308 END IF
1309
1310 IF (imconv==1) THEN
1311#include "lockon.inc"
1312 econtv = econtv + econvt
1313 econt = econt + econtt
1314#include "lockoff.inc"
1315 ENDIF
1316
1317 DO i=1,jlt
1318
1319 IF(pene(i) == zero)cycle
1320
1321 fx1(i)=-fxi(i)*hs1(i)
1322 fy1(i)=-fyi(i)*hs1(i)
1323 fz1(i)=-fzi(i)*hs1(i)
1324
1325 fx2(i)=-fxi(i)*hs2(i)
1326 fy2(i)=-fyi(i)*hs2(i)
1327 fz2(i)=-fzi(i)*hs2(i)
1328
1329 fx3(i)=fxi(i)*hm1(i)
1330 fy3(i)=fyi(i)*hm1(i)
1331 fz3(i)=fzi(i)*hm1(i)
1332
1333 fx4(i)=fxi(i)*hm2(i)
1334 fy4(i)=fyi(i)*hm2(i)
1335 fz4(i)=fzi(i)*hm2(i)
1336
1337 ENDDO
1338
1339 DO i=1,jlt
1340 stif(i) = two*stif(i)
1341 ENDDO
1342
1343
1344 IF(kdtint==1.OR.idtmins==2)THEN
1345 IF( (visc/=zero)
1346 . .AND.(ivis2==0.OR.ivis2==1))THEN
1347 DO i=1,jlt
1348 cx= c(i)*c(i)
1349
1350 IF(ms1(i)==zero)THEN
1351 k1(i) =zero
1352 c1(i) =zero
1353 ELSE
1354 k1(i)=kt(i)*abs(hs1(i))
1355 c1(i)=c(i)*abs(hs1(i))
1356 cx =four*c1(i)*c1(i)
1357 cy =eight*ms1(i)*k1(i)
1358 aux = sqrt(cx+cy)+two*c1(i)
1359 st1(i)= k1(i)*aux*aux/
max(cy,em30)
1360 cfi = cf(i)*abs(hs1(i))
1361 aux = two*cfi*cfi/
max(ms1(i),em20)
1362 IF(aux>st1(i))THEN
1363 k1(i) =zero
1364 c1(i) =cfi
1365 ENDIF
1366 ENDIF
1367
1368 IF(ms2(i)==zero)THEN
1369 k2(i) =zero
1370 c2(i) =zero
1371 ELSE
1372 k2(i)=kt(i)*abs(hs2(i))
1373 c2(i)=c(i)*abs(hs2(i))
1374 cx =four*c2(i)*c2(i)
1375 cy =eight*ms2(i)*k2(i)
1376 aux = sqrt(cx+cy)+two*c2(i)
1377 st2(i)= k2(i)*aux*aux/
max(cy,em30)
1378 cfi = cf(i)*abs(hs2(i))
1379 aux = two*cfi*cfi/
max(ms2(i),em20)
1380 IF(aux>st2(i))THEN
1381 k2(i) =zero
1382 c2(i) =cfi
1383 ENDIF
1384 ENDIF
1385
1386 IF(mm1(i)==zero)THEN
1387 k3(i) =zero
1388 c3(i) =zero
1389 ELSE
1390 k3(i)=kt(i)*abs(hm1(i))
1391 c3(i)=c(i)*abs(hm1(i))
1392 cx =four*c3(i)*c3(i)
1393 cy =eight*mm1(i)*k3(i)
1394 aux = sqrt(cx+cy)+two*c3(i)
1395 st3(i)= k3(i)*aux*aux/
max(cy,em30)
1396 cfi = cf(i)*abs(hm1(i))
1397 aux = two*cfi*cfi/
max(mm1(i),em20)
1398 IF(aux>st3(i))THEN
1399 k3(i) =zero
1400 c3(i) =cfi
1401 ENDIF
1402 ENDIF
1403
1404 IF(mm2(i)==zero)THEN
1405 k4(i) =zero
1406 c4(i) =zero
1407 ELSE
1408 k4(i)=kt(i)*abs(hm2(i))
1409 c4(i)=c(i)*abs(hm2(i))
1410 cx =four*c4(i)*c4(i)
1411 cy =eight*mm2(i)*k4(i)
1412 aux = sqrt(cx+cy)+two*c4(i)
1413 st4(i)= k4(i)*aux*aux/
max(cy,em30)
1414 cfi = cf(i)*abs(hm2(i))
1415 aux = two*cfi*cfi/
max(mm2(i),em20)
1416 IF(aux>st4(i))THEN
1417 k4(i) =zero
1418 c4(i) =cfi
1419 ENDIF
1420 ENDIF
1421 ENDDO
1422 ELSE
1423 DO i=1,jlt
1424 k1(i) =stif(i)*abs(hs1(i))
1425 c1(i) =zero
1426 k2(i) =stif(i)*abs(hs2(i))
1427 c2(i) =zero
1428 k3(i) =stif(i)*abs(hm1(i))
1429 c3(i) =zero
1430 k4(i) =stif(i)*abs(hm2(i))
1431 c4(i) =zero
1432 ENDDO
1433 ENDIF
1434 ENDIF
1435
1436 tagip(1:jlt) = 0
1437 IF(ninloadp > 0) THEN
1438 DO k = kloadpinter(nin)+1, kloadpinter(nin+1)
1439 pp = loadpinter(k)
1440 ppl = loadp_hyd_inter(pp)
1441 dgapload = dgaploadint(k)
1442 DO i=1,jlt
1443 gapp= gapve(i) + dgapload
1444 IF(pene(i) > zero .OR.dist(i) <= gapp) THEN
1445 tagip(i) = 1
1446 tagncont(ppl,m1(i)) = 1
1447 tagncont(ppl,m2(i)) = 1
1448 IF(cs_loc(i)<=nedge) THEN
1449
1450 tagncont(ppl,n1(i)) = 1
1451 tagncont(ppl,n2(i)) = 1
1452 ENDIF
1453 ENDIF
1454 ENDDO
1455 ENDDO
1456
1457 ENDIF
1458
1459
1460
1461
1462 IF(iparit==0)THEN
1463 IF(kdtint==0)THEN
1465 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1466 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,
1467 4 fz2 ,fx3 ,fy3 ,fz3 ,fx4 ,
1468 5 fy4 ,fz4 ,a ,stifn,stif ,
1469 6 nedge,nin ,jtask,pene ,ibm )
1470 ELSE
1471 CALL i25asse05(jlt ,cs_loc ,n1 ,n2 ,m1 ,
1472 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1473 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,
1474 4 fz2 ,fx3 ,fy3 ,fz3 ,fx4 ,
1475 5 fy4 ,fz4 ,a ,stifn,nedge,
1476 6 k1 ,k2 ,k3 ,k4 ,c1 ,
1477 7 c2 ,c3 ,c4 ,viscn,nin ,
1478 8 jtask ,pene ,ibm )
1479 END IF
1480 ELSE
1481 IF(kdtint==0)THEN
1482 CALL i25asse2(jlt ,cs_loc ,n1 ,n2 ,m1 ,
1483 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1484 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,
1485 4 fz2 ,fx3 ,fy3 ,fz3 ,fx4 ,
1486 5 fy4 ,fz4 ,fskyi ,isky ,niskyfie,
1487 6 stif ,nedge ,nin ,noint ,pene ,
1488 7 ibm ,edge_id,tagip )
1489 ELSE
1490 CALL i25asse25(jlt ,cs_loc ,n1 ,n2 ,m1 ,
1491 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1492 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,
1493 4 fz2 ,fx3 ,fy3 ,fz3 ,fx4 ,
1494 5 fy4 ,fz4 ,isky ,niskyfie,nedge ,
1495 6 k1 ,k2 ,k3 ,k4 ,c1 ,
1496 7 c2 ,c3 ,c4 ,nin , noint,
1497 8 pene ,ibm ,tagip )
1498 END IF
1499 END IF
1500
1501 IF(idtmins==2)THEN
1502 dti=dt2t
1503 IF(iparit==0)THEN
1504 CALL i25sms0e(jlt ,cs_loc ,n1 ,n2 ,m1 ,
1505 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1506 3 stif ,nin ,noint ,mskyi_sms ,iskyi_sms,
1507 4 nsms ,k1 ,k2 ,k3 ,k4 ,
1508 5 c1 ,c2 ,c3 ,c4 ,nedge ,
1509 6 ibm )
1510 ELSE
1511 CALL i25sms2e(jlt ,cs_loc ,n1 ,n2 ,m1 ,
1512 2 m2 ,hs1 ,hs2 ,hm1 ,hm2 ,
1513 3 stif ,nin ,noint ,mskyi_sms ,iskyi_sms,
1514 4 nsms ,k1 ,k2 ,k3 ,k4 ,
1515 5 c1 ,c2 ,c3 ,c4 ,nedge ,
1516 6 ibm, edge_id)
1517 END IF
1518 END IF
1519
1520
1521
1522 IF (nspmd>1) THEN
1523
1524#include "mic_lockon.inc"
1525 DO i = 1,jlt
1526 IF(cs_loc(i)>nedge)THEN
1527 ni = cs_loc(i)-nedge
1528
1529
1530
1531
1532 IF(pene(i) /= zero.OR.tagip(i)==1) THEN
1534 ENDIF
1535 ENDIF
1536 ENDDO
1537
1538#include "mic_lockoff.inc"
1539 ENDIF
1540
1541
1542 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT>0)THEN
1543#include "lockon.inc"
1544
1545 DO i=1,jlt
1546
1547 IF(pene(i)==zero) cycle
1548
1549 IF(cs_loc(i)<=nedge) THEN
1550 fcont(1,n1(i)) =fcont(1,n1(i)) + fx1(i)
1551 fcont(2,n1(i)) =fcont(2,n1(i)) + fy1(i)
1552 fcont(3,n1(i)) =fcont(3,n1(i)) + fz1(i)
1553 fcont(1,n2(i)) =fcont(1,n2(i)) + fx2(i)
1554 fcont(2,n2(i)) =fcont(2,n2(i)) + fy2(i)
1555 fcont(3,n2(i)) =fcont(3,n2(i)) + fz2(i)
1556 END IF
1557 fcont(1,m1(i)) =fcont(1,m1(i)) + fx3(i)
1558 fcont(2,m1(i)) =fcont(2,m1(i)) + fy3(i)
1559 fcont(3,m1(i)) =fcont(3,m1(i)) + fz3(i)
1560 fcont(1,m2(i)) =fcont(1,m2(i)) + fx4(i)
1561 fcont(2,m2(i)) =fcont(2,m2(i)) + fy4(i)
1562 fcont(3,m2(i)) =fcont(3,m2(i)) + fz4(i)
1563 ENDDO
1564#include "lockoff.inc"
1565 ENDIF
1566
1567
1568 IF(isecin>0)THEN
1569 k0=nstrf(25)
1570 IF(nstrf(1)+nstrf(2)/=0)THEN
1571 DO i=1,nsect
1572 nbinter=nstrf(k0+14)
1573 k1s=k0+30
1574 DO j=1,nbinter
1575 IF(nstrf(k1s)==noint)THEN
1576 IF(isecut/=0)THEN
1577#include "lockon.inc"
1578 DO k=1,jlt
1579
1580 IF(pene(k) == zero)cycle
1581
1582 IF(cs_loc(k)<=nedge) THEN
1583 IF(secfcum(4,n1(k),i)==1.)THEN
1584 secfcum(1,n1(k),i)=secfcum(1,n1(k),i)-fx1(k)
1585 secfcum(2,n1(k),i)=secfcum(2,n1(k),i)-fy1(k)
1586 secfcum(3,n1(k),i)=secfcum(3,n1(k),i)-fz1(k)
1587 ENDIF
1588 IF(secfcum(4,n2(k),i)==1.)THEN
1589 secfcum(1,n2(k),i)=secfcum(1,n2(k),i)-fx2(k)
1590 secfcum(2,n2(k),i)=secfcum(2,n2(k),i)-fy2(k)
1591 secfcum(3,n2(k),i)=secfcum(3,n2(k),i)-fz2(k)
1592 ENDIF
1593 END IF
1594 IF(secfcum(4,m1(k),i)==1.)THEN
1595 secfcum(1,m1(k),i)=secfcum(1,m1(k),i)-fx3(k)
1596 secfcum(2,m1(k),i)=secfcum(2,m1(k),i)-fy3(k)
1597 secfcum(3,m1(k),i)=secfcum(3,m1(k),i)-fz3(k)
1598 ENDIF
1599 IF(secfcum(4,m2(k),i)==1.)THEN
1600 secfcum(1,m2(k),i)=secfcum(1,m2(k),i)-fx4(k)
1601 secfcum(2,m2(k),i)=secfcum(2,m2(k),i)-fy4(k)
1602 secfcum(3,m2(k),i)=secfcum(3,m2(k),i)-fz4(k)
1603 ENDIF
1604 ENDDO
1605#include "lockoff.inc"
1606 ENDIF
1607
1608 ENDIF
1609 k1s=k1s+1
1610 ENDDO
1611 k0=nstrf(k0+24)
1612 ENDDO
1613 ENDIF
1614 ENDIF
1615
1616 RETURN
integer function bitget(i, n)
subroutine i25asse05(jlt, cs_loc, n1, n2, m1, m2, hs1, hs2, hm1, hm2, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, a, stifn, nedge, k1, k2, k3, k4, c1, c2, c3, c4, viscn, nin, jtask, pene, ibm)
subroutine i25asse2(jlt, cs_loc, n1, n2, m1, m2, hs1, hs2, hm1, hm2, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fskyi, isky, niskyfie, stif, nedge, nin, noint, pene, ibm, edge_id, tagip)
subroutine i25asse25(jlt, cs_loc, n1, n2, m1, m2, hs1, hs2, hm1, hm2, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, isky, niskyfie, nedge, k1, k2, k3, k4, c1, c2, c3, c4, nin, noint, pene, ibm, tagip)
subroutine i25asse0(jlt, cs_loc, n1, n2, m1, m2, hs1, hs2, hm1, hm2, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, a, stifn, stif, nedge, nin, jtask, pene, ibm)
subroutine i25sms0e(jlt, cs_loc, n1, n2, m1, m2, hs1, hs2, hm1, hm2, stif, nin, noint, mskyi_sms, iskyi_sms, nsms, k1, k2, k3, k4, c1, c2, c3, c4, nrts, ibm)
subroutine i25sms2e(jlt, cs_loc, n1, n2, m1, m2, hs1, hs2, hm1, hm2, stif, nin, noint, mskyi_sms, iskyi_sms, nsms, k1, k2, k3, k4, c1, c2, c3, c4, nrts, ibm, edge_id)
type(int_pointer), dimension(:), allocatable nisubsfie
type(int_pointer), dimension(:), allocatable inflg_subsfie
type(int_pointer), dimension(:), allocatable lisubsfie
type(int_pointer), dimension(:), allocatable addsubsfie
type(int_pointer), dimension(:), allocatable nsnfie
type(int_pointer), dimension(:), allocatable nsvfie