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