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