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