37
38
39
40#include "implicit_f.inc"
41
42
43
44#include "mvsiz_p.inc"
45
46
47
48#include "com01_c.inc"
49
50
51
52 INTEGER JFT,JLT
53 my_real,
DIMENSION(MVSIZ),
INTENT(IN) ::
54 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
56 . e1x, e2x, e3x, e1y, e2y, e3y ,e1z, e2z, e3z
57
58
59
60 INTEGER I
62 my_real,
DIMENSION(MVSIZ) :: x21,y21,z21,x31,y31,z31,x42,y42,z42,sum
63
64 DO i=jft,jlt
65 x21(i)=x2(i)-x1(i)
66 y21(i)=y2(i)-y1(i)
67 z21(i)=z2(i)-z1(i)
68 x31(i)=x3(i)-x1(i)
69 y31(i)=y3(i)-y1(i)
70 z31(i)=z3(i)-z1(i)
71 x42(i)=x4(i)-x2(i)
72 y42(i)=y4(i)-y2(i)
73 z42(i)=z4(i)-z2(i)
74 ENDDO
75
76 DO i=jft,jlt
77 e3x(i)=y31(i)*z42(i)-z31(i)*y42(i)
78 e3y(i)=z31(i)*x42(i)-x31(i)*z42(i)
79 e3z(i)=x31(i)*y42(i)-y31(i)*x42(i)
80 sum(i)=sqrt(e3x(i)*e3x(i)+e3y(i)*e3y(i)+e3z(i)*e3z(i))
82 ENDDO
83
84 IF (ishfram == 1) THEN
85
86 DO i=jft,jlt
87 e3x(i)=e3x(i)/sum(i)
88 e3y(i)=e3y(i)/sum(i)
89 e3z(i)=e3z(i)/sum(i)
90 ENDDO
91 DO i=jft,jlt
92 sum(i)= x21(i)*e3x(i)+y21(i)*e3y(i)+z21(i)*e3z(i)
93 e1x(i)= x21(i)-e3x(i)*sum(i)
94 e1y(i)= y21(i)-e3y(i)*sum(i)
95 e1z(i)= z21(i)-e3z(i)*sum(i)
96 ENDDO
97
98 DO i=jft,jlt
99 sum(i)=sqrt(e1x(i)*e1x(i)+e1y(i)*e1y(i)+e1z(i)*e1z(i))
100 e1x(i)=e1x(i)/sum(i)
101 e1y(i)=e1y(i)/sum(i)
102 e1z(i)=e1z(i)/sum(i)
103 ENDDO
104
105 DO i=jft,jlt
106 e2x(i)=e3y(i)*e1z(i)-e3z(i)*e1y(i)
107 e2y(i)=e3z(i)*e1x(i)-e3x(i)*e1z(i)
108 e2z(i)=e3x(i)*e1y(i)-e3y(i)*e1x(i)
109 ENDDO
110
111 DO i=jft,jlt
112 sum(i)=sqrt(e2x(i)*e2x(i)+e2y(i)*e2y(i)+e2z(i)*e2z(i))
113 e2x(i)=e2x(i)/sum(i)
114 e2y(i)=e2y(i)/sum(i)
115 e2z(i)=e2z(i)/sum(i)
116 ENDDO
117 ELSEIF (ishfram == 2) THEN
118
119 DO i=jft,jlt
120 e1x(i)= x2(i)+x3(i)-x1(i)-x4(i)
121 e1y(i)= y2(i)+y3(i)-y1(i)-y4(i)
122 e1z(i)= z2(i)+z3(i)-z1(i)-z4(i)
123
124 e2x(i)= x3(i)+x4(i)-x1(i)-x2(i)
125 e2y(i)= y3(i)+y4(i)-y1(i)-y2(i)
126 e2z(i)= z3(i)+z4(i)-z1(i)-z2(i)
127
128 e3x(i)=e1y(i)*e2z(i)-e1z(i)*e2y(i)
129 e3y(i)=e1z(i)*e2x(i)-e1x(i)*e2z(i)
130 e3z(i)=e1x(i)*e2y(i)-e1y(i)*e2x(i)
131 ENDDO
132
133 DO i=jft,jlt
134 suma=e2x(i)*e2x(i)+e2y(i)*e2y(i)+e2z(i)*e2z(i)
135 e1x(i) = e1x(i)*suma + e2y(i)*e3z(i)-e2z(i)*e3y(i)
136 e1y(i) = e1y(i)*suma + e2z(i)*e3x(i)-e2x(i)*e3z(i)
137 e1z(i) = e1z(i)*suma + e2x(i)*e3y(i)-e2y(i)*e3x(i)
138 ENDDO
139 DO i=jft,jlt
140 suma=e1x(i)*e1x(i)+e1y(i)*e1y(i)+e1z(i)*e1z(i)
141 suma=one/
max(sqrt(suma),em20)
142 e1x(i)=e1x(i)*suma
143 e1y(i)=e1y(i)*suma
144 e1z(i)=e1z(i)*suma
145 ENDDO
146
147 DO i=jft,jlt
148 suma=e3x(i)*e3x(i)+e3y(i)*e3y(i)+e3z(i)*e3z(i)
149 suma=one/
max(sqrt(suma),em20)
150 e3x(i)=e3x(i)*suma
151 e3y(i)=e3y(i)*suma
152 e3z(i)=e3z(i)*suma
153
154 e2x(i)=e3y(i)*e1z(i)-e3z(i)*e1y(i)
155 e2y(i)=e3z(i)*e1x(i)-e3x(i)*e1z(i)
156 e2z(i)=e3x(i)*e1y(i)-e3y(i)*e1x(i)
157 ENDDO
158
159 DO i=jft,jlt
160 suma=e2x(i)*e2x(i)+e2y(i)*e2y(i)+e2z(i)*e2z(i)
161 suma=one/
max(sqrt(suma),em20)
162 e2x(i)=e2x(i)*suma
163 e2y(i)=e2y(i)*suma
164 e2z(i)=e2z(i)*suma
165 ENDDO
166 ELSE
167
168 DO i=jft,jlt
169 e1x(i) = x2(i)+x3(i)-x1(i)-x4(i)
170 e1y(i) = y2(i)+y3(i)-y1(i)-y4(i)
171 e1z(i) = z2(i)+z3(i)-z1(i)-z4(i)
172
173 e2x(i) = x3(i)+x4(i)-x1(i)-x2(i)
174 e2y(i) = y3(i)+y4(i)-y1(i)-y2(i)
175 e2z(i) = z3(i)+z4(i)-z1(i)-z2(i)
176
177 e3x(i) = e1y(i)*e2z(i)-e1z(i)*e2y(i)
178 e3y(i) = e1z(i)*e2x(i)-e1x(i)*e2z(i)
179 e3z(i) = e1x(i)*e2y(i)-e1y(i)*e2x(i)
180 ENDDO
181
182 DO i=jft,jlt
183 suma = e3x(i)*e3x(i)+e3y(i)*e3y(i)+e3z(i)*e3z(i)
184 suma = one/
max(sqrt(suma),em20)
185 e3x(i) = e3x(i)*suma
186 e3y(i) = e3y(i)*suma
187 e3z(i) = e3z(i)*suma
188
189 s1 = e1x(i)*e1x(i)+e1y(i)*e1y(i)+e1z(i)*e1z(i)
190 s2 = e2x(i)*e2x(i)+e2y(i)*e2y(i)+e2z(i)*e2z(i)
191 suma = sqrt(s1/s2)
192 e1x(i) = e1x(i) + (e2y(i)*e3z(i)-e2z(i)*e3y(i))*suma
193 e1y(i) = e1y(i) + (e2z(i)*e3x(i)-e2x(i)*e3z(i))*suma
194 e1z(i) = e1z(i) + (e2x(i)*e3y(i)-e2y(i)*e3x(i))*suma
195
196 suma = e1x(i)*e1x(i)+e1y(i)*e1y(i)+e1z(i)*e1z(i)
197 suma = one/
max(sqrt(suma),em20)
198 e1x(i) = e1x(i)*suma
199 e1y(i) = e1y(i)*suma
200 e1z(i) = e1z(i)*suma
201
202 e2x(i) = e3y(i) * e1z(i) - e3z(i) * e1y(i)
203 e2y(i) = e3z(i) * e1x(i) - e3x(i) * e1z(i)
204 e2z(i) = e3x(i) * e1y(i) - e3y(i) * e1x(i)
205 ENDDO
206 ENDIF
207
208 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)