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