265
266
267
268#include "implicit_f.inc"
269
270
271
272#include "com01_c.inc"
273#include "com08_c.inc"
274
275
276
277 INTEGER NSN, IC, ICR, IFLAG
278 INTEGER NOD(*),WEIGHT(*)
280 . ms(*), in(*), x(3,*), a(3,*), ar(3,*), v(3,*), vr(3,*),
281 . xframe(*)
282 DOUBLE PRECISION FRL6(15,6)
283
284
285
286 INTEGER :: IC1, ICC, IC2, IC3, I, N, K
287 my_real :: mass, iner, ax, ay, az, vx, vy, vz,
288 . rx, ry, rz, sx, sy, sz,
289 . tx, ty, tz, ox, oy, oz, aa, mr2,vtmr2,atmr2,rx0, ry0, rz0,r,
290 . f1(nsn), f2(nsn), f3(nsn), f4(nsn), f5(nsn), f6(nsn), f7(nsn)
291
292
293 IF(iflag == 1)THEN
294
295
296 DO k = 1, 6
297 frl6(1,k) = zero
298 frl6(2,k) = zero
299 frl6(3,k) = zero
300 frl6(4,k) = zero
301 frl6(5,k) = zero
302 frl6(6,k) = zero
303 frl6(7,k) = zero
304 frl6(8,k) = zero
305 frl6(9,k) = zero
306 frl6(10,k) = zero
307 frl6(11,k) = zero
308 frl6(12,k) = zero
309 frl6(13,k) = zero
310 frl6(14,k) = zero
311 frl6(15,k) = zero
312 END DO
313
314 sx = xframe(1)
315 sy = xframe(2)
316 sz = xframe(3)
317 aa = one / sqrt(sx*sx + sy*sy + sz*sz)
318 sx = sx * aa
319 sy = sy * aa
320 sz = sz * aa
321 ox = xframe(10)
322 oy = xframe(11)
323 oz = xframe(12)
324
325
326
327
328
329
330 DO i=1,nsn
331 n = nod(i)
332 rx0 = x(1,n) + half * dt2 * v(1,n) - ox
333 ry0 = x(2,n) + half * dt2 * v(2,n) - oy
334 rz0 = x(3,n) + half * dt2 * v(3,n) - oz
335 tx = ry0 * sz - rz0 * sy
336 ty = rz0 * sx - rx0 * sz
337 tz = rx0 * sy - ry0 * sx
338 aa = one / sqrt(tx*tx + ty*ty + tz*tz)
339 tx = tx * aa
340 ty = ty * aa
341 tz = tz * aa
342 rx = sy * tz - sz * ty
343 ry = sz * tx - sx * tz
344 rz = sx * ty - sy * tx
345 r = rx * rx0 + ry * ry0 + rz * rz0
347
348
349
350 ax = a(1,n)
351 ay = a(2,n)
352 az = a(3,n)
353 a(1,n) = rx * ax + ry * ay + rz * az
354 a(2,n) = sx * ax + sy * ay + sz * az
355 a(3,n) = (tx * ax + ty * ay + tz * az) / r
356 ax = v(1,n)
357 ay = v(2,n)
358 az = v(3,n)
359 v(1,n) = rx * ax + ry * ay + rz * az
360 v(2,n) = sx * ax + sy * ay + sz * az
361 v(3,n) = (tx * ax + ty * ay + tz * az) / r
362 IF(iroddl/=0)THEN
363 ax = ar(1,n)
364 ay = ar(2,n)
365 az = ar(3,n)
366 ar(1,n) = rx * ax + ry * ay + rz * az
367 ar(2,n) = sx * ax + sy * ay + sz * az
368 ar(3,n) = tx * ax + ty * ay + tz * az
369 ax = vr(1,n)
370 ay = vr(2,n)
371 az = vr(3,n)
372 vr(1,n) = rx * ax + ry * ay + rz * az
373 vr(2,n) = sx * ax + sy * ay + sz * az
374 vr(3,n) = tx * ax + ty * ay + tz * az
375 ENDIF
376 IF(weight(n)==1) THEN
377 f1(i) = r*r*ms(n)
378 f2(i) = r*r*ms(n)*v(3,n)
379 f3(i) = r*r*ms(n)*a(3,n)
380
381
382
383 ELSE
384 f1(i) = zero
385 f2(i) = zero
386 f3(i) = zero
387 ENDIF
388 ENDDO
389
390 IF(ic/=0)THEN
391
392
393
397
398
399
400
401 ax =zero
402 ay =zero
403 vx =zero
404 vy =zero
405 mass=zero
406 DO i=1,nsn
407 n = nod(i)
408 IF(weight(n)==1) THEN
409 f1(i)=ms(n)
410 f2(i)=ms(n)*a(1,n)
411 f3(i)=ms(n)*a(2,n)
412 f4(i)=ms(n)*v(1,n)
413 f5(i)=ms(n)*v(2,n)
414 ELSE
415 f1(i)=zero
416 f2(i)=zero
417 f3(i)=zero
418 f4(i)=zero
419 f5(i)=zero
420 ENDIF
421 ENDDO
422
423
424
430 ENDIF
431
432
433
434 IF(iroddl/=0.AND.icr/=0)THEN
435
436
437
438
439
440
441
442
443
444
445
446 DO i=1,nsn
447 n = nod(i)
448 IF(weight(n)==1) THEN
449 f1(i)=in(n)
450 f2(i)=in(n)*ar(1,n)
451 f3(i)=in(n)*ar(2,n)
452 f4(i)=in(n)*ar(3,n)
453 f5(i)=in(n)*vr(1,n)
454 f6(i)=in(n)*vr(2,n)
455 f7(i)=in(n)*vr(3,n)
456 ENDIF
457 ENDDO
458
459
460
461
469 ENDIF
470
471 ELSEIF(iflag == 2)THEN
472 sx = xframe(1)
473 sy = xframe(2)
474 sz = xframe(3)
475 aa = one / sqrt(sx*sx + sy*sy + sz*sz)
476 sx = sx * aa
477 sy = sy * aa
478 sz = sz * aa
479 ox = xframe(10)
480 oy = xframe(11)
481 oz = xframe(12)
482
483 IF(ic/=0)THEN
484 ic1=ic/4
485 icc=ic-4*ic1
486 ic2=icc/2
487 ic3=icc-2*ic2
488
489 mr2 = frl6(1,1)+frl6(1,2)+frl6(1,3)+
490 + frl6(1,4)+frl6(1,5)+frl6(1,6)
491 vtmr2= frl6(2,1)+frl6(2,2)+frl6(2,3)+
492 + frl6(2,4)+frl6(2,5)+frl6(2,6)
493 atmr2= frl6(3,1)+frl6(3,2)+frl6(3,3)+
494 + frl6(3,4)+frl6(3,5)+frl6(3,6)
495 mass = frl6(4,1)+frl6(4,2)+frl6(4,3)+
496 + frl6(4,4)+frl6(4,5)+frl6(4,6)
497 ax = frl6(5,1)+frl6(5,2)+frl6(5,3)+
498 + frl6(5,4)+frl6(5,5)+frl6(5,6)
499 ay = frl6(6,1)+frl6(6,2)+frl6(6,3)+
500 + frl6(6,4)+frl6(6,5)+frl6(6,6)
501 vx = frl6(7,1)+frl6(7,2)+frl6(7,3)+
502 + frl6(7,4)+frl6(7,5)+frl6(7,6)
503 vy = frl6(8,1)+frl6(8,2)+frl6(8,3)+
504 + frl6(8,4)+frl6(8,5)+frl6(8,6)
505
506 ax= ax/mass
507 ay= ay/mass
508 az= atmr2/mr2
509 vx= vx/mass
510 vy= vy/mass
511 vz= vtmr2/mr2
512
513
514
515 DO i=1,nsn
516 n = nod(i)
517 a(1,n) =a(1,n)-ic1*(a(1,n)-ax)
518 a(2,n) =a(2,n)-ic2*(a(2,n)-ay)
519 a(3,n) =a(3,n)-ic3*(a(3,n)-az)
520
521 v(1,n) =v(1,n)-ic1*(v(1,n)-vx)
522 v(2,n) =v(2,n)-ic2*(v(2,n)-vy)
523 v(3,n) =v(3,n)-ic3*(v(3,n)-vz)
524
525 ENDDO
526 END IF
527
528
529
530 IF(iroddl/=0.AND.icr/=0)THEN
531 ic1=icr/4
532 icc=icr-4*ic1
533 ic2=icc/2
534 ic3=icc-2*ic2
535
536 iner= frl6(9,1)+frl6(9,2)+frl6(9,3)+
537 + frl6(9,4)+frl6(9,5)+frl6(9,6)
538 ax = frl6(10,1)+frl6(10,2)+frl6(10,3)+
539 + frl6(10,4)+frl6(10,5)+frl6(10,6)
540 ay = frl6(11,1)+frl6(11,2)+frl6(11,3)+
541 + frl6(11,4)+frl6(11,5)+frl6(11,6)
542 az = frl6(12,1)+frl6(12,2)+frl6(12,3)+
543 + frl6(12,4)+frl6(12,5)+frl6(12,6)
544 vx = frl6(13,1)+frl6(13,2)+frl6(13,3)+
545 + frl6(13,4)+frl6(13,5)+frl6(13,6)
546 vy = frl6(14,1)+frl6(14,2)+frl6(14,3)+
547 + frl6(14,4)+frl6(14,5)+frl6(14,6)
548 vz = frl6(15,1)+frl6(15,2)+frl6(15,3)+
549 + frl6(15,4)+frl6(15,5)+frl6(15,6)
550 IF(iner/=zero)THEN
551 ax=ax/iner
552 ay=ay/iner
553 az=az/iner
554 vx=vx/iner
555 vy=vy/iner
556 vz=vz/iner
557
558
559
560 DO i=1,nsn
561 n = nod(i)
562 ar(1,n) =ar(1,n)-ic1*(ar(1,n)-ax)
563 ar(2,n) =ar(2,n)-ic2*(ar(2,n)-ay)
564 ar(3,n) =ar(3,n)-ic3*(ar(3,n)-az)
565
566 vr(1,n) =vr(1,n)-ic1*(vr(1,n)-vx)
567 vr(2,n) =vr(2,n)-ic2*(vr(2,n)-vy)
568 vr(3,n) =vr(3,n)-ic3*(vr(3,n)-vz)
569
570 ENDDO
571 ENDIF
572 ENDIF
573
574
575
576 DO i=1,nsn
577 n = nod(i)
578 rx0 = x(1,n) + half * dt2 * v(1,n) - ox
579 ry0 = x(2,n) + half * dt2 * v(2,n) - oy
580 rz0 = x(3,n) + half * dt2 * v(3,n) - oz
581 tx = ry0 * sz - rz0 * sy
582 ty = rz0 * sx - rx0 * sz
583 tz = rx0 * sy - ry0 * sx
584 aa = one / sqrt(tx*tx + ty*ty + tz*tz)
585 tx = tx * aa
586 ty = ty * aa
587 tz = tz * aa
588 rx = sy * tz - sz * ty
589 ry = sz * tx - sx * tz
590 rz = sx * ty - sy * tx
591 r = rx * rx0 + ry * ry0 + rz * rz0
593
594
595
596 ax = rx * a(1,n) + sx * a(2,n) + tx * a(3,n) * r
597 ay = ry * a(1,n) + sy * a(2,n) + ty * a(3,n) * r
598 az = rz * a(1,n) + sz * a(2,n) + tz * a(3,n) * r
599 a(1,n) = ax
600 a(2,n) = ay
601 a(3,n) = az
602 ax = rx * v(1,n) + sx * v(2,n) + tx * v(3,n) * r
603 ay = ry * v(1,n) + sy * v(2,n) + ty * v(3,n) * r
604 az = rz * v(1,n) + sz * v(2,n) + tz * v(3,n) * r
605 v(1,n) = ax
606 v(2,n) = ay
607 v(3,n) = az
608 IF(iroddl/=0)THEN
609 ax = rx * ar(1,n) + sx * ar(2,n) + tx * ar(3,n)
610 ay = ry * ar(1,n) + sy * ar(2,n) + ty * ar(3,n)
611 az = rz * ar(1,n) + sz * ar(2,n) + tz * ar(3,n)
612 ar(1,n) = ax
613 ar(2,n) = ay
614 ar(3,n) = az
615 ax = rx * vr(1,n) + sx * vr(2,n) + tx * vr(3,n)
616 ay = ry * vr(1,n) + sy * vr(2,n) + ty * vr(3,n)
617 az = rz * vr(1,n) + sz * vr(2,n) + tz * vr(3,n)
618 vr(1,n) = ax
619 vr(2,n) = ay
620 vr(3,n) = az
621 ENDIF
622 ENDDO
623 END IF
624
625 RETURN
subroutine sum_6_float(jft, jlt, f, f6, n)