32 USE ebcs_mod
34
35
36
37#include "implicit_f.inc"
38
39
40
41#include "com08_c.inc"
42#include "scr11_c.inc"
43
44
45
46 INTEGER NSEG,NOD,ISEG(NSEG),(NOD),IRECT(4,NSEG)
47 my_real a(3,*),v(3,*),x(3,*),la(3,nod),ms(*),stifn(*),fv(*)
48 TYPE(t_ebcs_vel), INTENT(IN) :: EBCS
49 TYPE(t_segvar) :: SEGVAR
50
51
52
53 INTEGER I,IS,KSEG,N1,N2,N3,N4,NG1,NG2,NG3,NG4,N,
54 . IVX,IVY,IVZ,IRHO,IENER
56 . orient,rho,c,roc,fac,
57 . x13,y13,z13,x24,y24,z24,nx,ny,nz,s,
58 . roou,enou,vmx,vmy,vmz,fluxi,fluxo,p,dvx,dvy,dvz,ener,
59 . vx,vy,vz
60
61 ivx=ebcs%ivx
62 ivy=ebcs%ivy
63 ivz=ebcs%ivz
64 irho=ebcs%irho
65 iener=ebcs%iener
66 IF(ivx > 0)THEN
67 vx=ebcs%vx*fv(ivx)
68 ELSE
69 vx=ebcs%vx
70 ENDIF
71 IF(ivy > 0)THEN
72 vy=ebcs%vy*fv(ivy)
73 ELSE
74 vy=ebcs%vy
75 ENDIF
76 IF(ivz > 0)THEN
77 vz=ebcs%vz*fv(ivz)
78 ELSE
79 vz=ebcs%vz
80 ENDIF
81 IF(irho > 0)THEN
82 rho=ebcs%rho*fv(irho)
83 ELSE
84 rho=ebcs%rho
85 ENDIF
86 IF(iener > 0)THEN
87 ener=ebcs%ener*fv(iener)
88 ELSE
89 ener=ebcs%ener
90 ENDIF
91 c=ebcs%c
92 roc=rho*c
93
94
95
96 DO i=1,nod
97 la(1,i)=zero
98 la(2,i)=zero
99 la(3,i)=zero
100 ENDDO
101 DO is=1,nseg
102 kseg=abs(iseg(is))
103 orient=float(iseg(is)/kseg)
104 n1=irect(1,is)
105 n2=irect(2,is)
106 n3=irect(3,is)
107 n4=irect(4,is)
108 IF(n4==0 .OR. n4==n3) THEN
109 fac=one_over_6*orient
110 n4=n3
111 ELSE
112 fac=one_over_8*orient
113 ENDIF
114 ng1=liste(n1)
115 ng2=liste(n2)
116 ng3=liste(n3)
117 ng4=liste(n4)
118 x13=x(1,ng3)-x(1,ng1)
119 y13=x(2,ng3)-x(2,ng1)
120 z13=x(3,ng3)-x(3,ng1)
121 x24=x(1,ng4)-x(1,ng2)
122 y24=x(2,ng4)-x(2,ng2)
123 z24=x(3,ng4)-x(3,ng2)
124 nx=(y13*z24-z13*y24)*fac
125 ny=(z13*x24-x13*z24)*fac
126 nz=(x13*y24-y13*x24)*fac
127 la(1,n1)=la(1,n1)+nx
128 la(2,n1)=la(2,n1)+ny
129 la(3,n1)=la(3,n1)+nz
130 la(1,n2)=la(1,n2)+nx
131 la(2,n2)=la(2,n2)+ny
132 la(3,n2)=la(3,n2)+nz
133 la(1,n3)=la(1,n3)+nx
134 la(2,n3)=la(2,n3)+ny
135 la(3,n3)=la(3,n3)+nz
136 vmx=v(1,ng1)+v(1,ng2)+v(1,ng3)
137 vmy=v(2,ng1)+v(2,ng2)+v(2,ng3)
138 vmz=v(3,ng1)+v(3,ng2)+v(3,ng3)
139 IF(n4/=n3) THEN
140 la(1,n4)=la(1,n4)+nx
141 la(2,n4)=la(2,n4)+ny
142 la(3,n4)=la(3,n4)+nz
143 vmx=vmx+v(1,ng4)
144 vmy=vmy+v(2,ng4)
145 vmz=vmz+v(3,ng4)
146 ENDIF
147
148 roou = segvar%RHO(kseg)
149 enou = segvar%EINT(kseg)
150 fluxo=(vmx*nx+vmy*ny+vmz*nz)*dt1
151 fluxi=
min(fluxo,zero)
152 fluxo=
max(fluxo,zero)
153 dmf=dmf-fluxo*roou-fluxi*rho
154 def=def-fluxo*enou-fluxi*ener
155
156 segvar%RHO(kseg)=rho
157 segvar%EINT(kseg)=ener
158 ENDDO
159 IF(tt == 0)THEN
160 DO i=1,nod
161 n=liste(i)
162 v(1,n)=vx
163 v(2,n)=vy
164 v(3,n)=vz
165 ENDDO
166 ENDIF
167
168 DO i=1,nod
169 n=liste(i)
170 s=sqrt(la(1,i)**2+la(2,i)**2+la(3,i)**2)
171 dvx=v(1,n)-vx
172 dvy=v(2,n)-vy
173 dvz=v(3,n)-vz
174 p=roc*(dvx*la(1,i)+dvy*la(2,i)+dvz*la(3,i))/s
175 a(1,n)=a(1,n)-p*la(1,i)
176 a(2,n)=a(2,n)-p*la(2,i)
177 a(3,n)=a(3,n)-p*la(3,i)
178 stifn(n)=stifn(n)+(two*(s*roc)**2)/ms(n)
179 ENDDO
180
181 RETURN