xref: /phasta/phSolver/compressible/e3bdg.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)
1*59599516SKenneth E. Jansen        subroutine e3bdg (shp, 	   shg,    WdetJ,
2*59599516SKenneth E. Jansen     &		          A1, 	   A2, 	   A3,
3*59599516SKenneth E. Jansen     &                    A0,      bcool,  tau,
4*59599516SKenneth E. Jansen     &		          u1,      u2,     u3,
5*59599516SKenneth E. Jansen     &			  BDiagl,
6*59599516SKenneth E. Jansen     &                    rmu,     rlm2mu, con)
7*59599516SKenneth E. Jansen        include "common.h"
8*59599516SKenneth E. Jansen
9*59599516SKenneth E. Jansenc
10*59599516SKenneth E. Jansenc  passed arrays
11*59599516SKenneth E. Jansenc
12*59599516SKenneth E. Jansen        dimension BDiagl(npro,nshl,nflow,nflow),
13*59599516SKenneth E. Jansen     &            shp(npro,nshl),          bcool(npro),
14*59599516SKenneth E. Jansen     &            shg(npro,nshl,nsd),      WdetJ(npro),
15*59599516SKenneth E. Jansen     &            A1tauA0(npro,nflow,nflow), A2tauA0(npro,nflow,nflow),
16*59599516SKenneth E. Jansen     &            A3tauA0(npro,nflow,nflow), Atau(npro,nflow,nflow),
17*59599516SKenneth E. Jansen     &            A1(npro,nflow,nflow),      A2(npro,nflow,nflow),
18*59599516SKenneth E. Jansen     &            A3(npro,nflow,nflow),      A0(npro,nflow,nflow),
19*59599516SKenneth E. Jansen     &            tau(npro,5),             u1(npro),
20*59599516SKenneth E. Jansen     &            u2(npro),                u3(npro),
21*59599516SKenneth E. Jansen     &            rlm2mu(npro),
22*59599516SKenneth E. Jansen     &            rmu(npro),               con(npro)
23*59599516SKenneth E. Jansenc
24*59599516SKenneth E. Jansenc
25*59599516SKenneth E. Jansenc  passed work arrays for local variables
26*59599516SKenneth E. Jansenc
27*59599516SKenneth E. Jansen	dimension tmp1(npro),             tmp2(npro)
28*59599516SKenneth E. Jansen	dimension tmp(npro)
29*59599516SKenneth E. Jansen
30*59599516SKenneth E. Jansen	ttim(30) = ttim(30) - secs(0.0)
31*59599516SKenneth E. Jansen
32*59599516SKenneth E. Jansenc $$ Ex-E3conv
33*59599516SKenneth E. Jansenc
34*59599516SKenneth E. Jansenc.... calculate the contribution in non-integrated by part form
35*59599516SKenneth E. Jansenc.... add (N_b A_tilde_i N_b,i) to BDiagl
36*59599516SKenneth E. Jansenc
37*59599516SKenneth E. Jansen       do j = 1, nshl   ! May be worth eliminating zeros in A(prim) matrices
38*59599516SKenneth E. Jansen        tmp=shp(:,j)*WdetJ
39*59599516SKenneth E. Jansen            BDiagl(:,j,1,1) = BDiagl(:,j,1,1)
40*59599516SKenneth E. Jansen     &                          + tmp * (
41*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,1,1) +
42*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,1,1) +
43*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,1,1)
44*59599516SKenneth E. Jansen     &                                                            )
45*59599516SKenneth E. Jansen            BDiagl(:,j,1,2) = BDiagl(:,j,1,2)
46*59599516SKenneth E. Jansen     &                          + tmp * (
47*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,1,2)
48*59599516SKenneth E. Jansenc    &                                  +shg(:,j,2) * A2(:,1,2)
49*59599516SKenneth E. Jansenc    &                                  +shg(:,j,3) * A3(:,1,2)
50*59599516SKenneth E. Jansen     &                                                            )
51*59599516SKenneth E. Jansen            BDiagl(:,j,1,3) = BDiagl(:,j,1,3)
52*59599516SKenneth E. Jansen     &                          + tmp * (
53*59599516SKenneth E. Jansenc    &                                   shg(:,j,1) * A1(:,1,3)
54*59599516SKenneth E. Jansen     &                                  +shg(:,j,2) * A2(:,1,3)
55*59599516SKenneth E. Jansenc    &                                  +shg(:,j,3) * A3(:,1,3)
56*59599516SKenneth E. Jansen     &                                                            )
57*59599516SKenneth E. Jansen            BDiagl(:,j,1,4) = BDiagl(:,j,1,4)
58*59599516SKenneth E. Jansen     &                          + tmp * (
59*59599516SKenneth E. Jansenc    &                                   shg(:,j,1) * A1(:,1,4) +
60*59599516SKenneth E. Jansenc    &                                   shg(:,j,2) * A2(:,1,4) +
61*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,1,4)
62*59599516SKenneth E. Jansen     &                                                            )
63*59599516SKenneth E. Jansen            BDiagl(:,j,1,5) = BDiagl(:,j,1,5)
64*59599516SKenneth E. Jansen     &                          + tmp * (
65*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,1,5) +
66*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,1,5) +
67*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,1,5)
68*59599516SKenneth E. Jansen     &                                                            )
69*59599516SKenneth E. Jansen            BDiagl(:,j,2,1) = BDiagl(:,j,2,1)
70*59599516SKenneth E. Jansen     &                          + tmp * (
71*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,2,1) +
72*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,2,1) +
73*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,2,1)
74*59599516SKenneth E. Jansen     &                                                            )
75*59599516SKenneth E. Jansen            BDiagl(:,j,2,2) = BDiagl(:,j,2,2)
76*59599516SKenneth E. Jansen     &                          + tmp * (
77*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,2,2) +
78*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,2,2) +
79*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,2,2)
80*59599516SKenneth E. Jansen     &                                                            )
81*59599516SKenneth E. Jansen            BDiagl(:,j,2,3) = BDiagl(:,j,2,3)
82*59599516SKenneth E. Jansen     &                          + tmp * (
83*59599516SKenneth E. Jansenc    &                                   shg(:,j,1) * A1(:,2,3)
84*59599516SKenneth E. Jansen     &                                  +shg(:,j,2) * A2(:,2,3)
85*59599516SKenneth E. Jansenc    &                                  +shg(:,j,3) * A3(:,2,3)
86*59599516SKenneth E. Jansen     &                                                            )
87*59599516SKenneth E. Jansen            BDiagl(:,j,2,4) = BDiagl(:,j,2,4)
88*59599516SKenneth E. Jansen     &                          + tmp * (
89*59599516SKenneth E. Jansenc    &                                   shg(:,j,1) * A1(:,2,4) +
90*59599516SKenneth E. Jansenc    &                                   shg(:,j,2) * A2(:,2,4) +
91*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,2,4)
92*59599516SKenneth E. Jansen     &                                                            )
93*59599516SKenneth E. Jansen            BDiagl(:,j,2,5) = BDiagl(:,j,2,5)
94*59599516SKenneth E. Jansen     &                          + tmp * (
95*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,2,5) +
96*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,2,5) +
97*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,2,5)
98*59599516SKenneth E. Jansen     &                                                            )
99*59599516SKenneth E. Jansen            BDiagl(:,j,3,1) = BDiagl(:,j,3,1)
100*59599516SKenneth E. Jansen     &                          + tmp * (
101*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,3,1) +
102*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,3,1) +
103*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,3,1)
104*59599516SKenneth E. Jansen     &                                                            )
105*59599516SKenneth E. Jansen            BDiagl(:,j,3,2) = BDiagl(:,j,3,2)
106*59599516SKenneth E. Jansen     &                          + tmp * (
107*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,3,2)
108*59599516SKenneth E. Jansenc    &                                  +shg(:,j,2) * A2(:,3,2)
109*59599516SKenneth E. Jansenc    &                                  +shg(:,j,3) * A3(:,3,2)
110*59599516SKenneth E. Jansen     &                                                            )
111*59599516SKenneth E. Jansen            BDiagl(:,j,3,3) = BDiagl(:,j,3,3)
112*59599516SKenneth E. Jansen     &                          + tmp * (
113*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,3,3) +
114*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,3,3) +
115*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,3,3)
116*59599516SKenneth E. Jansen     &                                                            )
117*59599516SKenneth E. Jansen            BDiagl(:,j,3,4) = BDiagl(:,j,3,4)
118*59599516SKenneth E. Jansen     &                          + tmp * (
119*59599516SKenneth E. Jansenc    &                                   shg(:,j,1) * A1(:,3,4) +
120*59599516SKenneth E. Jansenc    &                                   shg(:,j,2) * A2(:,3,4) +
121*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,3,4)
122*59599516SKenneth E. Jansen     &                                                            )
123*59599516SKenneth E. Jansen            BDiagl(:,j,3,5) = BDiagl(:,j,3,5)
124*59599516SKenneth E. Jansen     &                          + tmp * (
125*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,3,5) +
126*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,3,5) +
127*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,3,5)
128*59599516SKenneth E. Jansen     &                                                            )
129*59599516SKenneth E. Jansen            BDiagl(:,j,4,1) = BDiagl(:,j,4,1)
130*59599516SKenneth E. Jansen     &                          + tmp * (
131*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,4,1) +
132*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,4,1) +
133*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,4,1)
134*59599516SKenneth E. Jansen     &                                                            )
135*59599516SKenneth E. Jansen            BDiagl(:,j,4,2) = BDiagl(:,j,4,2)
136*59599516SKenneth E. Jansen     &                          + tmp * (
137*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,4,2)
138*59599516SKenneth E. Jansenc    &                                  +shg(:,j,2) * A2(:,4,2)
139*59599516SKenneth E. Jansenc    &                                  +shg(:,j,3) * A3(:,4,2)
140*59599516SKenneth E. Jansen     &                                                            )
141*59599516SKenneth E. Jansen            BDiagl(:,j,4,3) = BDiagl(:,j,4,3)
142*59599516SKenneth E. Jansen     &                          + tmp * (
143*59599516SKenneth E. Jansenc    &                                   shg(:,j,1) * A1(:,4,3)
144*59599516SKenneth E. Jansen     &                                  +shg(:,j,2) * A2(:,4,3)
145*59599516SKenneth E. Jansenc    &                                  +shg(:,j,3) * A3(:,4,3)
146*59599516SKenneth E. Jansen     &                                                            )
147*59599516SKenneth E. Jansen            BDiagl(:,j,4,4) = BDiagl(:,j,4,4)
148*59599516SKenneth E. Jansen     &                          + tmp * (
149*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,4,4) +
150*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,4,4) +
151*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,4,4)
152*59599516SKenneth E. Jansen     &                                                            )
153*59599516SKenneth E. Jansen            BDiagl(:,j,4,5) = BDiagl(:,j,4,5)
154*59599516SKenneth E. Jansen     &                          + tmp * (
155*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,4,5) +
156*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,4,5) +
157*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,4,5)
158*59599516SKenneth E. Jansen     &                                                            )
159*59599516SKenneth E. Jansen            BDiagl(:,j,5,1) = BDiagl(:,j,5,1)
160*59599516SKenneth E. Jansen     &                          + tmp * (
161*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,5,1) +
162*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,5,1) +
163*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,5,1)
164*59599516SKenneth E. Jansen     &                                                            )
165*59599516SKenneth E. Jansen            BDiagl(:,j,5,2) = BDiagl(:,j,5,2)
166*59599516SKenneth E. Jansen     &                          + tmp * (
167*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,5,2) +
168*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,5,2) +
169*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,5,2)
170*59599516SKenneth E. Jansen     &                                                            )
171*59599516SKenneth E. Jansen            BDiagl(:,j,5,3) = BDiagl(:,j,5,3)
172*59599516SKenneth E. Jansen     &                          + tmp * (
173*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,5,3) +
174*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,5,3) +
175*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,5,3)
176*59599516SKenneth E. Jansen     &                                                            )
177*59599516SKenneth E. Jansen            BDiagl(:,j,5,4) = BDiagl(:,j,5,4)
178*59599516SKenneth E. Jansen     &                          + tmp * (
179*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,5,4) +
180*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,5,4) +
181*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,5,4)
182*59599516SKenneth E. Jansen     &                                                            )
183*59599516SKenneth E. Jansen            BDiagl(:,j,5,5) = BDiagl(:,j,5,5)
184*59599516SKenneth E. Jansen     &                          + tmp * (
185*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1(:,5,5) +
186*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2(:,5,5) +
187*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3(:,5,5)
188*59599516SKenneth E. Jansen     &                                                            )
189*59599516SKenneth E. Jansen       enddo
190*59599516SKenneth E. Jansen        if (ngauss .eq. 1 .and. nshl .eq. 4) then  ! Exact integ of mass term for tets
191*59599516SKenneth E. Jansen
192*59599516SKenneth E. Jansenc $$ Ex-e3mael
193*59599516SKenneth E. Jansenc
194*59599516SKenneth E. Jansenc.... calculate factor   V             * sum(column of M)
195*59599516SKenneth E. Jansenc                     =WdetJ/(3/4)Qwt(lcsyst,intp) * (2+1+1+1)/20
196*59599516SKenneth E. Jansenc
197*59599516SKenneth E. Jansen	tmp =  0.4d0 * WdetJ * Dtgl/Qwt(lcsyst,intp)/three*almi/gami/alfi
198*59599516SKenneth E. Jansenc
199*59599516SKenneth E. Jansen        do j=1,nshl  ! take advantage of zeros in A0(Prim)
200*59599516SKenneth E. Jansen	  BDiagl(:,j,2,2) = BDiagl(:,j,2,2) + tmp * A0(:,2,2)
201*59599516SKenneth E. Jansen	  BDiagl(:,j,3,3) = BDiagl(:,j,3,3) + tmp * A0(:,3,3)
202*59599516SKenneth E. Jansen	  BDiagl(:,j,4,4) = BDiagl(:,j,4,4) + tmp * A0(:,4,4)
203*59599516SKenneth E. Jansen	  BDiagl(:,j,1,5) = BDiagl(:,j,1,5) + tmp * A0(:,1,5)
204*59599516SKenneth E. Jansen	  BDiagl(:,j,2,5) = BDiagl(:,j,2,5) + tmp * A0(:,2,5)
205*59599516SKenneth E. Jansen	  BDiagl(:,j,3,5) = BDiagl(:,j,3,5) + tmp * A0(:,3,5)
206*59599516SKenneth E. Jansen	  BDiagl(:,j,4,5) = BDiagl(:,j,4,5) + tmp * A0(:,4,5)
207*59599516SKenneth E. Jansen	  BDiagl(:,j,1,1) = BDiagl(:,j,1,1) + tmp * A0(:,1,1)
208*59599516SKenneth E. Jansen	  BDiagl(:,j,2,1) = BDiagl(:,j,2,1) + tmp * A0(:,2,1)
209*59599516SKenneth E. Jansen	  BDiagl(:,j,3,1) = BDiagl(:,j,3,1) + tmp * A0(:,3,1)
210*59599516SKenneth E. Jansen	  BDiagl(:,j,4,1) = BDiagl(:,j,4,1) + tmp * A0(:,4,1)
211*59599516SKenneth E. Jansen	  BDiagl(:,j,5,1) = BDiagl(:,j,5,1) + tmp * A0(:,5,1)
212*59599516SKenneth E. Jansen	  BDiagl(:,j,5,2) = BDiagl(:,j,5,2) + tmp * A0(:,5,2)
213*59599516SKenneth E. Jansen	  BDiagl(:,j,5,3) = BDiagl(:,j,5,3) + tmp * A0(:,5,3)
214*59599516SKenneth E. Jansen	  BDiagl(:,j,5,4) = BDiagl(:,j,5,4) + tmp * A0(:,5,4)
215*59599516SKenneth E. Jansen	  BDiagl(:,j,5,5) = BDiagl(:,j,5,5) + tmp * A0(:,5,5)
216*59599516SKenneth E. Jansen	enddo
217*59599516SKenneth E. Jansen
218*59599516SKenneth E. Jansen        else
219*59599516SKenneth E. Jansen
220*59599516SKenneth E. Jansenc $$ Ex-e3mass
221*59599516SKenneth E. Jansenc
222*59599516SKenneth E. Jansen        ff = almi / gami / alfi
223*59599516SKenneth E. Jansen        tmp = WdetJ * (Dtgl * ff + bcool)
224*59599516SKenneth E. Jansenc
225*59599516SKenneth E. Jansen        do j = 1, nshl
226*59599516SKenneth E. Jansenc
227*59599516SKenneth E. Jansen          tmp2 = (shp(:,j)*shp(:,j)) * tmp
228*59599516SKenneth E. Jansenc
229*59599516SKenneth E. Jansen	  BDiagl(:,j,2,2) = BDiagl(:,j,2,2) + tmp2 * A0(:,2,2)
230*59599516SKenneth E. Jansen	  BDiagl(:,j,3,3) = BDiagl(:,j,3,3) + tmp2 * A0(:,3,3)
231*59599516SKenneth E. Jansen	  BDiagl(:,j,4,4) = BDiagl(:,j,4,4) + tmp2 * A0(:,4,4)
232*59599516SKenneth E. Jansen	  BDiagl(:,j,1,5) = BDiagl(:,j,1,5) + tmp2 * A0(:,1,5)
233*59599516SKenneth E. Jansen	  BDiagl(:,j,2,5) = BDiagl(:,j,2,5) + tmp2 * A0(:,2,5)
234*59599516SKenneth E. Jansen	  BDiagl(:,j,3,5) = BDiagl(:,j,3,5) + tmp2 * A0(:,3,5)
235*59599516SKenneth E. Jansen	  BDiagl(:,j,4,5) = BDiagl(:,j,4,5) + tmp2 * A0(:,4,5)
236*59599516SKenneth E. Jansen	  BDiagl(:,j,1,1) = BDiagl(:,j,1,1) + tmp2 * A0(:,1,1)
237*59599516SKenneth E. Jansen	  BDiagl(:,j,2,1) = BDiagl(:,j,2,1) + tmp2 * A0(:,2,1)
238*59599516SKenneth E. Jansen	  BDiagl(:,j,3,1) = BDiagl(:,j,3,1) + tmp2 * A0(:,3,1)
239*59599516SKenneth E. Jansen	  BDiagl(:,j,4,1) = BDiagl(:,j,4,1) + tmp2 * A0(:,4,1)
240*59599516SKenneth E. Jansen	  BDiagl(:,j,5,1) = BDiagl(:,j,5,1) + tmp2 * A0(:,5,1)
241*59599516SKenneth E. Jansen	  BDiagl(:,j,5,2) = BDiagl(:,j,5,2) + tmp2 * A0(:,5,2)
242*59599516SKenneth E. Jansen	  BDiagl(:,j,5,3) = BDiagl(:,j,5,3) + tmp2 * A0(:,5,3)
243*59599516SKenneth E. Jansen	  BDiagl(:,j,5,4) = BDiagl(:,j,5,4) + tmp2 * A0(:,5,4)
244*59599516SKenneth E. Jansen	  BDiagl(:,j,5,5) = BDiagl(:,j,5,5) + tmp2 * A0(:,5,5)
245*59599516SKenneth E. Jansenc
246*59599516SKenneth E. Jansen        enddo
247*59599516SKenneth E. Jansen
248*59599516SKenneth E. Jansen        endif
249*59599516SKenneth E. Jansen
250*59599516SKenneth E. Jansenc
251*59599516SKenneth E. Jansenc.... calculate (Atau) <-- (A_1 tau) (Recall that we are using a
252*59599516SKenneth E. Jansenc                                     diagonal tau here)
253*59599516SKenneth E. Jansenc
254*59599516SKenneth E. Jansen       if (ivart .ge. 2) then
255*59599516SKenneth E. Jansen       do i = 1, nflow
256*59599516SKenneth E. Jansen          Atau(:,i,1) = A1(:,i,1)*tau(:,1)
257*59599516SKenneth E. Jansen          Atau(:,i,2) = A1(:,i,2)*tau(:,2)
258*59599516SKenneth E. Jansen          Atau(:,i,3) = A1(:,i,3)*tau(:,2)
259*59599516SKenneth E. Jansen          Atau(:,i,4) = A1(:,i,4)*tau(:,2)
260*59599516SKenneth E. Jansen          Atau(:,i,5) = A1(:,i,5)*tau(:,3)
261*59599516SKenneth E. Jansen       enddo
262*59599516SKenneth E. Jansenc
263*59599516SKenneth E. Jansenc.... calculate (A_1 tau A_0) (for L.S. time term of EGmass)
264*59599516SKenneth E. Jansenc
265*59599516SKenneth E. Jansen       do j = 1, nflow
266*59599516SKenneth E. Jansen          do i = 1, nflow
267*59599516SKenneth E. Jansen             A1tauA0(:,i,j) =
268*59599516SKenneth E. Jansen     &            Atau(:,i,1)*A0(:,1,j) +
269*59599516SKenneth E. Jansen     &            Atau(:,i,2)*A0(:,2,j) +
270*59599516SKenneth E. Jansen     &            Atau(:,i,3)*A0(:,3,j) +
271*59599516SKenneth E. Jansen     &            Atau(:,i,4)*A0(:,4,j) +
272*59599516SKenneth E. Jansen     &            Atau(:,i,5)*A0(:,5,j)
273*59599516SKenneth E. Jansen          enddo
274*59599516SKenneth E. Jansen       enddo
275*59599516SKenneth E. Jansenc
276*59599516SKenneth E. Jansenc.... calculate (Atau) <-- (A_2 tau) (Recall that we are using a
277*59599516SKenneth E. Jansenc                                     diagonal tau here)
278*59599516SKenneth E. Jansenc
279*59599516SKenneth E. Jansen       do i = 1, nflow
280*59599516SKenneth E. Jansen          Atau(:,i,1) = A2(:,i,1)*tau(:,1)
281*59599516SKenneth E. Jansen          Atau(:,i,2) = A2(:,i,2)*tau(:,2)
282*59599516SKenneth E. Jansen          Atau(:,i,3) = A2(:,i,3)*tau(:,2)
283*59599516SKenneth E. Jansen          Atau(:,i,4) = A2(:,i,4)*tau(:,2)
284*59599516SKenneth E. Jansen          Atau(:,i,5) = A2(:,i,5)*tau(:,3)
285*59599516SKenneth E. Jansen       enddo
286*59599516SKenneth E. Jansenc
287*59599516SKenneth E. Jansenc.... calculate (A_2 tau A_0) (for L.S. time term of EGmass)
288*59599516SKenneth E. Jansenc
289*59599516SKenneth E. Jansen       do j = 1, nflow
290*59599516SKenneth E. Jansen          do i = 1, nflow
291*59599516SKenneth E. Jansen             A2tauA0(:,i,j) =
292*59599516SKenneth E. Jansen     &            Atau(:,i,1)*A0(:,1,j) +
293*59599516SKenneth E. Jansen     &            Atau(:,i,2)*A0(:,2,j) +
294*59599516SKenneth E. Jansen     &            Atau(:,i,3)*A0(:,3,j) +
295*59599516SKenneth E. Jansen     &            Atau(:,i,4)*A0(:,4,j) +
296*59599516SKenneth E. Jansen     &            Atau(:,i,5)*A0(:,5,j)
297*59599516SKenneth E. Jansen          enddo
298*59599516SKenneth E. Jansen       enddo
299*59599516SKenneth E. Jansenc
300*59599516SKenneth E. Jansen       do i = 1, nflow
301*59599516SKenneth E. Jansen          Atau(:,i,1) = A3(:,i,1)*tau(:,1)
302*59599516SKenneth E. Jansen          Atau(:,i,2) = A3(:,i,2)*tau(:,2)
303*59599516SKenneth E. Jansen          Atau(:,i,3) = A3(:,i,3)*tau(:,2)
304*59599516SKenneth E. Jansen          Atau(:,i,4) = A3(:,i,4)*tau(:,2)
305*59599516SKenneth E. Jansen          Atau(:,i,5) = A3(:,i,5)*tau(:,3)
306*59599516SKenneth E. Jansen       enddo
307*59599516SKenneth E. Jansenc
308*59599516SKenneth E. Jansenc.... calculate (A_3 tau A_0) (for L.S. time term of EGmass)
309*59599516SKenneth E. Jansenc
310*59599516SKenneth E. Jansen       do j = 1, nflow
311*59599516SKenneth E. Jansen          do i = 1, nflow
312*59599516SKenneth E. Jansen             A3tauA0(:,i,j) =
313*59599516SKenneth E. Jansen     &            Atau(:,i,1)*A0(:,1,j) +
314*59599516SKenneth E. Jansen     &            Atau(:,i,2)*A0(:,2,j) +
315*59599516SKenneth E. Jansen     &            Atau(:,i,3)*A0(:,3,j) +
316*59599516SKenneth E. Jansen     &            Atau(:,i,4)*A0(:,4,j) +
317*59599516SKenneth E. Jansen     &            Atau(:,i,5)*A0(:,5,j)
318*59599516SKenneth E. Jansen          enddo
319*59599516SKenneth E. Jansen       enddo
320*59599516SKenneth E. Jansenc
321*59599516SKenneth E. Jansenc
322*59599516SKenneth E. Jansenc.... add least squares time term to BDiagl
323*59599516SKenneth E. Jansenc
324*59599516SKenneth E. Jansenc
325*59599516SKenneth E. Jansenc.... loop through rows (nodes i)
326*59599516SKenneth E. Jansenc
327*59599516SKenneth E. Jansen       do i = 1, nshl
328*59599516SKenneth E. Jansenc
329*59599516SKenneth E. Jansenc.... loop through column nodes, add (N_a,i A_i tau N_b) to BDiagl
330*59599516SKenneth E. Jansenc
331*59599516SKenneth E. Jansen          tmp1 = shp(:,i) * WdetJ * almi/gami/alfi*dtgl
332*59599516SKenneth E. Jansen          do idof = 1, nflow
333*59599516SKenneth E. Jansen             do jdof = 1, nflow
334*59599516SKenneth E. Jansen                BDiagl(:,i,idof,jdof) = BDiagl(:,i,idof,jdof)
335*59599516SKenneth E. Jansen     &                                + tmp1*(
336*59599516SKenneth E. Jansen     &               shg(:,i,1) * A1tauA0(:,idof,jdof) +
337*59599516SKenneth E. Jansen     &               shg(:,i,2) * A2tauA0(:,idof,jdof) +
338*59599516SKenneth E. Jansen     &               shg(:,i,3) * A3tauA0(:,idof,jdof))
339*59599516SKenneth E. Jansen             enddo
340*59599516SKenneth E. Jansen          enddo
341*59599516SKenneth E. Jansen       enddo
342*59599516SKenneth E. Jansen       endif
343*59599516SKenneth E. Jansenc
344*59599516SKenneth E. Jansenc.... add contribution of stiffness to BDiagl
345*59599516SKenneth E. Jansenc
346*59599516SKenneth E. Jansenc....  we no longer build stiff so we have to get it in here.
347*59599516SKenneth E. Jansenc      recall that this is the (A_i tau A_j + K_{ij}) that
348*59599516SKenneth E. Jansenc      multiplies N_{a,i} N_{b,j} (Actually for BDiag a=b).
349*59599516SKenneth E. Jansenc
350*59599516SKenneth E. Jansenc...   we are through with A0 so use it now for the sub-blocks
351*59599516SKenneth E. Jansenc      of what used to be stiff
352*59599516SKenneth E. Jansenc
353*59599516SKenneth E. Jansenc
354*59599516SKenneth E. Jansen
355*59599516SKenneth E. Jansenc
356*59599516SKenneth E. Jansenc  start with 1 1 contribution
357*59599516SKenneth E. Jansenc
358*59599516SKenneth E. Jansenc
359*59599516SKenneth E. Jansenc  Note that I comment out lines where zeros appear (for Pvariables)
360*59599516SKenneth E. Jansenc
361*59599516SKenneth E. Jansen          if(ivart .ge. 2) then
362*59599516SKenneth E. Jansenc
363*59599516SKenneth E. Jansenc  row one
364*59599516SKenneth E. Jansenc
365*59599516SKenneth E. Jansen          A0(:,1,1)=
366*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,1,1) * A1(:,1,1)
367*59599516SKenneth E. Jansen     &           + tau(:,2) * (
368*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A1(:,2,1)
369*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A1(:,3,1)
370*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A1(:,4,1)
371*59599516SKenneth E. Jansen     &                                                  )
372*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,1,5) * A1(:,5,1)
373*59599516SKenneth E. Jansenc
374*59599516SKenneth E. Jansen          A0(:,1,2)=
375*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,1,1) * A1(:,1,2)
376*59599516SKenneth E. Jansen     &           + tau(:,2) * (
377*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A1(:,2,2)
378*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A1(:,3,2)
379*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A1(:,4,2)
380*59599516SKenneth E. Jansen     &                                                  )
381*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,1,5) * A1(:,5,2)
382*59599516SKenneth E. Jansenc
383*59599516SKenneth E. Jansen          A0(:,1,3)=
384*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,1,1) * A1(:,1,3)
385*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
386*59599516SKenneth E. Jansenc    &                         + A1(:,1,2) * A1(:,2,3)
387*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A1(:,3,3)
388*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A1(:,4,3)
389*59599516SKenneth E. Jansenc    &                                                  )
390*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,1,5) * A1(:,5,3)
391*59599516SKenneth E. Jansenc
392*59599516SKenneth E. Jansen          A0(:,1,4)=
393*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,1,1) * A1(:,1,4)
394*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
395*59599516SKenneth E. Jansenc    &                         + A1(:,1,2) * A1(:,2,4)
396*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A1(:,3,4)
397*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A1(:,4,4)
398*59599516SKenneth E. Jansenc    &                                                  )
399*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,1,5) * A1(:,5,4)
400*59599516SKenneth E. Jansenc
401*59599516SKenneth E. Jansen          A0(:,1,5)=
402*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,1,1) * A1(:,1,5)
403*59599516SKenneth E. Jansen     &           + tau(:,2) * (
404*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A1(:,2,5)
405*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A1(:,3,5)
406*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A1(:,4,5)
407*59599516SKenneth E. Jansen     &                                                  )
408*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,1,5) * A1(:,5,5)
409*59599516SKenneth E. Jansenc
410*59599516SKenneth E. Jansenc  row two
411*59599516SKenneth E. Jansenc
412*59599516SKenneth E. Jansen          A0(:,2,1)=
413*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,2,1) * A1(:,1,1)
414*59599516SKenneth E. Jansen     &           + tau(:,2) * (
415*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A1(:,2,1)
416*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A1(:,3,1)
417*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A1(:,4,1)
418*59599516SKenneth E. Jansen     &                                                  )
419*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,2,5) * A1(:,5,1)
420*59599516SKenneth E. Jansenc
421*59599516SKenneth E. Jansen          A0(:,2,2)=
422*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,2,1) * A1(:,1,2)
423*59599516SKenneth E. Jansen     &           + tau(:,2) * (
424*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A1(:,2,2)
425*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A1(:,3,2)
426*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A1(:,4,2)
427*59599516SKenneth E. Jansen     &                                                  )
428*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,2,5) * A1(:,5,2)
429*59599516SKenneth E. Jansenc
430*59599516SKenneth E. Jansen          A0(:,2,3)=
431*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,2,1) * A1(:,1,3)
432*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
433*59599516SKenneth E. Jansenc    &                         + A1(:,2,2) * A1(:,2,3)
434*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A1(:,3,3)
435*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A1(:,4,3)
436*59599516SKenneth E. Jansenc    &                                                  )
437*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,2,5) * A1(:,5,3)
438*59599516SKenneth E. Jansenc
439*59599516SKenneth E. Jansen          A0(:,2,4)=
440*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,2,1) * A1(:,1,4)
441*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
442*59599516SKenneth E. Jansenc    &                         + A1(:,2,2) * A1(:,2,4)
443*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A1(:,3,4)
444*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A1(:,4,4)
445*59599516SKenneth E. Jansenc    &                                                  )
446*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,2,5) * A1(:,5,4)
447*59599516SKenneth E. Jansenc
448*59599516SKenneth E. Jansen          A0(:,2,5)=
449*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,2,1) * A1(:,1,5)
450*59599516SKenneth E. Jansen     &           + tau(:,2) * (
451*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A1(:,2,5)
452*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A1(:,3,5)
453*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A1(:,4,5)
454*59599516SKenneth E. Jansen     &                                                  )
455*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,2,5) * A1(:,5,5)
456*59599516SKenneth E. Jansenc
457*59599516SKenneth E. Jansenc  row three
458*59599516SKenneth E. Jansenc
459*59599516SKenneth E. Jansen          A0(:,3,1)=
460*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,3,1) * A1(:,1,1)
461*59599516SKenneth E. Jansen     &           + tau(:,2) * (
462*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A1(:,2,1)
463*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A1(:,3,1)
464*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A1(:,4,1)
465*59599516SKenneth E. Jansen     &                                                  )
466*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,3,5) * A1(:,5,1)
467*59599516SKenneth E. Jansenc
468*59599516SKenneth E. Jansen          A0(:,3,2)=
469*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,3,1) * A1(:,1,2)
470*59599516SKenneth E. Jansen     &           + tau(:,2) * (
471*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A1(:,2,2)
472*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A1(:,3,2)
473*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A1(:,4,2)
474*59599516SKenneth E. Jansen     &                                                  )
475*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,3,5) * A1(:,5,2)
476*59599516SKenneth E. Jansenc
477*59599516SKenneth E. Jansen          A0(:,3,3)=
478*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,3,1) * A1(:,1,3)
479*59599516SKenneth E. Jansen     &           + tau(:,2) * (
480*59599516SKenneth E. Jansenc    &                         + A1(:,3,2) * A1(:,2,3)
481*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A1(:,3,3)
482*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A1(:,4,3)
483*59599516SKenneth E. Jansen     &                                                  )
484*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,3,5) * A1(:,5,3)
485*59599516SKenneth E. Jansenc
486*59599516SKenneth E. Jansen          A0(:,3,4)=
487*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,3,1) * A1(:,1,4)
488*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
489*59599516SKenneth E. Jansenc    &                         + A1(:,3,2) * A1(:,2,4)
490*59599516SKenneth E. Jansenc    &                         + A1(:,3,3) * A1(:,3,4)
491*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A1(:,4,4)
492*59599516SKenneth E. Jansenc    &                                                  )
493*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,3,5) * A1(:,5,4)
494*59599516SKenneth E. Jansenc
495*59599516SKenneth E. Jansen          A0(:,3,5)=
496*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,3,1) * A1(:,1,5)
497*59599516SKenneth E. Jansen     &           + tau(:,2) * (
498*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A1(:,2,5)
499*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A1(:,3,5)
500*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A1(:,4,5)
501*59599516SKenneth E. Jansen     &                                                  )
502*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,3,5) * A1(:,5,5)
503*59599516SKenneth E. Jansenc
504*59599516SKenneth E. Jansenc  row four
505*59599516SKenneth E. Jansenc
506*59599516SKenneth E. Jansen          A0(:,4,1)=
507*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,4,1) * A1(:,1,1)
508*59599516SKenneth E. Jansen     &           + tau(:,2) * (
509*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A1(:,2,1)
510*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A1(:,3,1)
511*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A1(:,4,1)
512*59599516SKenneth E. Jansen     &                                                  )
513*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,4,5) * A1(:,5,1)
514*59599516SKenneth E. Jansenc
515*59599516SKenneth E. Jansen          A0(:,4,2)=
516*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,4,1) * A1(:,1,2)
517*59599516SKenneth E. Jansen     &           + tau(:,2) * (
518*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A1(:,2,2)
519*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A1(:,3,2)
520*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A1(:,4,2)
521*59599516SKenneth E. Jansen     &                                                  )
522*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,4,5) * A1(:,5,2)
523*59599516SKenneth E. Jansenc
524*59599516SKenneth E. Jansen          A0(:,4,3)=
525*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,4,1) * A1(:,1,3)
526*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
527*59599516SKenneth E. Jansenc    &                         + A1(:,4,2) * A1(:,2,3)
528*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A1(:,3,3)
529*59599516SKenneth E. Jansenc    &                         + A1(:,4,4) * A1(:,4,3)
530*59599516SKenneth E. Jansenc    &                                                  )
531*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,4,5) * A1(:,5,3)
532*59599516SKenneth E. Jansenc
533*59599516SKenneth E. Jansen          A0(:,4,4)=
534*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,4,1) * A1(:,1,4)
535*59599516SKenneth E. Jansen     &           + tau(:,2) * (
536*59599516SKenneth E. Jansenc    &                         + A1(:,4,2) * A1(:,2,4)
537*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A1(:,3,4)
538*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A1(:,4,4)
539*59599516SKenneth E. Jansen     &                                                  )
540*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,4,5) * A1(:,5,4)
541*59599516SKenneth E. Jansenc
542*59599516SKenneth E. Jansen          A0(:,4,5)=
543*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,4,1) * A1(:,1,5)
544*59599516SKenneth E. Jansen     &           + tau(:,2) * (
545*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A1(:,2,5)
546*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A1(:,3,5)
547*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A1(:,4,5)
548*59599516SKenneth E. Jansen     &                                                  )
549*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,4,5) * A1(:,5,5)
550*59599516SKenneth E. Jansenc
551*59599516SKenneth E. Jansenc  row two
552*59599516SKenneth E. Jansenc
553*59599516SKenneth E. Jansen          A0(:,5,1)=
554*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,5,1) * A1(:,1,1)
555*59599516SKenneth E. Jansen     &           + tau(:,2) * (
556*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A1(:,2,1)
557*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A1(:,3,1)
558*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A1(:,4,1)
559*59599516SKenneth E. Jansen     &                                                  )
560*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,5,5) * A1(:,5,1)
561*59599516SKenneth E. Jansenc
562*59599516SKenneth E. Jansen          A0(:,5,2)=
563*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,5,1) * A1(:,1,2)
564*59599516SKenneth E. Jansen     &           + tau(:,2) * (
565*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A1(:,2,2)
566*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A1(:,3,2)
567*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A1(:,4,2)
568*59599516SKenneth E. Jansen     &                                                  )
569*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,5,5) * A1(:,5,2)
570*59599516SKenneth E. Jansenc
571*59599516SKenneth E. Jansen          A0(:,5,3)=
572*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,5,1) * A1(:,1,3)
573*59599516SKenneth E. Jansen     &           + tau(:,2) * (
574*59599516SKenneth E. Jansenc    &                         + A1(:,5,2) * A1(:,2,3)
575*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A1(:,3,3)
576*59599516SKenneth E. Jansenc    &                         + A1(:,5,4) * A1(:,4,3)
577*59599516SKenneth E. Jansen     &                                                  )
578*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,5,5) * A1(:,5,3)
579*59599516SKenneth E. Jansenc
580*59599516SKenneth E. Jansen          A0(:,5,4)=
581*59599516SKenneth E. Jansenc    &             tau(:,1) *    A1(:,5,1) * A1(:,1,4)
582*59599516SKenneth E. Jansen     &           + tau(:,2) * (
583*59599516SKenneth E. Jansenc    &                         + A1(:,5,2) * A1(:,2,4)
584*59599516SKenneth E. Jansenc    &                         + A1(:,5,3) * A1(:,3,4)
585*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A1(:,4,4)
586*59599516SKenneth E. Jansen     &                                                  )
587*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,5,5) * A1(:,5,4)
588*59599516SKenneth E. Jansenc
589*59599516SKenneth E. Jansen          A0(:,5,5)=
590*59599516SKenneth E. Jansen     &             tau(:,1) *    A1(:,5,1) * A1(:,1,5)
591*59599516SKenneth E. Jansen     &           + tau(:,2) * (
592*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A1(:,2,5)
593*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A1(:,3,5)
594*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A1(:,4,5)
595*59599516SKenneth E. Jansen     &                                                  )
596*59599516SKenneth E. Jansen     &           + tau(:,3) *    A1(:,5,5) * A1(:,5,5)
597*59599516SKenneth E. Jansenc
598*59599516SKenneth E. Jansen         else
599*59599516SKenneth E. Jansen          A0=zero
600*59599516SKenneth E. Jansen         endif
601*59599516SKenneth E. Jansenc
602*59599516SKenneth E. Jansen         if (Navier .eq. 1) then
603*59599516SKenneth E. Jansen           A0(:,2,2) = A0(:,2,2) + rlm2mu !rlm2mu
604*59599516SKenneth E. Jansen           A0(:,3,3) = A0(:,3,3) + rmu !rmu
605*59599516SKenneth E. Jansen           A0(:,4,4) = A0(:,4,4) + rmu !rmu
606*59599516SKenneth E. Jansen           A0(:,5,2) = A0(:,5,2) + rlm2mu * u1
607*59599516SKenneth E. Jansen           A0(:,5,3) = A0(:,5,3) + rmu * u2
608*59599516SKenneth E. Jansen           A0(:,5,4) = A0(:,5,4) + rmu * u3
609*59599516SKenneth E. Jansen           A0(:,5,5) = A0(:,5,5) + con !con
610*59599516SKenneth E. Jansen        endif
611*59599516SKenneth E. Jansen        do j = 1, nshl
612*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,1) * shg(:,j,1)
613*59599516SKenneth E. Jansen           do i=1,nflow
614*59599516SKenneth E. Jansen           do k=1,nflow
615*59599516SKenneth E. Jansen             BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k)
616*59599516SKenneth E. Jansen           enddo
617*59599516SKenneth E. Jansen           enddo
618*59599516SKenneth E. Jansen        enddo
619*59599516SKenneth E. Jansenc
620*59599516SKenneth E. Jansenc  start with 2 2 contribution
621*59599516SKenneth E. Jansenc
622*59599516SKenneth E. Jansen          if(ivart.ge.2) then
623*59599516SKenneth E. Jansenc
624*59599516SKenneth E. Jansenc  row one
625*59599516SKenneth E. Jansenc
626*59599516SKenneth E. Jansen          A0(:,1,1)=
627*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,1,1) * A2(:,1,1)
628*59599516SKenneth E. Jansen     &           + tau(:,2) * (
629*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A2(:,2,1)
630*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A2(:,3,1)
631*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A2(:,4,1)
632*59599516SKenneth E. Jansen     &                                                  )
633*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,1,5) * A2(:,5,1)
634*59599516SKenneth E. Jansenc
635*59599516SKenneth E. Jansen          A0(:,1,2)=
636*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,1,1) * A2(:,1,2)
637*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
638*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A2(:,2,2)
639*59599516SKenneth E. Jansenc    &                         + A2(:,1,3) * A2(:,3,2)
640*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A2(:,4,2)
641*59599516SKenneth E. Jansenc    &                                                  )
642*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,1,5) * A2(:,5,2)
643*59599516SKenneth E. Jansenc
644*59599516SKenneth E. Jansen          A0(:,1,3)=
645*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,1,1) * A2(:,1,3)
646*59599516SKenneth E. Jansen     &           + tau(:,2) * (
647*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A2(:,2,3)
648*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A2(:,3,3)
649*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A2(:,4,3)
650*59599516SKenneth E. Jansen     &                                                  )
651*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,1,5) * A2(:,5,3)
652*59599516SKenneth E. Jansenc
653*59599516SKenneth E. Jansen          A0(:,1,4)=
654*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,1,1) * A2(:,1,4)
655*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
656*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A2(:,2,4)
657*59599516SKenneth E. Jansenc    &                         + A2(:,1,3) * A2(:,3,4)
658*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A2(:,4,4)
659*59599516SKenneth E. Jansenc    &                                                  )
660*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,1,5) * A2(:,5,4)
661*59599516SKenneth E. Jansenc
662*59599516SKenneth E. Jansen          A0(:,1,5)=
663*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,1,1) * A2(:,1,5)
664*59599516SKenneth E. Jansen     &           + tau(:,2) * (
665*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A2(:,2,5)
666*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A2(:,3,5)
667*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A2(:,4,5)
668*59599516SKenneth E. Jansen     &                                                  )
669*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,1,5) * A2(:,5,5)
670*59599516SKenneth E. Jansenc
671*59599516SKenneth E. Jansenc  row two
672*59599516SKenneth E. Jansenc
673*59599516SKenneth E. Jansen          A0(:,2,1)=
674*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,2,1) * A2(:,1,1)
675*59599516SKenneth E. Jansen     &           + tau(:,2) * (
676*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A2(:,2,1)
677*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A2(:,3,1)
678*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A2(:,4,1)
679*59599516SKenneth E. Jansen     &                                                  )
680*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,2,5) * A2(:,5,1)
681*59599516SKenneth E. Jansenc
682*59599516SKenneth E. Jansen          A0(:,2,2)=
683*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,2,1) * A2(:,1,2)
684*59599516SKenneth E. Jansen     &           + tau(:,2) * (
685*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A2(:,2,2)
686*59599516SKenneth E. Jansenc    &                         + A2(:,2,3) * A2(:,3,2)
687*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A2(:,4,2)
688*59599516SKenneth E. Jansen     &                                                  )
689*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,2,5) * A2(:,5,2)
690*59599516SKenneth E. Jansenc
691*59599516SKenneth E. Jansen          A0(:,2,3)=
692*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,2,1) * A2(:,1,3)
693*59599516SKenneth E. Jansen     &           + tau(:,2) * (
694*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A2(:,2,3)
695*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A2(:,3,3)
696*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A2(:,4,3)
697*59599516SKenneth E. Jansen     &                                                  )
698*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,2,5) * A2(:,5,3)
699*59599516SKenneth E. Jansenc
700*59599516SKenneth E. Jansen          A0(:,2,4)=
701*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,2,1) * A2(:,1,4)
702*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
703*59599516SKenneth E. Jansenc    &                         + A2(:,2,2) * A2(:,2,4)
704*59599516SKenneth E. Jansenc    &                         + A2(:,2,3) * A2(:,3,4)
705*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A2(:,4,4)
706*59599516SKenneth E. Jansenc    &                                                  )
707*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,2,5) * A2(:,5,4)
708*59599516SKenneth E. Jansenc
709*59599516SKenneth E. Jansen          A0(:,2,5)=
710*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,2,1) * A2(:,1,5)
711*59599516SKenneth E. Jansen     &           + tau(:,2) * (
712*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A2(:,2,5)
713*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A2(:,3,5)
714*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A2(:,4,5)
715*59599516SKenneth E. Jansen     &                                                  )
716*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,2,5) * A2(:,5,5)
717*59599516SKenneth E. Jansenc
718*59599516SKenneth E. Jansenc  row three
719*59599516SKenneth E. Jansenc
720*59599516SKenneth E. Jansen          A0(:,3,1)=
721*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,3,1) * A2(:,1,1)
722*59599516SKenneth E. Jansen     &           + tau(:,2) * (
723*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A2(:,2,1)
724*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A2(:,3,1)
725*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A2(:,4,1)
726*59599516SKenneth E. Jansen     &                                                  )
727*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,3,5) * A2(:,5,1)
728*59599516SKenneth E. Jansenc
729*59599516SKenneth E. Jansen          A0(:,3,2)=
730*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,3,1) * A2(:,1,2)
731*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
732*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A2(:,2,2)
733*59599516SKenneth E. Jansenc    &                         + A2(:,3,3) * A2(:,3,2)
734*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A2(:,4,2)
735*59599516SKenneth E. Jansenc    &                                                  )
736*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,3,5) * A2(:,5,2)
737*59599516SKenneth E. Jansenc
738*59599516SKenneth E. Jansen          A0(:,3,3)=
739*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,3,1) * A2(:,1,3)
740*59599516SKenneth E. Jansen     &           + tau(:,2) * (
741*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A2(:,2,3)
742*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A2(:,3,3)
743*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A2(:,4,3)
744*59599516SKenneth E. Jansen     &                                                  )
745*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,3,5) * A2(:,5,3)
746*59599516SKenneth E. Jansenc
747*59599516SKenneth E. Jansen          A0(:,3,4)=
748*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,3,1) * A2(:,1,4)
749*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
750*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A2(:,2,4)
751*59599516SKenneth E. Jansenc    &                         + A2(:,3,3) * A2(:,3,4)
752*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A2(:,4,4)
753*59599516SKenneth E. Jansenc    &                                                  )
754*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,3,5) * A2(:,5,4)
755*59599516SKenneth E. Jansenc
756*59599516SKenneth E. Jansen          A0(:,3,5)=
757*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,3,1) * A2(:,1,5)
758*59599516SKenneth E. Jansen     &           + tau(:,2) * (
759*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A2(:,2,5)
760*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A2(:,3,5)
761*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A2(:,4,5)
762*59599516SKenneth E. Jansen     &                                                  )
763*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,3,5) * A2(:,5,5)
764*59599516SKenneth E. Jansenc
765*59599516SKenneth E. Jansenc  row four
766*59599516SKenneth E. Jansenc
767*59599516SKenneth E. Jansen          A0(:,4,1)=
768*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,4,1) * A2(:,1,1)
769*59599516SKenneth E. Jansen     &           + tau(:,2) * (
770*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A2(:,2,1)
771*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A2(:,3,1)
772*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A2(:,4,1)
773*59599516SKenneth E. Jansen     &                                                  )
774*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,4,5) * A2(:,5,1)
775*59599516SKenneth E. Jansenc
776*59599516SKenneth E. Jansen          A0(:,4,2)=
777*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,4,1) * A2(:,1,2)
778*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
779*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A2(:,2,2)
780*59599516SKenneth E. Jansenc    &                         + A2(:,4,3) * A2(:,3,2)
781*59599516SKenneth E. Jansenc    &                         + A2(:,4,4) * A2(:,4,2)
782*59599516SKenneth E. Jansenc    &                                                  )
783*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,4,5) * A2(:,5,2)
784*59599516SKenneth E. Jansenc
785*59599516SKenneth E. Jansen          A0(:,4,3)=
786*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,4,1) * A2(:,1,3)
787*59599516SKenneth E. Jansen     &           + tau(:,2) * (
788*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A2(:,2,3)
789*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A2(:,3,3)
790*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A2(:,4,3)
791*59599516SKenneth E. Jansen     &                                                  )
792*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,4,5) * A2(:,5,3)
793*59599516SKenneth E. Jansenc
794*59599516SKenneth E. Jansen          A0(:,4,4)=
795*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,4,1) * A2(:,1,4)
796*59599516SKenneth E. Jansen     &           + tau(:,2) * (
797*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A2(:,2,4)
798*59599516SKenneth E. Jansenc    &                         + A2(:,4,3) * A2(:,3,4)
799*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A2(:,4,4)
800*59599516SKenneth E. Jansen     &                                                  )
801*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,4,5) * A2(:,5,4)
802*59599516SKenneth E. Jansenc
803*59599516SKenneth E. Jansen          A0(:,4,5)=
804*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,4,1) * A2(:,1,5)
805*59599516SKenneth E. Jansen     &           + tau(:,2) * (
806*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A2(:,2,5)
807*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A2(:,3,5)
808*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A2(:,4,5)
809*59599516SKenneth E. Jansen     &                                                  )
810*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,4,5) * A2(:,5,5)
811*59599516SKenneth E. Jansenc
812*59599516SKenneth E. Jansenc  row two
813*59599516SKenneth E. Jansenc
814*59599516SKenneth E. Jansen          A0(:,5,1)=
815*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,5,1) * A2(:,1,1)
816*59599516SKenneth E. Jansen     &           + tau(:,2) * (
817*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A2(:,2,1)
818*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A2(:,3,1)
819*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A2(:,4,1)
820*59599516SKenneth E. Jansen     &                                                  )
821*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,5,5) * A2(:,5,1)
822*59599516SKenneth E. Jansenc
823*59599516SKenneth E. Jansen          A0(:,5,2)=
824*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,5,1) * A2(:,1,2)
825*59599516SKenneth E. Jansen     &           + tau(:,2) * (
826*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A2(:,2,2)
827*59599516SKenneth E. Jansenc    &                         + A2(:,5,3) * A2(:,3,2)
828*59599516SKenneth E. Jansenc    &                         + A2(:,5,4) * A2(:,4,2)
829*59599516SKenneth E. Jansen     &                                                  )
830*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,5,5) * A2(:,5,2)
831*59599516SKenneth E. Jansenc
832*59599516SKenneth E. Jansen          A0(:,5,3)=
833*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,5,1) * A2(:,1,3)
834*59599516SKenneth E. Jansen     &           + tau(:,2) * (
835*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A2(:,2,3)
836*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A2(:,3,3)
837*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A2(:,4,3)
838*59599516SKenneth E. Jansen     &                                                  )
839*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,5,5) * A2(:,5,3)
840*59599516SKenneth E. Jansenc
841*59599516SKenneth E. Jansen          A0(:,5,4)=
842*59599516SKenneth E. Jansenc    &             tau(:,1) *    A2(:,5,1) * A2(:,1,4)
843*59599516SKenneth E. Jansen     &           + tau(:,2) * (
844*59599516SKenneth E. Jansenc    &                         + A2(:,5,2) * A2(:,2,4)
845*59599516SKenneth E. Jansenc    &                         + A2(:,5,3) * A2(:,3,4)
846*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A2(:,4,4)
847*59599516SKenneth E. Jansen     &                                                  )
848*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,5,5) * A2(:,5,4)
849*59599516SKenneth E. Jansenc
850*59599516SKenneth E. Jansen          A0(:,5,5)=
851*59599516SKenneth E. Jansen     &             tau(:,1) *    A2(:,5,1) * A2(:,1,5)
852*59599516SKenneth E. Jansen     &           + tau(:,2) * (
853*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A2(:,2,5)
854*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A2(:,3,5)
855*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A2(:,4,5)
856*59599516SKenneth E. Jansen     &                                                  )
857*59599516SKenneth E. Jansen     &           + tau(:,3) *    A2(:,5,5) * A2(:,5,5)
858*59599516SKenneth E. Jansen          else
859*59599516SKenneth E. Jansen             A0 = zero
860*59599516SKenneth E. Jansen          endif
861*59599516SKenneth E. Jansenc
862*59599516SKenneth E. Jansen          if (Navier .eq. 1) then
863*59599516SKenneth E. Jansen           A0(:,2,2) = A0(:,2,2) + rmu !rmu
864*59599516SKenneth E. Jansen           A0(:,3,3) = A0(:,3,3) + rlm2mu !rlm2mu
865*59599516SKenneth E. Jansen           A0(:,4,4) = A0(:,4,4) + rmu !rmu
866*59599516SKenneth E. Jansen           A0(:,5,2) = A0(:,5,2) + rmu * u1
867*59599516SKenneth E. Jansen           A0(:,5,3) = A0(:,5,3) + rlm2mu * u2
868*59599516SKenneth E. Jansen           A0(:,5,4) = A0(:,5,4) + rmu * u3
869*59599516SKenneth E. Jansen           A0(:,5,5) = A0(:,5,5) + con !con
870*59599516SKenneth E. Jansen        endif
871*59599516SKenneth E. Jansenc
872*59599516SKenneth E. Jansen        do j = 1, nshl
873*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,2) * shg(:,j,2)
874*59599516SKenneth E. Jansen           do i=1,nflow
875*59599516SKenneth E. Jansen           do k=1,nflow
876*59599516SKenneth E. Jansen             BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k)
877*59599516SKenneth E. Jansen           enddo
878*59599516SKenneth E. Jansen           enddo
879*59599516SKenneth E. Jansen        enddo
880*59599516SKenneth E. Jansenc
881*59599516SKenneth E. Jansenc  start with 3 3 contribution
882*59599516SKenneth E. Jansenc
883*59599516SKenneth E. Jansen       if(ivart.ge.2) then
884*59599516SKenneth E. Jansenc
885*59599516SKenneth E. Jansenc  row one
886*59599516SKenneth E. Jansenc
887*59599516SKenneth E. Jansen          A0(:,1,1)=
888*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,1,1) * A3(:,1,1)
889*59599516SKenneth E. Jansen     &           + tau(:,2) * (
890*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A3(:,2,1)
891*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A3(:,3,1)
892*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A3(:,4,1)
893*59599516SKenneth E. Jansen     &                                                  )
894*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,1,5) * A3(:,5,1)
895*59599516SKenneth E. Jansenc
896*59599516SKenneth E. Jansen          A0(:,1,2)=
897*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,1,1) * A3(:,1,2)
898*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
899*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A3(:,2,2)
900*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A3(:,3,2)
901*59599516SKenneth E. Jansenc    &                         + A3(:,1,4) * A3(:,4,2)
902*59599516SKenneth E. Jansenc    &                                                  )
903*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,1,5) * A3(:,5,2)
904*59599516SKenneth E. Jansenc
905*59599516SKenneth E. Jansen          A0(:,1,3)=
906*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,1,1) * A3(:,1,3)
907*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
908*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A3(:,2,3)
909*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A3(:,3,3)
910*59599516SKenneth E. Jansenc    &                         + A3(:,1,4) * A3(:,4,3)
911*59599516SKenneth E. Jansenc    &                                                  )
912*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,1,5) * A3(:,5,3)
913*59599516SKenneth E. Jansenc
914*59599516SKenneth E. Jansen          A0(:,1,4)=
915*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,1,1) * A3(:,1,4)
916*59599516SKenneth E. Jansen     &           + tau(:,2) * (
917*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A3(:,2,4)
918*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A3(:,3,4)
919*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A3(:,4,4)
920*59599516SKenneth E. Jansen     &                                                  )
921*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,1,5) * A3(:,5,4)
922*59599516SKenneth E. Jansenc
923*59599516SKenneth E. Jansen          A0(:,1,5)=
924*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,1,1) * A3(:,1,5)
925*59599516SKenneth E. Jansen     &           + tau(:,2) * (
926*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A3(:,2,5)
927*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A3(:,3,5)
928*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A3(:,4,5)
929*59599516SKenneth E. Jansen     &                                                  )
930*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,1,5) * A3(:,5,5)
931*59599516SKenneth E. Jansenc
932*59599516SKenneth E. Jansenc  row two
933*59599516SKenneth E. Jansenc
934*59599516SKenneth E. Jansen          A0(:,2,1)=
935*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,2,1) * A3(:,1,1)
936*59599516SKenneth E. Jansen     &           + tau(:,2) * (
937*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A3(:,2,1)
938*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A3(:,3,1)
939*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A3(:,4,1)
940*59599516SKenneth E. Jansen     &                                                  )
941*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,2,5) * A3(:,5,1)
942*59599516SKenneth E. Jansenc
943*59599516SKenneth E. Jansen          A0(:,2,2)=
944*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,2,1) * A3(:,1,2)
945*59599516SKenneth E. Jansen     &           + tau(:,2) * (
946*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A3(:,2,2)
947*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A3(:,3,2)
948*59599516SKenneth E. Jansenc    &                         + A3(:,2,4) * A3(:,4,2)
949*59599516SKenneth E. Jansen     &                                                  )
950*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,2,5) * A3(:,5,2)
951*59599516SKenneth E. Jansenc
952*59599516SKenneth E. Jansen          A0(:,2,3)=
953*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,2,1) * A3(:,1,3)
954*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
955*59599516SKenneth E. Jansenc    &                         + A3(:,2,2) * A3(:,2,3)
956*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A3(:,3,3)
957*59599516SKenneth E. Jansenc    &                         + A3(:,2,4) * A3(:,4,3)
958*59599516SKenneth E. Jansenc    &                                                  )
959*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,2,5) * A3(:,5,3)
960*59599516SKenneth E. Jansenc
961*59599516SKenneth E. Jansen          A0(:,2,4)=
962*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,2,1) * A3(:,1,4)
963*59599516SKenneth E. Jansen     &           + tau(:,2) * (
964*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A3(:,2,4)
965*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A3(:,3,4)
966*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A3(:,4,4)
967*59599516SKenneth E. Jansen     &                                                  )
968*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,2,5) * A3(:,5,4)
969*59599516SKenneth E. Jansenc
970*59599516SKenneth E. Jansen          A0(:,2,5)=
971*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,2,1) * A3(:,1,5)
972*59599516SKenneth E. Jansen     &           + tau(:,2) * (
973*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A3(:,2,5)
974*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A3(:,3,5)
975*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A3(:,4,5)
976*59599516SKenneth E. Jansen     &                                                  )
977*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,2,5) * A3(:,5,5)
978*59599516SKenneth E. Jansenc
979*59599516SKenneth E. Jansenc  row three
980*59599516SKenneth E. Jansenc
981*59599516SKenneth E. Jansen          A0(:,3,1)=
982*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,3,1) * A3(:,1,1)
983*59599516SKenneth E. Jansen     &           + tau(:,2) * (
984*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A3(:,2,1)
985*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A3(:,3,1)
986*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A3(:,4,1)
987*59599516SKenneth E. Jansen     &                                                  )
988*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,3,5) * A3(:,5,1)
989*59599516SKenneth E. Jansenc
990*59599516SKenneth E. Jansen          A0(:,3,2)=
991*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,3,1) * A3(:,1,2)
992*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
993*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A3(:,2,2)
994*59599516SKenneth E. Jansenc    &                         + A3(:,3,3) * A3(:,3,2)
995*59599516SKenneth E. Jansenc    &                         + A3(:,3,4) * A3(:,4,2)
996*59599516SKenneth E. Jansenc    &                                                  )
997*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,3,5) * A3(:,5,2)
998*59599516SKenneth E. Jansenc
999*59599516SKenneth E. Jansen          A0(:,3,3)=
1000*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,3,1) * A3(:,1,3)
1001*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1002*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A3(:,2,3)
1003*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A3(:,3,3)
1004*59599516SKenneth E. Jansenc    &                         + A3(:,3,4) * A3(:,4,3)
1005*59599516SKenneth E. Jansen     &                                                  )
1006*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,3,5) * A3(:,5,3)
1007*59599516SKenneth E. Jansenc
1008*59599516SKenneth E. Jansen          A0(:,3,4)=
1009*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,3,1) * A3(:,1,4)
1010*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1011*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A3(:,2,4)
1012*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A3(:,3,4)
1013*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A3(:,4,4)
1014*59599516SKenneth E. Jansen     &                                                  )
1015*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,3,5) * A3(:,5,4)
1016*59599516SKenneth E. Jansenc
1017*59599516SKenneth E. Jansen          A0(:,3,5)=
1018*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,3,1) * A3(:,1,5)
1019*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1020*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A3(:,2,5)
1021*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A3(:,3,5)
1022*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A3(:,4,5)
1023*59599516SKenneth E. Jansen     &                                                  )
1024*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,3,5) * A3(:,5,5)
1025*59599516SKenneth E. Jansenc
1026*59599516SKenneth E. Jansenc  row four
1027*59599516SKenneth E. Jansenc
1028*59599516SKenneth E. Jansen          A0(:,4,1)=
1029*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,4,1) * A3(:,1,1)
1030*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1031*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A3(:,2,1)
1032*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A3(:,3,1)
1033*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A3(:,4,1)
1034*59599516SKenneth E. Jansen     &                                                  )
1035*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,4,5) * A3(:,5,1)
1036*59599516SKenneth E. Jansenc
1037*59599516SKenneth E. Jansen          A0(:,4,2)=
1038*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,4,1) * A3(:,1,2)
1039*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
1040*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A3(:,2,2)
1041*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A3(:,3,2)
1042*59599516SKenneth E. Jansenc    &                         + A3(:,4,4) * A3(:,4,2)
1043*59599516SKenneth E. Jansenc    &                                                  )
1044*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,4,5) * A3(:,5,2)
1045*59599516SKenneth E. Jansenc
1046*59599516SKenneth E. Jansen          A0(:,4,3)=
1047*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,4,1) * A3(:,1,3)
1048*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
1049*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A3(:,2,3)
1050*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A3(:,3,3)
1051*59599516SKenneth E. Jansenc    &                         + A3(:,4,4) * A3(:,4,3)
1052*59599516SKenneth E. Jansenc    &                                                  )
1053*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,4,5) * A3(:,5,3)
1054*59599516SKenneth E. Jansenc
1055*59599516SKenneth E. Jansen          A0(:,4,4)=
1056*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,4,1) * A3(:,1,4)
1057*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1058*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A3(:,2,4)
1059*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A3(:,3,4)
1060*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A3(:,4,4)
1061*59599516SKenneth E. Jansen     &                                                  )
1062*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,4,5) * A3(:,5,4)
1063*59599516SKenneth E. Jansenc
1064*59599516SKenneth E. Jansen          A0(:,4,5)=
1065*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,4,1) * A3(:,1,5)
1066*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1067*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A3(:,2,5)
1068*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A3(:,3,5)
1069*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A3(:,4,5)
1070*59599516SKenneth E. Jansen     &                                                  )
1071*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,4,5) * A3(:,5,5)
1072*59599516SKenneth E. Jansenc
1073*59599516SKenneth E. Jansenc  row two
1074*59599516SKenneth E. Jansenc
1075*59599516SKenneth E. Jansen          A0(:,5,1)=
1076*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,5,1) * A3(:,1,1)
1077*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1078*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A3(:,2,1)
1079*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A3(:,3,1)
1080*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A3(:,4,1)
1081*59599516SKenneth E. Jansen     &                                                  )
1082*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,5,5) * A3(:,5,1)
1083*59599516SKenneth E. Jansenc
1084*59599516SKenneth E. Jansen          A0(:,5,2)=
1085*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,5,1) * A3(:,1,2)
1086*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1087*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A3(:,2,2)
1088*59599516SKenneth E. Jansenc    &                         + A3(:,5,3) * A3(:,3,2)
1089*59599516SKenneth E. Jansenc    &                         + A3(:,5,4) * A3(:,4,2)
1090*59599516SKenneth E. Jansen     &                                                  )
1091*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,5,5) * A3(:,5,2)
1092*59599516SKenneth E. Jansenc
1093*59599516SKenneth E. Jansen          A0(:,5,3)=
1094*59599516SKenneth E. Jansenc    &             tau(:,1) *    A3(:,5,1) * A3(:,1,3)
1095*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1096*59599516SKenneth E. Jansenc    &                         + A3(:,5,2) * A3(:,2,3)
1097*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A3(:,3,3)
1098*59599516SKenneth E. Jansenc    &                         + A3(:,5,4) * A3(:,4,3)
1099*59599516SKenneth E. Jansen     &                                                  )
1100*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,5,5) * A3(:,5,3)
1101*59599516SKenneth E. Jansenc
1102*59599516SKenneth E. Jansen          A0(:,5,4)=
1103*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,5,1) * A3(:,1,4)
1104*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1105*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A3(:,2,4)
1106*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A3(:,3,4)
1107*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A3(:,4,4)
1108*59599516SKenneth E. Jansen     &                                                  )
1109*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,5,5) * A3(:,5,4)
1110*59599516SKenneth E. Jansenc
1111*59599516SKenneth E. Jansen          A0(:,5,5)=
1112*59599516SKenneth E. Jansen     &             tau(:,1) *    A3(:,5,1) * A3(:,1,5)
1113*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1114*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A3(:,2,5)
1115*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A3(:,3,5)
1116*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A3(:,4,5)
1117*59599516SKenneth E. Jansen     &                                                  )
1118*59599516SKenneth E. Jansen     &           + tau(:,3) *    A3(:,5,5) * A3(:,5,5)
1119*59599516SKenneth E. Jansen          else
1120*59599516SKenneth E. Jansen             A0 = zero
1121*59599516SKenneth E. Jansen          endif
1122*59599516SKenneth E. Jansenc
1123*59599516SKenneth E. Jansen          if (Navier .eq. 1) then
1124*59599516SKenneth E. Jansen           A0(:,2,2) = A0(:,2,2) + rmu !rmu
1125*59599516SKenneth E. Jansen           A0(:,3,3) = A0(:,3,3) + rmu !rmu
1126*59599516SKenneth E. Jansen           A0(:,4,4) = A0(:,4,4) + rlm2mu !rlm2mu
1127*59599516SKenneth E. Jansen           A0(:,5,2) = A0(:,5,2) + rmu * u1
1128*59599516SKenneth E. Jansen           A0(:,5,3) = A0(:,5,3) + rmu * u2
1129*59599516SKenneth E. Jansen           A0(:,5,4) = A0(:,5,4) + rlm2mu * u3
1130*59599516SKenneth E. Jansen           A0(:,5,5) = A0(:,5,5) + con !con
1131*59599516SKenneth E. Jansen          endif
1132*59599516SKenneth E. Jansenc
1133*59599516SKenneth E. Jansen        do j = 1, nshl
1134*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,3) * shg(:,j,3)
1135*59599516SKenneth E. Jansen           do i=1,nflow
1136*59599516SKenneth E. Jansen           do k=1,nflow
1137*59599516SKenneth E. Jansen             BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k)
1138*59599516SKenneth E. Jansen           enddo
1139*59599516SKenneth E. Jansen           enddo
1140*59599516SKenneth E. Jansen        enddo
1141*59599516SKenneth E. Jansenc
1142*59599516SKenneth E. Jansenc now for the 1 2  plus 2 1
1143*59599516SKenneth E. Jansenc
1144*59599516SKenneth E. Jansen          if(ivart.ge.2) then
1145*59599516SKenneth E. Jansenc
1146*59599516SKenneth E. Jansenc  row one
1147*59599516SKenneth E. Jansenc
1148*59599516SKenneth E. Jansen          A0(:,1,1)=
1149*59599516SKenneth E. Jansen     &             tau(:,1) * (
1150*59599516SKenneth E. Jansen     &                         + A1(:,1,1) * A2(:,1,1)
1151*59599516SKenneth E. Jansen     &                         + A2(:,1,1) * A1(:,1,1)
1152*59599516SKenneth E. Jansen     &		                                        )
1153*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1154*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A2(:,2,1)
1155*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A2(:,3,1)
1156*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A2(:,4,1)
1157*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A1(:,2,1)
1158*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A1(:,3,1)
1159*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A1(:,4,1)
1160*59599516SKenneth E. Jansen     &                                                  )
1161*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1162*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A2(:,5,1)
1163*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A1(:,5,1)
1164*59599516SKenneth E. Jansen     &		                                        )
1165*59599516SKenneth E. Jansenc
1166*59599516SKenneth E. Jansen          A0(:,1,2)=
1167*59599516SKenneth E. Jansen     &             tau(:,1) * (
1168*59599516SKenneth E. Jansenc    &                         + A1(:,1,1) * A2(:,1,2)
1169*59599516SKenneth E. Jansen     &                         + A2(:,1,1) * A1(:,1,2)
1170*59599516SKenneth E. Jansen     &		                                        )
1171*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1172*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A2(:,2,2)
1173*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A2(:,3,2)
1174*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A2(:,4,2)
1175*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A1(:,2,2)
1176*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A1(:,3,2)
1177*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A1(:,4,2)
1178*59599516SKenneth E. Jansen     &                                                  )
1179*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1180*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A2(:,5,2)
1181*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A1(:,5,2)
1182*59599516SKenneth E. Jansen     &		                                        )
1183*59599516SKenneth E. Jansenc
1184*59599516SKenneth E. Jansen          A0(:,1,3)=
1185*59599516SKenneth E. Jansen     &             tau(:,1) * (
1186*59599516SKenneth E. Jansen     &                         + A1(:,1,1) * A2(:,1,3)
1187*59599516SKenneth E. Jansenc    &                         + A2(:,1,1) * A1(:,1,3)
1188*59599516SKenneth E. Jansen     &		                                        )
1189*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1190*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A2(:,2,3)
1191*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A2(:,3,3)
1192*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A2(:,4,3)
1193*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A1(:,2,3)
1194*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A1(:,3,3)
1195*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A1(:,4,3)
1196*59599516SKenneth E. Jansen     &                                                  )
1197*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1198*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A2(:,5,3)
1199*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A1(:,5,3)
1200*59599516SKenneth E. Jansen     &		                                        )
1201*59599516SKenneth E. Jansenc
1202*59599516SKenneth E. Jansen          A0(:,1,4)=
1203*59599516SKenneth E. Jansenc    &             tau(:,1) * (
1204*59599516SKenneth E. Jansenc    &                         + A1(:,1,1) * A2(:,1,4)
1205*59599516SKenneth E. Jansenc    &                         + A2(:,1,1) * A1(:,1,4)
1206*59599516SKenneth E. Jansenc    &		                                        )
1207*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
1208*59599516SKenneth E. Jansenc    &                         + A1(:,1,2) * A2(:,2,4)
1209*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A2(:,3,4)
1210*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A2(:,4,4)
1211*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A1(:,2,4)
1212*59599516SKenneth E. Jansenc    &                         + A2(:,1,3) * A1(:,3,4)
1213*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A1(:,4,4)
1214*59599516SKenneth E. Jansenc    &                                                  )
1215*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1216*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A2(:,5,4)
1217*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A1(:,5,4)
1218*59599516SKenneth E. Jansen     &		                                        )
1219*59599516SKenneth E. Jansenc
1220*59599516SKenneth E. Jansen          A0(:,1,5)=
1221*59599516SKenneth E. Jansen     &             tau(:,1) * (
1222*59599516SKenneth E. Jansen     &                         + A1(:,1,1) * A2(:,1,5)
1223*59599516SKenneth E. Jansen     &                         + A2(:,1,1) * A1(:,1,5)
1224*59599516SKenneth E. Jansen     &		                                        )
1225*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1226*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A2(:,2,5)
1227*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A2(:,3,5)
1228*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A2(:,4,5)
1229*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A1(:,2,5)
1230*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A1(:,3,5)
1231*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A1(:,4,5)
1232*59599516SKenneth E. Jansen     &                                                  )
1233*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1234*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A2(:,5,5)
1235*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A1(:,5,5)
1236*59599516SKenneth E. Jansen     &		                                        )
1237*59599516SKenneth E. Jansenc
1238*59599516SKenneth E. Jansenc  row two
1239*59599516SKenneth E. Jansenc
1240*59599516SKenneth E. Jansen          A0(:,2,1)=
1241*59599516SKenneth E. Jansen     &             tau(:,1) * (
1242*59599516SKenneth E. Jansen     &                         + A1(:,2,1) * A2(:,1,1)
1243*59599516SKenneth E. Jansen     &                         + A2(:,2,1) * A1(:,1,1)
1244*59599516SKenneth E. Jansen     &		                                        )
1245*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1246*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A2(:,2,1)
1247*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A2(:,3,1)
1248*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A2(:,4,1)
1249*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A1(:,2,1)
1250*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A1(:,3,1)
1251*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A1(:,4,1)
1252*59599516SKenneth E. Jansen     &                                                  )
1253*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1254*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A2(:,5,1)
1255*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A1(:,5,1)
1256*59599516SKenneth E. Jansen     &		                                        )
1257*59599516SKenneth E. Jansenc
1258*59599516SKenneth E. Jansen          A0(:,2,2)=
1259*59599516SKenneth E. Jansen     &             tau(:,1) * (
1260*59599516SKenneth E. Jansenc    &                         + A1(:,2,1) * A2(:,1,2)
1261*59599516SKenneth E. Jansen     &                         + A2(:,2,1) * A1(:,1,2)
1262*59599516SKenneth E. Jansen     &		                                        )
1263*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1264*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A2(:,2,2)
1265*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A2(:,3,2)
1266*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A2(:,4,2)
1267*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A1(:,2,2)
1268*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A1(:,3,2)
1269*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A1(:,4,2)
1270*59599516SKenneth E. Jansen     &                                                  )
1271*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1272*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A2(:,5,2)
1273*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A1(:,5,2)
1274*59599516SKenneth E. Jansen     &		                                        )
1275*59599516SKenneth E. Jansenc
1276*59599516SKenneth E. Jansen          A0(:,2,3)=
1277*59599516SKenneth E. Jansen     &             tau(:,1) * (
1278*59599516SKenneth E. Jansen     &                         + A1(:,2,1) * A2(:,1,3)
1279*59599516SKenneth E. Jansenc    &                         + A2(:,2,1) * A1(:,1,3)
1280*59599516SKenneth E. Jansen     &		                                        )
1281*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1282*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A2(:,2,3)
1283*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A2(:,3,3)
1284*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A2(:,4,3)
1285*59599516SKenneth E. Jansenc    &                         + A2(:,2,2) * A1(:,2,3)
1286*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A1(:,3,3)
1287*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A1(:,4,3)
1288*59599516SKenneth E. Jansen     &                                                  )
1289*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1290*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A2(:,5,3)
1291*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A1(:,5,3)
1292*59599516SKenneth E. Jansen     &		                                        )
1293*59599516SKenneth E. Jansenc
1294*59599516SKenneth E. Jansen          A0(:,2,4)=
1295*59599516SKenneth E. Jansenc    &             tau(:,1) * (
1296*59599516SKenneth E. Jansenc    &                         + A1(:,2,1) * A2(:,1,4)
1297*59599516SKenneth E. Jansenc    &                         + A2(:,2,1) * A1(:,1,4)
1298*59599516SKenneth E. Jansenc    &		                                        )
1299*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
1300*59599516SKenneth E. Jansenc    &                         + A1(:,2,2) * A2(:,2,4)
1301*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A2(:,3,4)
1302*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A2(:,4,4)
1303*59599516SKenneth E. Jansenc    &                         + A2(:,2,2) * A1(:,2,4)
1304*59599516SKenneth E. Jansenc    &                         + A2(:,2,3) * A1(:,3,4)
1305*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A1(:,4,4)
1306*59599516SKenneth E. Jansenc    &                                                  )
1307*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1308*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A2(:,5,4)
1309*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A1(:,5,4)
1310*59599516SKenneth E. Jansen     &		                                        )
1311*59599516SKenneth E. Jansenc
1312*59599516SKenneth E. Jansen          A0(:,2,5)=
1313*59599516SKenneth E. Jansen     &             tau(:,1) * (
1314*59599516SKenneth E. Jansen     &                         + A1(:,2,1) * A2(:,1,5)
1315*59599516SKenneth E. Jansen     &                         + A2(:,2,1) * A1(:,1,5)
1316*59599516SKenneth E. Jansen     &		                                        )
1317*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1318*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A2(:,2,5)
1319*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A2(:,3,5)
1320*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A2(:,4,5)
1321*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A1(:,2,5)
1322*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A1(:,3,5)
1323*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A1(:,4,5)
1324*59599516SKenneth E. Jansen     &                                                  )
1325*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1326*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A2(:,5,5)
1327*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A1(:,5,5)
1328*59599516SKenneth E. Jansen     &		                                        )
1329*59599516SKenneth E. Jansenc
1330*59599516SKenneth E. Jansenc  row three
1331*59599516SKenneth E. Jansenc
1332*59599516SKenneth E. Jansen          A0(:,3,1)=
1333*59599516SKenneth E. Jansen     &             tau(:,1) * (
1334*59599516SKenneth E. Jansen     &                         + A1(:,3,1) * A2(:,1,1)
1335*59599516SKenneth E. Jansen     &                         + A2(:,3,1) * A1(:,1,1)
1336*59599516SKenneth E. Jansen     &		                                        )
1337*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1338*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A2(:,2,1)
1339*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A2(:,3,1)
1340*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A2(:,4,1)
1341*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A1(:,2,1)
1342*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A1(:,3,1)
1343*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A1(:,4,1)
1344*59599516SKenneth E. Jansen     &                                                  )
1345*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1346*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A2(:,5,1)
1347*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A1(:,5,1)
1348*59599516SKenneth E. Jansen     &		                                        )
1349*59599516SKenneth E. Jansenc
1350*59599516SKenneth E. Jansen          A0(:,3,2)=
1351*59599516SKenneth E. Jansen     &             tau(:,1) * (
1352*59599516SKenneth E. Jansenc    &                         + A1(:,3,1) * A2(:,1,2)
1353*59599516SKenneth E. Jansen     &                         + A2(:,3,1) * A1(:,1,2)
1354*59599516SKenneth E. Jansen     &		                                        )
1355*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1356*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A2(:,2,2)
1357*59599516SKenneth E. Jansenc    &                         + A1(:,3,3) * A2(:,3,2)
1358*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A2(:,4,2)
1359*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A1(:,2,2)
1360*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A1(:,3,2)
1361*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A1(:,4,2)
1362*59599516SKenneth E. Jansen     &                                                  )
1363*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1364*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A2(:,5,2)
1365*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A1(:,5,2)
1366*59599516SKenneth E. Jansen     &		                                        )
1367*59599516SKenneth E. Jansenc
1368*59599516SKenneth E. Jansen          A0(:,3,3)=
1369*59599516SKenneth E. Jansen     &             tau(:,1) * (
1370*59599516SKenneth E. Jansen     &                         + A1(:,3,1) * A2(:,1,3)
1371*59599516SKenneth E. Jansenc    &                         + A2(:,3,1) * A1(:,1,3)
1372*59599516SKenneth E. Jansen     &		                                        )
1373*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1374*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A2(:,2,3)
1375*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A2(:,3,3)
1376*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A2(:,4,3)
1377*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A1(:,2,3)
1378*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A1(:,3,3)
1379*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A1(:,4,3)
1380*59599516SKenneth E. Jansen     &                                                  )
1381*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1382*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A2(:,5,3)
1383*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A1(:,5,3)
1384*59599516SKenneth E. Jansen     &		                                        )
1385*59599516SKenneth E. Jansenc
1386*59599516SKenneth E. Jansen          A0(:,3,4)=
1387*59599516SKenneth E. Jansenc    &             tau(:,1) * (
1388*59599516SKenneth E. Jansenc    &                         + A1(:,3,1) * A2(:,1,4)
1389*59599516SKenneth E. Jansenc    &                         + A2(:,3,1) * A1(:,1,4)
1390*59599516SKenneth E. Jansenc    &		                                        )
1391*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
1392*59599516SKenneth E. Jansenc    &                         + A1(:,3,2) * A2(:,2,4)
1393*59599516SKenneth E. Jansenc    &                         + A1(:,3,3) * A2(:,3,4)
1394*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A2(:,4,4)
1395*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A1(:,2,4)
1396*59599516SKenneth E. Jansenc    &                         + A2(:,3,3) * A1(:,3,4)
1397*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A1(:,4,4)
1398*59599516SKenneth E. Jansenc    &                                                  )
1399*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1400*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A2(:,5,4)
1401*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A1(:,5,4)
1402*59599516SKenneth E. Jansen     &		                                        )
1403*59599516SKenneth E. Jansenc
1404*59599516SKenneth E. Jansen          A0(:,3,5)=
1405*59599516SKenneth E. Jansen     &             tau(:,1) * (
1406*59599516SKenneth E. Jansen     &                         + A1(:,3,1) * A2(:,1,5)
1407*59599516SKenneth E. Jansen     &                         + A2(:,3,1) * A1(:,1,5)
1408*59599516SKenneth E. Jansen     &		                                        )
1409*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1410*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A2(:,2,5)
1411*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A2(:,3,5)
1412*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A2(:,4,5)
1413*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A1(:,2,5)
1414*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A1(:,3,5)
1415*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A1(:,4,5)
1416*59599516SKenneth E. Jansen     &                                                  )
1417*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1418*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A2(:,5,5)
1419*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A1(:,5,5)
1420*59599516SKenneth E. Jansen     &		                                        )
1421*59599516SKenneth E. Jansenc
1422*59599516SKenneth E. Jansenc  row four
1423*59599516SKenneth E. Jansenc
1424*59599516SKenneth E. Jansen          A0(:,4,1)=
1425*59599516SKenneth E. Jansen     &             tau(:,1) * (
1426*59599516SKenneth E. Jansen     &                         + A1(:,4,1) * A2(:,1,1)
1427*59599516SKenneth E. Jansen     &                         + A2(:,4,1) * A1(:,1,1)
1428*59599516SKenneth E. Jansen     &		                                        )
1429*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1430*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A2(:,2,1)
1431*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A2(:,3,1)
1432*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A2(:,4,1)
1433*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A1(:,2,1)
1434*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A1(:,3,1)
1435*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A1(:,4,1)
1436*59599516SKenneth E. Jansen     &                                                  )
1437*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1438*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A2(:,5,1)
1439*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A1(:,5,1)
1440*59599516SKenneth E. Jansen     &		                                        )
1441*59599516SKenneth E. Jansenc
1442*59599516SKenneth E. Jansen          A0(:,4,2)=
1443*59599516SKenneth E. Jansen     &             tau(:,1) * (
1444*59599516SKenneth E. Jansenc    &                         + A1(:,4,1) * A2(:,1,2)
1445*59599516SKenneth E. Jansen     &                         + A2(:,4,1) * A1(:,1,2)
1446*59599516SKenneth E. Jansen     &		                                        )
1447*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1448*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A2(:,2,2)
1449*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A2(:,3,2)
1450*59599516SKenneth E. Jansenc    &                         + A1(:,4,4) * A2(:,4,2)
1451*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A1(:,2,2)
1452*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A1(:,3,2)
1453*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A1(:,4,2)
1454*59599516SKenneth E. Jansen     &                                                  )
1455*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1456*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A2(:,5,2)
1457*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A1(:,5,2)
1458*59599516SKenneth E. Jansen     &		                                        )
1459*59599516SKenneth E. Jansenc
1460*59599516SKenneth E. Jansen          A0(:,4,3)=
1461*59599516SKenneth E. Jansen     &             tau(:,1) * (
1462*59599516SKenneth E. Jansen     &                         + A1(:,4,1) * A2(:,1,3)
1463*59599516SKenneth E. Jansenc    &                         + A2(:,4,1) * A1(:,1,3)
1464*59599516SKenneth E. Jansen     &		                                        )
1465*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1466*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A2(:,2,3)
1467*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A2(:,3,3)
1468*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A2(:,4,3)
1469*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A1(:,2,3)
1470*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A1(:,3,3)
1471*59599516SKenneth E. Jansenc    &                         + A2(:,4,4) * A1(:,4,3)
1472*59599516SKenneth E. Jansen     &                                                  )
1473*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1474*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A2(:,5,3)
1475*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A1(:,5,3)
1476*59599516SKenneth E. Jansen     &		                                        )
1477*59599516SKenneth E. Jansenc
1478*59599516SKenneth E. Jansen          A0(:,4,4)=
1479*59599516SKenneth E. Jansenc    &             tau(:,1) * (
1480*59599516SKenneth E. Jansenc    &                         + A1(:,4,1) * A2(:,1,4)
1481*59599516SKenneth E. Jansenc    &                         + A2(:,4,1) * A1(:,1,4)
1482*59599516SKenneth E. Jansenc    &		                                        )
1483*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1484*59599516SKenneth E. Jansenc    &                         + A1(:,4,2) * A2(:,2,4)
1485*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A2(:,3,4)
1486*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A2(:,4,4)
1487*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A1(:,2,4)
1488*59599516SKenneth E. Jansenc    &                         + A2(:,4,3) * A1(:,3,4)
1489*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A1(:,4,4)
1490*59599516SKenneth E. Jansen     &                                                  )
1491*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1492*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A2(:,5,4)
1493*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A1(:,5,4)
1494*59599516SKenneth E. Jansen     &		                                        )
1495*59599516SKenneth E. Jansenc
1496*59599516SKenneth E. Jansen          A0(:,4,5)=
1497*59599516SKenneth E. Jansen     &             tau(:,1) * (
1498*59599516SKenneth E. Jansen     &                         + A1(:,4,1) * A2(:,1,5)
1499*59599516SKenneth E. Jansen     &                         + A2(:,4,1) * A1(:,1,5)
1500*59599516SKenneth E. Jansen     &		                                        )
1501*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1502*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A2(:,2,5)
1503*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A2(:,3,5)
1504*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A2(:,4,5)
1505*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A1(:,2,5)
1506*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A1(:,3,5)
1507*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A1(:,4,5)
1508*59599516SKenneth E. Jansen     &                                                  )
1509*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1510*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A2(:,5,5)
1511*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A1(:,5,5)
1512*59599516SKenneth E. Jansen     &		                                        )
1513*59599516SKenneth E. Jansenc
1514*59599516SKenneth E. Jansenc  row five
1515*59599516SKenneth E. Jansenc
1516*59599516SKenneth E. Jansen          A0(:,5,1)=
1517*59599516SKenneth E. Jansen     &             tau(:,1) * (
1518*59599516SKenneth E. Jansen     &                         + A1(:,5,1) * A2(:,1,1)
1519*59599516SKenneth E. Jansen     &                         + A2(:,5,1) * A1(:,1,1)
1520*59599516SKenneth E. Jansen     &		                                        )
1521*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1522*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A2(:,2,1)
1523*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A2(:,3,1)
1524*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A2(:,4,1)
1525*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A1(:,2,1)
1526*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A1(:,3,1)
1527*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A1(:,4,1)
1528*59599516SKenneth E. Jansen     &                                                  )
1529*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1530*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A2(:,5,1)
1531*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A1(:,5,1)
1532*59599516SKenneth E. Jansen     &		                                        )
1533*59599516SKenneth E. Jansenc
1534*59599516SKenneth E. Jansen          A0(:,5,2)=
1535*59599516SKenneth E. Jansen     &             tau(:,1) * (
1536*59599516SKenneth E. Jansenc    &                         + A1(:,5,1) * A2(:,1,2)
1537*59599516SKenneth E. Jansen     &                         + A2(:,5,1) * A1(:,1,2)
1538*59599516SKenneth E. Jansen     &		                                        )
1539*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1540*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A2(:,2,2)
1541*59599516SKenneth E. Jansenc    &                         + A1(:,5,3) * A2(:,3,2)
1542*59599516SKenneth E. Jansenc    &                         + A1(:,5,4) * A2(:,4,2)
1543*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A1(:,2,2)
1544*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A1(:,3,2)
1545*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A1(:,4,2)
1546*59599516SKenneth E. Jansen     &                                                  )
1547*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1548*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A2(:,5,2)
1549*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A1(:,5,2)
1550*59599516SKenneth E. Jansen     &		                                        )
1551*59599516SKenneth E. Jansenc
1552*59599516SKenneth E. Jansen          A0(:,5,3)=
1553*59599516SKenneth E. Jansen     &             tau(:,1) * (
1554*59599516SKenneth E. Jansen     &                         + A1(:,5,1) * A2(:,1,3)
1555*59599516SKenneth E. Jansenc    &                         + A2(:,5,1) * A1(:,1,3)
1556*59599516SKenneth E. Jansen     &		                                        )
1557*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1558*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A2(:,2,3)
1559*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A2(:,3,3)
1560*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A2(:,4,3)
1561*59599516SKenneth E. Jansenc    &                         + A2(:,5,2) * A1(:,2,3)
1562*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A1(:,3,3)
1563*59599516SKenneth E. Jansenc    &                         + A2(:,5,4) * A1(:,4,3)
1564*59599516SKenneth E. Jansen     &                                                  )
1565*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1566*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A2(:,5,3)
1567*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A1(:,5,3)
1568*59599516SKenneth E. Jansen     &		                                        )
1569*59599516SKenneth E. Jansenc
1570*59599516SKenneth E. Jansen          A0(:,5,4)=
1571*59599516SKenneth E. Jansenc    &             tau(:,1) * (
1572*59599516SKenneth E. Jansenc    &                         + A1(:,5,1) * A2(:,1,4)
1573*59599516SKenneth E. Jansenc    &                         + A2(:,5,1) * A1(:,1,4)
1574*59599516SKenneth E. Jansenc    &		                                        )
1575*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1576*59599516SKenneth E. Jansenc    &                         + A1(:,5,2) * A2(:,2,4)
1577*59599516SKenneth E. Jansenc    &                         + A1(:,5,3) * A2(:,3,4)
1578*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A2(:,4,4)
1579*59599516SKenneth E. Jansenc    &                         + A2(:,5,2) * A1(:,2,4)
1580*59599516SKenneth E. Jansenc    &                         + A2(:,5,3) * A1(:,3,4)
1581*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A1(:,4,4)
1582*59599516SKenneth E. Jansen     &                                                  )
1583*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1584*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A2(:,5,4)
1585*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A1(:,5,4)
1586*59599516SKenneth E. Jansen     &		                                        )
1587*59599516SKenneth E. Jansenc
1588*59599516SKenneth E. Jansen          A0(:,5,5)=
1589*59599516SKenneth E. Jansen     &             tau(:,1) * (
1590*59599516SKenneth E. Jansen     &                         + A1(:,5,1) * A2(:,1,5)
1591*59599516SKenneth E. Jansen     &                         + A2(:,5,1) * A1(:,1,5)
1592*59599516SKenneth E. Jansen     &		                                        )
1593*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1594*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A2(:,2,5)
1595*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A2(:,3,5)
1596*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A2(:,4,5)
1597*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A1(:,2,5)
1598*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A1(:,3,5)
1599*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A1(:,4,5)
1600*59599516SKenneth E. Jansen     &                                                  )
1601*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1602*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A2(:,5,5)
1603*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A1(:,5,5)
1604*59599516SKenneth E. Jansen     &		                                        )
1605*59599516SKenneth E. Jansen        else
1606*59599516SKenneth E. Jansen           A0 = zero
1607*59599516SKenneth E. Jansen        endif
1608*59599516SKenneth E. Jansenc
1609*59599516SKenneth E. Jansen        if (Navier .eq. 1) then
1610*59599516SKenneth E. Jansen           A0(:,2,3) = A0(:,2,3) + rlm2mu - rmu
1611*59599516SKenneth E. Jansen           A0(:,3,2) = A0(:,3,2) + rlm2mu - rmu
1612*59599516SKenneth E. Jansen           A0(:,5,2) = A0(:,5,2) + (rlm2mu - rmu) * u2
1613*59599516SKenneth E. Jansen           A0(:,5,3) = A0(:,5,3) + (rlm2mu- rmu) * u1
1614*59599516SKenneth E. Jansen        endif
1615*59599516SKenneth E. Jansenc
1616*59599516SKenneth E. Jansen        do j = 1, nshl
1617*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,1) * shg(:,j,2)
1618*59599516SKenneth E. Jansen           do i=1,nflow
1619*59599516SKenneth E. Jansen           do k=1,nflow
1620*59599516SKenneth E. Jansen             BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k)
1621*59599516SKenneth E. Jansen           enddo
1622*59599516SKenneth E. Jansen           enddo
1623*59599516SKenneth E. Jansen         enddo
1624*59599516SKenneth E. Jansenc
1625*59599516SKenneth E. Jansenc now for the 1 3  plus 3 1
1626*59599516SKenneth E. Jansenc
1627*59599516SKenneth E. Jansen          if(ivart.ge.2) then
1628*59599516SKenneth E. Jansenc
1629*59599516SKenneth E. Jansenc  row one
1630*59599516SKenneth E. Jansenc
1631*59599516SKenneth E. Jansen          A0(:,1,1)=
1632*59599516SKenneth E. Jansen     &             tau(:,1) * (
1633*59599516SKenneth E. Jansen     &                         + A1(:,1,1) * A3(:,1,1)
1634*59599516SKenneth E. Jansen     &                         + A3(:,1,1) * A1(:,1,1)
1635*59599516SKenneth E. Jansen     &		                                        )
1636*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1637*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A3(:,2,1)
1638*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A3(:,3,1)
1639*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A3(:,4,1)
1640*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A1(:,2,1)
1641*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A1(:,3,1)
1642*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A1(:,4,1)
1643*59599516SKenneth E. Jansen     &                                                  )
1644*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1645*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A3(:,5,1)
1646*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A1(:,5,1)
1647*59599516SKenneth E. Jansen     &		                                        )
1648*59599516SKenneth E. Jansenc
1649*59599516SKenneth E. Jansen          A0(:,1,2)=
1650*59599516SKenneth E. Jansen     &             tau(:,1) * (
1651*59599516SKenneth E. Jansenc    &                         + A1(:,1,1) * A3(:,1,2)
1652*59599516SKenneth E. Jansen     &                         + A3(:,1,1) * A1(:,1,2)
1653*59599516SKenneth E. Jansen     &		                                        )
1654*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1655*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A3(:,2,2)
1656*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A3(:,3,2)
1657*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A3(:,4,2)
1658*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A1(:,2,2)
1659*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A1(:,3,2)
1660*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A1(:,4,2)
1661*59599516SKenneth E. Jansen     &                                                  )
1662*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1663*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A3(:,5,2)
1664*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A1(:,5,2)
1665*59599516SKenneth E. Jansen     &		                                        )
1666*59599516SKenneth E. Jansenc
1667*59599516SKenneth E. Jansen          A0(:,1,3)=
1668*59599516SKenneth E. Jansenc    &             tau(:,1) * (
1669*59599516SKenneth E. Jansenc    &                         + A1(:,1,1) * A3(:,1,3)
1670*59599516SKenneth E. Jansenc    &                         + A3(:,1,1) * A1(:,1,3)
1671*59599516SKenneth E. Jansenc    &		                                        )
1672*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
1673*59599516SKenneth E. Jansenc    &                         + A1(:,1,2) * A3(:,2,3)
1674*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A3(:,3,3)
1675*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A3(:,4,3)
1676*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A1(:,2,3)
1677*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A1(:,3,3)
1678*59599516SKenneth E. Jansenc    &                         + A3(:,1,4) * A1(:,4,3)
1679*59599516SKenneth E. Jansenc    &                                                  )
1680*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1681*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A3(:,5,3)
1682*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A1(:,5,3)
1683*59599516SKenneth E. Jansen     &		                                        )
1684*59599516SKenneth E. Jansenc
1685*59599516SKenneth E. Jansen          A0(:,1,4)=
1686*59599516SKenneth E. Jansen     &             tau(:,1) * (
1687*59599516SKenneth E. Jansen     &                         + A1(:,1,1) * A3(:,1,4)
1688*59599516SKenneth E. Jansenc    &                         + A3(:,1,1) * A1(:,1,4)
1689*59599516SKenneth E. Jansen     &		                                        )
1690*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1691*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A3(:,2,4)
1692*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A3(:,3,4)
1693*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A3(:,4,4)
1694*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A1(:,2,4)
1695*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A1(:,3,4)
1696*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A1(:,4,4)
1697*59599516SKenneth E. Jansen     &                                                  )
1698*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1699*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A3(:,5,4)
1700*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A1(:,5,4)
1701*59599516SKenneth E. Jansen     &		                                        )
1702*59599516SKenneth E. Jansenc
1703*59599516SKenneth E. Jansen          A0(:,1,5)=
1704*59599516SKenneth E. Jansen     &             tau(:,1) * (
1705*59599516SKenneth E. Jansen     &                         + A1(:,1,1) * A3(:,1,5)
1706*59599516SKenneth E. Jansen     &                         + A3(:,1,1) * A1(:,1,5)
1707*59599516SKenneth E. Jansen     &		                                        )
1708*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1709*59599516SKenneth E. Jansen     &                         + A1(:,1,2) * A3(:,2,5)
1710*59599516SKenneth E. Jansenc    &                         + A1(:,1,3) * A3(:,3,5)
1711*59599516SKenneth E. Jansenc    &                         + A1(:,1,4) * A3(:,4,5)
1712*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A1(:,2,5)
1713*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A1(:,3,5)
1714*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A1(:,4,5)
1715*59599516SKenneth E. Jansen     &                                                  )
1716*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1717*59599516SKenneth E. Jansen     &		               + A1(:,1,5) * A3(:,5,5)
1718*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A1(:,5,5)
1719*59599516SKenneth E. Jansen     &		                                        )
1720*59599516SKenneth E. Jansenc
1721*59599516SKenneth E. Jansenc  row two
1722*59599516SKenneth E. Jansenc
1723*59599516SKenneth E. Jansen          A0(:,2,1)=
1724*59599516SKenneth E. Jansen     &             tau(:,1) * (
1725*59599516SKenneth E. Jansen     &                         + A1(:,2,1) * A3(:,1,1)
1726*59599516SKenneth E. Jansen     &                         + A3(:,2,1) * A1(:,1,1)
1727*59599516SKenneth E. Jansen     &		                                        )
1728*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1729*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A3(:,2,1)
1730*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A3(:,3,1)
1731*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A3(:,4,1)
1732*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A1(:,2,1)
1733*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A1(:,3,1)
1734*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A1(:,4,1)
1735*59599516SKenneth E. Jansen     &                                                  )
1736*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1737*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A3(:,5,1)
1738*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A1(:,5,1)
1739*59599516SKenneth E. Jansen     &		                                        )
1740*59599516SKenneth E. Jansenc
1741*59599516SKenneth E. Jansen          A0(:,2,2)=
1742*59599516SKenneth E. Jansen     &             tau(:,1) * (
1743*59599516SKenneth E. Jansenc    &                         + A1(:,2,1) * A3(:,1,2)
1744*59599516SKenneth E. Jansen     &                         + A3(:,2,1) * A1(:,1,2)
1745*59599516SKenneth E. Jansen     &		                                        )
1746*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1747*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A3(:,2,2)
1748*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A3(:,3,2)
1749*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A3(:,4,2)
1750*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A1(:,2,2)
1751*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A1(:,3,2)
1752*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A1(:,4,2)
1753*59599516SKenneth E. Jansen     &                                                  )
1754*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1755*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A3(:,5,2)
1756*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A1(:,5,2)
1757*59599516SKenneth E. Jansen     &		                                        )
1758*59599516SKenneth E. Jansenc
1759*59599516SKenneth E. Jansen          A0(:,2,3)=
1760*59599516SKenneth E. Jansenc    &             tau(:,1) * (
1761*59599516SKenneth E. Jansenc    &                         + A1(:,2,1) * A3(:,1,3)
1762*59599516SKenneth E. Jansenc    &                         + A3(:,2,1) * A1(:,1,3)
1763*59599516SKenneth E. Jansenc    &		                                        )
1764*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
1765*59599516SKenneth E. Jansenc    &                         + A1(:,2,2) * A3(:,2,3)
1766*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A3(:,3,3)
1767*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A3(:,4,3)
1768*59599516SKenneth E. Jansenc    &                         + A3(:,2,2) * A1(:,2,3)
1769*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A1(:,3,3)
1770*59599516SKenneth E. Jansenc    &                         + A3(:,2,4) * A1(:,4,3)
1771*59599516SKenneth E. Jansenc    &                                                  )
1772*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1773*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A3(:,5,3)
1774*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A1(:,5,3)
1775*59599516SKenneth E. Jansen     &		                                        )
1776*59599516SKenneth E. Jansenc
1777*59599516SKenneth E. Jansen          A0(:,2,4)=
1778*59599516SKenneth E. Jansen     &             tau(:,1) * (
1779*59599516SKenneth E. Jansen     &                         + A1(:,2,1) * A3(:,1,4)
1780*59599516SKenneth E. Jansenc    &                         + A3(:,2,1) * A1(:,1,4)
1781*59599516SKenneth E. Jansen     &		                                        )
1782*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1783*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A3(:,2,4)
1784*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A3(:,3,4)
1785*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A3(:,4,4)
1786*59599516SKenneth E. Jansenc    &                         + A3(:,2,2) * A1(:,2,4)
1787*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A1(:,3,4)
1788*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A1(:,4,4)
1789*59599516SKenneth E. Jansen     &                                                  )
1790*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1791*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A3(:,5,4)
1792*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A1(:,5,4)
1793*59599516SKenneth E. Jansen     &		                                        )
1794*59599516SKenneth E. Jansenc
1795*59599516SKenneth E. Jansen          A0(:,2,5)=
1796*59599516SKenneth E. Jansen     &             tau(:,1) * (
1797*59599516SKenneth E. Jansen     &                         + A1(:,2,1) * A3(:,1,5)
1798*59599516SKenneth E. Jansen     &                         + A3(:,2,1) * A1(:,1,5)
1799*59599516SKenneth E. Jansen     &		                                        )
1800*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1801*59599516SKenneth E. Jansen     &                         + A1(:,2,2) * A3(:,2,5)
1802*59599516SKenneth E. Jansenc    &                         + A1(:,2,3) * A3(:,3,5)
1803*59599516SKenneth E. Jansenc    &                         + A1(:,2,4) * A3(:,4,5)
1804*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A1(:,2,5)
1805*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A1(:,3,5)
1806*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A1(:,4,5)
1807*59599516SKenneth E. Jansen     &                                                  )
1808*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1809*59599516SKenneth E. Jansen     &		               + A1(:,2,5) * A3(:,5,5)
1810*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A1(:,5,5)
1811*59599516SKenneth E. Jansen     &		                                        )
1812*59599516SKenneth E. Jansenc
1813*59599516SKenneth E. Jansenc  row three
1814*59599516SKenneth E. Jansenc
1815*59599516SKenneth E. Jansen          A0(:,3,1)=
1816*59599516SKenneth E. Jansen     &             tau(:,1) * (
1817*59599516SKenneth E. Jansen     &                         + A1(:,3,1) * A3(:,1,1)
1818*59599516SKenneth E. Jansen     &                         + A3(:,3,1) * A1(:,1,1)
1819*59599516SKenneth E. Jansen     &		                                        )
1820*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1821*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A3(:,2,1)
1822*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A3(:,3,1)
1823*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A3(:,4,1)
1824*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A1(:,2,1)
1825*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A1(:,3,1)
1826*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A1(:,4,1)
1827*59599516SKenneth E. Jansen     &                                                  )
1828*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1829*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A3(:,5,1)
1830*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A1(:,5,1)
1831*59599516SKenneth E. Jansen     &		                                        )
1832*59599516SKenneth E. Jansenc
1833*59599516SKenneth E. Jansen          A0(:,3,2)=
1834*59599516SKenneth E. Jansen     &             tau(:,1) * (
1835*59599516SKenneth E. Jansenc    &                         + A1(:,3,1) * A3(:,1,2)
1836*59599516SKenneth E. Jansen     &                         + A3(:,3,1) * A1(:,1,2)
1837*59599516SKenneth E. Jansen     &		                                        )
1838*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1839*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A3(:,2,2)
1840*59599516SKenneth E. Jansenc    &                         + A1(:,3,3) * A3(:,3,2)
1841*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A3(:,4,2)
1842*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A1(:,2,2)
1843*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A1(:,3,2)
1844*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A1(:,4,2)
1845*59599516SKenneth E. Jansen     &                                                  )
1846*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1847*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A3(:,5,2)
1848*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A1(:,5,2)
1849*59599516SKenneth E. Jansen     &		                                        )
1850*59599516SKenneth E. Jansenc
1851*59599516SKenneth E. Jansen          A0(:,3,3)=
1852*59599516SKenneth E. Jansenc    &             tau(:,1) * (
1853*59599516SKenneth E. Jansenc    &                         + A1(:,3,1) * A3(:,1,3)
1854*59599516SKenneth E. Jansenc    &                         + A3(:,3,1) * A1(:,1,3)
1855*59599516SKenneth E. Jansenc    &		                                        )
1856*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1857*59599516SKenneth E. Jansenc    &                         + A1(:,3,2) * A3(:,2,3)
1858*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A3(:,3,3)
1859*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A3(:,4,3)
1860*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A1(:,2,3)
1861*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A1(:,3,3)
1862*59599516SKenneth E. Jansenc    &                         + A3(:,3,4) * A1(:,4,3)
1863*59599516SKenneth E. Jansen     &                                                  )
1864*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1865*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A3(:,5,3)
1866*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A1(:,5,3)
1867*59599516SKenneth E. Jansen     &		                                        )
1868*59599516SKenneth E. Jansenc
1869*59599516SKenneth E. Jansen          A0(:,3,4)=
1870*59599516SKenneth E. Jansen     &             tau(:,1) * (
1871*59599516SKenneth E. Jansen     &                         + A1(:,3,1) * A3(:,1,4)
1872*59599516SKenneth E. Jansenc    &                         + A3(:,3,1) * A1(:,1,4)
1873*59599516SKenneth E. Jansen     &		                                        )
1874*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1875*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A3(:,2,4)
1876*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A3(:,3,4)
1877*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A3(:,4,4)
1878*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A1(:,2,4)
1879*59599516SKenneth E. Jansenc    &                         + A3(:,3,3) * A1(:,3,4)
1880*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A1(:,4,4)
1881*59599516SKenneth E. Jansen     &                                                  )
1882*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1883*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A3(:,5,4)
1884*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A1(:,5,4)
1885*59599516SKenneth E. Jansen     &		                                        )
1886*59599516SKenneth E. Jansenc
1887*59599516SKenneth E. Jansen          A0(:,3,5)=
1888*59599516SKenneth E. Jansen     &             tau(:,1) * (
1889*59599516SKenneth E. Jansen     &                         + A1(:,3,1) * A3(:,1,5)
1890*59599516SKenneth E. Jansen     &                         + A3(:,3,1) * A1(:,1,5)
1891*59599516SKenneth E. Jansen     &		                                        )
1892*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1893*59599516SKenneth E. Jansen     &                         + A1(:,3,2) * A3(:,2,5)
1894*59599516SKenneth E. Jansen     &                         + A1(:,3,3) * A3(:,3,5)
1895*59599516SKenneth E. Jansenc    &                         + A1(:,3,4) * A3(:,4,5)
1896*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A1(:,2,5)
1897*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A1(:,3,5)
1898*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A1(:,4,5)
1899*59599516SKenneth E. Jansen     &                                                  )
1900*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1901*59599516SKenneth E. Jansen     &		               + A1(:,3,5) * A3(:,5,5)
1902*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A1(:,5,5)
1903*59599516SKenneth E. Jansen     &		                                        )
1904*59599516SKenneth E. Jansenc
1905*59599516SKenneth E. Jansenc  row four
1906*59599516SKenneth E. Jansenc
1907*59599516SKenneth E. Jansen          A0(:,4,1)=
1908*59599516SKenneth E. Jansen     &             tau(:,1) * (
1909*59599516SKenneth E. Jansen     &                         + A1(:,4,1) * A3(:,1,1)
1910*59599516SKenneth E. Jansen     &                         + A3(:,4,1) * A1(:,1,1)
1911*59599516SKenneth E. Jansen     &		                                        )
1912*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1913*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A3(:,2,1)
1914*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A3(:,3,1)
1915*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A3(:,4,1)
1916*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A1(:,2,1)
1917*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A1(:,3,1)
1918*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A1(:,4,1)
1919*59599516SKenneth E. Jansen     &                                                  )
1920*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1921*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A3(:,5,1)
1922*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A1(:,5,1)
1923*59599516SKenneth E. Jansen     &		                                        )
1924*59599516SKenneth E. Jansenc
1925*59599516SKenneth E. Jansen          A0(:,4,2)=
1926*59599516SKenneth E. Jansen     &             tau(:,1) * (
1927*59599516SKenneth E. Jansenc    &                         + A1(:,4,1) * A3(:,1,2)
1928*59599516SKenneth E. Jansen     &                         + A3(:,4,1) * A1(:,1,2)
1929*59599516SKenneth E. Jansen     &		                                        )
1930*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1931*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A3(:,2,2)
1932*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A3(:,3,2)
1933*59599516SKenneth E. Jansenc    &                         + A1(:,4,4) * A3(:,4,2)
1934*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A1(:,2,2)
1935*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A1(:,3,2)
1936*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A1(:,4,2)
1937*59599516SKenneth E. Jansen     &                                                  )
1938*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1939*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A3(:,5,2)
1940*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A1(:,5,2)
1941*59599516SKenneth E. Jansen     &		                                        )
1942*59599516SKenneth E. Jansenc
1943*59599516SKenneth E. Jansen          A0(:,4,3)=
1944*59599516SKenneth E. Jansenc    &             tau(:,1) * (
1945*59599516SKenneth E. Jansenc    &                         + A1(:,4,1) * A3(:,1,3)
1946*59599516SKenneth E. Jansenc    &                         + A3(:,4,1) * A1(:,1,3)
1947*59599516SKenneth E. Jansenc    &		                                        )
1948*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
1949*59599516SKenneth E. Jansenc    &                         + A1(:,4,2) * A3(:,2,3)
1950*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A3(:,3,3)
1951*59599516SKenneth E. Jansenc    &                         + A1(:,4,4) * A3(:,4,3)
1952*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A1(:,2,3)
1953*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A1(:,3,3)
1954*59599516SKenneth E. Jansenc    &                         + A3(:,4,4) * A1(:,4,3)
1955*59599516SKenneth E. Jansenc    &                                                  )
1956*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1957*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A3(:,5,3)
1958*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A1(:,5,3)
1959*59599516SKenneth E. Jansen     &		                                        )
1960*59599516SKenneth E. Jansenc
1961*59599516SKenneth E. Jansen          A0(:,4,4)=
1962*59599516SKenneth E. Jansen     &             tau(:,1) * (
1963*59599516SKenneth E. Jansen     &                         + A1(:,4,1) * A3(:,1,4)
1964*59599516SKenneth E. Jansenc    &                         + A3(:,4,1) * A1(:,1,4)
1965*59599516SKenneth E. Jansen     &		                                        )
1966*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1967*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A3(:,2,4)
1968*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A3(:,3,4)
1969*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A3(:,4,4)
1970*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A1(:,2,4)
1971*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A1(:,3,4)
1972*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A1(:,4,4)
1973*59599516SKenneth E. Jansen     &                                                  )
1974*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1975*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A3(:,5,4)
1976*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A1(:,5,4)
1977*59599516SKenneth E. Jansen     &		                                        )
1978*59599516SKenneth E. Jansenc
1979*59599516SKenneth E. Jansen          A0(:,4,5)=
1980*59599516SKenneth E. Jansen     &             tau(:,1) * (
1981*59599516SKenneth E. Jansen     &                         + A1(:,4,1) * A3(:,1,5)
1982*59599516SKenneth E. Jansen     &                         + A3(:,4,1) * A1(:,1,5)
1983*59599516SKenneth E. Jansen     &		                                        )
1984*59599516SKenneth E. Jansen     &           + tau(:,2) * (
1985*59599516SKenneth E. Jansen     &                         + A1(:,4,2) * A3(:,2,5)
1986*59599516SKenneth E. Jansenc    &                         + A1(:,4,3) * A3(:,3,5)
1987*59599516SKenneth E. Jansen     &                         + A1(:,4,4) * A3(:,4,5)
1988*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A1(:,2,5)
1989*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A1(:,3,5)
1990*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A1(:,4,5)
1991*59599516SKenneth E. Jansen     &                                                  )
1992*59599516SKenneth E. Jansen     &           + tau(:,3) * (
1993*59599516SKenneth E. Jansen     &		               + A1(:,4,5) * A3(:,5,5)
1994*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A1(:,5,5)
1995*59599516SKenneth E. Jansen     &		                                        )
1996*59599516SKenneth E. Jansenc
1997*59599516SKenneth E. Jansenc  row five
1998*59599516SKenneth E. Jansenc
1999*59599516SKenneth E. Jansen          A0(:,5,1)=
2000*59599516SKenneth E. Jansen     &             tau(:,1) * (
2001*59599516SKenneth E. Jansen     &                         + A1(:,5,1) * A3(:,1,1)
2002*59599516SKenneth E. Jansen     &                         + A3(:,5,1) * A1(:,1,1)
2003*59599516SKenneth E. Jansen     &		                                        )
2004*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2005*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A3(:,2,1)
2006*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A3(:,3,1)
2007*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A3(:,4,1)
2008*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A1(:,2,1)
2009*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A1(:,3,1)
2010*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A1(:,4,1)
2011*59599516SKenneth E. Jansen     &                                                  )
2012*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2013*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A3(:,5,1)
2014*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A1(:,5,1)
2015*59599516SKenneth E. Jansen     &		                                        )
2016*59599516SKenneth E. Jansenc
2017*59599516SKenneth E. Jansen          A0(:,5,2)=
2018*59599516SKenneth E. Jansen     &             tau(:,1) * (
2019*59599516SKenneth E. Jansenc    &                         + A1(:,5,1) * A3(:,1,2)
2020*59599516SKenneth E. Jansen     &                         + A3(:,5,1) * A1(:,1,2)
2021*59599516SKenneth E. Jansen     &		                                        )
2022*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2023*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A3(:,2,2)
2024*59599516SKenneth E. Jansenc    &                         + A1(:,5,3) * A3(:,3,2)
2025*59599516SKenneth E. Jansenc    &                         + A1(:,5,4) * A3(:,4,2)
2026*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A1(:,2,2)
2027*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A1(:,3,2)
2028*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A1(:,4,2)
2029*59599516SKenneth E. Jansen     &                                                  )
2030*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2031*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A3(:,5,2)
2032*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A1(:,5,2)
2033*59599516SKenneth E. Jansen     &		                                        )
2034*59599516SKenneth E. Jansenc
2035*59599516SKenneth E. Jansen          A0(:,5,3)=
2036*59599516SKenneth E. Jansenc    &             tau(:,1) * (
2037*59599516SKenneth E. Jansenc    &                         + A1(:,5,1) * A3(:,1,3)
2038*59599516SKenneth E. Jansenc    &                         + A3(:,5,1) * A1(:,1,3)
2039*59599516SKenneth E. Jansenc    &		                                        )
2040*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2041*59599516SKenneth E. Jansenc    &                         + A1(:,5,2) * A3(:,2,3)
2042*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A3(:,3,3)
2043*59599516SKenneth E. Jansenc    &                         + A1(:,5,4) * A3(:,4,3)
2044*59599516SKenneth E. Jansenc    &                         + A3(:,5,2) * A1(:,2,3)
2045*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A1(:,3,3)
2046*59599516SKenneth E. Jansenc    &                         + A3(:,5,4) * A1(:,4,3)
2047*59599516SKenneth E. Jansen     &                                                  )
2048*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2049*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A3(:,5,3)
2050*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A1(:,5,3)
2051*59599516SKenneth E. Jansen     &		                                        )
2052*59599516SKenneth E. Jansenc
2053*59599516SKenneth E. Jansen          A0(:,5,4)=
2054*59599516SKenneth E. Jansen     &             tau(:,1) * (
2055*59599516SKenneth E. Jansen     &                         + A1(:,5,1) * A3(:,1,4)
2056*59599516SKenneth E. Jansenc    &                         + A3(:,5,1) * A1(:,1,4)
2057*59599516SKenneth E. Jansen     &		                                        )
2058*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2059*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A3(:,2,4)
2060*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A3(:,3,4)
2061*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A3(:,4,4)
2062*59599516SKenneth E. Jansenc    &                         + A3(:,5,2) * A1(:,2,4)
2063*59599516SKenneth E. Jansenc    &                         + A3(:,5,3) * A1(:,3,4)
2064*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A1(:,4,4)
2065*59599516SKenneth E. Jansen     &                                                  )
2066*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2067*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A3(:,5,4)
2068*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A1(:,5,4)
2069*59599516SKenneth E. Jansen     &		                                        )
2070*59599516SKenneth E. Jansenc
2071*59599516SKenneth E. Jansen          A0(:,5,5)=
2072*59599516SKenneth E. Jansen     &             tau(:,1) * (
2073*59599516SKenneth E. Jansen     &                         + A1(:,5,1) * A3(:,1,5)
2074*59599516SKenneth E. Jansen     &                         + A3(:,5,1) * A1(:,1,5)
2075*59599516SKenneth E. Jansen     &		                                        )
2076*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2077*59599516SKenneth E. Jansen     &                         + A1(:,5,2) * A3(:,2,5)
2078*59599516SKenneth E. Jansen     &                         + A1(:,5,3) * A3(:,3,5)
2079*59599516SKenneth E. Jansen     &                         + A1(:,5,4) * A3(:,4,5)
2080*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A1(:,2,5)
2081*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A1(:,3,5)
2082*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A1(:,4,5)
2083*59599516SKenneth E. Jansen     &                                                  )
2084*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2085*59599516SKenneth E. Jansen     &		               + A1(:,5,5) * A3(:,5,5)
2086*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A1(:,5,5)
2087*59599516SKenneth E. Jansen     &		                                        )
2088*59599516SKenneth E. Jansen        else
2089*59599516SKenneth E. Jansen           A0 = zero
2090*59599516SKenneth E. Jansen        endif
2091*59599516SKenneth E. Jansenc
2092*59599516SKenneth E. Jansen        if (Navier .eq. 1) then
2093*59599516SKenneth E. Jansen           A0(:,2,4) = A0(:,2,4) + rlm2mu - rmu
2094*59599516SKenneth E. Jansen           A0(:,4,2) = A0(:,4,2) + rlm2mu - rmu
2095*59599516SKenneth E. Jansen           A0(:,5,2) = A0(:,5,2) + (rlm2mu-rmu) * u3
2096*59599516SKenneth E. Jansen           A0(:,5,4) = A0(:,5,4) + (rlm2mu-rmu) * u1
2097*59599516SKenneth E. Jansen        endif
2098*59599516SKenneth E. Jansenc
2099*59599516SKenneth E. Jansen        do j = 1, nshl
2100*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,1) * shg(:,j,3)
2101*59599516SKenneth E. Jansen           do i=1,nflow
2102*59599516SKenneth E. Jansen           do k=1,nflow
2103*59599516SKenneth E. Jansen             BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k)
2104*59599516SKenneth E. Jansen           enddo
2105*59599516SKenneth E. Jansen           enddo
2106*59599516SKenneth E. Jansen         enddo
2107*59599516SKenneth E. Jansenc
2108*59599516SKenneth E. Jansenc now for the 2 3  plus 3 2
2109*59599516SKenneth E. Jansenc
2110*59599516SKenneth E. Jansen          if(ivart.ge.2) then
2111*59599516SKenneth E. Jansenc
2112*59599516SKenneth E. Jansenc  row one
2113*59599516SKenneth E. Jansenc
2114*59599516SKenneth E. Jansen          A0(:,1,1)=
2115*59599516SKenneth E. Jansen     &             tau(:,1) * (
2116*59599516SKenneth E. Jansen     &                         + A2(:,1,1) * A3(:,1,1)
2117*59599516SKenneth E. Jansen     &                         + A3(:,1,1) * A2(:,1,1)
2118*59599516SKenneth E. Jansen     &		                                        )
2119*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2120*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A3(:,2,1)
2121*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A3(:,3,1)
2122*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A3(:,4,1)
2123*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A2(:,2,1)
2124*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A2(:,3,1)
2125*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A2(:,4,1)
2126*59599516SKenneth E. Jansen     &                                                  )
2127*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2128*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A3(:,5,1)
2129*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A2(:,5,1)
2130*59599516SKenneth E. Jansen     &		                                        )
2131*59599516SKenneth E. Jansenc
2132*59599516SKenneth E. Jansen          A0(:,1,2)=
2133*59599516SKenneth E. Jansenc    &             tau(:,1) * (
2134*59599516SKenneth E. Jansenc    &                         + A2(:,1,1) * A3(:,1,2)
2135*59599516SKenneth E. Jansenc    &                         + A3(:,1,1) * A2(:,1,2)
2136*59599516SKenneth E. Jansenc    &		                                        )
2137*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
2138*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A3(:,2,2)
2139*59599516SKenneth E. Jansenc    &                         + A2(:,1,3) * A3(:,3,2)
2140*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A3(:,4,2)
2141*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A2(:,2,2)
2142*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A2(:,3,2)
2143*59599516SKenneth E. Jansenc    &                         + A3(:,1,4) * A2(:,4,2)
2144*59599516SKenneth E. Jansenc    &                                                  )
2145*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2146*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A3(:,5,2)
2147*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A2(:,5,2)
2148*59599516SKenneth E. Jansen     &		                                        )
2149*59599516SKenneth E. Jansenc
2150*59599516SKenneth E. Jansen          A0(:,1,3)=
2151*59599516SKenneth E. Jansen     &             tau(:,1) * (
2152*59599516SKenneth E. Jansenc    &                         + A2(:,1,1) * A3(:,1,3)
2153*59599516SKenneth E. Jansen     &                         + A3(:,1,1) * A2(:,1,3)
2154*59599516SKenneth E. Jansen     &		                                        )
2155*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2156*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A3(:,2,3)
2157*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A3(:,3,3)
2158*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A3(:,4,3)
2159*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A2(:,2,3)
2160*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A2(:,3,3)
2161*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A2(:,4,3)
2162*59599516SKenneth E. Jansen     &                                                  )
2163*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2164*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A3(:,5,3)
2165*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A2(:,5,3)
2166*59599516SKenneth E. Jansen     &		                                        )
2167*59599516SKenneth E. Jansenc
2168*59599516SKenneth E. Jansen          A0(:,1,4)=
2169*59599516SKenneth E. Jansen     &             tau(:,1) * (
2170*59599516SKenneth E. Jansen     &                         + A2(:,1,1) * A3(:,1,4)
2171*59599516SKenneth E. Jansenc    &                         + A3(:,1,1) * A2(:,1,4)
2172*59599516SKenneth E. Jansen     &		                                        )
2173*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2174*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A3(:,2,4)
2175*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A3(:,3,4)
2176*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A3(:,4,4)
2177*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A2(:,2,4)
2178*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A2(:,3,4)
2179*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A2(:,4,4)
2180*59599516SKenneth E. Jansen     &                                                  )
2181*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2182*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A3(:,5,4)
2183*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A2(:,5,4)
2184*59599516SKenneth E. Jansen     &		                                        )
2185*59599516SKenneth E. Jansenc
2186*59599516SKenneth E. Jansen          A0(:,1,5)=
2187*59599516SKenneth E. Jansen     &             tau(:,1) * (
2188*59599516SKenneth E. Jansen     &                         + A2(:,1,1) * A3(:,1,5)
2189*59599516SKenneth E. Jansen     &                         + A3(:,1,1) * A2(:,1,5)
2190*59599516SKenneth E. Jansen     &		                                        )
2191*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2192*59599516SKenneth E. Jansenc    &                         + A2(:,1,2) * A3(:,2,5)
2193*59599516SKenneth E. Jansen     &                         + A2(:,1,3) * A3(:,3,5)
2194*59599516SKenneth E. Jansenc    &                         + A2(:,1,4) * A3(:,4,5)
2195*59599516SKenneth E. Jansenc    &                         + A3(:,1,2) * A2(:,2,5)
2196*59599516SKenneth E. Jansenc    &                         + A3(:,1,3) * A2(:,3,5)
2197*59599516SKenneth E. Jansen     &                         + A3(:,1,4) * A2(:,4,5)
2198*59599516SKenneth E. Jansen     &                                                  )
2199*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2200*59599516SKenneth E. Jansen     &		               + A2(:,1,5) * A3(:,5,5)
2201*59599516SKenneth E. Jansen     &		               + A3(:,1,5) * A2(:,5,5)
2202*59599516SKenneth E. Jansen     &		                                        )
2203*59599516SKenneth E. Jansenc
2204*59599516SKenneth E. Jansenc  row two
2205*59599516SKenneth E. Jansenc
2206*59599516SKenneth E. Jansen          A0(:,2,1)=
2207*59599516SKenneth E. Jansen     &             tau(:,1) * (
2208*59599516SKenneth E. Jansen     &                         + A2(:,2,1) * A3(:,1,1)
2209*59599516SKenneth E. Jansen     &                         + A3(:,2,1) * A2(:,1,1)
2210*59599516SKenneth E. Jansen     &		                                        )
2211*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2212*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A3(:,2,1)
2213*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A3(:,3,1)
2214*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A3(:,4,1)
2215*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A2(:,2,1)
2216*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A2(:,3,1)
2217*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A2(:,4,1)
2218*59599516SKenneth E. Jansen     &                                                  )
2219*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2220*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A3(:,5,1)
2221*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A2(:,5,1)
2222*59599516SKenneth E. Jansen     &		                                        )
2223*59599516SKenneth E. Jansenc
2224*59599516SKenneth E. Jansen          A0(:,2,2)=
2225*59599516SKenneth E. Jansenc    &             tau(:,1) * (
2226*59599516SKenneth E. Jansenc    &                         + A2(:,2,1) * A3(:,1,2)
2227*59599516SKenneth E. Jansenc    &                         + A3(:,2,1) * A2(:,1,2)
2228*59599516SKenneth E. Jansenc    &		                                        )
2229*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2230*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A3(:,2,2)
2231*59599516SKenneth E. Jansenc    &                         + A2(:,2,3) * A3(:,3,2)
2232*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A3(:,4,2)
2233*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A2(:,2,2)
2234*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A2(:,3,2)
2235*59599516SKenneth E. Jansenc    &                         + A3(:,2,4) * A2(:,4,2)
2236*59599516SKenneth E. Jansen     &                                                  )
2237*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2238*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A3(:,5,2)
2239*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A2(:,5,2)
2240*59599516SKenneth E. Jansen     &		                                        )
2241*59599516SKenneth E. Jansenc
2242*59599516SKenneth E. Jansen          A0(:,2,3)=
2243*59599516SKenneth E. Jansen     &             tau(:,1) * (
2244*59599516SKenneth E. Jansenc    &                         + A2(:,2,1) * A3(:,1,3)
2245*59599516SKenneth E. Jansen     &                         + A3(:,2,1) * A2(:,1,3)
2246*59599516SKenneth E. Jansen     &		                                        )
2247*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2248*59599516SKenneth E. Jansenc    &                         + A2(:,2,2) * A3(:,2,3)
2249*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A3(:,3,3)
2250*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A3(:,4,3)
2251*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A2(:,2,3)
2252*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A2(:,3,3)
2253*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A2(:,4,3)
2254*59599516SKenneth E. Jansen     &                                                  )
2255*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2256*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A3(:,5,3)
2257*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A2(:,5,3)
2258*59599516SKenneth E. Jansen     &		                                        )
2259*59599516SKenneth E. Jansenc
2260*59599516SKenneth E. Jansen          A0(:,2,4)=
2261*59599516SKenneth E. Jansen     &             tau(:,1) * (
2262*59599516SKenneth E. Jansen     &                         + A2(:,2,1) * A3(:,1,4)
2263*59599516SKenneth E. Jansenc    &                         + A3(:,2,1) * A2(:,1,4)
2264*59599516SKenneth E. Jansen     &		                                        )
2265*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2266*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A3(:,2,4)
2267*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A3(:,3,4)
2268*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A3(:,4,4)
2269*59599516SKenneth E. Jansenc    &                         + A3(:,2,2) * A2(:,2,4)
2270*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A2(:,3,4)
2271*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A2(:,4,4)
2272*59599516SKenneth E. Jansen     &                                                  )
2273*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2274*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A3(:,5,4)
2275*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A2(:,5,4)
2276*59599516SKenneth E. Jansen     &		                                        )
2277*59599516SKenneth E. Jansenc
2278*59599516SKenneth E. Jansen          A0(:,2,5)=
2279*59599516SKenneth E. Jansen     &             tau(:,1) * (
2280*59599516SKenneth E. Jansen     &                         + A2(:,2,1) * A3(:,1,5)
2281*59599516SKenneth E. Jansen     &                         + A3(:,2,1) * A2(:,1,5)
2282*59599516SKenneth E. Jansen     &		                                        )
2283*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2284*59599516SKenneth E. Jansen     &                         + A2(:,2,2) * A3(:,2,5)
2285*59599516SKenneth E. Jansen     &                         + A2(:,2,3) * A3(:,3,5)
2286*59599516SKenneth E. Jansenc    &                         + A2(:,2,4) * A3(:,4,5)
2287*59599516SKenneth E. Jansen     &                         + A3(:,2,2) * A2(:,2,5)
2288*59599516SKenneth E. Jansenc    &                         + A3(:,2,3) * A2(:,3,5)
2289*59599516SKenneth E. Jansen     &                         + A3(:,2,4) * A2(:,4,5)
2290*59599516SKenneth E. Jansen     &                                                  )
2291*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2292*59599516SKenneth E. Jansen     &		               + A2(:,2,5) * A3(:,5,5)
2293*59599516SKenneth E. Jansen     &		               + A3(:,2,5) * A2(:,5,5)
2294*59599516SKenneth E. Jansen     &		                                        )
2295*59599516SKenneth E. Jansenc
2296*59599516SKenneth E. Jansenc  row three
2297*59599516SKenneth E. Jansenc
2298*59599516SKenneth E. Jansen          A0(:,3,1)=
2299*59599516SKenneth E. Jansen     &             tau(:,1) * (
2300*59599516SKenneth E. Jansen     &                         + A2(:,3,1) * A3(:,1,1)
2301*59599516SKenneth E. Jansen     &                         + A3(:,3,1) * A2(:,1,1)
2302*59599516SKenneth E. Jansen     &		                                        )
2303*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2304*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A3(:,2,1)
2305*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A3(:,3,1)
2306*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A3(:,4,1)
2307*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A2(:,2,1)
2308*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A2(:,3,1)
2309*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A2(:,4,1)
2310*59599516SKenneth E. Jansen     &                                                  )
2311*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2312*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A3(:,5,1)
2313*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A2(:,5,1)
2314*59599516SKenneth E. Jansen     &		                                        )
2315*59599516SKenneth E. Jansenc
2316*59599516SKenneth E. Jansen          A0(:,3,2)=
2317*59599516SKenneth E. Jansenc    &             tau(:,1) * (
2318*59599516SKenneth E. Jansenc    &                         + A2(:,3,1) * A3(:,1,2)
2319*59599516SKenneth E. Jansenc    &                         + A3(:,3,1) * A2(:,1,2)
2320*59599516SKenneth E. Jansenc    &		                                        )
2321*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
2322*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A3(:,2,2)
2323*59599516SKenneth E. Jansenc    &                         + A2(:,3,3) * A3(:,3,2)
2324*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A3(:,4,2)
2325*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A2(:,2,2)
2326*59599516SKenneth E. Jansenc    &                         + A3(:,3,3) * A2(:,3,2)
2327*59599516SKenneth E. Jansenc    &                         + A3(:,3,4) * A2(:,4,2)
2328*59599516SKenneth E. Jansenc    &                                                  )
2329*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2330*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A3(:,5,2)
2331*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A2(:,5,2)
2332*59599516SKenneth E. Jansen     &		                                        )
2333*59599516SKenneth E. Jansenc
2334*59599516SKenneth E. Jansen          A0(:,3,3)=
2335*59599516SKenneth E. Jansen     &             tau(:,1) * (
2336*59599516SKenneth E. Jansenc    &                         + A2(:,3,1) * A3(:,1,3)
2337*59599516SKenneth E. Jansen     &                         + A3(:,3,1) * A2(:,1,3)
2338*59599516SKenneth E. Jansen     &		                                        )
2339*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2340*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A3(:,2,3)
2341*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A3(:,3,3)
2342*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A3(:,4,3)
2343*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A2(:,2,3)
2344*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A2(:,3,3)
2345*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A2(:,4,3)
2346*59599516SKenneth E. Jansen     &                                                  )
2347*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2348*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A3(:,5,3)
2349*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A2(:,5,3)
2350*59599516SKenneth E. Jansen     &		                                        )
2351*59599516SKenneth E. Jansenc
2352*59599516SKenneth E. Jansen          A0(:,3,4)=
2353*59599516SKenneth E. Jansen     &             tau(:,1) * (
2354*59599516SKenneth E. Jansen     &                         + A2(:,3,1) * A3(:,1,4)
2355*59599516SKenneth E. Jansenc    &                         + A3(:,3,1) * A2(:,1,4)
2356*59599516SKenneth E. Jansen     &		                                        )
2357*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2358*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A3(:,2,4)
2359*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A3(:,3,4)
2360*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A3(:,4,4)
2361*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A2(:,2,4)
2362*59599516SKenneth E. Jansenc    &                         + A3(:,3,3) * A2(:,3,4)
2363*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A2(:,4,4)
2364*59599516SKenneth E. Jansen     &                                                  )
2365*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2366*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A3(:,5,4)
2367*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A2(:,5,4)
2368*59599516SKenneth E. Jansen     &		                                        )
2369*59599516SKenneth E. Jansenc
2370*59599516SKenneth E. Jansen          A0(:,3,5)=
2371*59599516SKenneth E. Jansen     &             tau(:,1) * (
2372*59599516SKenneth E. Jansen     &                         + A2(:,3,1) * A3(:,1,5)
2373*59599516SKenneth E. Jansen     &                         + A3(:,3,1) * A2(:,1,5)
2374*59599516SKenneth E. Jansen     &		                                        )
2375*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2376*59599516SKenneth E. Jansenc    &                         + A2(:,3,2) * A3(:,2,5)
2377*59599516SKenneth E. Jansen     &                         + A2(:,3,3) * A3(:,3,5)
2378*59599516SKenneth E. Jansenc    &                         + A2(:,3,4) * A3(:,4,5)
2379*59599516SKenneth E. Jansenc    &                         + A3(:,3,2) * A2(:,2,5)
2380*59599516SKenneth E. Jansen     &                         + A3(:,3,3) * A2(:,3,5)
2381*59599516SKenneth E. Jansen     &                         + A3(:,3,4) * A2(:,4,5)
2382*59599516SKenneth E. Jansen     &                                                  )
2383*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2384*59599516SKenneth E. Jansen     &		               + A2(:,3,5) * A3(:,5,5)
2385*59599516SKenneth E. Jansen     &		               + A3(:,3,5) * A2(:,5,5)
2386*59599516SKenneth E. Jansen     &		                                        )
2387*59599516SKenneth E. Jansenc
2388*59599516SKenneth E. Jansenc  row four
2389*59599516SKenneth E. Jansenc
2390*59599516SKenneth E. Jansen          A0(:,4,1)=
2391*59599516SKenneth E. Jansen     &             tau(:,1) * (
2392*59599516SKenneth E. Jansen     &                         + A2(:,4,1) * A3(:,1,1)
2393*59599516SKenneth E. Jansen     &                         + A3(:,4,1) * A2(:,1,1)
2394*59599516SKenneth E. Jansen     &		                                        )
2395*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2396*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A3(:,2,1)
2397*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A3(:,3,1)
2398*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A3(:,4,1)
2399*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A2(:,2,1)
2400*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A2(:,3,1)
2401*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A2(:,4,1)
2402*59599516SKenneth E. Jansen     &                                                  )
2403*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2404*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A3(:,5,1)
2405*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A2(:,5,1)
2406*59599516SKenneth E. Jansen     &		                                        )
2407*59599516SKenneth E. Jansenc
2408*59599516SKenneth E. Jansen          A0(:,4,2)=
2409*59599516SKenneth E. Jansenc    &             tau(:,1) * (
2410*59599516SKenneth E. Jansenc    &                         + A2(:,4,1) * A3(:,1,2)
2411*59599516SKenneth E. Jansenc    &                         + A3(:,4,1) * A2(:,1,2)
2412*59599516SKenneth E. Jansenc    &		                                        )
2413*59599516SKenneth E. Jansenc    &           + tau(:,2) * (
2414*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A3(:,2,2)
2415*59599516SKenneth E. Jansenc    &                         + A2(:,4,3) * A3(:,3,2)
2416*59599516SKenneth E. Jansenc    &                         + A2(:,4,4) * A3(:,4,2)
2417*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A2(:,2,2)
2418*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A2(:,3,2)
2419*59599516SKenneth E. Jansenc    &                         + A3(:,4,4) * A2(:,4,2)
2420*59599516SKenneth E. Jansenc    &                                                  )
2421*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2422*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A3(:,5,2)
2423*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A2(:,5,2)
2424*59599516SKenneth E. Jansen     &		                                        )
2425*59599516SKenneth E. Jansenc
2426*59599516SKenneth E. Jansen          A0(:,4,3)=
2427*59599516SKenneth E. Jansen     &             tau(:,1) * (
2428*59599516SKenneth E. Jansenc    &                         + A2(:,4,1) * A3(:,1,3)
2429*59599516SKenneth E. Jansen     &                         + A3(:,4,1) * A2(:,1,3)
2430*59599516SKenneth E. Jansen     &		                                        )
2431*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2432*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A3(:,2,3)
2433*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A3(:,3,3)
2434*59599516SKenneth E. Jansenc    &                         + A2(:,4,4) * A3(:,4,3)
2435*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A2(:,2,3)
2436*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A2(:,3,3)
2437*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A2(:,4,3)
2438*59599516SKenneth E. Jansen     &                                                  )
2439*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2440*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A3(:,5,3)
2441*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A2(:,5,3)
2442*59599516SKenneth E. Jansen     &		                                        )
2443*59599516SKenneth E. Jansenc
2444*59599516SKenneth E. Jansen          A0(:,4,4)=
2445*59599516SKenneth E. Jansen     &             tau(:,1) * (
2446*59599516SKenneth E. Jansen     &                         + A2(:,4,1) * A3(:,1,4)
2447*59599516SKenneth E. Jansenc    &                         + A3(:,4,1) * A2(:,1,4)
2448*59599516SKenneth E. Jansen     &		                                        )
2449*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2450*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A3(:,2,4)
2451*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A3(:,3,4)
2452*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A3(:,4,4)
2453*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A2(:,2,4)
2454*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A2(:,3,4)
2455*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A2(:,4,4)
2456*59599516SKenneth E. Jansen     &                                                  )
2457*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2458*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A3(:,5,4)
2459*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A2(:,5,4)
2460*59599516SKenneth E. Jansen     &		                                        )
2461*59599516SKenneth E. Jansenc
2462*59599516SKenneth E. Jansen          A0(:,4,5)=
2463*59599516SKenneth E. Jansen     &             tau(:,1) * (
2464*59599516SKenneth E. Jansen     &                         + A2(:,4,1) * A3(:,1,5)
2465*59599516SKenneth E. Jansen     &                         + A3(:,4,1) * A2(:,1,5)
2466*59599516SKenneth E. Jansen     &		                                        )
2467*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2468*59599516SKenneth E. Jansenc    &                         + A2(:,4,2) * A3(:,2,5)
2469*59599516SKenneth E. Jansen     &                         + A2(:,4,3) * A3(:,3,5)
2470*59599516SKenneth E. Jansen     &                         + A2(:,4,4) * A3(:,4,5)
2471*59599516SKenneth E. Jansenc    &                         + A3(:,4,2) * A2(:,2,5)
2472*59599516SKenneth E. Jansenc    &                         + A3(:,4,3) * A2(:,3,5)
2473*59599516SKenneth E. Jansen     &                         + A3(:,4,4) * A2(:,4,5)
2474*59599516SKenneth E. Jansen     &                                                  )
2475*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2476*59599516SKenneth E. Jansen     &		               + A2(:,4,5) * A3(:,5,5)
2477*59599516SKenneth E. Jansen     &		               + A3(:,4,5) * A2(:,5,5)
2478*59599516SKenneth E. Jansen     &		                                        )
2479*59599516SKenneth E. Jansenc
2480*59599516SKenneth E. Jansenc  row five
2481*59599516SKenneth E. Jansenc
2482*59599516SKenneth E. Jansen          A0(:,5,1)=
2483*59599516SKenneth E. Jansen     &             tau(:,1) * (
2484*59599516SKenneth E. Jansen     &                         + A2(:,5,1) * A3(:,1,1)
2485*59599516SKenneth E. Jansen     &                         + A3(:,5,1) * A2(:,1,1)
2486*59599516SKenneth E. Jansen     &		                                        )
2487*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2488*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A3(:,2,1)
2489*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A3(:,3,1)
2490*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A3(:,4,1)
2491*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A2(:,2,1)
2492*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A2(:,3,1)
2493*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A2(:,4,1)
2494*59599516SKenneth E. Jansen     &                                                  )
2495*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2496*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A3(:,5,1)
2497*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A2(:,5,1)
2498*59599516SKenneth E. Jansen     &		                                        )
2499*59599516SKenneth E. Jansenc
2500*59599516SKenneth E. Jansen          A0(:,5,2)=
2501*59599516SKenneth E. Jansenc    &             tau(:,1) * (
2502*59599516SKenneth E. Jansenc    &                         + A2(:,5,1) * A3(:,1,2)
2503*59599516SKenneth E. Jansenc    &                         + A3(:,5,1) * A2(:,1,2)
2504*59599516SKenneth E. Jansenc    &		                                        )
2505*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2506*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A3(:,2,2)
2507*59599516SKenneth E. Jansenc    &                         + A2(:,5,3) * A3(:,3,2)
2508*59599516SKenneth E. Jansenc    &                         + A2(:,5,4) * A3(:,4,2)
2509*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A2(:,2,2)
2510*59599516SKenneth E. Jansenc    &                         + A3(:,5,3) * A2(:,3,2)
2511*59599516SKenneth E. Jansenc    &                         + A3(:,5,4) * A2(:,4,2)
2512*59599516SKenneth E. Jansen     &                                                  )
2513*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2514*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A3(:,5,2)
2515*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A2(:,5,2)
2516*59599516SKenneth E. Jansen     &		                                        )
2517*59599516SKenneth E. Jansenc
2518*59599516SKenneth E. Jansen          A0(:,5,3)=
2519*59599516SKenneth E. Jansen     &             tau(:,1) * (
2520*59599516SKenneth E. Jansenc    &                         + A2(:,5,1) * A3(:,1,3)
2521*59599516SKenneth E. Jansen     &                         + A3(:,5,1) * A2(:,1,3)
2522*59599516SKenneth E. Jansen     &		                                        )
2523*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2524*59599516SKenneth E. Jansenc    &                         + A2(:,5,2) * A3(:,2,3)
2525*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A3(:,3,3)
2526*59599516SKenneth E. Jansenc    &                         + A2(:,5,4) * A3(:,4,3)
2527*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A2(:,2,3)
2528*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A2(:,3,3)
2529*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A2(:,4,3)
2530*59599516SKenneth E. Jansen     &                                                  )
2531*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2532*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A3(:,5,3)
2533*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A2(:,5,3)
2534*59599516SKenneth E. Jansen     &		                                        )
2535*59599516SKenneth E. Jansenc
2536*59599516SKenneth E. Jansen          A0(:,5,4)=
2537*59599516SKenneth E. Jansen     &             tau(:,1) * (
2538*59599516SKenneth E. Jansen     &                         + A2(:,5,1) * A3(:,1,4)
2539*59599516SKenneth E. Jansenc    &                         + A3(:,5,1) * A2(:,1,4)
2540*59599516SKenneth E. Jansen     &		                                        )
2541*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2542*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A3(:,2,4)
2543*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A3(:,3,4)
2544*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A3(:,4,4)
2545*59599516SKenneth E. Jansenc    &                         + A3(:,5,2) * A2(:,2,4)
2546*59599516SKenneth E. Jansenc    &                         + A3(:,5,3) * A2(:,3,4)
2547*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A2(:,4,4)
2548*59599516SKenneth E. Jansen     &                                                  )
2549*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2550*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A3(:,5,4)
2551*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A2(:,5,4)
2552*59599516SKenneth E. Jansen     &		                                        )
2553*59599516SKenneth E. Jansenc
2554*59599516SKenneth E. Jansen          A0(:,5,5)=
2555*59599516SKenneth E. Jansen     &             tau(:,1) * (
2556*59599516SKenneth E. Jansen     &                         + A2(:,5,1) * A3(:,1,5)
2557*59599516SKenneth E. Jansen     &                         + A3(:,5,1) * A2(:,1,5)
2558*59599516SKenneth E. Jansen     &		                                        )
2559*59599516SKenneth E. Jansen     &           + tau(:,2) * (
2560*59599516SKenneth E. Jansen     &                         + A2(:,5,2) * A3(:,2,5)
2561*59599516SKenneth E. Jansen     &                         + A2(:,5,3) * A3(:,3,5)
2562*59599516SKenneth E. Jansen     &                         + A2(:,5,4) * A3(:,4,5)
2563*59599516SKenneth E. Jansen     &                         + A3(:,5,2) * A2(:,2,5)
2564*59599516SKenneth E. Jansen     &                         + A3(:,5,3) * A2(:,3,5)
2565*59599516SKenneth E. Jansen     &                         + A3(:,5,4) * A2(:,4,5)
2566*59599516SKenneth E. Jansen     &                                                  )
2567*59599516SKenneth E. Jansen     &           + tau(:,3) * (
2568*59599516SKenneth E. Jansen     &		               + A2(:,5,5) * A3(:,5,5)
2569*59599516SKenneth E. Jansen     &		               + A3(:,5,5) * A2(:,5,5)
2570*59599516SKenneth E. Jansen     &		                                        )
2571*59599516SKenneth E. Jansen        else
2572*59599516SKenneth E. Jansen           A0 = zero
2573*59599516SKenneth E. Jansen        endif
2574*59599516SKenneth E. Jansenc
2575*59599516SKenneth E. Jansen        if (Navier .eq. 1) then
2576*59599516SKenneth E. Jansen           A0(:,3,4) = A0(:,3,4) + rlm2mu - rmu
2577*59599516SKenneth E. Jansen           A0(:,4,3) = A0(:,4,3) + rlm2mu - rmu
2578*59599516SKenneth E. Jansen           A0(:,5,3) = A0(:,5,3) + (rlm2mu-rmu) * u3
2579*59599516SKenneth E. Jansen           A0(:,5,4) = A0(:,5,4) + (rlm2mu-rmu) * u2
2580*59599516SKenneth E. Jansen        endif
2581*59599516SKenneth E. Jansen
2582*59599516SKenneth E. Jansen        do j = 1, nshl
2583*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,2) * shg(:,j,3)
2584*59599516SKenneth E. Jansen           do i=1,nflow
2585*59599516SKenneth E. Jansen           do k=1,nflow
2586*59599516SKenneth E. Jansen             BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k)
2587*59599516SKenneth E. Jansen           enddo
2588*59599516SKenneth E. Jansen           enddo
2589*59599516SKenneth E. Jansen        enddo
2590*59599516SKenneth E. Jansen
2591*59599516SKenneth E. Jansen	ttim(30) = ttim(30) + tmr()
2592*59599516SKenneth E. Jansen
2593*59599516SKenneth E. Jansen        return
2594*59599516SKenneth E. Jansen        end
2595*59599516SKenneth E. Jansen
2596*59599516SKenneth E. Jansen
2597*59599516SKenneth E. Jansen        subroutine e3bdgSclr  (	shp, 	shg, 	WdetJ,
2598*59599516SKenneth E. Jansen     &				A1t, 	A2t, 	A3t, 	A0t,
2599*59599516SKenneth E. Jansen     &                        	taut,
2600*59599516SKenneth E. Jansen     &   			Sclr,
2601*59599516SKenneth E. Jansen     &				u1,  	u2,  	u3,
2602*59599516SKenneth E. Jansen     &				rmu,    BDiagtl)
2603*59599516SKenneth E. Jansen
2604*59599516SKenneth E. Jansen        include "common.h"
2605*59599516SKenneth E. Jansen
2606*59599516SKenneth E. Jansenc
2607*59599516SKenneth E. Jansenc  passed arrays
2608*59599516SKenneth E. Jansenc
2609*59599516SKenneth E. Jansen        dimension BDiagtl(npro,nshl),
2610*59599516SKenneth E. Jansen     &            shp(npro,nshl),      Sclr(npro),
2611*59599516SKenneth E. Jansen     &            shg(npro,nshl,nsd),  WdetJ(npro),
2612*59599516SKenneth E. Jansen     &            A1tautA0(npro),      A2tautA0(npro),
2613*59599516SKenneth E. Jansen     &            A3tautA0(npro),      Ataut(npro),
2614*59599516SKenneth E. Jansen     &            A1t(npro),           A2t(npro),
2615*59599516SKenneth E. Jansen     &            A3t(npro),           A0t(npro),
2616*59599516SKenneth E. Jansen     &            taut(npro),          u1(npro),
2617*59599516SKenneth E. Jansen     &            u2(npro),            u3(npro),
2618*59599516SKenneth E. Jansen     &            rmu(npro)
2619*59599516SKenneth E. Jansenc
2620*59599516SKenneth E. Jansenc
2621*59599516SKenneth E. Jansenc  passed work arrays for local variables
2622*59599516SKenneth E. Jansenc
2623*59599516SKenneth E. Jansen	dimension tmp1(npro),             tmp2(npro)
2624*59599516SKenneth E. Jansen	dimension tmp(npro)
2625*59599516SKenneth E. Jansen
2626*59599516SKenneth E. Jansen	ttim(30) = ttim(30) - tmr()
2627*59599516SKenneth E. Jansen
2628*59599516SKenneth E. Jansenc $$ Ex-E3conv
2629*59599516SKenneth E. Jansenc
2630*59599516SKenneth E. Jansenc.... calculate the contribution in non-integrated by part form
2631*59599516SKenneth E. Jansenc.... add (N_b A_tilde_i N_b,i) to BDiagl
2632*59599516SKenneth E. Jansenc
2633*59599516SKenneth E. Jansen       do j = 1, nshl   ! May be worth eliminating zeros in A(prim) matrices
2634*59599516SKenneth E. Jansen        tmp=shp(:,j)*WdetJ
2635*59599516SKenneth E. Jansen            BDiagtl(:,j) = BDiagtl(:,j)
2636*59599516SKenneth E. Jansen     &                          + tmp * (
2637*59599516SKenneth E. Jansen     &                                   shg(:,j,1) * A1t(:) +
2638*59599516SKenneth E. Jansen     &                                   shg(:,j,2) * A2t(:) +
2639*59599516SKenneth E. Jansen     &                                   shg(:,j,3) * A3t(:)
2640*59599516SKenneth E. Jansen     &                                                            )
2641*59599516SKenneth E. Jansen
2642*59599516SKenneth E. Jansen       enddo
2643*59599516SKenneth E. Jansen        if (ngauss .eq. 1 .and. nshl .eq. 4) then  ! Exact integ of mass term for tets
2644*59599516SKenneth E. Jansen
2645*59599516SKenneth E. Jansenc $$ Ex-e3mael
2646*59599516SKenneth E. Jansenc
2647*59599516SKenneth E. Jansenc.... calculate factor   V             * sum(column of M)
2648*59599516SKenneth E. Jansenc                     =WdetJ/(3/4)Qwt(lcsyst,intp)* (2+1+1+1)/20
2649*59599516SKenneth E. Jansenc
2650*59599516SKenneth E. Jansen	tmp =  0.4d0 * WdetJ * Dtgl/Qwt(lcsyst,intp)/three*almi/gami/alfi
2651*59599516SKenneth E. Jansenc
2652*59599516SKenneth E. Jansen        do j=1,nshl  ! take advantage of zeros in A0t(Prim)
2653*59599516SKenneth E. Jansen	  BDiagtl(:,j) = BDiagtl(:,j) + tmp * A0t(:)
2654*59599516SKenneth E. Jansen	enddo
2655*59599516SKenneth E. Jansen
2656*59599516SKenneth E. Jansen        else
2657*59599516SKenneth E. Jansen
2658*59599516SKenneth E. Jansenc $$ Ex-e3mass
2659*59599516SKenneth E. Jansenc
2660*59599516SKenneth E. Jansen        ff = almi / gami / alfi
2661*59599516SKenneth E. Jansen        tmp = WdetJ * Dtgl * ff
2662*59599516SKenneth E. Jansenc
2663*59599516SKenneth E. Jansen        do j = 1, nshl
2664*59599516SKenneth E. Jansenc
2665*59599516SKenneth E. Jansen          tmp2 = (shp(:,j)*shp(:,j)) * tmp
2666*59599516SKenneth E. Jansenc
2667*59599516SKenneth E. Jansen	  BDiagtl(:,j) = BDiagtl(:,j) + tmp2 * A0t(:)
2668*59599516SKenneth E. Jansen
2669*59599516SKenneth E. Jansen        enddo
2670*59599516SKenneth E. Jansen
2671*59599516SKenneth E. Jansen        endif
2672*59599516SKenneth E. Jansen
2673*59599516SKenneth E. Jansenc
2674*59599516SKenneth E. Jansenc.... calculate (Atau) <-- (A_1 tau)
2675*59599516SKenneth E. Jansenc
2676*59599516SKenneth E. Jansenc
2677*59599516SKenneth E. Jansen       if (ivart .ge. 2) then
2678*59599516SKenneth E. Jansen       do i = 1, nflow
2679*59599516SKenneth E. Jansen          Ataut(:) = A1t(:)*taut(:)
2680*59599516SKenneth E. Jansen       enddo
2681*59599516SKenneth E. Jansenc
2682*59599516SKenneth E. Jansenc.... calculate (A_1 tau A_0) (for L.S. time term of EGmass)
2683*59599516SKenneth E. Jansenc
2684*59599516SKenneth E. Jansen             A1tautA0(:) = Ataut(:)*A0t(:)
2685*59599516SKenneth E. Jansenc
2686*59599516SKenneth E. Jansenc.... calculate (Atau) <-- (A_2 tau)
2687*59599516SKenneth E. Jansenc
2688*59599516SKenneth E. Jansenc
2689*59599516SKenneth E. Jansen          Ataut(:) = A2t(:)*taut(:)
2690*59599516SKenneth E. Jansen
2691*59599516SKenneth E. Jansenc
2692*59599516SKenneth E. Jansenc.... calculate (A_2 tau A_0) (for L.S. time term of EGmass)
2693*59599516SKenneth E. Jansenc
2694*59599516SKenneth E. Jansen
2695*59599516SKenneth E. Jansen             A2tautA0(:) = Ataut(:)*A0t(:)
2696*59599516SKenneth E. Jansenc
2697*59599516SKenneth E. Jansenc.... calculate (Atau) <-- (A_3 tau)
2698*59599516SKenneth E. Jansenc
2699*59599516SKenneth E. Jansenc
2700*59599516SKenneth E. Jansen
2701*59599516SKenneth E. Jansen          Ataut(:) = A3t(:)*taut(:)
2702*59599516SKenneth E. Jansenc
2703*59599516SKenneth E. Jansenc.... calculate (A_3 tau A_0) (for L.S. time term of EGmass)
2704*59599516SKenneth E. Jansenc
2705*59599516SKenneth E. Jansen
2706*59599516SKenneth E. Jansen             A3tautA0(:) = Ataut(:)*A0t(:)
2707*59599516SKenneth E. Jansen
2708*59599516SKenneth E. Jansen
2709*59599516SKenneth E. Jansenc.... add least squares time term to BDiagl
2710*59599516SKenneth E. Jansenc
2711*59599516SKenneth E. Jansenc
2712*59599516SKenneth E. Jansenc.... loop through rows (nodes i)
2713*59599516SKenneth E. Jansenc
2714*59599516SKenneth E. Jansen       do i = 1, nshl
2715*59599516SKenneth E. Jansenc
2716*59599516SKenneth E. Jansenc.... loop through column nodes, add (N_a,i A_i tau N_b) to BDiagl
2717*59599516SKenneth E. Jansenc
2718*59599516SKenneth E. Jansen          tmp1 = shp(:,i) * WdetJ * almi/gami/alfi*dtgl
2719*59599516SKenneth E. Jansen
2720*59599516SKenneth E. Jansen                BDiagtl(:,i) = BDiagtl(:,i)  + tmp1*(
2721*59599516SKenneth E. Jansen     &               shg(:,i,1) * A1tautA0(:) +
2722*59599516SKenneth E. Jansen     &               shg(:,i,2) * A2tautA0(:) +
2723*59599516SKenneth E. Jansen     &               shg(:,i,3) * A3tautA0(:))
2724*59599516SKenneth E. Jansen
2725*59599516SKenneth E. Jansen       enddo
2726*59599516SKenneth E. Jansen       endif
2727*59599516SKenneth E. Jansen
2728*59599516SKenneth E. Jansenc
2729*59599516SKenneth E. Jansenc $$ Ex-e3wmlt
2730*59599516SKenneth E. Jansenc
2731*59599516SKenneth E. Jansenc.... add contribution of stiffness to BDiagl
2732*59599516SKenneth E. Jansenc
2733*59599516SKenneth E. Jansenc....  we no longer build stiff so we have to get it in here.
2734*59599516SKenneth E. Jansenc      recall that this is the (A_i tau A_j + K_{ij}) that
2735*59599516SKenneth E. Jansenc      multiplies N_{a,i} N_{b,j} (Actually for BDiag a=b).
2736*59599516SKenneth E. Jansenc
2737*59599516SKenneth E. Jansenc...   we are through with A0 so use it now for the sub-blocks
2738*59599516SKenneth E. Jansenc      of what used to be stiff
2739*59599516SKenneth E. Jansenc
2740*59599516SKenneth E. Jansenc
2741*59599516SKenneth E. Jansen
2742*59599516SKenneth E. Jansenc
2743*59599516SKenneth E. Jansenc  start with 1 1 contribution
2744*59599516SKenneth E. Jansenc
2745*59599516SKenneth E. Jansen        if(ivart .ge. 2) then
2746*59599516SKenneth E. Jansen
2747*59599516SKenneth E. Jansen
2748*59599516SKenneth E. Jansen          A0t(:)= A1t(:)*taut(:) *A1t(:)
2749*59599516SKenneth E. Jansen        else
2750*59599516SKenneth E. Jansen          A0t=zero
2751*59599516SKenneth E. Jansen        endif
2752*59599516SKenneth E. Jansenc
2753*59599516SKenneth E. Jansen        if (Navier .eq. 1) then
2754*59599516SKenneth E. Jansenc.... K11
2755*59599516SKenneth E. Jansenc
2756*59599516SKenneth E. Jansen           A0t(:) = A0t(:) +1.5*(rmu+Sclr)
2757*59599516SKenneth E. Jansen
2758*59599516SKenneth E. Jansen        endif
2759*59599516SKenneth E. Jansen        do j = 1, nshl
2760*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,1) * shg(:,j,1)
2761*59599516SKenneth E. Jansen
2762*59599516SKenneth E. Jansen             BDiagtl(:,j) = BDiagtl(:,j) + tmp * A0t(:)
2763*59599516SKenneth E. Jansen        enddo
2764*59599516SKenneth E. Jansenc
2765*59599516SKenneth E. Jansenc  start with 2 2 contribution
2766*59599516SKenneth E. Jansenc
2767*59599516SKenneth E. Jansen        if(ivart .ge. 2) then
2768*59599516SKenneth E. Jansen
2769*59599516SKenneth E. Jansen
2770*59599516SKenneth E. Jansen          A0t(:)= A2t(:)*taut(:) *A2t(:)
2771*59599516SKenneth E. Jansen        else
2772*59599516SKenneth E. Jansen          A0t=zero
2773*59599516SKenneth E. Jansen        endif
2774*59599516SKenneth E. Jansenc
2775*59599516SKenneth E. Jansen        if (Navier .eq. 1) then
2776*59599516SKenneth E. Jansenc.... K22
2777*59599516SKenneth E. Jansenc
2778*59599516SKenneth E. Jansen           A0t(:) = A0t(:) +1.5*(rmu+Sclr)
2779*59599516SKenneth E. Jansen
2780*59599516SKenneth E. Jansen        endif
2781*59599516SKenneth E. Jansen        do j = 1, nshl
2782*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,2) * shg(:,j,2)
2783*59599516SKenneth E. Jansen
2784*59599516SKenneth E. Jansen             BDiagtl(:,j) = BDiagtl(:,j) + tmp * A0t(:)
2785*59599516SKenneth E. Jansen        enddo
2786*59599516SKenneth E. Jansen
2787*59599516SKenneth E. Jansenc
2788*59599516SKenneth E. Jansenc  start with 3 3 contribution
2789*59599516SKenneth E. Jansenc
2790*59599516SKenneth E. Jansen        if(ivart .ge. 2) then
2791*59599516SKenneth E. Jansen
2792*59599516SKenneth E. Jansen
2793*59599516SKenneth E. Jansen          A0t(:)= A3t(:)*taut(:) *A3t(:)
2794*59599516SKenneth E. Jansen        else
2795*59599516SKenneth E. Jansen          A0t=zero
2796*59599516SKenneth E. Jansen        endif
2797*59599516SKenneth E. Jansenc
2798*59599516SKenneth E. Jansen        if (Navier .eq. 1) then
2799*59599516SKenneth E. Jansenc.... K33
2800*59599516SKenneth E. Jansenc
2801*59599516SKenneth E. Jansen           A0t(:) = A0t(:) +1.5*(rmu+Sclr)
2802*59599516SKenneth E. Jansen
2803*59599516SKenneth E. Jansen        endif
2804*59599516SKenneth E. Jansen        do j = 1, nshl
2805*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,3) * shg(:,j,3)
2806*59599516SKenneth E. Jansen
2807*59599516SKenneth E. Jansen             BDiagtl(:,j) = BDiagtl(:,j) + tmp * A0t(:)
2808*59599516SKenneth E. Jansen        enddo
2809*59599516SKenneth E. Jansen
2810*59599516SKenneth E. Jansen
2811*59599516SKenneth E. Jansen
2812*59599516SKenneth E. Jansenc
2813*59599516SKenneth E. Jansenc now for the 1 2  plus 2 1
2814*59599516SKenneth E. Jansenc
2815*59599516SKenneth E. Jansen        if(ivart.ge.2) then
2816*59599516SKenneth E. Jansen
2817*59599516SKenneth E. Jansen          A0t(:)= taut(:) * two * A1t(:) * A2t(:)
2818*59599516SKenneth E. Jansen
2819*59599516SKenneth E. Jansen        else
2820*59599516SKenneth E. Jansen           A0t = zero
2821*59599516SKenneth E. Jansen        endif
2822*59599516SKenneth E. Jansen
2823*59599516SKenneth E. Jansen        do j = 1, nshl
2824*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,1) * shg(:,j,2)
2825*59599516SKenneth E. Jansen
2826*59599516SKenneth E. Jansen             BDiagtl(:,j) = BDiagtl(:,j) + tmp * A0t(:)
2827*59599516SKenneth E. Jansen        enddo
2828*59599516SKenneth E. Jansenc
2829*59599516SKenneth E. Jansenc now for the 1 3  plus 3 1
2830*59599516SKenneth E. Jansenc
2831*59599516SKenneth E. Jansen        if(ivart.ge.2) then
2832*59599516SKenneth E. Jansen          A0t(:)= taut(:) * two * A1t(:) * A3t(:)
2833*59599516SKenneth E. Jansen
2834*59599516SKenneth E. Jansen        else
2835*59599516SKenneth E. Jansen           A0t = zero
2836*59599516SKenneth E. Jansen        endif
2837*59599516SKenneth E. Jansen
2838*59599516SKenneth E. Jansen        do j = 1, nshl
2839*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,1) * shg(:,j,3)
2840*59599516SKenneth E. Jansen
2841*59599516SKenneth E. Jansen             BDiagtl(:,j) = BDiagtl(:,j) + tmp * A0t(:)
2842*59599516SKenneth E. Jansen        enddo
2843*59599516SKenneth E. Jansenc
2844*59599516SKenneth E. Jansenc now for the 2 3  plus 3 2
2845*59599516SKenneth E. Jansenc
2846*59599516SKenneth E. Jansen        if(ivart.ge.2) then
2847*59599516SKenneth E. Jansen
2848*59599516SKenneth E. Jansen          A0t(:)= taut(:) * two * A2t(:) * A3t(:)
2849*59599516SKenneth E. Jansen
2850*59599516SKenneth E. Jansen        else
2851*59599516SKenneth E. Jansen           A0t = zero
2852*59599516SKenneth E. Jansen        endif
2853*59599516SKenneth E. Jansen
2854*59599516SKenneth E. Jansen        do j = 1, nshl
2855*59599516SKenneth E. Jansen          tmp = WdetJ * shg(:,j,2) * shg(:,j,3)
2856*59599516SKenneth E. Jansen
2857*59599516SKenneth E. Jansen             BDiagtl(:,j) = BDiagtl(:,j) + tmp * A0t(:)
2858*59599516SKenneth E. Jansen        enddo
2859*59599516SKenneth E. Jansen
2860*59599516SKenneth E. Jansen
2861*59599516SKenneth E. Jansen	ttim(30) = ttim(30) + tmr()
2862*59599516SKenneth E. Jansen
2863*59599516SKenneth E. Jansen        return
2864*59599516SKenneth E. Jansen        end
2865