33
34
35
37
38
39
40#include "implicit_f.inc"
41
42
43
44#include "mvsiz_p.inc"
45
46
47
49 . pm(npropm,*), ms(*), vol(*),wma(*), msq(*), partsav(20, *),
50 . x(3,*),v(3,*)
51 INTEGER MAT(*),IPART(*), IX1(*), IX2(*), IX3(*), IX4(*)
52
53
54
55#include "vect01_c.inc"
56#include "param_c.inc"
57
58
59
60 INTEGER I,IP ,I1 ,I2 ,I3 ,I4
62 . yy ,zz ,yz ,
63 . mass(mvsiz)
64
65 DO i=lft,llt
66 mass(i)=pm(89,mat(i))*vol(i)*fourth
67 ENDDO
68
69
70 DO i=lft,llt
71 msq(i)=mass(i)
72 ENDDO
73
74 IF(jale > 0 .AND.
ale%GRID%NWALE == 4)
THEN
75 DO i=lft,llt
76 wma(ix1(i))=wma(ix1(i))+one
77 wma(ix2(i))=wma(ix2(i))+one
78 wma(ix3(i))=wma(ix3(i))+one
79 wma(ix4(i))=wma(ix4(i))+one
80 ENDDO
81 ENDIF
82
83 DO i=lft,llt
84 ip = ipart(i)
85 partsav(1,ip)=partsav(1,ip) + four*mass(i)
86
87 i1 = ix1(i)
88 i2 = ix2(i)
89 i3 = ix3(i)
90 i4 = ix4(i)
91 yy = (x(2,i1)*x(2,i1)+x(2,i2)*x(2,i2)
92 . +x(2,i3)*x(2,i3)+x(2,i4)*x(2,i4))
93 yz = (x(2,i1)*x(3,i1)+x(2,i2)*x(3,i2)
94 . +x(2,i3)*x(3,i3)+x(2,i4)*x(3,i4))
95 zz = (x(3,i1)*x(3,i1)+x(3,i2)*x(3,i2)
96 . +x(3,i3)*x(3,i3)+x(3,i4)*x(3,i4))
97
98 partsav(5,ip) =partsav(5,ip) + mass(i) * (yy+zz)
99 partsav(6,ip) =partsav(6,ip) + mass(i) * zz
100 partsav(7,ip) =partsav(7,ip) + mass(i) * yy
101 partsav(9,ip) =partsav(9,ip) - mass(i) * yz
102
103 partsav(12,ip)=partsav(12,ip) + mass(i)*
104 . (v(2,i1)+v(2,i2)+v(2,i3)+v(2,i4))
105 partsav(13,ip)=partsav(13,ip) + mass(i)*
106 . (v(3,i1)+v(3,i2)+v(3,i3)+v(3,i4))
107 partsav(14,ip)=partsav(14,ip) + half * mass(i) *
108 . (v(2,i1)*v(2,i1)+v(3,i1)*v(3,i1)
109 . +v(2,i2)*v(2,i2)+v(3,i2)*v(3,i2)
110 . +v(2,i3)*v(2,i3)+v(3,i3)*v(3,i3)
111 . +v(2,i4)*v(2,i4)+v(3,i4)*v(3,i4))
112
113 ENDDO
114
115 RETURN
116