41
42
43
44 USE sensor_mod
45
46
47
48
49
50
51
52
53
54
55
56#include "implicit_f.inc"
57
58
59
60#include "mvsiz_p.inc"
61#include "param_c.inc"
62#include "com08_c.inc"
63#include "com04_c.inc"
64#include "com01_c.inc"
65#include "scr05_c.inc"
66
67
68
69 INTEGER ,INTENT(IN) :: NSENSOR
70 INTEGER JFT, JLT, IOUT, NUVAR, IPROP, IXR(NIXR,*),IGTYP,
71 . ISENS,NC1(*),NC2(*)
72
74 . rot1(3,mvsiz),rot2(3,mvsiz),rby(*),
75 . dx(*), dy(*), dz(*), rx(*), ry(*), rz(*),vr(3,*)
76 DOUBLE PRECISION XDP(3,*),XL(MVSIZ,3)
77 TYPE () ,DIMENSION(NSENSOR) :: SENSOR_TAB
78
79
80
81 INTEGER I, J, K, JTYP, IERR, NRB, SKFLG,
82 . IDSK1,IDSK2,ISK1,ISK2,IP1,IP2,
83 . N1,N2,N3,USENS,INSENS,
84 . ISENS_OLD,ISENS_ACT
86 . u(lskew),v(lskew),u0(lskew),v0(lskew),ex(lskew),
87 . r1(3),r2(3),rm(3),rl1(3),rl2(3),t(3),
88 . x1(lskew),x2(lskew),q(lskew),a0(lskew),b0(lskew),
89 . a(lskew),b(lskew),exi(lskew),
90 . get_u_geo,nr,dt(3),dex(lskew),exprec(lskew)
91 DOUBLE PRECISION X21(MVSIZ),Y21(MVSIZ),Z21(MVSIZ)
92
93
94
95 INTEGER GET_U_SKEW
97
98 DO i=jft,jlt
99 nc1(i)= ixr(2,i)
100 nc2(i)= ixr(3,i)
101 ENDDO
102
103 ierr = 0
104 isens = 0
105 isens_act = 0
106 jtyp = nint(get_u_geo(1,iprop))
107 IF (igtyp==33) THEN
108
109 idsk1 = nint(get_u_geo(2,iprop))
110 idsk2 = nint(get_u_geo(3,iprop))
111 skflg = nint(get_u_geo(14,iprop))
114 ELSE
115 skflg = 0
116
117 usens = nint(get_u_geo(2,iprop))
118 isens_act = 0
119 IF (usens>0) THEN
120 isens_old = nint(uvar(16,jft))
121 DO k=1,nsensor
122 IF(usens==sensor_tab(k)%SENS_ID) insens=k
123 ENDDO
124
125 IF (tt>sensor_tab(insens)%TSTART) THEN
126 isens = 1
127 uvar(16,jft) = isens
128 ENDIF
129 IF (isens/=isens_old) isens_act = 1
130 ENDIF
131 ENDIF
132
133
134
135 DO i=jft,jlt
136
137
138
139 IF ((ncycle==0).AND.(tt==0)) THEN
140 IF (igtyp==33) THEN
141 DO j=1,9
142 a0(j)= uvar(3+j,i)
143 END DO
144 ELSE
145 rx(i)= uvar(7,i)
146 ry(i)= uvar(8,i)
147 rz(i)= uvar(9,i)
148 ENDIF
149 ENDIF
150
151
152
153 IF ((igtyp==45).AND.(isens_act==1)) THEN
154 uvar(4,i) = dx(i)
155 uvar(5,i) = dy(i)
156 uvar(6,i) = dz(i)
157 uvar(7,i) = rx(i)
158 uvar(8,i) = ry(i)
159 uvar(9,i) = rz(i)
160 ENDIF
161
162
163
164 IF (jtyp==5) THEN
165
166
167
168 IF (igtyp==45) THEN
169
170 dt(1)= vr(1,nc1(i))*dt1
171 dt(2)= vr(2,nc1(i))*dt1
172 dt(3)= vr(3,nc1(i))*dt1
173 u(1)=uvar(10,i) - uvar(11,i)*dt(3)+uvar(12,i)*dt(2)
174 u(2)=uvar(11,i) - uvar(12,i)*dt(1)+uvar(10,i)*dt(3)
175 u(3)=uvar(12,i) - uvar(10,i)*dt(2)+uvar(11,i)*dt(1)
176 nr =sqrt(u(1)*u(1)+u(2)*u(2)+u(3)*u(3))
177 IF (nr>0) THEN
178 u(1)=u(1)/nr
179 u(2)=u(2)/nr
180 u(3)=u(3)/nr
181 ENDIF
182
183 uvar(10,i) = u(1)
184 uvar(11,i) = u(2)
185 uvar(12,i) = u(3)
186
187
188 dt(1)= vr(1,nc2(i))*dt1
189 dt(2)= vr(2,nc2(i))*dt1
190 dt(3)= vr(3,nc2(i))*dt1
191 v(1)=uvar(13,i) - uvar(14,i)*dt(3)+uvar(15,i)*dt(2)
192 v(2)=uvar(14,i) - uvar(15,i)*dt(1)+uvar(13,i)*dt(3)
193 v(3)=uvar(15,i) - uvar(13,i)*dt(2)+uvar(14,i)*dt(1)
194 nr =sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
195 IF (nr>0) THEN
196 v(1)=v(1)/nr
197 v(2)=v(2)/nr
198 v(3)=v(3)/nr
199 ENDIF
200
201 uvar(13,i) = v(1)
202 uvar(14,i) = v(2)
203 uvar(15,i) = v(3)
204 ENDIF
205
206 ex(1) = u(2)*v(3) - u(3)*v(2)
207 ex(2) = u(3)*v(1) - u(1)*v(3)
208 ex(3) = u(1)*v(2) - u(2)*v(1)
209 nx = sqrt(ex(1)*ex(1)+ex(2)*ex(2)+ex(3)*ex(3))
210 ex(1) = ex(1) / nx
211 ex(2) = ex(2) / nx
212 ex(3) = ex(3) / nx
213 ex(4) = u(1)
214 ex(5) = u(2)
215 ex(6) = u(3)
216 ex(7) = v(1)
217 ex(8) = v(2)
218 ex(9) = v(3)
219
221
222
223 x21(i) = (vr(1,nc2(i))-vr(1,nc1(i)))*dt1
224 y21(i) = (vr(2,nc2(i))-vr(2,nc1(i)))*dt1
225 z21(i) = (vr(3,nc2(i))-vr(3,nc1(i)))*dt1
226
227 rm(1) = rx(i)+exi(1)*x21(i)+exi(4)*y21(i)+exi(7)*z21(i)
228 rm(2) = ry(i)+exi(2)*x21(i)+exi(5)*y21(i)+exi(8)*z21(i)
229 rm(3) = rz(i)+exi(3)*x21(i)+exi(6)*y21(i)+exi(9)*z21(i)
230
231 r2(1) = 0.5*rm(1)
232 r2(2) = 0.5*rm(2)
233 r2(3) = 0.5*rm(3)
234 r1(1) = -0.5*rm(1)
235 r1(2) = -0.5*rm(2)
236 r1(3) = -0.5*rm(3)
237
238 x21(i) = x(1,nc2(i))-x(1,nc1(i))
239 y21(i) = x(2,nc2(i))-x(2,nc1(i))
240 z21(i) = x(3,nc2(i))-x(3,nc1(i))
241 xl(i,1)=exi(1)*x21(i)+exi(4)*y21(i)+exi(7)*z21(i)
242 xl(i,2)=exi(2)*x21(i)+exi(5)*y21(i)+exi(8)*z21(i)
243 xl(i,3)=exi(3)*x21(i)+exi(6)*y21(i)+exi(9)*z21(i)
244
245
246 ELSE
247
248 IF (skflg==1) THEN
249
250
251
252
253 DO j=1,9
254 ex(j) = u(j)
255 ENDDO
256
257 ELSE
258
259
260
261
262
263 IF ((ncycle==0).AND.(tt==0).AND.(igtyp==33)) THEN
265 DO j=1,9
266 exprec(j) = a(j)
267 END DO
268 ELSE
269 DO j=1,9
270 exprec(j) = uvar(21+j,i)
271 END DO
272 ENDIF
273
274
275 x21(i) = half*(vr(1,nc2(i))+vr(1,nc1(i)))*dt1
276 y21(i) = half*(vr(2,nc2(i))+vr(2,nc1(i)))*dt1
277 z21(i) = half*(vr(3,nc2(i))+vr(3,nc1(i)))*dt1
278 dt(1)=exprec(1)*x21(i)+exprec(2)*y21(i)+exprec(3)*z21(i)
279 dt(2)=exprec(4)*x21(i)+exprec(5)*y21(i)+exprec(6)*z21(i)
280 dt(3)=exprec(7)*x21(i)+exprec(8)*y21(i)+exprec(9)*z21(i)
281
282 nr =sqrt(dt(1)*dt(1)+dt(2)*dt(2)+dt(3)*dt(3))
283 IF (nr>0) THEN
284 dt(1)=dt(1)/nr
285 dt(2)=dt(2)/nr
286 dt(3)=dt(3)/nr
287 ENDIF
288 co = cos(nr)
289 si = sin(nr)
290 CALL qrot33(dex, dt, co, si)
291
292
294
295
296 ENDIF
297
298
299 x21(i) = (vr(1,nc2(i))-vr(1,nc1(i)))*dt1
300 y21(i) = (vr(2,nc2(i))-vr(2,nc1(i)))*dt1
301 z21(i) = (vr(3,nc2(i))-vr(3,nc1(i)))*dt1
302
303 rm(1) = rx(i)+ex(1)*x21(i)+ex(2)*y21(i)+ex(3)*z21(i)
304 rm(2) = ry(i)+ex(4)*x21(i)+ex(5)*y21(i)+ex(6)*z21(i)
305 rm(3) = rz(i)+ex(7)*x21(i)+ex(8)*y21(i)+ex(9)*z21(i)
306
307 r2(1) = 0.5*rm(1)
308 r2(2) = 0.5*rm(2)
309 r2(3) = 0.5*rm(3)
310 r1(1) = -0.5*rm(1)
311 r1(2) = -0.5*rm(2)
312 r1(3) = -0.5*rm(3)
313
314 IF (iresp == 1) THEN
315
316 x21(i) = xdp(1,nc2(i))-xdp(1,nc1(i))
317 y21(i) = xdp(2,nc2(i))-xdp(2,nc1(i))
318 z21(i) = xdp(3,nc2(i))-xdp(3,nc1(i))
319 ELSE
320
321 x21(i) = x(1,nc2(i))-x(1,nc1(i))
322 y21(i) = x(2,nc2(i))-x(2,nc1(i))
323 z21(i) = x(3,nc2(i))-x(3,nc1(i))
324 ENDIF
325
326 xl(i,1)=ex(1)*x21(i)+ex(2)*y21(i)+ex(3)*z21(i)
327 xl(i,2)=ex(4)*x21(i)+ex(5)*y21(i)+ex(6)*z21(i)
328 xl(i,3)=ex(7)*x21(i)+ex(8)*y21(i)+ex(9)*z21(i)
329
330
331 ENDIF
332
333
334 rot1(1,i) = r1(1)
335 rot1(2,i) = r1(2)
336 rot1(3,i) = r1(3)
337 rot2(1,i) = r2(1)
338 rot2(2,i) = r2(2)
339 rot2(3,i) = r2(3)
340
341 DO j=1,9
342 uvar(21+j,i) = ex(j)
343 END DO
344 END DO
345
346
347 RETURN
subroutine qrot33(skew, t, c, s)
subroutine prod_ab(a, b, x)
integer function get_u_skew(idskw, n1, n2, n3, v)