OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
emomt2.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| emomt2 ../engine/source/ale/euler2d/emomt2.F
25!||--- called by ------------------------------------------------------
26!|| qforc2 ../engine/source/elements/solid_2d/quad/qforc2.F
27!||--- calls -----------------------------------------------------
28!|| upwind ../engine/source/elements/solid/solide/upwind.F
29!||--- uses -----------------------------------------------------
30!|| ale_mod ../common_source/modules/ale/ale_mod.f
31!||====================================================================
32 SUBROUTINE emomt2(PM,RHO,
33 . Y1,Y2,Y3,Y4,Z1,Z2,Z3,Z4,
34 . VY1,VY2,VY3,VY4,VZ1,VZ2,VZ3,VZ4,
35 . T11, T12, T13, T14, T21, T22, T23 , T24,
36 . PY1, PY2, PZ1, PZ2, AIRE,
37 . DYY, DZZ, DYZ, DZY,VDY,VDZ,
38 . DELTAX,VIS,
39 . MAT,
40 . QMV,BUFMAT,IPM)
41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE ale_mod
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C G l o b a l P a r a m e t e r s
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "com04_c.inc"
57#include "com08_c.inc"
58#include "vect01_c.inc"
59#include "param_c.inc"
60C-----------------------------------------------
61C D u m m y A r g u m e n t s
62C-----------------------------------------------
63 my_real pm(npropm,nummat), rho(*),
64 . y1(*),y2(*),y3(*),y4(*),z1(*),z2(*),z3(*),z4(*),
65 . vy1(*), vy2(*), vy3(*), vy4(*), vz1(*), vz2(*), vz3(*), vz4(*),
66 . t11(*), t12(*), t13(*), t14(*), t21(*), t22(*), t23(*), t24(*),
67 . py1(*), py2(*), pz1(*), pz2(*), aire(*),
68 . dyy(*), dzz(*), dyz(*), dzy(*),vdy(*),vdz(*),
69 . deltax(*),vis(*),qmv(8,*),bufmat(*)
70 INTEGER,INTENT(IN) :: MAT(*)
71 INTEGER,INTENT(IN)::IPM(NPROPMI,NUMMAT)
72C-----------------------------------------------
73C L o c a l V a r i a b l e s
74C-----------------------------------------------
75 my_real
76 . Y12(MVSIZ),Z12(MVSIZ),Y23(MVSIZ),Z23(MVSIZ),
77 . Y34(MVSIZ),Z34(MVSIZ),Y41(MVSIZ),Z41(MVSIZ),
78 . y24(mvsiz),z24(mvsiz),y13(mvsiz),z13(mvsiz),
79 . dyyp(mvsiz), dzzp(mvsiz), dyzp(mvsiz), dzyp(mvsiz),
80 . f1(mvsiz), f2(mvsiz),a1(mvsiz),a2(mvsiz),a3(mvsiz),gam(mvsiz),
81 . f1p(mvsiz),f2p(mvsiz),
82 . s(3,mvsiz),t(3,mvsiz),
83 . tr(mvsiz),fac,dt1h,
84 . pd,ps,pp,pc
85 INTEGER I
86C-----------------------------------------------
87 my_real AAA,DM,DMY,DMZ
88 INTEGER IADBUF,IFLG
89C-----------------------------------------------
90C S o u r c e L i n e s
91C-----------------------------------------------
92C
93C-----------------------------------------------
94C Reynolds transportation theorem
95C-----------------------------------------------
96C DMi = 0.5*DT1*QMV(i,I) incoming mass from face i
97 IF(mtn == 51 .AND. ale%UPWIND%UPWM /= 3)THEN
98 iadbuf = ipm(27,mat(1))
99 iflg = nint(bufmat(31+iadbuf-1))
100C
101 IF(iflg > 1)RETURN
102C
103 IF(ale%UPWIND%UPWM == 0.)THEN
104 DO i=lft,llt
105 gam(i)= pm(15,mat(i))
106 ENDDO
107 ELSE
108 DO i=lft,llt
109 gam(i)= ale%UPWIND%CUPWM
110 ENDDO
111 ENDIF
112C
113 DO i=lft,llt
114 aaa = qmv(1,i)+qmv(2,i)+qmv(3,i)+qmv(4,i)
115 aaa = fourth * aaa
116 qmv(1,i) = fourth * (qmv(1,i) - aaa)
117 qmv(2,i) = fourth * (qmv(2,i) - aaa)
118 qmv(3,i) = fourth * (qmv(3,i) - aaa)
119 qmv(4,i) = fourth * (qmv(4,i) - aaa)
120 dmy = zero
121 dmz = zero
122 dm = qmv(4,i)+qmv(1,i)
123 dmy = dmy + vy1(i)*dm
124 dmz = dmz + vz1(i)*dm
125 dm = qmv(1,i)+qmv(2,i)
126 dmy = dmy + vy2(i)*dm
127 dmz = dmz + vz2(i)*dm
128 dm = qmv(2,i)+qmv(3,i)
129 dmy = dmy + vy3(i)*dm
130 dmz = dmz + vz3(i)*dm
131 dm = qmv(3,i)+qmv(4,i)
132 dmy = dmy + vy4(i)*dm
133 dmz = dmz + vz4(i)*dm
134C
135 dmy = -fourth * dmy
136 dmz = -fourth * dmz
137C
138 a1(i) = py1(i)*vdy(i)+pz1(i)*vdz(i)
139 a2(i) = py2(i)*vdy(i)+pz2(i)*vdz(i)
140 a1(i) = sign(gam(i),a1(i))
141 a2(i) = sign(gam(i),a2(i))
142C
143 t11(i) = t11(i)+(one+a1(i))*dmy
144 t12(i) = t12(i)+(one+a2(i))*dmy
145 t13(i) = t13(i)+(one-a1(i))*dmy
146 t14(i) = t14(i)+(one-a2(i))*dmy
147C
148 t21(i) = t21(i)+(one+a1(i))*dmz
149 t22(i) = t22(i)+(one+a2(i))*dmz
150 t23(i) = t23(i)+(one-a1(i))*dmz
151 t24(i) = t24(i)+(one-a2(i))*dmz
152 ENDDO
153 RETURN
154 ENDIF
155 IF(ale%UPWIND%UPWM == 0.AND.mtn == 11)RETURN
156 IF(ale%UPWIND%UPWM > 1)THEN
157 pc=fourth
158 pd=one_over_16
159 ps=one_over_16
160 pp=one_over_16
161 dt1h=half*dt1*ale%UPWIND%CUPWM
162 DO i=lft,llt
163 y13(i)=y3(i)-y1(i)
164 z13(i)=z3(i)-z1(i)
165 y24(i)=y4(i)-y2(i)
166 z24(i)=z4(i)-z2(i)
167 y12(i)=y2(i)-y1(i)
168 z12(i)=z2(i)-z1(i)
169 y23(i)=y3(i)-y2(i)
170 z23(i)=z3(i)-z2(i)
171 y34(i)=y4(i)-y3(i)
172 z34(i)=z4(i)-z3(i)
173 y41(i)=y1(i)-y4(i)
174 z41(i)=z1(i)-z4(i)
175 ENDDO
176 DO i=lft,llt
177 s(2,i)=half*(y12(i)-y34(i))
178 s(3,i)=half*(z12(i)-z34(i))
179 t(2,i)=half*(-y41(i)+y23(i))
180 t(3,i)=half*(-z41(i)+z23(i))
181 ENDDO
182 ENDIF
183C
184C FORCE AT CENTROID
185C
186 IF(ale%UPWIND%UPWM < 2.OR.mtn == 11)THEN
187 IF(mtn == 11.AND.ale%UPWIND%UPWM > 1)THEN
188 CALL upwind(
189 1 rho, vis, vdy, vdy,
190 2 vdz, s, s, t,
191 3 deltax, gam, llt)
192 DO i=lft,llt
193 fac=four*gam(i)/aire(i)
194 a1(i) = fac*(py1(i)*vdy(i)+pz1(i)*vdz(i))
195 a2(i) = fac*(py2(i)*vdy(i)+pz2(i)*vdz(i))
196 ENDDO
197 ELSE
198 DO i=lft,llt
199 gam(i)= pm(15,mat(i))
200 a1(i) = py1(i)*vdy(i)+pz1(i)*vdz(i)
201 a2(i) = py2(i)*vdy(i)+pz2(i)*vdz(i)
202 a1(i) = sign(gam(i),a1(i))
203 a2(i) = sign(gam(i),a2(i))
204 ENDDO
205 ENDIF
206 DO i=lft,llt
207 fac = fourth*rho(i)*aire(i)
208 f1(i) = (vdy(i)*dyy(i)+vdz(i)*dyz(i))*fac
209 f2(i) = (vdy(i)*dzy(i)+vdz(i)*dzz(i))*fac
210 ENDDO
211 DO i=lft,llt
212 t11(i) = t11(i)+(one+a1(i))*f1(i)
213 t12(i) = t12(i)+(one+a2(i))*f1(i)
214 t13(i) = t13(i)+(one-a1(i))*f1(i)
215 t14(i) = t14(i)+(one-a2(i))*f1(i)
216 t21(i) = t21(i)+(one+a1(i))*f2(i)
217 t22(i) = t22(i)+(one+a2(i))*f2(i)
218 t23(i) = t23(i)+(one-a1(i))*f2(i)
219 t24(i) = t24(i)+(one-a2(i))*f2(i)
220 ENDDO
221C
222 ELSE
223C
224C TRANSPORTATION FORCE WHICH MINIMIZE HOURGLASS ISSUES
225C EVALUATE <PHI,UJ*DUI/DXJ> DIRECTLY TO NODES
226C DERI PHI1
227C P1Y1=-Z24/TR1 P1Z1= Y24/TR1 P1Y2=-Z23/TR2 P1Z2= Y23/TR2
228C P1Y4=-Z34/TR4 P1Z4= Y34/TR4
229C DERI PHI2
230C P2Y1=-Z41/TR1 P2Z1= Y41/TR1 P2Y2= Z13/TR2 P2Z2=-Y13/TR2
231C P2Y3=-Z34/TR3 P2Z3= Y34/TR3
232C DERI PHI3
233C P3Y2=-Z12/TR2 P3Z2= Y12/TR2
234C P3Y3= Z24/TR3 P3Z3=-Y24/TR3 P3Y4=-Z41/TR4 P3Z4= Y41/TR4
235C DERIVATE DE PHI4
236C P4Y1=-Z12/TR1 P4Z1= Y12/TR1
237C P4Y3=-Z23/TR3 P4Z3= Y23/TR3 P4Y4=-Z13/TR4 P4Z4= Y13/TR4
238C
239C DERIV POINT 1
240C P1Y1=-Z24/TR1 P1Z1= Y24/TR1 P2Y1=-Z41/TR1 P2Z1= Y41/TR1
241C P4Y1=-Z12/TR1 P4Z1= Y12/TR1
242C DERIV POINT 2
243C P1Y2=-Z23/TR2 P1Z2= Y23/TR2 P2Y2= Z13/TR2 P2Z2=-Y13/TR2
244C P3Y2=-Z12/TR2 P3Z2= Y12/TR2
245C DERIV POINT 3
246C P2Y3=-Z34/TR3 P2Z3= Y34/TR3
247C P3Y3= Z24/TR3 P3Z3=-Y24/TR3 P4Y3=-Z23/TR3 P4Z3= Y23/TR3
248C DERIV POINT 4
249C P1Y4=-Z34/TR4 P1Z4= Y34/TR4
250C P3Y4=-Z41/TR4 P3Z4= Y41/TR4 P4Y4=-Z13/TR4 P4Z4= Y13/TR4
251C
252C NODE 1
253 CALL upwind(
254 1 rho, vis, vy1, vy1,
255 2 vz1, s, s, t,
256 3 deltax, gam, llt)
257 DO i=lft,llt
258 tr(i)=(y41(i)*z12(i)-y12(i)*z41(i))
259 dyyp(i)=(-z24(i)*vy1(i)-z41(i)*vy2(i)-z12(i)*vy4(i))
260 dyzp(i)=( y24(i)*vy1(i)+y41(i)*vy2(i)+y12(i)*vy4(i))
261 dzzp(i)=( y24(i)*vz1(i)+y41(i)*vz2(i)+y12(i)*vz4(i))
262 dzyp(i)=(-z24(i)*vz1(i)-z41(i)*vz2(i)-z12(i)*vz4(i))
263 f1p(i)=rho(i)*(vy1(i)*dyyp(i)+vz1(i)*dyzp(i))
264 f2p(i)=rho(i)*(vy1(i)*dzyp(i)+vz1(i)*dzzp(i))
265 fac=pc*gam(i)/tr(i)
266 a1(i)=fac*(-z24(i)*vy1(i)+y24(i)*vz1(i))
267 a2(i)=fac*(-z41(i)*vy1(i)+y41(i)*vz1(i))
268 a3(i)=fac*(-z12(i)*vy1(i)+y12(i)*vz1(i))
269C
270 t11(i) = t11(i)+(pp+a1(i))*f1p(i)
271 t12(i) = t12(i)+(ps+a2(i))*f1p(i)
272 t13(i) = t13(i)+pd*f1p(i)
273 t14(i) = t14(i)+(ps+a3(i))*f1p(i)
274 t21(i) = t21(i)+(pp+a1(i))*f2p(i)
275 t22(i) = t22(i)+(ps+a2(i))*f2p(i)
276 t23(i) = t23(i)+pd*f2p(i)
277 t24(i) = t24(i)+(ps+a3(i))*f2p(i)
278 ENDDO
279C
280 IF(ale%UPWIND%UPWM == 2)THEN
281 DO i=lft,llt
282 fac=-dt1h/tr(i)
283 f1(i)=fac*(-f1p(i)*dzzp(i)+f2p(i)*dyzp(i))
284 f2(i)=fac*( f1p(i)*dzyp(i)-f2p(i)*dyyp(i))
285 t11(i) = t11(i)+pp*f1(i)
286 t12(i) = t12(i)+ps*f1(i)
287 t13(i) = t13(i)+pd*f1(i)
288 t14(i) = t14(i)+ps*f1(i)
289 t21(i) = t21(i)+pp*f2(i)
290 t22(i) = t22(i)+ps*f2(i)
291 t23(i) = t23(i)+pd*f2(i)
292 t24(i) = t24(i)+ps*f2(i)
293 ENDDO
294 ENDIF
295C NODE 2
296 CALL upwind(
297 1 rho, vis, vy2, vy2,
298 2 vz2, s, s, t,
299 3 deltax, gam, llt)
300 DO i=lft,llt
301 tr(i)=(y12(i)*z23(i)-y23(i)*z12(i))
302 dyyp(i)=(-z23(i)*vy1(i)+z13(i)*vy2(i)-z12(i)*vy3(i))
303 dyzp(i)=( y23(i)*vy1(i)-y13(i)*vy2(i)+y12(i)*vy3(i))
304 dzzp(i)=( y23(i)*vz1(i)-y13(i)*vz2(i)+y12(i)*vz3(i))
305 dzyp(i)=(-z23(i)*vz1(i)+z13(i)*vz2(i)-z12(i)*vz3(i))
306 f1p(i)=rho(i)*(vy2(i)*dyyp(i)+vz2(i)*dyzp(i))
307 f2p(i)=rho(i)*(vy2(i)*dzyp(i)+vz2(i)*dzzp(i))
308 fac=pc*gam(i)/tr(i)
309 a1(i)=fac*(-z23(i)*vy2(i)+y23(i)*vz2(i))
310 a2(i)=fac*( z13(i)*vy2(i)-y13(i)*vz2(i))
311 a3(i)=fac*(-z12(i)*vy2(i)+y12(i)*vz2(i))
312 t11(i) = t11(i)+(ps+a1(i))*f1p(i)
313 t12(i) = t12(i)+(pp+a2(i))*f1p(i)
314 t13(i) = t13(i)+(ps+a3(i))*f1p(i)
315 t14(i) = t14(i)+pd*f1p(i)
316 t21(i) = t21(i)+(ps+a1(i))*f2p(i)
317 t22(i) = t22(i)+(pp+a2(i))*f2p(i)
318 t23(i) = t23(i)+(ps+a3(i))*f2p(i)
319 t24(i) = t24(i)+pd*f2p(i)
320 ENDDO
321 IF(ale%UPWIND%UPWM == 2)THEN
322 DO i=lft,llt
323 fac=-dt1h/tr(i)
324 f1(i)=fac*(-f1p(i)*dzzp(i)+f2p(i)*dyzp(i))
325 f2(i)=fac*( f1p(i)*dzyp(i)-f2p(i)*dyyp(i))
326 t11(i) = t11(i)+ps*f1(i)
327 t12(i) = t12(i)+pp*f1(i)
328 t13(i) = t13(i)+ps*f1(i)
329 t14(i) = t14(i)+pd*f1(i)
330 t21(i) = t21(i)+ps*f2(i)
331 t22(i) = t22(i)+pp*f2(i)
332 t23(i) = t23(i)+ps*f2(i)
333 t24(i) = t24(i)+pd*f2(i)
334 ENDDO
335 ENDIF
336C NODE 3
337 CALL upwind(
338 1 rho, vis, vy3, vy3,
339 2 vz3, s, s, t,
340 3 deltax, gam, llt)
341 DO i=lft,llt
342 tr(i)=(y23(i)*z34(i)-y34(i)*z23(i))
343 dyyp(i)=(-z34(i)*vy2(i)+z24(i)*vy3(i)-z23(i)*vy4(i))
344 dyzp(i)=( y34(i)*vy2(i)-y24(i)*vy3(i)+y23(i)*vy4(i))
345 dzzp(i)=( y34(i)*vz2(i)-y24(i)*vz3(i)+y23(i)*vz4(i))
346 dzyp(i)=(-z34(i)*vz2(i)+z24(i)*vz3(i)-z23(i)*vz4(i))
347 f1p(i)=rho(i)*(vy3(i)*dyyp(i)+vz3(i)*dyzp(i))
348 f2p(i)=rho(i)*(vy3(i)*dzyp(i)+vz3(i)*dzzp(i))
349 fac=pc*gam(i)/tr(i)
350 a1(i)=fac*(-z34(i)*vy3(i)+y34(i)*vz3(i))
351 a2(i)=fac*( z24(i)*vy3(i)-y24(i)*vz3(i))
352 a3(i)=fac*(-z23(i)*vy3(i)+y23(i)*vz3(i))
353 t11(i) = t11(i)+pd*f1p(i)
354 t12(i) = t12(i)+(ps+a1(i))*f1p(i)
355 t13(i) = t13(i)+(pp+a2(i))*f1p(i)
356 t14(i) = t14(i)+(ps+a3(i))*f1p(i)
357 t21(i) = t21(i)+pd*f2p(i)
358 t22(i) = t22(i)+(ps+a1(i))*f2p(i)
359 t23(i) = t23(i)+(pp+a2(i))*f2p(i)
360 t24(i) = t24(i)+(ps+a3(i))*f2p(i)
361 ENDDO
362 IF(ale%UPWIND%UPWM == 2)THEN
363 DO i=lft,llt
364 fac=-dt1h/tr(i)
365 f1(i)=fac*(-f1p(i)*dzzp(i)+f2p(i)*dyzp(i))
366 f2(i)=fac*( f1p(i)*dzyp(i)-f2p(i)*dyyp(i))
367 t11(i) = t11(i)+pd*f1(i)
368 t12(i) = t12(i)+ps*f1(i)
369 t13(i) = t13(i)+pp*f1(i)
370 t14(i) = t14(i)+ps*f1(i)
371 t21(i) = t21(i)+pd*f2(i)
372 t22(i) = t22(i)+ps*f2(i)
373 t23(i) = t23(i)+pp*f2(i)
374 t24(i) = t24(i)+ps*f2(i)
375 ENDDO
376 ENDIF
377C NODE 4
378 CALL upwind(
379 1 rho, vis, vy4, vy4,
380 2 vz4, s, s, t,
381 3 deltax, gam, llt)
382 DO i=lft,llt
383 tr(i)=(y34(i)*z41(i)-y41(i)*z34(i))
384 dyyp(i)=(-z34(i)*vy1(i)-z41(i)*vy3(i)-z13(i)*vy4(i))
385 dyzp(i)=( y34(i)*vy1(i)+y41(i)*vy3(i)+y13(i)*vy4(i))
386 dzzp(i)=( y34(i)*vz1(i)+y41(i)*vz3(i)+y13(i)*vz4(i))
387 dzyp(i)=(-z34(i)*vz1(i)-z41(i)*vz3(i)-z13(i)*vz4(i))
388 f1p(i)=rho(i)*(vy4(i)*dyyp(i)+vz4(i)*dyzp(i))
389 f2p(i)=rho(i)*(vy4(i)*dzyp(i)+vz4(i)*dzzp(i))
390 fac=pc*gam(i)/tr(i)
391 a1(i)=fac*(-z34(i)*vy4(i)+y34(i)*vz4(i))
392 a2(i)=fac*(-z41(i)*vy4(i)+y41(i)*vz4(i))
393 a3(i)=fac*(-z13(i)*vy4(i)+y13(i)*vz4(i))
394 t11(i) = t11(i)+(ps+a1(i))*f1p(i)
395 t12(i) = t12(i)+pd*f1p(i)
396 t13(i) = t13(i)+(ps+a2(i))*f1p(i)
397 t14(i) = t14(i)+(pp+a3(i))*f1p(i)
398 t21(i) = t21(i)+(ps+a1(i))*f2p(i)
399 t22(i) = t22(i)+pd*f2p(i)
400 t23(i) = t23(i)+(ps+a2(i))*f2p(i)
401 t24(i) = t24(i)+(pp+a3(i))*f2p(i)
402 ENDDO
403 IF(ale%UPWIND%UPWM == 2)THEN
404 DO i=lft,llt
405 fac=-dt1h/tr(i)
406 f1(i)=fac*(-f1p(i)*dzzp(i)+f2p(i)*dyzp(i))
407 f2(i)=fac*( f1p(i)*dzyp(i)-f2p(i)*dyyp(i))
408 t11(i) = t11(i)+ps*f1(i)
409 t12(i) = t12(i)+pd*f1(i)
410 t13(i) = t13(i)+ps*f1(i)
411 t14(i) = t14(i)+pp*f1(i)
412 t21(i) = t21(i)+ps*f2(i)
413 t22(i) = t22(i)+pd*f2(i)
414 t23(i) = t23(i)+ps*f2(i)
415 t24(i) = t24(i)+pp*f2(i)
416 ENDDO
417 ENDIF
418 ENDIF
419C-----------------------------------------------
420 RETURN
421 END
my_real
#define my_real
Definition cppsort.cpp:32
emomt2
subroutine emomt2(pm, rho, y1, y2, y3, y4, z1, z2, z3, z4, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, t11, t12, t13, t14, t21, t22, t23, t24, py1, py2, pz1, pz2, aire, dyy, dzz, dyz, dzy, vdy, vdz, deltax, vis, mat, qmv, bufmat, ipm)
Definition emomt2.F:41
ale_mod
Definition ale_mod.F:149
ale_mod::ale
type(ale_) ale
Definition ale_mod.F:249
upwind
subroutine upwind(rho, vis, vdx, vdy, vdz, r, s, t, deltax, gam, nel)
Definition upwind.F:35