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