35
36
37
39
40
41
42#include "implicit_f.inc"
43
44
45
46#include "com01_c.inc"
47#include "com04_c.inc"
48#include "scr08_a_c.inc"
49#include "param_c.inc"
50#include "warn_c.inc"
51#include "tabsiz_c.inc"
52
53
54
55 INTEGER IRECT(4,*), MSR(*), NSV(*), IRTL(*), LCODE(*), ISKEW(*), ITAB(*)
56 my_real skew(lskew,*), a(sa), e(*), msm(*), crst(2,*), v(sv),nor(3,*)
57 INTEGER, INTENT(IN) :: NSN,NMN
58
59
60
61 INTEGER JBC(3), NIR, I, J, I3, J3, I2, J2, I1, J1, ISK, LCOD, II, L, JJ, NN, LK, IBC
62 my_real xn(3), yn(3), zn(3), h(4), n1, n2, n3, ss, tt, ax, ay, az,
63 . vx, vy, vz, amn, vmn, amod, vmod, bvz, baz, bvx, bvy, bax, bay,
64 . a11, a12, a13, a21, a22, a23, a31, a32, a33, det
65
66
67
68 nir=2
69 IF(n2d == 0)nir=4
70
71
72
73 DO i=1,nmn
74 j=msr(i)
75 i3=3*i
76 j3=3*j
77 i2=i3-1
78 j2=j3-1
79 i1=i2-1
80 j1=j2-1
81 IF(msm(i) > zero)THEN
82 a(j1)=a(j1)+e(i1)/msm(i)
83 a(j2)=a(j2)+e(i2)/msm(i)
84 a(j3)=a(j3)+e(i3)/msm(i)
85 ENDIF
86 isk = iskew(j)
87 lcod = lcode(j)
88 IF(lcod>0)
CALL bcs2(a(j1),skew(1,isk),isk,lcod)
89 ENDDO
90
91
92
93 DO ii=1,nsn
94 i=nsv(ii)
95 l=irtl(ii)
96 DO jj=1,nir
97 nn=irect(jj,l)
98 ix(jj)=msr(nn)
99 ENDDO
100 ss=crst(1,ii)
101 tt=crst(2,ii)
102 n1=nor(1,ii)
103 n2=nor(2,ii)
104 n3=nor(3,ii)
106 i3=3*i
107 i2=i3-1
108 i1=i2-1
109 ax=zero
110 ay=zero
111 az=zero
112 vx=zero
113 vy=zero
114 vz=zero
115 DO jj=1,nir
116 j3=3*ix(jj)
117 j2=j3-1
118 j1=j2-1
119 ax=ax+a(j1)*h(jj)
120 ay=ay+a(j2)*h(jj)
121 az=az+a(j3)*h(jj)
122 vx=vx+v(j1)*h(jj)
123 vy=vy+v(j2)*h(jj)
124 vz=vz+v(j3)*h(jj)
125 ENDDO
126
127 amn = n1*ax+n2*ay+n3*az
128 vmn = n1*vx+n2*vy+n3*vz
129 amod= amn-n1*a(i1)-n2*a(i2)-n3*a(i3)
130 vmod= vmn-n1*v(i1)-n2*v(i2)-n3*v(i3)
131
132 lcod=lcode(i)
133
134 IF(lcod /= 0)THEN
135
136
137
138 xn(1)=n1
139 yn(1)=n2
140 zn(1)=n3
141
142 jbc(1:3) = 0
143 IF(lcode(i) > 0) THEN
144 jbc(3) = iand(lcode(i), 1)
145 jbc(2) = iand(lcode(i), 2)
146 jbc(1) = iand(lcode(i), 4)
147 ENDIF
148
149 lk=iskew(i)
150 ibc=2
151
152 IF(jbc(1) /= 0)THEN
153 xn(ibc)=skew(1,lk)
154 yn(ibc)=skew(2,lk)
155 zn(ibc)=skew(3,lk)
156 ibc=ibc+1
157 ENDIF
158 IF(jbc(2) /= 0)THEN
159 xn(ibc)=skew(4,lk)
160 yn(ibc)=skew(5,lk)
161 zn(ibc)=skew(6,lk)
162 ibc=ibc+1
163 ENDIF
164 IF(jbc(3) /= 0)THEN
165 IF(ibc == 4)THEN
166
167 CALL ancmsg(msgid=11,anmode=aninfo,i1=itab(i))
168 ierr=1
169 EXIT
170 ELSE
171 xn(ibc)=skew(7,lk)
172 yn(ibc)=skew(8,lk)
173 zn(ibc)=skew(9,lk)
174 ibc=ibc+1
175 ENDIF
176 ENDIF
177 IF(ibc == 3)THEN
178
179
180
181 xn(3)=yn(1)*zn(2)-zn(1)*yn(2)
182 yn(3)=zn(1)*xn(2)-xn(1)*zn(2)
183 zn(3)=xn(1)*yn(2)-yn(1)*xn(2)
184 bvz=v(i1)*xn(3)+v(i2)*yn(3)+v(i3)*zn(3)
185 baz=a(i1)*xn(3)+a(i2)*yn(3)+a(i3)*zn(3)
186 ELSE
187
188
189
190 bvz=zero
191 baz=zero
192 ENDIF
193
194 bvx=vmn
195 bvy=zero
196 bax=amn
197 bay=zero
198
199 a11=yn(2)*zn(3)-zn(2)*yn(3)
200 a12=zn(2)*xn(3)-xn(2)*zn(3)
201 a13=xn(2)*yn(3)-yn(2)*xn(3)
202 a21=yn(3)*zn(1)-zn(3)*yn(1)
203 a22=zn(3)*xn(1)-xn(3)*zn(1)
204 a23=xn(3)*yn(1)-yn(3)*xn(1)
205 a31=yn(1)*zn(2)-zn(1)*yn(2)
206 a32=zn(1)*xn(2)-xn(1)*zn(2)
207 a33=xn(1)*yn(2)-yn(1)*xn(2)
208 det=xn(1)*a11+yn(1)*a12+zn(1)*a13
209
210
211
212
213 IF(det /= zero)THEN
214 v(i1)=(a11*bvx+a21*bvy+a31*bvz)/det
215 v(i2)=(a12*bvx+a22*bvy+a32*bvz)/det
216 v(i3)=(a13*bvx+a23*bvy+a33*bvz)/det
217 a(i1)=(a11*bax+a21*bay+a31*baz)/det
218 a(i2)=(a12*bax+a22*bay+a32*baz)/det
219 a(i3)=(a13*bax+a23*bay+a33*baz)/det
220 ENDIF
221
222 ELSEIF(lcod == 0)THEN
223
224
225
226 a(i1)=a(i1)+amod*n1
227 a(i2)=a(i2)+amod*n2
228 a(i3)=a(i3)+amod*n3
229 v(i1)=v(i1)+vmod*n1
230 v(i2)=v(i2)+vmod*n2
231 v(i3)=v(i3)+vmod*n3
232 ENDIF
233
234 ENDDO
235
236 RETURN
subroutine bcs2(a, b, j, k)
subroutine shapeh(h, s, t)
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)