41
42
43
45 use damping_vref_compute_dampa_mod
46
47
48
49#include "implicit_f.inc"
50#include "comlock.inc"
51
52
53
54#include "com01_c.inc"
55#include "com04_c.inc"
56#include "com06_c.inc"
57#include "com08_c.inc"
58#include "sms_c.inc"
59#include "param_c.inc"
60
61
62
63 INTEGER DIM,ITASK,WEIGHT(*),ICONTACT(*),I_DAMP_RDOF_TAB(*),TAGSLV_RBY(*)
64 INTEGER ,INTENT(IN) :: NDAMP_VREL,ID_DAMP_VREL(NDAMP_VREL)
65 INTEGER ,INTENT(IN) :: FR_DAMP_VREL(NSPMD+2,NDAMP_VREL),IPARIT,ISPMD
66 my_real v(3,*), vr(3,*), a(3,*), ar(3,*) ,ms(*), in(*),dampr(nrdamp,*), damp(dim,*), skew(lskew,*)
67 my_real,
DIMENSION(:,:),
ALLOCATABLE :: vskw,askw,dampskw
68 DOUBLE PRECISION,INTENT(INOUT):: WFEXT
69
70 TYPE (GROUP_) , DIMENSION(NGRNOD) :: IGRNOD
71
72
73
74 LOGICAL :: NEED_COMPUTE
75 INTEGER I,J,N,ND,IGR,ISK,FL_VREL,ID_RBY,FL_FREQ_RANGE,ITYPE
76 my_real,
DIMENSION(3) :: damp_a,damp_b
77 my_real,
DIMENSION(3) :: da_,dw_local,omega_,betasdt_
78 my_real :: aa,dw,da,betasdt,omega,factb,dampt,d_tstart,d_tstop,
79 . v_ref(3,ndamp),a_ref(3,ndamp),v_refskw(3),a_refskw(3),
80 . dam(3),darm(3),dfm,get_u_func
81 EXTERNAL get_u_func
82
83
84
85 ALLOCATE(vskw(3,numnod))
86 ALLOCATE(askw(3,numnod))
87 ALLOCATE(dampskw(3,numnod))
88
89 IF(idtmins==2.OR.idtmins_int/=0)GOTO 1000
90
91
92
93 IF (ndamp_vrel > 0) THEN
95 . weight,v_ref,a_ref,ndamp_vrel,id_damp_vrel,
96 . fr_damp_vrel,ngrnod,numnod,ndamp,nrdamp,
97 . nspmd,iparit,ispmd)
98 ENDIF
99
100 dw = zero
101 DO nd=1,ndamp
102 igr = nint(dampr(2,nd))
103 isk = nint(dampr(15,nd))
104 factb = one
105
106 dampt =
min(dt1,dt2)*factb
107 d_tstart = dampr(17,nd)
108 d_tstop = dampr(18,nd)
109 itype = nint(dampr(21,nd))
110 fl_vrel = 0
111 fl_freq_range = 0
112 SELECT CASE (itype)
113 CASE(2)
114 fl_vrel = 1
115 CASE(3)
116 fl_freq_range = 1
117 END SELECT
118 id_rby = nint(dampr(25,nd))
119
120 IF (tt>=d_tstart .AND. tt<=d_tstop) THEN
121
122 IF(isk<=1)THEN
123
124
125
126
127 IF ((fl_vrel == 0).AND.(fl_freq_range == 0)) THEN
128
129 IF (dampr(19,nd)==0) THEN
130
131
132 damp_a(1) = dampr(3,nd)
133 damp_b(1) = dampr(4,nd)
134 betasdt_(1)= -
min(damp_b(1),dampt)*dt1/
max(dt1*dt1,em30)
135 omega_(1) = one/ (one + half * damp_a(1) * dt1)
136
137 damp_a(2) = dampr(5,nd)
138 damp_b(2) = dampr(6,nd)
139 betasdt_(2)= -
min(damp_b(2),dampt)*dt1/
max(dt1*dt1,em30)
140 omega_(2) = one/ (one + half * damp_a(2) * dt1)
141
142 damp_a(3) = dampr(7,nd)
143 damp_b(3) = dampr(8,nd)
144 betasdt_(3)= -
min(damp_b(3),dampt)*dt1/
max(dt1*dt1,em30)
145 omega_(3) = one/ (one + half * damp_a(3) * dt1)
146
147
148 DO n=1,igrnod(igr)%NENTITY
149 i=igrnod(igr)%ENTITY(n)
150 IF(tagslv_rby(i)==0) THEN
151
152 da_(1) = a(1,i) - damp_a(1)*v(1,i) - betasdt_(1) *(a(1,i) - damp(1,i))
153 da_(1) = da_(1) * omega_(1) - a(1,i)
154 damp(1,i) = a(1,i)
155 a(1,i) = a(1,i) + da_(1)
156
157 da_(2) = a(2,i) - damp_a(2)*v(2,i) - betasdt_(2) *(a(2,i) - damp(2,i))
158 da_(2) = da_(2) * omega_(2) - a(2,i)
159 damp(2,i) = a(2,i)
160 a(2,i) = a(2,i) + da_(2)
161
162 da_(3) = a(3,i) - damp_a(3)*v(3,i) - betasdt_(3) *(a(3,i) - damp(3,i))
163 da_(3) = da_(3) * omega_(3) - a(3,i)
164 damp(3,i) = a(3,i)
165 a(3,i) = a(3,i) + da_(3)
166
167 dw_local(1) = da_(1)*(v(1,i)+half*a(1,i)*dt1)
168 dw_local(2) = da_(2)*(v(2,i)+half*a(2,i)*dt1)
169 dw_local(3) = da_(3)*(v(3,i)+half*a(3,i)*dt1)
170 dw =dw+ms(i)*(dw_local(1)+dw_local(2)+dw_local(3))*dt12*weight(i)
171 ENDIF
172 ENDDO
173 ENDIF
174
175 IF(iroddl/=0)THEN
176
177 damp_a(1) = dampr(9,nd)
178 damp_b(1) = dampr(10,nd)
179 betasdt_(1)= -
min(damp_b(1),dampt)*dt1/
max(dt1*dt1,em30)
180 omega_(1) = one/ (one + half * damp_a(1) * dt1)
181
182 damp_a(2) = dampr(11,nd)
183 damp_b(2) = dampr(12,nd)
184 betasdt_(2)= -
min(damp_b(2),dampt)*dt1/
max(dt1*dt1,em30)
185 omega_(2) = one/ (one + half * damp_a(2) * dt1)
186
187 damp_a(3) = dampr(13,nd)
188 damp_b(3) = dampr(14,nd)
189 betasdt_(3)= -
min(damp_b(3),dampt)*dt1/
max(dt1*dt1,em30)
190 omega_(3) = one/ (one + half * damp_a(3) * dt1)
191
192
193 DO n=1,igrnod(igr)%NENTITY
194 i=igrnod(igr)%ENTITY(n)
195 IF(tagslv_rby(i)==0) THEN
196
197 need_compute = .true.
198 IF (dampr(19,nd)/=0) THEN
199 IF (icontact(i)/=0) i_damp_rdof_tab(i) = dampr(19,nd)
200 IF (i_damp_rdof_tab(i)==0) need_compute = .false.
201 ENDIF
202
203 IF(need_compute) THEN
204
205 da_(1) = ar(1,i) - damp_a(1)*vr(1,i) - betasdt_(1) *(ar(1,i)-damp(4,i))
206 da_(1) = da_(1) * omega_(1) - ar(1,i)
207 damp(4,i) = ar(1,i)
208 ar(1,i) = ar(1,i) + da_(1)
209
210 da_(2) = ar(2,i) - damp_a(2)*vr(2,i) - betasdt_(2) *(ar(2,i)-damp(5,i))
211 da_(2) = da_(2) * omega_(2) - ar(2,i)
212 damp(5,i) = ar(2,i)
213 ar(2,i) = ar(2,i) + da_(2)
214
215 da_(3) = ar(3,i) - damp_a(3)*vr(3,i) - betasdt_(3) *(ar(3,i)-damp(6,i))
216 da_(3) = da_(3) * omega_(3) - ar(3,i)
217 damp(6,i) = ar(3,i)
218 ar(3,i) = ar(3,i) + da_(3)
219
220 dw_local(1) = da_(1)*(vr(1,i)+half*ar(1,i)*dt1)
221 dw_local(2) = da_(2)*(vr(2,i)+half*ar(2,i)*dt1)
222 dw_local(3) = da_(3)*(vr(3,i)+half*ar(3,i)*dt1)
223 dw = dw+in(i)*(dw_local(1)+dw_local(2)+dw_local(3)) *dt12*weight(i)
224
225 IF (dampr(19,nd)/=0) i_damp_rdof_tab(i)=i_damp_rdof_tab(i)-1
226 ENDIF
227 ENDIF
228 ENDDO
229 END IF
230
231 ELSEIF ((fl_vrel == 1).AND.(id_rby == 0)) THEN
232
233
234
235 call damping_vref_compute_dampa(nd,ndamp,nrdamp,dampr,dt1,tt,damp_a)
236
237 omega_(1) = one/ (one + half * damp_a(1) * dt1)
238 omega_(2) = one/ (one + half * damp_a(2) * dt1)
239 omega_(3) = one/ (one + half * damp_a(3) * dt1)
240
241
242 DO n=1,igrnod(igr)%NENTITY
243 i=igrnod(igr)%ENTITY(n)
244 IF(tagslv_rby(i)==0) THEN
245 da_(1) = a(1,i)-damp_a(1)*(v(1,i)-v_ref(1,nd))
246 . +half*damp_a(1)*dt1*a_ref(1,nd)
247 da_(1) = da_(1) * omega_(1) - a(1,i)
248 damp(1,i) = a(1,i)
249 a(1,i) = a(1,i) + da_(1)
250 dw =dw+ms(i)*da_(1)*(v(1,i)+half*a(1,i)*dt1)*dt12*weight(i)
251 ENDIF
252 ENDDO
253
254
255 DO n=1,igrnod(igr)%NENTITY
256 i=igrnod(igr)%ENTITY(n)
257 IF(tagslv_rby(i)==0) THEN
258 da_(2) = a(2,i)-damp_a(2)*(v(2,i)-v_ref(2,nd))
259 . +half*damp_a(2)*dt1*a_ref(2,nd)
260 da_(2) = da_(2) * omega_(2) - a(2,i)
261 damp(2,i) = a(2,i)
262 a(2,i) = a(2,i) + da_(2)
263 dw =dw+ms(i)*da_(2)*(v(2,i)+half*a(2,i)*dt1)*dt12*weight(i)
264 ENDIF
265 ENDDO
266
267 DO n=1,igrnod(igr)%NENTITY
268 i=igrnod(igr)%ENTITY(n)
269 IF(tagslv_rby(i)==0) THEN
270 da_(3) = a(3,i)-damp_a(3)*(v(3,i)-v_ref(3,nd))
271 . +half*damp_a(3)*dt1*a_ref(3,nd)
272 da_(3) = da_(3) * omega_(3) - a(3,i)
273 damp(3,i) = a(3,i)
274 a(3,i) = a(3,i) + da_(3)
275 dw =dw+ms(i)*da_(3)*(v(3,i)+half*a(3,i)*dt1)*dt12*weight(i)
276 ENDIF
277 ENDDO
278
279 ENDIF
280
281 ELSE
282
283
284
285
286 IF ((fl_vrel == 0).AND.(fl_freq_range == 0)) THEN
287
288 IF (dampr(19,nd)==0) THEN
289
290#include "vectorize.inc"
291 DO n=1,igrnod(igr)%NENTITY
292 i=igrnod(igr)%ENTITY(n)
293 vskw(1,i)= skew(1,isk)*v(1,i)+skew(2,isk)*v(2,i)+skew(3,isk)*v(3,i)
294 vskw(2,i)= skew(4,isk)*v(1,i)+skew(5,isk)*v(2,i)+skew(6,isk)*v(3,i)
295 vskw(3,i)= skew(7,isk)*v(1,i)+skew(8,isk)*v(2,i)+skew(9,isk)*v(3,i)
296 askw(1,i)= skew(1,isk)*a(1,i)+skew(2,isk)*a(2,i)+skew(3,isk)*a(3,i)
297 askw(2,i)= skew(4,isk)*a(1,i)+skew(5,isk)*a(2,i)+skew(6,isk)*a(3,i)
298 askw(3,i)= skew(7,isk)*a
299 dampskw(1,i)= skew(1,isk)*damp(1,i)+skew(2,isk)*damp(2,i)+skew(3,isk)*damp(3,i)
300 dampskw(2,i)= skew(4,isk)*damp(1,i)+skew(5,isk)*damp(2,i)+skew(6,isk)*damp(3,i)
301 dampskw(3,i)= skew(7,isk)*damp(1,i)+skew(8,isk)*damp(2,i)+skew(9,isk)*damp(3,i)
302 END DO
303 damp_a(1) = dampr(3,nd)
304 damp_b(1) = dampr(4,nd)
305 betasdt_(1)= -
min(damp_b(1),dampt)*dt1/
max(dt1*dt1,em30)
306 omega_(1) = one/ (one + half * damp_a(1) * dt1)
307
308 damp_a(2) = dampr(5,nd)
309 damp_b(2) = dampr(6,nd)
310 betasdt_(2)= -
min(damp_b(2),dampt)*dt1/
max(dt1*dt1,em30)
311 omega_(2) = one/ (one + half * damp_a(2) * dt1)
312
313 damp_a(3) = dampr(7,nd)
314 damp_b(3) = dampr(8,nd)
315 betasdt_(3)= -
min(damp_b(3),dampt)*dt1/
max(dt1*dt1,em30)
316 omega_(3) = one/ (one + half * damp_a(3) * dt1)
317
318 DO n=1,igrnod(igr)%NENTITY
319 i=igrnod(igr)%ENTITY(n)
320 IF(tagslv_rby(i)==0) THEN
321 da_(1) = askw(1,i) - damp_a(1)*vskw(1,i) - betasdt_(1) *(askw(1,i) - dampskw(1,i))
322 da_(1) = da_(1) * omega_(1) - askw(1,i)
323 dampskw(1,i) = askw(1,i)
324 askw(1,i) = askw(1,i) + da_(1)
325 dw_local(1)= da_(1)*(vskw(1,i)+half*askw(1,i)*dt1)
326
327 da_(2) = askw(2,i) - damp_a(2)*vskw(2,i)- betasdt_(2) *(askw(2,i) - dampskw(2,i))
328 da_(2) = da_(2) * omega_(2) - askw(2,i)
329 dampskw(2,i) = askw(2,i)
330 askw(2,i) = askw(2,i) + da_(2)
331 dw_local(2)= da_(2)*(vskw(2,i)+half*askw(2,i)*dt1)
332
333 da_(3) = askw(3,i) - damp_a(3)*vskw(3,i) - betasdt_(3) *(askw(3,i) - dampskw(3,i))
334 da_(3) = da_(3) * omega_(3) - askw(3,i)
335 dampskw(3,i) = askw(3,i)
336 askw(3,i) = askw(3,i) + da_(3)
337 dw_local(3)= da_(3)*(vskw(3,i)+half*askw(3,i)*dt1)
338
339 dw =dw +ms(i)*(dw_local(1)+dw_local(2)+dw_local(3))*dt12*weight(i)
340 ENDIF
341 ENDDO
342
343#include "vectorize.inc"
344 DO n=1,igrnod(igr)%NENTITY
345 i=igrnod(igr)%ENTITY(n)
346 a(1,i)= skew(1,isk)*askw(1,i)
347 . +skew(4,isk)*askw(2,i)
348 . +skew(7,isk)*askw(3,i)
349 a(2,i)= skew(2,isk)*askw(1,i)
350 . +skew(5,isk)*askw(2,i)
351 . +skew(8,isk)*askw(3,i)
352 a(3,i)= skew(3,isk)*askw(1,i)
353 . +skew(6,isk)*askw(2,i)
354 . +skew(9,isk)*askw(3,i)
355 damp(1,i)= skew(1,isk)*dampskw(1,i)
356 . +skew(4,isk)*dampskw(2,i)
357 . +skew(7,isk)*dampskw(3,i)
358 damp(2,i)= skew(2,isk)*dampskw(1,i)
359 . +skew(5,isk)*dampskw(2,i)
360 . +skew(8,isk)*dampskw(3,i)
361 damp(3,i)= skew(3,isk)*dampskw(1,i)
362 . +skew(6,isk)*dampskw(2,i)
363 . +skew(9,isk)*dampskw(3,i)
364 END DO
365 ENDIF
366
367 IF(iroddl/=0)THEN
368#include "vectorize.inc"
369 DO n=1,igrnod(igr)%NENTITY
370 i=igrnod(igr)%ENTITY(n)
371 vskw(1,i)= skew(1,isk)*vr(1,i)
372 . +skew(2,isk)*vr(2,i)
373 . +skew(3,isk)*vr(3,i)
374 vskw(2,i)= skew(4,isk)*vr(1,i)
375 . +skew(5,isk)*vr(2,i)
376 . +skew(6,isk)*vr(3,i)
377 vskw(3,i)= skew(7,isk)*vr(1,i)
378 . +skew(8,isk)*vr(2,i)
379 . +skew(9,isk)*vr(3,i)
380 askw(1,i)= skew(1,isk)*ar(1,i)
381 . +skew(2,isk)*ar(2,i)
382 . +skew(3,isk)*ar(3,i)
383 askw(2,i)= skew(4,isk)*ar(1,i)
384 . +skew(5,isk)*ar(2,i)
385 . +skew(6,isk)*ar(3,i)
386 askw(3,i)= skew(7,isk)*ar(1,i)
387 . +skew(8,isk)*ar(2,i)
388 . +skew(9,isk)*ar(3,i)
389 dampskw(1,i)= skew(1,isk)*damp(4,i)
390 . +skew(2,isk)*damp(5,i)
391 . +skew(3,isk)*damp(6,i)
392 dampskw(2,i)= skew(4,isk)*damp(4,i)
393 . +skew(5,isk)*damp(5,i)
394 . +skew(6,isk)*damp(6,i)
395 dampskw(3,i)= skew(7,isk)*damp(4,i)
396 . +skew(8,isk)*damp(5,i)
397 . +skew(9,isk)*damp(6,i)
398 END DO
399
400 damp_a(1) = dampr(9,nd)
401 damp_b(1) = dampr(10,nd)
402 betasdt_(1)= -
min(damp_b(1),dampt)*dt1/
max(dt1*dt1,em30)
403 omega_(1) = one/ (one + half * damp_a(1) * dt1)
404
405 damp_a(2) = dampr(11,nd)
406 damp_b(2) = dampr(12,nd)
407 betasdt_(2)= -
min(damp_b(2),dampt)*dt1/
max(dt1*dt1,em30)
408 omega_(2) = one/ (one + half * damp_a(2) * dt1)
409
410 damp_a(3) = dampr(13,nd)
411 damp_b(3) = dampr(14,nd)
412 betasdt_(3)= -
min(damp_b(3),dampt)*dt1/
max(dt1*dt1,em30)
413 omega_(3) = one/ (one + half * damp_a(3) * dt1)
414
415 DO n=1,igrnod(igr)%NENTITY
416 i=igrnod(igr)%ENTITY(n)
417 IF(tagslv_rby(i)==0) THEN
418
419 need_compute = .true.
420 IF (dampr(19,nd)/=0) THEN
421 IF (icontact(i)/=0) i_damp_rdof_tab(i) = dampr(19,nd)
422 IF (i_damp_rdof_tab(i)==0) need_compute=.false.
423 ENDIF
424
425 IF(need_compute) THEN
426 da_(1) = askw(1,i) - damp_a(1)*vskw(1,i)
427 . - betasdt_(1) *(askw(1,i) - dampskw(1,i))
428 da_(1) = da_(1) * omega_(1) - askw(1,i)
429 dampskw(1,i) = askw(1,i)
430 askw(1,i) = askw(1,i) + da_(1)
431 dw_local(1)=da_(1)*(vskw(1,i)+half*askw(1,i)*dt1)
432
433 da_(2) = askw(2,i) - damp_a(2)*vskw(2,i)
434 . - betasdt_(2) *(askw(2,i) - dampskw(2,i))
435 da_(2) = da_(2) * omega_(2) - askw(2,i)
436 dampskw(2,i) = askw(2,i)
437 askw(2,i) = askw(2,i) + da_(2)
438 dw_local(2) =da_(2)*(vskw(2,i)+half*askw(2,i)*dt1)
439
440 da_(3) = askw(3,i) - damp_a(3)*vskw(3,i)
441 . - betasdt_(3) *(askw(3,i) - dampskw(3,i))
442 da_(3) = da_(3) * omega_(3) - askw(3,i)
443 dampskw(3,i) = askw(3,i)
444 askw(3,i) = askw(3,i) + da_(3)
445 dw_local(3)= da_(3)*(vskw(3,i)+half*askw(3,i)*dt1)
446
447 dw =dw +in(i)*(dw_local(1)+dw_local(2)+dw_local(3))*dt12*weight(i)
448 IF (dampr(19,nd)/=0) i_damp_rdof_tab(i)=i_damp_rdof_tab(i)-1
449 ENDIF
450 ENDIF
451 ENDDO
452 ENDIF
453#include "vectorize.inc"
454 DO n=1,igrnod(igr)%NENTITY
455 i=igrnod(igr)%ENTITY(n)
456 ar(1,i)= skew(1,isk)*askw(1,i)
457 . +skew(4,isk)*askw(2,i)
458 . +skew(7,isk)*askw(3,i)
459 ar(2,i)= skew(2,isk)*askw(1,i)
460 . +skew(5,isk)*askw(2,i)
461 . +skew(8,isk)*askw(3,i)
462 ar(3,i)= skew(3,isk)*askw(1,i)
463 . +skew(6,isk)*askw(2,i)
464 . +skew(9,isk)*askw(3,i)
465 damp(4,i)= skew(1,isk)*dampskw(1,i)
466 . +skew(4,isk)*dampskw(2,i)
467 . +skew(7,isk)*dampskw(3,i)
468 damp(5,i)= skew(2,isk)*dampskw(1,i)
469 . +skew(5,isk)*dampskw(2,i)
470 . +skew(8,isk)*dampskw(3,i)
471 damp(6,i)= skew(3,isk)*dampskw(1,i)
472 . +skew(6,isk)*dampskw(2,i)
473 . +skew(9,isk)*dampskw(3,i)
474 END DO
475
476 ELSEIF ((fl_vrel == 1).AND.(id_rby == 0)) THEN
477
478
479#include "vectorize.inc"
480 DO n=1,igrnod(igr)%NENTITY
481 i=igrnod(igr)%ENTITY(n)
482 vskw(1,i)= skew(1,isk)*v(1,i)+skew(2,isk)*v(2,i)+skew(3,isk)*v(3,i)
483 vskw(2,i)= skew(4,isk)*v(1,i)+skew(5,isk)*v(2,i)+skew(6,isk)*v(3,i)
484 vskw(3,i)= skew(7,isk)*v(1,i)+skew(8,isk)*v(2,i)+skew(9,isk)*v(3,i)
485 askw(1,i)= skew(1,isk)*a(1,i)+skew(2,isk)*a(2,i)+skew(3,isk)*a(3,i)
486 askw(2,i)= skew(4,isk)*a(1,i)+skew(5,isk)*a(2,i)+skew(6,isk)*a(3,i)
487 askw(3,i)= skew(7,isk)*a(1,i)+skew(8,isk)*a(2,i)+skew(9,isk)*a(3,i)
488 END DO
489
490
491
492 v_refskw(1)= skew(1,isk)*v_ref(1,nd)+skew(2,isk)*v_ref(2,nd)
493 . +skew(3,isk)*v_ref(3,nd)
494 v_refskw(2)= skew(4,isk)*v_ref(1,nd)+skew(5,isk)*v_ref(2,nd)
495 . +skew(6,isk)*v_ref(3,nd)
496 v_refskw(3)= skew(7,isk)*v_ref(1,nd)+skew(8,isk)*v_ref(2,nd)
497 . +skew(9,isk)*v_ref(3,nd)
498 a_refskw(1)= skew(1,isk)*a_ref(1,nd)+skew(2,isk)*a_ref(2,nd)
499 . +skew(3,isk)*a_ref(3,nd)
500 a_refskw(2)= skew(4,isk)*a_ref(1,nd)+skew(5,isk)*a_ref(2,nd)
501 . +skew(6,isk)*a_ref(3,nd)
502 a_refskw(3)= skew(7,isk)*a_ref(1,nd)+skew(8,isk)*a_ref(2,nd)
503 . +skew(9,isk)*a_ref(3,nd)
504
505
506 call damping_vref_compute_dampa(nd,ndamp,nrdamp,dampr,dt1,tt,damp_a)
507
508 omega_(1) = one/ (one + half * damp_a(1) * dt1)
509 omega_(2) = one/ (one + half * damp_a(2) * dt1)
510 omega_(3) = one/ (one + half * damp_a(3) * dt1)
511
512 DO n=1,igrnod(igr)%NENTITY
513 i=igrnod(igr)%ENTITY(n)
514 IF(tagslv_rby(i)==0) THEN
515 da_(1) = askw(1,i)-damp_a(1)*(vskw(1,i)-v_refskw(1))
516 . +half*damp_a(1)*dt1*a_refskw(1)
517 da_(1) = da_(1) * omega_(1) - askw(1,i)
518 askw(1,i) = askw(1,i) + da_(1)
519 dw_local(1)=da_(1)*(vskw(1,i)+half*askw(1,i)*dt1)
520
521 da_(2) = askw(2,i)-damp_a(2)*(vskw(2,i)-v_refskw(2))
522 . +half*damp_a(2)*dt1*a_refskw(2)
523 da_(2) = da_(2) * omega_(2) - askw(2,i)
524 askw(2,i) = askw(2,i) + da_(2)
525 dw_local(2)=da_(2)*(vskw(2,i)+half*askw(2,i)*dt1)
526
527 da_(3) = askw(3,i)-damp_a(3)*(vskw(3,i)-v_refskw(3))
528 . +half*damp_a(3)*dt1*a_refskw(3)
529 da_(3) = da_(3) * omega_(3) - askw(3,i)
530 askw(3,i) = askw(3,i) + da_(3)
531 dw_local(3)=da_(3)*(vskw(3,i)+half*askw(3,i)*dt1)
532
533 dw =dw +ms(i)*(dw_local(1)+dw_local(2)+dw_local(3))*dt12*weight(i)
534 ENDIF
535 ENDDO
536
537#include "vectorize.inc"
538 DO n=1,igrnod(igr)%NENTITY
539 i=igrnod(igr)%ENTITY(n)
540 a(1,i)= skew(1,isk)*askw(1,i)+skew(4,isk)*askw(2,i)
541 . +skew(7,isk)*askw(3,i)
542 a(2,i)= skew(2,isk)*askw(1,i)+skew(5,isk)*askw(2,i)
543 . +skew(8,isk)*askw(3,i)
544 a(3,i)= skew(3,isk)*askw(1,i)+skew(6,isk)*askw(2,i)
545 . +skew(9,isk)*askw(3,i)
546 END DO
547
548 ENDIF
549
550 END IF
551
552 ENDIF
553 ENDDO
554
555#include "lockon.inc"
556 wfext = wfext + dw
557#include "lockoff.inc"
558 RETURN
559
560
561
562 1000 CONTINUE
563 dw = zero
564 DO nd=1,ndamp
565 igr = nint(dampr(2,nd))
566 isk = nint(dampr(15,nd))
567 factb = dampr(16,nd)
568 dampt =
min(dt1,dt2)*factb
569 d_tstart = dampr(17,nd)
570 d_tstop = dampr(18,nd)
571 fl_freq_range = nint(dampr(31,nd))
572 IF ((tt>=d_tstart .AND. tt<=d_tstop).AND.(fl_freq_range==0)) THEN
573 IF(isk<=1)THEN
574
575 IF(iroddl/=0)THEN
576 dampa = dampr(9,nd)
577 dampb = dampr(10,nd)
578 betasdt= -
min(dampb,dampt)*dt1/
max(dt1*dt1,em30)
579 omega = one/ (one + half * dampa * dt1)
580
581 DO n=1,igrnod(igr)%NENTITY
582 i=igrnod(igr)%ENTITY(n)
583 IF(tagslv_rby(i)/=0) cycle
584
585 IF (dampr(19,nd)/=0) THEN
586 IF (icontact(i)/=0) i_damp_rdof_tab(i) = dampr(19,nd)
587 IF (i_damp_rdof_tab(i)==0) GOTO 351
588 ENDIF
589
590 da = ar(1,i) - dampa*vr(1,i)
591 . - betasdt *(ar(1,i)-damp(4,i))
592 da = da * omega - ar(1,i)
593 damp(4,i) = ar(1,i)
594 ar(1,i) = ar(1,i) + da
595 dw = dw+in(i)*da*(vr(1,i)+half*ar(1,i)*dt1)*dt12*weight(i)
596351 CONTINUE
597 ENDDO
598 dampa = dampr(11,nd)
599 dampb = dampr(12,nd)
600 betasdt= -
min(dampb,dampt)*dt1/
max(dt1*dt1,em30)
601 omega = one/ (one + half * dampa * dt1)
602
603 DO n=1,igrnod(igr)%NENTITY
604 i=igrnod(igr)%ENTITY(n)
605 IF(tagslv_rby(i)/=0) cycle
606
607 IF (dampr(19,nd)/=0) THEN
608 IF (i_damp_rdof_tab(i)==0) GOTO 451
609 ENDIF
610
611 da = ar(2,i) - dampa*vr(2,i)
612 . - betasdt *(ar(2,i)-damp(5,i))
613 da = da * omega - ar(2,i)
614 damp(5,i) = ar(2,i)
615 ar(2,i) = ar(2,i) + da
616 dw = dw+in(i)*da*(vr(2,i)+half*ar(2,i)*dt1)*dt12*weight(i)
617451 CONTINUE
618 ENDDO
619 dampa = dampr(13,nd)
620 dampb = dampr(14,nd)
621 betasdt= -
min(dampb,dampt)*dt1/
max(dt1*dt1,em30)
622 omega = one/ (one + half * dampa * dt1)
623
624 DO n=1,igrnod(igr)%NENTITY
625 i=igrnod(igr)%ENTITY(n)
626 IF(tagslv_rby(i)/=0) cycle
627
628 IF (dampr(19,nd)/=0) THEN
629 IF (i_damp_rdof_tab(i)==0) GOTO 551
630 i_damp_rdof_tab(i)=i_damp_rdof_tab(i)-1
631 ENDIF
632
633 da = ar(3,i) - dampa*vr(3,i)
634 . - betasdt *(ar(3,i)-damp(6,i))
635 da = da * omega - ar(3,i)
636 damp(6,i) = ar(3,i)
637 ar(3,i) = ar(3,i) + da
638 dw = dw+in(i)*da*(vr(3,i)+half*ar(3,i)*dt1)*dt12*weight(i)
639551 CONTINUE
640 ENDDO
641 END IF
642 ELSE
643 IF(iroddl/=0)THEN
644#include "vectorize.inc"
645 DO n=1,igrnod(igr)%NENTITY
646 i=igrnod(igr)%ENTITY(n)
647 vskw(1,i)= skew(1,isk)*vr(1,i)
648 . +skew(2,isk)*vr(2,i)
649 . +skew(3,isk)*vr(3,i)
650 vskw(2,i)= skew(4,isk)*vr(1,i)
651 . +skew(5,isk)*vr(2,i)
652 . +skew(6,isk)*vr(3,i)
653 vskw(3,i)= skew(7,isk)*vr(1,i)
654 . +skew(8,isk)*vr(2,i)
655 . +skew(9,isk)*vr(3,i)
656 askw(1,i)= skew(1,isk)*ar(1,i)
657 . +skew(2,isk)*ar(2,i)
658 . +skew(3,isk)*ar(3,i)
659 askw(2,i)= skew(4,isk)*ar(1,i)
660 . +skew(5,isk)*ar(2,i)
661 . +skew(6,isk)*ar(3,i)
662 askw(3,i)= skew(7,isk)*ar(1,i)
663 . +skew(8,isk)*ar(2,i)
664 . +skew(9,isk)*ar(3,i)
665 dampskw(1,i)= skew(1,isk)*damp(4,i)
666 . +skew(2,isk)*damp(5,i)
667 . +skew(3,isk)*damp(6,i)
668 dampskw(2,i)= skew(4,isk)*damp(4,i)
669 . +skew(5,isk)*damp(5,i)
670 . +skew(6,isk)*damp(6,i)
671 dampskw(3,i)= skew(7,isk)*damp(4,i)
672 . +skew(8,isk)*damp(5,i)
673 . +skew(9,isk)*damp(6,i)
674 END DO
675 dampa = dampr(9,nd)
676 dampb = dampr(10,nd)
677 betasdt= -
min(dampb,dampt)*dt1/
max(dt1*dt1,em30)
678 omega = one/ (one + half * dampa * dt1)
679
680 DO n=1,igrnod(igr)%NENTITY
681 i=igrnod(igr)%ENTITY(n)
682 IF(tagslv_rby(i)/=0) cycle
683
684 IF (dampr(19,nd)/=0) THEN
685 IF (icontact(i)/=0) i_damp_rdof_tab(i) = dampr(19,nd)
686 IF (i_damp_rdof_tab(i)==0) GOTO 751
687 ENDIF
688
689 da = askw(1,i) - dampa*vskw(1,i)
690 . - betasdt *(askw(1,i) - dampskw(1,i))
691 da = da * omega - askw(1,i)
692 dampskw(1,i) = askw(1,i)
693 askw(1,i) = askw(1,i) + da
694 dw =dw
695 . +in(i)*da*(vskw(1,i)+half*askw(1,i)*dt1)*dt12*weight(i)
696751 CONTINUE
697 ENDDO
698 dampa = dampr(11,nd)
699 dampb = dampr(12,nd)
700 betasdt= -
min(dampb,dampt)*dt1/
max(dt1*dt1,em30)
701 omega = one/ (one + half * dampa * dt1)
702
703 DO n=1,igrnod(igr)%NENTITY
704 i=igrnod(igr)%ENTITY(n)
705 IF(tagslv_rby(i)/=0) cycle
706
707 IF (dampr(19,nd)/=0) THEN
708 IF (i_damp_rdof_tab(i)==0) GOTO 851
709 ENDIF
710
711 da = askw(2,i) - dampa*vskw(2,i)
712 . - betasdt *(askw(2,i) - dampskw(2,i))
713 da = da * omega - askw(2,i)
714 dampskw(2,i) = askw(2,i)
715 askw(2,i) = askw(2,i) + da
716 dw =dw
717 . +in(i)*da*(vskw(2,i)+half*askw(2,i)*dt1)*dt12*weight(i)
718851 CONTINUE
719 ENDDO
720 dampa = dampr(13,nd)
721 dampb = dampr(14,nd)
722 betasdt= -
min(dampb,dampt)*dt1/
max(dt1*dt1,em30)
723 omega = one/ (one + half * dampa * dt1)
724
725 DO n=1,igrnod(igr)%NENTITY
726 i=igrnod(igr)%ENTITY(n)
727 IF(tagslv_rby(i)/=0) cycle
728
729 IF (dampr(19,nd)/=0) THEN
730 IF (i_damp_rdof_tab(i)==0) GOTO 951
731 i_damp_rdof_tab(i)=i_damp_rdof_tab(i)-1
732 ENDIF
733
734 da = askw(3,i) - dampa*vskw(3,i)
735 . - betasdt *(askw(3,i) - dampskw(3,i))
736 da = da * omega - askw(3,i)
737 dampskw(3,i) = askw(3,i)
738 askw(3,i) = askw(3,i) + da
739 dw =dw
740 . +in(i)*da*(vskw(3,i)+half*askw(3,i)*dt1)*dt12*weight(i)
741951 CONTINUE
742 ENDDO
743#include "vectorize.inc"
744 DO n=1,igrnod(igr)%NENTITY
745 i=igrnod(igr)%ENTITY(n)
746 IF(tagslv_rby(i)/=0) cycle
747 ar(1,i)= skew(1,isk)*askw(1,i)
748 . +skew(4,isk
749 . +skew(7,isk)*askw(3,i)
750 ar(2,i)= skew(2,isk)*askw(1,i)
751 . +skew(5,isk)*askw(2,i)
752 . +skew(8,isk)*askw(3,i)
753 ar(3,i)= skew(3,isk)*askw(1,i)
754 . +skew(6,isk)*askw(2,i)
755 . +skew(9,isk)*askw(3,i)
756 damp(4,i)= skew(1,isk)*dampskw(1,i)
757 . +skew(4,isk)*dampskw(2,i)
758 . +skew(7,isk)*dampskw(3,i)
759 damp(5,i)= skew(2,isk)*dampskw(1,i)
760 . +skew(5,isk)*dampskw(2,i)
761 . +skew(8,isk)*dampskw(3,i)
762 damp(6,i)= skew(3,isk)*dampskw(1,i)
763 . +skew(6,isk)*dampskw(2,i)
764 . +skew(9,isk)*dampskw(3,i)
765 END DO
766 END IF
767 END IF
768 ENDIF
769 ENDDO
770
771#include "lockon.inc"
772 wfext = wfext + dw
773#include "lockoff.inc"
774 DEALLOCATE(vskw)
775 DEALLOCATE(askw)
776 DEALLOCATE(dampskw)
777
778
779 RETURN
subroutine damping_vref(igrnod, v, a, ms, dampr, weight, v_ref, a_ref, ndamp_vrel, id_damp_vrel, fr_damp_vrel, ngrnod, numnod, ndamp, nrdamp, nspmd, iparit, ispmd)