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