50
51
52
53#include "implicit_f.inc"
54
55
56
57#include "mvsiz_p.inc"
58
59
60
61
62
63
64 INTEGER, INTENT(INOUT) :: LFT
65 INTEGER, INTENT(INOUT) :: LLT
66 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: x1,x2,x3,x4,xi
67 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: y1,y2,y3,y4,yi
68 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: z1,z2,z3,z4,zi
69 my_real,
DIMENSION(MVSIZ),
INTENT(IN) ::
70 my_real,
DIMENSION(MVSIZ),
INTENT(OUT) :: ssc,ttc,
area
71 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: n1,n2,n3
72 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: x0,y0,z0
73 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: xx1,xx2,xx3,xx4
74 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: yy1,yy2,yy3,yy4
75 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: zz1,zz2,zz3,zz4
76 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: xi1,xi2,xi3
77 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: yi1,yi2,yi3,yi4
78 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: zi1,zi2,zi3,zi4
79 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT)
80 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: yn1,yn2,yn3,yn4
81 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: zn1,zn2,zn3,zn4
82
83
84
85 INTEGER I
87 . a12(mvsiz), a23(mvsiz), a34(mvsiz), a41(mvsiz),
88 . b12(mvsiz), b23(mvsiz), b34(mvsiz), b41(mvsiz),
89 . ab1(mvsiz), ab2(mvsiz), an
90
91 DO 100 i=lft,llt
92 x0(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
93 y0(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
94 z0(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i))
95
96 xx1(i) = x1(i)-x0(i)
97 xx2(i) = x2(i)-x0(i)
98 xx3(i) = x3(i)-x0(i)
99 xx4(i) = x4(i)-x0(i)
100 yy1(i) = y1(i)-y0(i)
101 yy2(i) = y2(i)-y0(i)
102 yy3(i) = y3(i)-y0(i)
103 yy4(i) = y4(i)-y0(i)
104 zz1(i) = z1(i)-z0(i)
105 zz2(i) = z2(i)-z0(i)
106 zz3(i) = z3(i)-z0(i)
107 zz4(i) = z4(i)-z0(i)
108
109 xi1(i) = x1(i)-xi(i)
110 xi2(i) = x2(i)-xi(i)
111 xi3(i) = x3(i)-xi(i)
112 xi4(i) = x4(i)-xi(i)
113 yi1(i) = y1(i)-yi(i)
114 yi2(i) = y2(i)-yi(i)
115 yi3(i) = y3(i)-yi(i)
116 yi4(i) = y4(i)-yi(i)
117 zi1(i) = z1(i)-zi(i)
118 zi2(i) = z2(i)-zi(i)
119 zi3(i) = z3(i)-zi(i)
120 zi4(i) = z4(i)-zi(i)
121 100 CONTINUE
122
123 DO 120 i=lft,llt
124 xn1(i) = yy1(i)*zz2(i) - yy2(i)*zz1(i)
125 yn1(i) = zz1(i)*xx2(i) - zz2(i)*xx1(i)
126 zn1(i) = xx1(i)*yy2(i) - xx2(i)*yy1(i)
127 n1(i)=xn1(i)
128 n2(i)=yn1(i)
129 n3(i)=zn1(i)
130 120 CONTINUE
131
132 DO 140 i=lft,llt
133 xn2(i) = yy2(i)*zz3(i) - yy3(i)*zz2(i)
134 yn2(i) = zz2(i)*xx3(i) - zz3(i)*xx2(i)
135 zn2(i) = xx2(i)*yy3(i) - xx3(i)*yy2(i)
136 n1(i)=n1(i)+xn2(i)
137 n2(i)=n2(i)+yn2(i)
138 n3(i)=n3(i)+zn2(i)
139 140 CONTINUE
140
141 DO 160 i=lft,llt
142 xn3(i) = yy3(i)*zz4(i) - yy4(i)*zz3(i)
143 yn3(i) = zz3(i)*xx4(i) - zz4(i)*xx3(i)
144 zn3(i) = xx3(i)*yy4(i) - xx4(i)*yy3(i)
145 n1(i)=n1(i)+xn3(i)
146 n2(i)=n2(i)+yn3(i)
147 n3(i)=n3(i)+zn3(i)
148 160 CONTINUE
149
150 DO 180 i=lft,llt
151 xn4(i) = yy4(i)*zz1(i) - yy1(i)*zz4(i)
152 yn4(i) = zz4(i)*xx1(i) - zz1(i)*xx4(i)
153 zn4(i) = xx4(i)*yy1(i) - xx1(i)*yy4(i)
154 n1(i)=n1(i)+xn4(i)
155 n2(i)=n2(i)+yn4(i)
156 n3(i)=n3(i)+zn4(i)
157 180 CONTINUE
158
159 DO 200 i=lft,llt
160 an=
max(em20,sqrt(n1(i)*n1(i)+n2
161 n1(i)=n1(i)/an
162 n2(i)=n2(i)/an
163 n3(i)=n3(i)/an
165 200 CONTINUE
166
167 DO 210 i=lft,llt
168 x0(i)=(n1(i)*xn1(i)+n2(i)*yn1(i)+n3(i)*zn1(i))
169 y0(i)=(n1(i)*xn2(i)+n2(i)*yn2(i)+n3(i)*zn2(i))
170 z0(i)=(n1(i)*xn3(i)+n2(i)*yn3(i)+n3(i)*zn3(i))
171 xx1(i)=(n1(i)*xn4(i)+n2(i)*yn4(i)+n3(i)*zn4(i))
172 210 CONTINUE
173
174 DO 220 i=lft,llt
175 xn1(i) = yi1(i)*zi2(i) - yi2(i)*zi1(i)
176 yn1(i) = zi1(i)*xi2(i) - zi2(i)*xi1(i)
177 zn1(i) = xi1(i)*yi2(i) - xi2(i)*yi1(i)
178 yy1(i)=(n1(i)*xn1(i)+n2(i)*yn1(i)+n3(i)*zn1(i))
179 220 CONTINUE
180
181 DO 240 i=lft,llt
182 xn2(i) = yi2(i)*zi3(i) - yi3(i)*zi2(i)
183 yn2(i) = zi2(i)*xi3(i) - zi3(i)*xi2(i)
184 zn2(i) = xi2(i)*yi3(i) - xi3(i)*yi2(i)
185 zz1(i)=(n1(i)*xn2(i)+n2(i)*yn2(i)+n3(i)*zn2(i))
186 240 CONTINUE
187
188 DO 260 i=lft,llt
189 xn3(i) = yi3(i)*zi4(i) - yi4(i)*zi3(i)
190 yn3(i) = zi3(i)*xi4(i) - zi4(i)*xi3(i)
191 zn3(i) = xi3(i)*yi4(i) - xi4(i)*yi3(i)
192 xx2(i)=(n1(i)*xn3(i)+n2(i)*yn3(i)+n3(i)*zn3(i))
193 260 CONTINUE
194
195 DO 280 i=lft,llt
196 xn4(i) = yi4(i)*zi1(i) - yi1(i)*zi4(i)
197 yn4(i) = zi4(i)*xi1(i) - zi1(i)*xi4(i)
198 zn4(i) = xi4(i)*yi1(i) - xi1(i)*yi4(i)
199 yy2(i)=(n1(i)*xn4(i)+n2(i)*yn4(i)+n3(i)*zn4(i))
200 280 CONTINUE
201
202 DO 300 i=lft,llt
203 zz2(i)=y0(i)*yy2(i)
204 xx3(i)=zz1(i)*xx1(i)
205 300 CONTINUE
206
207 DO 320 i=lft,llt
208 IF(xface(i)==zero)GOTO 320
209 IF(zz2(i)+xx3(i)/=zero)THEN
210 ssc(i)=(zz2(i)-xx3(i))/(zz2(i)+xx3(i))
211 ELSE
212 ssc(i)=zero
213 ENDIF
214 IF(z0(i)/=zero)THEN
215 zz2(i)=yy1(i)*z0(i)
216 xx3(i)=xx2(i)*x0(i)
217 IF(zz2(i)+xx3(i)/=zero)THEN
218 ttc(i)=(zz2(i)-xx3(i))/(zz2(i)+xx3(i))
219 ELSE
220 ttc(i)=zero
221 ENDIF
222 ELSE
223 ttc(i)=(yy1(i)
224 ENDIF
225 320 CONTINUE
226
227 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)