30
31
32
33#include "implicit_f.inc"
34
35
36
37 INTEGER NME
39 . cr(6,*), sr(*), dt1, dt2, fxbrpm(*), fxbvit(*)
40
41
42
43 INTEGER I,II
45 . v1(3),v2(3),v3(3),v4(3),v5(3),v6(3)
47 . d1,d2,d3,d4,d5,d6
48
49
50
51 v1(1)=fxbrpm(2)
52 v1(2)=fxbrpm(5)
53 v1(3)=fxbrpm(8)
54 v2(1)=fxbrpm(3)
55 v2(2)=fxbrpm(6)
56 v2(3)=fxbrpm(9)
57 v3(1)=fxbrpm(4)
58 v3(2)=fxbrpm(7)
59 v3(3)=fxbrpm(10)
60 v4(1)=v2(1)-v1(1)
61 v4(2)=v2(2)-v1(2)
62 v4(3)=v2(3)-v1(3)
63 v5(1)=v3(1)-v1(1)
64 v5(2)=v3(2)-v1(2)
65 v5(3)=v3(3)-v1(3)
66 v6(1)=v3(1)-v2(1)
67 v6(2)=v3(2)-v2(2)
68 v6(3)=v3(3)-v2(3)
69
70 d1=zero
71 d2=zero
72 d3=zero
73 d4=zero
74 d5=zero
75 d6=zero
76 DO i=1,3
77 d1=d1+v1(i)**2
78 d2=d2+v2(i)**2
79 d3=d3+v3(i)**2
80 d4=d4+v4(i)**2
81 d5=d5+v5(i)**2
82 d6=d6+v6(i)**2
83 END DO
84
85
86
87 DO i=1,6
88 DO ii=1,nme
89 cr(i,ii)=zero
90 ENDDO
91 ENDDO
92
93 cr(1,10)=-v1(1)-dt2*(fxbvit(1)-fxbvit(10))
94 cr(1,11)=-v1(2)-dt2*(fxbvit(2)-fxbvit(11))
95 cr(1,12)=-v1(3)-dt2*(fxbvit(3)-fxbvit(12))
96 cr(1,1)=-cr(1,10)
97 cr(1,2)=-cr(1,11)
98 cr(1,3)=-cr(1,12)
99
100 cr(2,10)=-v2(1)-dt2*(fxbvit(4)-fxbvit(10))
101 cr(2,11)=-v2(2)-dt2*(fxbvit(5)-fxbvit(11))
102 cr(2,12)=-v2(3)-dt2*(fxbvit(6)-fxbvit(12))
103 cr(2,4)=-cr(2,10)
104 cr(2,5)=-cr(2,11)
105 cr(2,6)=-cr(2,12)
106
107 cr(3,10)=-v3(1)-dt2*(fxbvit(7)-fxbvit(10))
108 cr(3,11)=-v3(2)-dt2*(fxbvit(8)-fxbvit(11))
109 cr(3,12)=-v3(3)-dt2*(fxbvit(9)-fxbvit(12))
110 cr(3,7)=-cr(3,10)
111 cr(3,8)=-cr(3,11)
112 cr(3,9)=-cr(3,12)
113
114 cr(4,1)=-v4(1)-dt2*(fxbvit(4)-fxbvit(1))
115 cr(4,2)=-v4(2)-dt2*(fxbvit(5)-fxbvit(2))
116 cr(4,3)=-v4(3)-dt2*(fxbvit(6)-fxbvit(3))
117 cr(4,4)=-cr(4,1)
118 cr(4,5)=-cr(4,2)
119 cr(4,6)=-cr(4,3)
120
121 cr(5,1)=-v5(1)-dt2*(fxbvit(7)-fxbvit(1))
122 cr(5,2)=-v5(2)-dt2*(fxbvit(8)-fxbvit(2))
123 cr(5,3)=-v5(3)-dt2*(fxbvit(9)-fxbvit(3))
124 cr(5,7)=-cr(5,1)
125 cr(5,8)=-cr(5,2)
126 cr(5,9)=-cr(5,3)
127
128 cr(6,4)=-v6(1)-dt2*(fxbvit(7)-fxbvit(4))
129 cr(6,5)=-v6(2)-dt2*(fxbvit(8)-fxbvit(5))
130 cr(6,6)=-v6(3)-dt2*(fxbvit(9)-fxbvit(6))
131 cr(6,7)=-cr(6,4)
132 cr(6,8)=-cr(6,5)
133 cr(6,9)=-cr(6,6)
134
135
136
137 DO i=1,6
138 sr(i)=zero
139 ENDDO
140 DO i=1,3
141 sr(1)=sr(1)-one/(dt1+dt2)*(two*v1(i)+
142 . dt2*(fxbvit(i)-fxbvit(9+i)))
143 . *(fxbvit(i)-fxbvit(9+i))
144 sr(2)=sr(2)-one/(dt1+dt2)*(two*v2(i)+
145 . dt2*(fxbvit(3+i)-fxbvit(9+i)))
146 . *(fxbvit(3+i)-fxbvit(9+i))
147 sr(3)=sr(3)-one/(dt1+dt2)*(two*v3(i)+
148 . dt2*(fxbvit(6+i)-fxbvit(9+i)))
149 . *(fxbvit(6+i)-fxbvit(9+i))
150 sr(4)=sr(4)-one/(dt1+dt2)*(two*v4(i)+
151 . dt2*(fxbvit(3+i)-fxbvit(i)))
152 . *(fxbvit(3+i)-fxbvit(i))
153 sr(5)=sr(5)-one/(dt1+dt2)*(two*v5(i)+
154 . dt2*(fxbvit(6+i)-fxbvit(i)))
155 . *(fxbvit(6+i)-fxbvit(i))
156 sr(6)=sr(6)-one/(dt1+dt2)*(two*v6(i)+
157 . dt2*(fxbvit(6+i)-fxbvit(3+i)))
158 . *(fxbvit(6+i)-fxbvit(3+i))
159 END DO
160 sr(1)=sr(1)+one/(dt2*(dt1+dt2))*(one-d1)
161 sr(2)=sr(2)+one/(dt2*(dt1+dt2))*(one-d2)
162 sr(3)=sr(3)+one/(dt2*(dt1+dt2))*(one-d3)
163 sr(4)=sr(4)+one/(dt2*(dt1+dt2))*(two-d4)
164 sr(5)=sr(5)+one/(dt2*(dt1+dt2))*(two-d5)
165 sr(6)=sr(6)+one/(dt2*(dt1+dt2))*(two-d6)
166
167 RETURN