37
39
40
41
42#include "implicit_f.inc"
43
44
45
46 INTEGER NOD(*), NBY(*), ISKEW(*),ITAB(*), WEIGHT(*),
47 . WEIGHT_MD(*)
48
50 . v(3,*), vr(3,*), x(3,*), rby(*), skew(lskew,*), fs(*),
51 . a(3,*),ar(3,*),in(*),ms(*),enrot_t,encin_t,xmass_t,
52 . xmomt_t,ymomt_t,zmomt_t,encin2_t,enrot2_t,ms_2d(*)
53
54
55
56#include "com08_c.inc"
57#include "parit_c.inc"
58#include "impl1_c.inc"
59#include "param_c.inc"
60#include "com01_c.inc"
61
62
63
64 INTEGER M, NSN, ICODR, ISK, I, N,ISENS
65
67 . vi(3),vg(3),dt05,mas,wewe2
68
69 double precision
70 . enrott,encint,xmasst,xmomtt,ymomtt,zmomtt,
71 . encin2t,enrot2t
72
73 m =nby(1)
74
75 enrott=zero
76 encint=zero
77 xmasst=zero
78 xmomtt=zero
79 ymomtt=zero
80 zmomtt=zero
81 encin2t=zero
82 enrot2t=zero
83
84
85
86 IF (m<0) RETURN
87 nsn =nby(2)
88
89
90
91 IF(impl_s>0.AND.idyna==0)THEN
92 IF(isens==0)THEN
93 vg(1)=vr(1,m)
94 vg(2)=vr(2,m)
95 vg(3)=vr(3,m)
96 vi(1)=rby(1)*vg(1)+rby(2)*vg(2)+rby(3)*vg(3)
97 vi(2)=rby(4)*vg(1)+rby(5)*vg(2)+rby(6)*vg(3)
98 vi(3)=rby(7)*vg(1)+rby(8)*vg(2)+rby(9)*vg(3)
99
100 enrott= - ( vg(1)*vg(1)
101 . + vg(2)*vg(2)
102 . + vg(3)*vg(3))*in(m)*weight_md(m)
103 . + ( rby(10)*vi(1)*vi(1)
104 . + rby(11)*vi(2)*vi(2)
105 . + rby(12)*vi(3)*vi(3))*weight_md(m)
106 encint=zero
107 xmasst=zero
108 xmomtt=zero
109 ymomtt=zero
110 zmomtt=zero
111 encin2t=zero
112 enrot2t=zero
113
114 IF(n2d==0) THEN
115 IF (nsn>=10.OR.iparit>0) THEN
116 DO i=1,nsn
117 n = nod(i)
118 mas=ms(n)*weight_md(n)
119 wewe2 = (1-weight_md(n))*weight(n)
120 vg(1)=v(1,n)
121 vg(2)=v(2,n)
122 vg(3)=v(3,n)
123 encint=encint - ( vg(1)*vg(1)
124 . + vg(2)*vg(2)
125 . + vg(3)*vg(3))*mas
126 encin2t=encin2t - ( vg(1)*vg(1)
127 . + vg(2)*vg(2)
128 . + vg(3)*vg(3))*ms(n)*wewe2
129 xmomtt=xmomtt-vg(1)*mas
130 ymomtt=ymomtt-vg(2)*mas
131 zmomtt=zmomtt-vg(3)*mas
132 vg(1)=vr(1,n)
133 vg(2)=vr(2,n)
134 vg(3)=vr(3,n)
135 enrott=enrott - ( vg(1)*vg(1)
136 . + vg(2)*vg(2)
137 . + vg(3)*vg(3))*in(n)*weight_md(n)
138 enrot2t=enrot2t - ( vg(1)*vg(1)
139 . + vg(2)*vg(2)
140 . + vg(3)*vg(3))*in(n)*wewe2
141 xmasst=xmasst-mas
142 ENDDO
143 ELSE
144 DO i=1,nsn
145 n = nod(i)
146 mas=ms(n)*weight_md(n)
147 wewe2 = (1-weight_md(n))*weight(n)
148 vg(1)=v(1,n)
149 vg(2)=v(2,n)
150 vg(3)=v(3,n)
151 encint=encint - ( vg(1)*vg(1)
152 . + vg(2)*vg(2)
153 . + vg(3)*vg(3))*mas
154 encin2t=encin2t - ( vg(1)*vg(1)
155 . + vg(2)*vg(2)
156 . + vg(3)*vg(3))*ms(n)*wewe2
157 xmomtt=xmomtt-vg(1)*mas
158 ymomtt=ymomtt-vg(2)*mas
159 zmomtt=zmomtt-vg(3)*mas
160 vg(1)=vr(1,n)
161 vg(2)=vr(2,n)
162 vg(3)=vr(3,n)
163 enrott=enrott - ( vg(1)*vg(1)
164 . + vg(2)*vg(2)
165 . + vg(3)*vg(3))*in(n)*weight_md(n)
166 enrot2t=enrot2t - ( vg(1)*vg(1)
167 . + vg(2)*vg(2)
168 . + vg(3)*vg(3))*in(n)*wewe2
169 xmasst=xmasst-mas
170
171 ENDDO
172 ENDIF
173
174 ELSE
175
176 IF (nsn>=10.OR.iparit>0) THEN
177 DO i=1,nsn
178 n = nod(i)
179 mas=ms_2d(n)*weight_md(n)
180 wewe2 = (1-weight_md(n))*weight(n)
181 vg(1)=v(1,n)
182 vg(2)=v(2,n)
183 vg(3)=v(3,n)
184 encint=encint - ( vg(1)*vg(1)
185 . + vg(2)*vg(2)
186 . + vg(3)*vg(3))*mas
187 encin2t=encin2t - ( vg(1)*vg(1)
188 . + vg(2)*vg(2)
189 . + vg(3)*vg(3))*ms_2d(n)*wewe2
190 xmomtt=xmomtt-vg(1)*mas
191 ymomtt=ymomtt-vg(2)*mas
192 zmomtt=zmomtt-vg(3)*mas
193 vg(1)=vr(1,n)
194 vg(2)=vr(2,n)
195 vg(3)=vr(3,n)
196 enrott=enrott - ( vg(1)*vg(1)
197 . + vg(2)*vg(2)
198 . + vg(3)*vg(3))*in(n)*weight_md(n)
199 enrot2t=enrot2t - ( vg(1)*vg(1)
200 . + vg(2)*vg(2)
201 . + vg(3)*vg(3))*in(n)*wewe2
202 xmasst=xmasst-mas
203 ENDDO
204 ELSE
205 DO i=1,nsn
206 n = nod(i)
207 mas=ms_2d(n)*weight_md(n)
208 wewe2 = (1-weight_md(n))*weight(n)
209 vg(1)=v(1,n)
210 vg(2)=v(2,n)
211 vg(3)=v(3,n)
212 encint=encint - ( vg(1)*vg(1)
213 . + vg(2)*vg(2)
214 . + vg(3)*vg(3))*mas
215 encin2t=encin2t - ( vg(1)*vg(1)
216 . + vg(2)*vg(2)
217 . + vg(3)*vg(3))*ms_2d(n)*wewe2
218 xmomtt=xmomtt-vg(1)*mas
219 ymomtt=ymomtt-vg(2)*mas
220 zmomtt=zmomtt-vg(3)*mas
221 vg(1)=vr(1,n)
222 vg(2)=vr(2,n)
223 vg(3)=vr(3,n)
224 enrott=enrott - ( vg(1)*vg(1)
225 . + vg(2)*vg(2)
226 . + vg(3)*vg(3))*in(n)*weight_md(n)
227 enrot2t=enrot2t - ( vg(1)*vg(1)
228 . + vg(2)*vg(2)
229 . + vg(3)*vg(3))*in(n)*wewe2
230 xmasst=xmasst-mas
231
232 ENDDO
233 ENDIF
234
235 ENDIF
236 ELSE
237
238 IF(n2d==0) THEN
239 mas=ms(m)*weight_md(m)
240 ELSE
241 mas=ms_2d(m)*weight_md(m)
242 ENDIF
243 vg(1)=vr(1,m)
244 vg(2)=vr(2,m)
245 vg(3)=vr(3,m)
246
247 enrott= - ( vg(1)*vg(1)
248 . + vg(2)*vg(2)
249 . + vg(3)*vg(3))*in(m)*weight_md(m)
250
251
252 vg(1)=v(1,m)
253 vg(2)=v(2,m)
254 vg(3)=v(3,m)
255 encint= - ( vg(1)*vg(1)
256 . + vg(2)*vg(2)
257 . + vg(3)*vg(3))*mas
258 xmasst=-mas
259 xmomtt=-vg(1)*mas
260 ymomtt=-vg(2)*mas
261 zmomtt=-vg(3)*mas
262
263 ENDIF
264 ELSE
265 IF(isens==0)THEN
266 dt05 = half*dt1
267 IF(idyna>0) dt05=(dy_g-one)*dt1
268 vg(1)=vr(1,m)+ar(1,m)*dt05
269 vg(2)=vr(2,m)+ar(2,m)*dt05
270 vg(3)=vr(3,m)+ar(3,m)*dt05
271 vi(1)=rby(1)*vg(1)+rby(2)*vg(2)+rby(3)*vg(3)
272 vi(2)=rby(4)*vg(1)+rby(5)*vg(2)+rby(6)*vg(3)
273 vi(3)=rby(7)*vg(1)+rby(8)*vg(2)+rby(9)*vg(3)
274
275 enrott= - ( vg(1)*vg(1)
276 . + vg(2)*vg(2)
277 . + vg(3)*vg(3))*in(m)*weight_md(m)
278 . + ( rby(10)*vi(1)*vi(1)
279 . + rby(11)*vi(2)*vi(2)
280 . + rby(12)*vi(3)*vi(3))*weight_md(m)
281 encint=zero
282 xmasst=zero
283 xmomtt=zero
284 ymomtt=zero
285 zmomtt=zero
286 encin2t=zero
287 enrot2t=zero
288
289
290 IF(n2d==0) THEN
291 IF (nsn>=10.OR.iparit>0) THEN
292 DO i=1,nsn
293
294 n = nod(i)
295
296 mas=ms(n)*weight_md(n)
297 wewe2 = (1-weight_md(n))*weight(n)
298 vg(1)=v(1,n)+a(1,n)*dt05
299 vg(2)=v(2,n)+a(2,n)*dt05
300 vg(3)=v(3,n)+a(3,n)*dt05
301 encint=encint - ( vg(1)*vg(1)
302 . + vg(2)*vg(2)
303 . + vg(3)*vg(3))*mas
304 encin2t=encin2t - ( vg(1)*vg(1)
305 . + vg(2)*vg(2)
306 . + vg(3)*vg(3))*ms(n)*wewe2
307 xmomtt=xmomtt-vg(1)*mas
308 ymomtt=ymomtt-vg(2)*mas
309 zmomtt=zmomtt-vg(3)*mas
310 vg(1)=vr(1,n)+ar(1,n)*dt05
311 vg(2)=vr(2,n)+ar(2,n)*dt05
312 vg(3)=vr(3,n)+ar(3,n)*dt05
313 enrott=enrott - ( vg(1)*vg(1)
314 . + vg(2)*vg(2)
315 . + vg(3)*vg(3))*in(n)*weight_md(n)
316 enrot2t=enrot2t - ( vg(1)*vg(1)
317 . + vg(2)*vg(2)
318 . + vg(3)*vg(3))*in(n)*wewe2
319 xmasst=xmasst-mas
320
321 ENDDO
322 ELSE
323 DO i=1,nsn
324
325 n = nod(i)
326
327 mas=ms(n)*weight_md(n)
328 wewe2 = (1-weight_md(n))*weight(n)
329 vg(1)=v(1,n)+a(1,n)*dt05
330 vg(2)=v(2,n)+a(2,n)*dt05
331 vg(3)=v(3,n)+a(3,n)*dt05
332 encint=encint - ( vg(1)*vg(1)
333 . + vg(2)*vg(2)
334 . + vg(3)*vg(3))*mas
335 encin2t=encin2t - ( vg(1)*vg(1)
336 . + vg(2)*vg(2)
337 . + vg(3)*vg(3))*ms(n)*wewe2
338 xmomtt=xmomtt-vg(1)*mas
339 ymomtt=ymomtt-vg(2)*mas
340 zmomtt=zmomtt-vg(3)*mas
341 vg(1)=vr(1,n)+ar(1,n)*dt05
342 vg(2)=vr(2,n)+ar(2,n)*dt05
343 vg(3)=vr(3,n)+ar(3,n)*dt05
344 enrott=enrott - ( vg(1)*vg(1)
345 . + vg(2)*vg(2)
346 . + vg(3)*vg(3))*in(n)*weight_md(n)
347 enrot2t=enrot2t - ( vg(1)*vg(1)
348 . + vg(2)*vg(2)
349 . + vg(3)*vg(3))*in(n)*wewe2
350 xmasst=xmasst-mas
351
352 ENDDO
353 ENDIF
354
355 ELSE
356
357 IF (nsn>=10.OR.iparit>0) THEN
358 DO i=1,nsn
359
360 n = nod(i)
361
362 mas=ms_2d(n)*weight_md(n)
363 wewe2 = (1-weight_md(n))*weight(n)
364 vg(1)=v(1,n)+a(1,n)*dt05
365 vg(2)=v(2,n)+a(2,n)*dt05
366 vg(3)=v(3,n)+a(3,n)*dt05
367 encint=encint - ( vg(1)*vg(1)
368 . + vg(2)*vg(2)
369 . + vg(3)*vg(3))*mas
370 encin2t=encin2t - ( vg(1)*vg(1)
371 . + vg(2)*vg(2)
372 . + vg(3)*vg(3))*ms_2d(n)*wewe2
373 xmomtt=xmomtt-vg(1)*mas
374 ymomtt=ymomtt-vg(2)*mas
375 zmomtt=zmomtt-vg(3)*mas
376 vg(1)=vr(1,n)+ar(1,n)*dt05
377 vg(2)=vr(2,n)+ar(2,n)*dt05
378 vg(3)=vr(3,n)+ar(3,n)*dt05
379 enrott=enrott - ( vg(1)*vg(1)
380 . + vg(2)*vg(2)
381 . + vg(3)*vg(3))*in(n)*weight_md(n)
382 enrot2t=enrot2t - ( vg(1)*vg(1)
383 . + vg(2)*vg(2)
384 . + vg(3)*vg(3))*in(n)*wewe2
385 xmasst=xmasst-mas
386
387 ENDDO
388 ELSE
389 DO i=1,nsn
390
391 n = nod(i)
392
393 mas=ms_2d(n)*weight_md(n)
394 wewe2 = (1-weight_md(n))*weight(n)
395 vg(1)=v(1,n)+a(1,n)*dt05
396 vg(2)=v(2,n)+a(2,n)*dt05
397 vg(3)=v(3,n)+a(3,n)*dt05
398 encint=encint - ( vg(1)*vg(1)
399 . + vg(2)*vg(2)
400 . + vg(3)*vg(3))*mas
401 encin2t=encin2t - ( vg(1)*vg(1)
402 . + vg(2)*vg(2)
403 . + vg(3)*vg(3))*ms_2d(n)*wewe2
404 xmomtt=xmomtt-vg(1)*mas
405 ymomtt=ymomtt-vg(2)*mas
406 zmomtt=zmomtt-vg(3)*mas
407 vg(1)=vr(1,n)+ar(1,n)*dt05
408 vg(2)=vr(2,n)+ar(2,n)*dt05
409 vg(3)=vr(3,n)+ar(3,n)*dt05
410 enrott=enrott - ( vg(1)*vg(1)
411 . + vg(2)*vg(2)
412 . + vg(3)*vg(3))*in(n)*weight_md(n)
413 enrot2t=enrot2t - ( vg(1)*vg(1)
414 . + vg(2)*vg(2)
415 . + vg(3)*vg(3))*in(n)*wewe2
416 xmasst=xmasst-mas
417
418 ENDDO
419 ENDIF
420
421 ENDIF
422 ELSE
423
424 dt05 = half*dt1
425 IF(idyna>0) dt05=(dy_g-one)*dt1
426 vg(1)=vr(1,m)+ar(1,m)*dt05
427 vg(2)=vr(2,m)+ar(2,m)*dt05
428 vg(3)=vr(3,m)+ar(3,m)*dt05
429
430 enrott= - ( vg(1)*vg(1)
431 . + vg(2)*vg(2)
432 . + vg(3)*vg(3))*in(m)*weight_md(m)
433
434 IF(n2d==0) THEN
435 mas=ms(m)*weight_md(m)
436 ELSE
437 mas=ms_2d(m)*weight_md(m)
438 ENDIF
439 vg(1)=v(1,m)+a(1,m)*dt05
440 vg(2)=v(2,m)+a(2,m)*dt05
441 vg(3)=v(3,m)+a(3,m)*dt05
442 encint= - ( vg(1)*vg(1)
443 . + vg(2)*vg(2)
444 . + vg(3)*vg(3))*mas
445 xmasst=-mas
446 xmomtt=-vg(1)*mas
447 ymomtt=-vg(2)*mas
448 zmomtt=-vg(3)*mas
449 ENDIF
450 ENDIF
451
452 enrot_t=enrot_t + enrott*half
453 encin_t=encin_t + encint*half
454 enrot2_t=enrot2_t + enrot2t*half
455 encin2_t=encin2_t + encin2t*half
456 xmass_t=xmass_t + xmasst
457 xmomt_t=xmomt_t + xmomtt
458 ymomt_t=ymomt_t + ymomtt
459 zmomt_t=zmomt_t + zmomtt
460
461 RETURN