634
635
636
637#include "implicit_f.inc"
638
639
640
641#include "com08_c.inc"
642
643
644
645 INTEGER NSN, IC, IFLAG, IDOWN, PMAIN, TAG_LNK
646 INTEGER NOD(*),WEIGHT(*), ITAB(*)
647
649 . ms(*), x(3,*), a(3,*), v(3,*),
650 . xframe(*), frl(4)
651 DOUBLE PRECISION FRL6(5,6)
652
653
654
655 INTEGER IC1, ICC, IC2, IC3, I, N, K
656
658 . mass, ax, ay, az,
659 . rx, ry, rz, sx, sy, sz,
660 . tx, ty, tz, ox, oy, oz, aa, r2, atr2,rx0, ry0, rz0, r,
661 . mr2, atmr2,
662 . f1(nsn), f2(nsn), f3(nsn), f4(nsn), f5(nsn)
663
664 SELECT CASE(idown)
665
666 CASE(0)
667
668
669 sx = xframe(1)
670 sy = xframe(2)
671 sz = xframe(3)
672 aa = one / sqrt(sx*sx + sy*sy + sz*sz)
673 sx = sx * aa
674 sy = sy * aa
675 sz = sz * aa
676 ox = xframe(10)
677 oy = xframe(11)
678 oz = xframe(12)
679
680 IF(iflag == 1)THEN
681
682
683 DO k = 1, 6
684 frl6(1,k) = zero
685 frl6(2,k) = zero
686 frl6(3,k) = zero
687 frl6(4,k) = zero
688 END DO
689
690
691
692
693
694
695 DO i=1,nsn
696 n = nod(i)
697 rx0 = x(1,n) + half * dt2 * v(1,n) - ox
698 ry0 = x(2,n) + half * dt2 * v(2,n) - oy
699 rz0 = x(3,n) + half * dt2 * v(3,n) - oz
700 tx = ry0 * sz - rz0 * sy
701 ty = rz0 * sx - rx0 * sz
702 tz = rx0 * sy - ry0 * sx
703 aa = one / sqrt(tx*tx + ty*ty + tz*tz)
704 tx = tx * aa
705 ty = ty * aa
706 tz = tz * aa
707 rx = sy * tz - sz * ty
708 ry = sz * tx - sx * tz
709 rz = sx * ty - sy * tx
710 r = rx * rx0 + ry * ry0 + rz * rz0
712
713
714
715 ax = a(1,n)
716 ay = a(2,n)
717 az = a(3,n)
718 a(1,n) = rx * ax + ry * ay + rz * az
719 a(2,n) = sx * ax + sy * ay + sz * az
720 a(3,n) = (tx * ax + ty * ay + tz * az) / r
721 IF(weight(n)==1) THEN
722 f1(i) = r*r
723 f2(i) = r*r*a(3,n)
724
725
726 ELSE
727 f1(i) = zero
728 f2(i) = zero
729 ENDIF
730 ENDDO
731
732
733
736
737
738
739
740 ax =zero
741 ay =zero
742 mass=zero
743 DO i=1,nsn
744 n = nod(i)
745 IF(weight(n)==1) THEN
746 f3(i)=a(1,n)
747 f4(i)=a(2,n)
748 ELSE
749 f3(i)=zero
750 f4(i)=zero
751 ENDIF
752 ENDDO
753
754
755
758
759 ELSEIF(iflag == 2)THEN
760
761 ic1=ic/4
762 icc=ic-4*ic1
763 ic2=icc/2
764 ic3=icc-2*ic2
765
766 r2 = frl6(1,1)+frl6(1,2)+frl6(1,3)+
767 + frl6(1,4)+frl6(1,5)+frl6(1,6)
768 atr2 = frl6(2,1)+frl6(2,2)+frl6(2,3)+
769 + frl6(2,4)+frl6(2,5)+frl6(2,6)
770 ax = frl6(3,1)+frl6(3,2)+frl6(3,3)+
771 + frl6(3,4)+frl6(3,5)+frl6(3,6)
772 ay = frl6(4,1)+frl6(4,2)+frl6(4,3)+
773 + frl6(4,4)+frl6(4,5)+frl6(4,6)
774
775 az= atr2/r2
776 DO i=1,nsn
777 n = nod(i)
778
779
780 a(1,n) =a(1,n)-ic1*(a(1,n)-ax)
781 a(2,n) =a(2,n)-ic2*(a(2,n)-ay)
782 a(3,n) =a(3,n)-ic3*(a(3,n)-az)
783 END DO
784
785
786
787 DO i=1,nsn
788 n = nod(i)
789 rx0 = x(1,n) + half * dt2 * v(1,n) - ox
790 ry0 = x(2,n) + half * dt2 * v(2,n) - oy
791 rz0 = x(3,n) + half * dt2 * v(3,n) - oz
792 tx = ry0 * sz - rz0 * sy
793 ty = rz0 * sx - rx0 * sz
794 tz = rx0 * sy - ry0 * sx
795 aa = one / sqrt(tx*tx + ty*ty + tz*tz)
796 tx = tx * aa
797 ty = ty * aa
798 tz = tz * aa
799 rx = sy * tz - sz * ty
800 ry = sz * tx - sx * tz
801 rz = sx * ty - sy * tx
802 r = rx * rx0 + ry * ry0 + rz * rz0
804
805
806
807 ax = rx * a(1,n) + sx * a(2,n) + tx * a(3,n) * r
808 ay = ry * a(1,n) + sy * a(2,n) + ty * a(3,n) * r
809 az = rz * a(1,n) + sz * a(2,n) + tz * a(3,n) * r
810 a(1,n) = ax
811 a(2,n) = ay
812 a(3,n) = az
813 ENDDO
814
815 END IF
816 CASE(1)
817
818
819
820 sx = xframe(1)
821 sy = xframe(2)
822 sz = xframe(3)
823 aa = one / sqrt(sx*sx + sy*sy + sz*sz)
824 sx = sx * aa
825 sy = sy * aa
826 sz = sz * aa
827 ox = xframe(10)
828 oy = xframe(11)
829 oz = xframe(12)
830
831
832
833 IF(iflag == 1)THEN
834
835
836 DO k = 1, 6
837 frl6(1,k) = zero
838 frl6(2,k) = zero
839 frl6(3,k) = zero
840 frl6(4,k) = zero
841 frl6(5,k) = zero
842 END DO
843
844
845
846
847
848
849 DO i=1,nsn
850 n = nod(i)
851 rx0 = x(1,n) + half * dt2 * v(1,n) - ox
852 ry0 = x(2,n) + half * dt2 * v(2,n) - oy
853 rz0 = x(3,n) + half * dt2 * v(3,n) - oz
854 tx = ry0 * sz - rz0 * sy
855 ty = rz0 * sx - rx0 * sz
856 tz = rx0 * sy - ry0 * sx
857 aa = one / sqrt(tx*tx + ty*ty + tz*tz)
858 tx = tx * aa
859 ty = ty * aa
860 tz = tz * aa
861 rx = sy * tz - sz * ty
862 ry = sz * tx - sx * tz
863 rz = sx * ty - sy * tx
864 r = rx * rx0 + ry * ry0 + rz * rz0
866
867
868
869 ax = a(1,n)
870 ay = a(2,n)
871 az = a(3,n)
872 a(1,n) = rx * ax + ry * ay + rz * az
873 a(2,n) = sx * ax + sy * ay + sz * az
874 a(3,n) = (tx * ax + ty * ay + tz * az) / r
875 IF(weight(n)==1) THEN
876 f1(i) = r*r*ms(n)
877 f2(i) = r*r*ms(n)*a(3,n)
878
879
880 ELSE
881 f1(i) = zero
882 f2(i) = zero
883 ENDIF
884 ENDDO
885
886
887
890
891
892
893
894 ax =zero
895 ay =zero
896 mass=zero
897 DO i=1,nsn
898 n = nod(i)
899 IF(weight(n)==1) THEN
900 f3(i)=a(1,n)
901 f4(i)=a(2,n)
902 f5(i)=ms(n)
903 ELSE
904 f3(i)=zero
905 f4(i)=zero
906 f5(i)=zero
907 ENDIF
908 ENDDO
909
910
911
915
916
917
918
919 ic1=ic/4
920 icc=ic-4*ic1
921 ic2=icc/2
922 ic3=icc-2*ic2
923
924 mr2 = frl6(1,1)+frl6(1,2)+frl6(1,3)+
925 + frl6(1,4)+frl6(1,5)+frl6(1,6)
926 atmr2= frl6(2,1)+frl6(2,2)+frl6(2,3)+
927 + frl6(2,4)+frl6(2,5)+frl6(2,6)
928 ax = frl6(3,1)+frl6(3,2)+frl6(3,3)+
929 + frl6(3,4)+frl6(3,5)+frl6(3,6)
930 ay = frl6(4,1)+frl6(4,2)+frl6(4,3)+
931 + frl6(4,4)+frl6(4,5)+frl6(4,6)
932 mass = frl6(5,1)+frl6(5,2)+frl6(5,3)+
933 + frl6(5,4)+frl6(5,5)+frl6(5,6)
934
935 ax= ax/mass
936 ay= ay/mass
937 az= atmr2/mr2
938
939
940
941 DO i=1,nsn
942 n = nod(i)
943 a(1,n) =a(1,n)-ic1*(a(1,n)-ax)
944 a(2,n) =a(2,n)-ic2*(a(2,n)-ay)
945 a(3,n) =a(3,n)-ic3*(a(3,n)-az)
946 ENDDO
947
948
949
950 DO i=1,nsn
951 n = nod(i)
952 rx0 = x(1,n) + half * dt2 * v(1,n
953 ry0 = x(2,n) + half * dt2 * v(2,n) - oy
954 rz0 = x(3,n) + half * dt2 * v(3,n) - oz
955 tx = ry0 * sz - rz0 * sy
956 ty = rz0 * sx - rx0 * sz
957 tz = rx0 * sy - ry0 * sx
958 aa = one / sqrt(tx*tx + ty*ty + tz*tz)
959 tx = tx * aa
960 ty = ty * aa
961 tz = tz * aa
962 rx = sy * tz - sz * ty
963 ry = sz * tx - sx * tz
964 rz = sx * ty - sy * tx
965 r = rx * rx0 + ry * ry0 + rz * rz0
967
968
969
970 ax = rx * a(1,n) + sx * a(2,n) + tx * a(3,n) * r
971 ay = ry * a(1,n) + sy * a(2,n) + ty * a(3,n) * r
972 az = rz * a(1,n) + sz * a(2,n) + tz * a(3,n) * r
973 a(1,n) = ax
974 a(2,n) = ay
975 a(3,n) = az
976 END DO
977 END IF
978
979 END SELECT
980
981 RETURN