1 subroutine e3bdg_nd (shp, shg, WdetJ, 2 & A1, A2, A3, 3 & A0, bcool, PTau, 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), Atau1(npro,nflow,nflow), 17 & A1(npro,nflow,nflow), A2(npro,nflow,nflow), 18 & A3(npro,nflow,nflow), A0(npro,nflow,nflow), 19 & PTau(npro,5,5), u1(npro), 20 & u2(npro), u3(npro), 21 & rlm2mu(npro), Atau2(npro,nflow,nflow), 22 & rmu(npro), con(npro), 23 & Atau3(npro,nflow,nflow) 24c 25c 26c passed work arrays for local variables 27c 28 dimension tmp1(npro), tmp2(npro) 29 dimension tmp(npro) 30 31 ttim(30) = ttim(30) - secs(0.0) 32 33c $$ Ex-E3conv 34c 35c.... calculate the contribution in non-integrated by part form 36c.... add (N_b A_tilde_i N_b,i) to BDiagl 37c 38 do j = 1, nshl ! May be worth eliminating zeros in A(prim) matrices 39 tmp=shp(:,j)*WdetJ 40 BDiagl(:,j,1,1) = BDiagl(:,j,1,1) 41 & + tmp * ( 42 & shg(:,j,1) * A1(:,1,1) + 43 & shg(:,j,2) * A2(:,1,1) + 44 & shg(:,j,3) * A3(:,1,1) 45 & ) 46 BDiagl(:,j,1,2) = BDiagl(:,j,1,2) 47 & + tmp * ( 48 & shg(:,j,1) * A1(:,1,2) 49c & +shg(:,j,2) * A2(:,1,2) 50c & +shg(:,j,3) * A3(:,1,2) 51 & ) 52 BDiagl(:,j,1,3) = BDiagl(:,j,1,3) 53 & + tmp * ( 54c & shg(:,j,1) * A1(:,1,3) 55 & +shg(:,j,2) * A2(:,1,3) 56c & +shg(:,j,3) * A3(:,1,3) 57 & ) 58 BDiagl(:,j,1,4) = BDiagl(:,j,1,4) 59 & + tmp * ( 60c & shg(:,j,1) * A1(:,1,4) + 61c & shg(:,j,2) * A2(:,1,4) + 62 & shg(:,j,3) * A3(:,1,4) 63 & ) 64 BDiagl(:,j,1,5) = BDiagl(:,j,1,5) 65 & + tmp * ( 66 & shg(:,j,1) * A1(:,1,5) + 67 & shg(:,j,2) * A2(:,1,5) + 68 & shg(:,j,3) * A3(:,1,5) 69 & ) 70 BDiagl(:,j,2,1) = BDiagl(:,j,2,1) 71 & + tmp * ( 72 & shg(:,j,1) * A1(:,2,1) + 73 & shg(:,j,2) * A2(:,2,1) + 74 & shg(:,j,3) * A3(:,2,1) 75 & ) 76 BDiagl(:,j,2,2) = BDiagl(:,j,2,2) 77 & + tmp * ( 78 & shg(:,j,1) * A1(:,2,2) + 79 & shg(:,j,2) * A2(:,2,2) + 80 & shg(:,j,3) * A3(:,2,2) 81 & ) 82 BDiagl(:,j,2,3) = BDiagl(:,j,2,3) 83 & + tmp * ( 84c & shg(:,j,1) * A1(:,2,3) 85 & +shg(:,j,2) * A2(:,2,3) 86c & +shg(:,j,3) * A3(:,2,3) 87 & ) 88 BDiagl(:,j,2,4) = BDiagl(:,j,2,4) 89 & + tmp * ( 90c & shg(:,j,1) * A1(:,2,4) + 91c & shg(:,j,2) * A2(:,2,4) + 92 & shg(:,j,3) * A3(:,2,4) 93 & ) 94 BDiagl(:,j,2,5) = BDiagl(:,j,2,5) 95 & + tmp * ( 96 & shg(:,j,1) * A1(:,2,5) + 97 & shg(:,j,2) * A2(:,2,5) + 98 & shg(:,j,3) * A3(:,2,5) 99 & ) 100 BDiagl(:,j,3,1) = BDiagl(:,j,3,1) 101 & + tmp * ( 102 & shg(:,j,1) * A1(:,3,1) + 103 & shg(:,j,2) * A2(:,3,1) + 104 & shg(:,j,3) * A3(:,3,1) 105 & ) 106 BDiagl(:,j,3,2) = BDiagl(:,j,3,2) 107 & + tmp * ( 108 & shg(:,j,1) * A1(:,3,2) 109c & +shg(:,j,2) * A2(:,3,2) 110c & +shg(:,j,3) * A3(:,3,2) 111 & ) 112 BDiagl(:,j,3,3) = BDiagl(:,j,3,3) 113 & + tmp * ( 114 & shg(:,j,1) * A1(:,3,3) + 115 & shg(:,j,2) * A2(:,3,3) + 116 & shg(:,j,3) * A3(:,3,3) 117 & ) 118 BDiagl(:,j,3,4) = BDiagl(:,j,3,4) 119 & + tmp * ( 120c & shg(:,j,1) * A1(:,3,4) + 121c & shg(:,j,2) * A2(:,3,4) + 122 & shg(:,j,3) * A3(:,3,4) 123 & ) 124 BDiagl(:,j,3,5) = BDiagl(:,j,3,5) 125 & + tmp * ( 126 & shg(:,j,1) * A1(:,3,5) + 127 & shg(:,j,2) * A2(:,3,5) + 128 & shg(:,j,3) * A3(:,3,5) 129 & ) 130 BDiagl(:,j,4,1) = BDiagl(:,j,4,1) 131 & + tmp * ( 132 & shg(:,j,1) * A1(:,4,1) + 133 & shg(:,j,2) * A2(:,4,1) + 134 & shg(:,j,3) * A3(:,4,1) 135 & ) 136 BDiagl(:,j,4,2) = BDiagl(:,j,4,2) 137 & + tmp * ( 138 & shg(:,j,1) * A1(:,4,2) 139c & +shg(:,j,2) * A2(:,4,2) 140c & +shg(:,j,3) * A3(:,4,2) 141 & ) 142 BDiagl(:,j,4,3) = BDiagl(:,j,4,3) 143 & + tmp * ( 144c & shg(:,j,1) * A1(:,4,3) 145 & +shg(:,j,2) * A2(:,4,3) 146c & +shg(:,j,3) * A3(:,4,3) 147 & ) 148 BDiagl(:,j,4,4) = BDiagl(:,j,4,4) 149 & + tmp * ( 150 & shg(:,j,1) * A1(:,4,4) + 151 & shg(:,j,2) * A2(:,4,4) + 152 & shg(:,j,3) * A3(:,4,4) 153 & ) 154 BDiagl(:,j,4,5) = BDiagl(:,j,4,5) 155 & + tmp * ( 156 & shg(:,j,1) * A1(:,4,5) + 157 & shg(:,j,2) * A2(:,4,5) + 158 & shg(:,j,3) * A3(:,4,5) 159 & ) 160 BDiagl(:,j,5,1) = BDiagl(:,j,5,1) 161 & + tmp * ( 162 & shg(:,j,1) * A1(:,5,1) + 163 & shg(:,j,2) * A2(:,5,1) + 164 & shg(:,j,3) * A3(:,5,1) 165 & ) 166 BDiagl(:,j,5,2) = BDiagl(:,j,5,2) 167 & + tmp * ( 168 & shg(:,j,1) * A1(:,5,2) + 169 & shg(:,j,2) * A2(:,5,2) + 170 & shg(:,j,3) * A3(:,5,2) 171 & ) 172 BDiagl(:,j,5,3) = BDiagl(:,j,5,3) 173 & + tmp * ( 174 & shg(:,j,1) * A1(:,5,3) + 175 & shg(:,j,2) * A2(:,5,3) + 176 & shg(:,j,3) * A3(:,5,3) 177 & ) 178 BDiagl(:,j,5,4) = BDiagl(:,j,5,4) 179 & + tmp * ( 180 & shg(:,j,1) * A1(:,5,4) + 181 & shg(:,j,2) * A2(:,5,4) + 182 & shg(:,j,3) * A3(:,5,4) 183 & ) 184 BDiagl(:,j,5,5) = BDiagl(:,j,5,5) 185 & + tmp * ( 186 & shg(:,j,1) * A1(:,5,5) + 187 & shg(:,j,2) * A2(:,5,5) + 188 & shg(:,j,3) * A3(:,5,5) 189 & ) 190 enddo 191 if (ngauss .eq. 1 .and. nshl .eq. 4) then ! Exact integ of mass term for tets 192 193c $$ Ex-e3mael 194c 195c.... calculate factor V * sum(column of M) 196c =WdetJ/(3/4)Qwt(lcsyst,intp) * (2+1+1+1)/20 197c 198 tmp = 0.4d0 * WdetJ * Dtgl/Qwt(lcsyst,intp)/three*almi/gami/alfi 199c 200 do j=1,nshl ! take advantage of zeros in A0(Prim) 201 BDiagl(:,j,2,2) = BDiagl(:,j,2,2) + tmp * A0(:,2,2) 202 BDiagl(:,j,3,3) = BDiagl(:,j,3,3) + tmp * A0(:,3,3) 203 BDiagl(:,j,4,4) = BDiagl(:,j,4,4) + tmp * A0(:,4,4) 204 BDiagl(:,j,1,5) = BDiagl(:,j,1,5) + tmp * A0(:,1,5) 205 BDiagl(:,j,2,5) = BDiagl(:,j,2,5) + tmp * A0(:,2,5) 206 BDiagl(:,j,3,5) = BDiagl(:,j,3,5) + tmp * A0(:,3,5) 207 BDiagl(:,j,4,5) = BDiagl(:,j,4,5) + tmp * A0(:,4,5) 208 BDiagl(:,j,1,1) = BDiagl(:,j,1,1) + tmp * A0(:,1,1) 209 BDiagl(:,j,2,1) = BDiagl(:,j,2,1) + tmp * A0(:,2,1) 210 BDiagl(:,j,3,1) = BDiagl(:,j,3,1) + tmp * A0(:,3,1) 211 BDiagl(:,j,4,1) = BDiagl(:,j,4,1) + tmp * A0(:,4,1) 212 BDiagl(:,j,5,1) = BDiagl(:,j,5,1) + tmp * A0(:,5,1) 213 BDiagl(:,j,5,2) = BDiagl(:,j,5,2) + tmp * A0(:,5,2) 214 BDiagl(:,j,5,3) = BDiagl(:,j,5,3) + tmp * A0(:,5,3) 215 BDiagl(:,j,5,4) = BDiagl(:,j,5,4) + tmp * A0(:,5,4) 216 BDiagl(:,j,5,5) = BDiagl(:,j,5,5) + tmp * A0(:,5,5) 217 enddo 218 219 else 220 221c $$ Ex-e3mass 222c 223 ff = almi / gami / alfi 224 tmp = WdetJ * (Dtgl * ff + bcool) 225c 226 do j = 1, nshl 227c 228 tmp2 = (shp(:,j)*shp(:,j)) * tmp 229c 230 BDiagl(:,j,2,2) = BDiagl(:,j,2,2) + tmp2 * A0(:,2,2) 231 BDiagl(:,j,3,3) = BDiagl(:,j,3,3) + tmp2 * A0(:,3,3) 232 BDiagl(:,j,4,4) = BDiagl(:,j,4,4) + tmp2 * A0(:,4,4) 233 BDiagl(:,j,1,5) = BDiagl(:,j,1,5) + tmp2 * A0(:,1,5) 234 BDiagl(:,j,2,5) = BDiagl(:,j,2,5) + tmp2 * A0(:,2,5) 235 BDiagl(:,j,3,5) = BDiagl(:,j,3,5) + tmp2 * A0(:,3,5) 236 BDiagl(:,j,4,5) = BDiagl(:,j,4,5) + tmp2 * A0(:,4,5) 237 BDiagl(:,j,1,1) = BDiagl(:,j,1,1) + tmp2 * A0(:,1,1) 238 BDiagl(:,j,2,1) = BDiagl(:,j,2,1) + tmp2 * A0(:,2,1) 239 BDiagl(:,j,3,1) = BDiagl(:,j,3,1) + tmp2 * A0(:,3,1) 240 BDiagl(:,j,4,1) = BDiagl(:,j,4,1) + tmp2 * A0(:,4,1) 241 BDiagl(:,j,5,1) = BDiagl(:,j,5,1) + tmp2 * A0(:,5,1) 242 BDiagl(:,j,5,2) = BDiagl(:,j,5,2) + tmp2 * A0(:,5,2) 243 BDiagl(:,j,5,3) = BDiagl(:,j,5,3) + tmp2 * A0(:,5,3) 244 BDiagl(:,j,5,4) = BDiagl(:,j,5,4) + tmp2 * A0(:,5,4) 245 BDiagl(:,j,5,5) = BDiagl(:,j,5,5) + tmp2 * A0(:,5,5) 246c 247 enddo 248 249 endif 250 251c 252c.... calculate (Atau) <-- (A_1 tau) (Recall that we are using a 253c non-diagonal tau here) 254c 255 if (ivart .ge. 2) then 256 Atau1 = zero 257 do i = 1, nflow 258 do j = 1, nflow 259 do k = 1, nflow 260 Atau1(:,i,j) = Atau1(:,i,j) + A1(:,i,k)*PTau(:,k,j) 261 enddo 262 enddo 263 enddo 264c 265c.... calculate (A_1 tau A_0) (for L.S. time term of EGmass) 266c 267 do j = 1, nflow 268 do i = 1, nflow 269 A1tauA0(:,i,j) = 270 & Atau1(:,i,1)*A0(:,1,j) + 271 & Atau1(:,i,2)*A0(:,2,j) + 272 & Atau1(:,i,3)*A0(:,3,j) + 273 & Atau1(:,i,4)*A0(:,4,j) + 274 & Atau1(:,i,5)*A0(:,5,j) 275 enddo 276 enddo 277c 278c.... calculate (Atau) <-- (A_2 tau) (Recall that we are using a 279c non-diagonal tau here) 280c 281 Atau2 = zero 282 do i = 1, nflow 283 do j = 1, nflow 284 do k = 1, nflow 285 Atau2(:,i,j) = Atau2(:,i,j) + A2(:,i,k)*PTau(:,k,j) 286 enddo 287 enddo 288 enddo 289 290c 291c.... calculate (A_2 tau A_0) (for L.S. time term of EGmass) 292c 293 do j = 1, nflow 294 do i = 1, nflow 295 A2tauA0(:,i,j) = 296 & Atau2(:,i,1)*A0(:,1,j) + 297 & Atau2(:,i,2)*A0(:,2,j) + 298 & Atau2(:,i,3)*A0(:,3,j) + 299 & Atau2(:,i,4)*A0(:,4,j) + 300 & Atau2(:,i,5)*A0(:,5,j) 301 enddo 302 enddo 303c 304c.... calculate (Atau) <-- (A_3 tau) (Recall that we are using a 305c non-diagonal tau here) 306 307 Atau3 = zero 308 do i = 1, nflow 309 do j = 1, nflow 310 do k = 1, nflow 311 Atau3(:,i,j) = Atau3(:,i,j) + A3(:,i,k)*PTau(:,k,j) 312 enddo 313 enddo 314 enddo 315 316 317c 318c.... calculate (A_3 tau A_0) (for L.S. time term of EGmass) 319c 320 do j = 1, nflow 321 do i = 1, nflow 322 A3tauA0(:,i,j) = 323 & Atau3(:,i,1)*A0(:,1,j) + 324 & Atau3(:,i,2)*A0(:,2,j) + 325 & Atau3(:,i,3)*A0(:,3,j) + 326 & Atau3(:,i,4)*A0(:,4,j) + 327 & Atau3(:,i,5)*A0(:,5,j) 328 enddo 329 enddo 330c 331c 332c.... add least squares time term to BDiagl 333c 334c 335c.... loop through rows (nodes i) 336c 337 do i = 1, nshl 338c 339c.... loop through column nodes, add (N_a,i A_i tau N_b) to BDiagl 340c 341 tmp1 = shp(:,i) * WdetJ * almi/gami/alfi*dtgl 342 do idof = 1, nflow 343 do jdof = 1, nflow 344 BDiagl(:,i,idof,jdof) = BDiagl(:,i,idof,jdof) 345 & + tmp1*( 346 & shg(:,i,1) * A1tauA0(:,idof,jdof) + 347 & shg(:,i,2) * A2tauA0(:,idof,jdof) + 348 & shg(:,i,3) * A3tauA0(:,idof,jdof)) 349 enddo 350 enddo 351 enddo 352 endif 353c 354c.... add contribution of stiffness to BDiagl 355c 356c.... we no longer build stiff so we have to get it in here. 357c recall that this is the (A_i tau A_j + K_{ij}) that 358c multiplies N_{a,i} N_{b,j} (Actually for BDiag a=b). 359c 360c... we are through with A0 so use it now for the sub-blocks 361c of what used to be stiff 362c 363c 364 365c 366c start with 1 1 contribution 367c 368c 369c Note that I comment out lines where zeros appear (for Pvariables) 370c 371 if(ivart .ge. 2) then 372c 373c row one 374c 375 A0(:,1,1) = 376 & Atau1(:,1,1) * A1(:,1,1) 377 & + Atau1(:,1,2) * A1(:,2,1) 378 & + Atau1(:,1,3) * A1(:,3,1) 379 & + Atau1(:,1,4) * A1(:,4,1) 380 & + Atau1(:,1,5) * A1(:,5,1) 381c 382 A0(:,1,2) = 383 & Atau1(:,1,1) * A1(:,1,2) 384 & + Atau1(:,1,2) * A1(:,2,2) 385 & + Atau1(:,1,3) * A1(:,3,2) 386 & + Atau1(:,1,4) * A1(:,4,2) 387 & + Atau1(:,1,5) * A1(:,5,2) 388c 389 A0(:,1,3) = 390c & Atau1(:,1,1) * A1(:,1,3) 391c & + Atau1(:,1,2) * A1(:,2,3) 392 & Atau1(:,1,3) * A1(:,3,3) 393c & + Atau1(:,1,4) * A1(:,4,3) 394 & + Atau1(:,1,5) * A1(:,5,3) 395c 396 A0(:,1,4)= 397c & Atau1(:,1,1) * A1(:,1,4) 398c & + Atau1(:,1,2) * A1(:,2,4) 399c & + Atau1(:,1,3) * A1(:,3,4) 400 & Atau1(:,1,4) * A1(:,4,4) 401 & + Atau1(:,1,5) * A1(:,5,4) 402c 403 A0(:,1,5)= 404 & Atau1(:,1,1) * A1(:,1,5) 405 & + Atau1(:,1,2) * A1(:,2,5) 406 & + Atau1(:,1,3) * A1(:,3,5) 407 & + Atau1(:,1,4) * A1(:,4,5) 408 & + Atau1(:,1,5) * A1(:,5,5) 409c 410c row two 411c 412 A0(:,2,1)= 413 & Atau1(:,2,1) * A1(:,1,1) 414 & + Atau1(:,2,2) * A1(:,2,1) 415 & + Atau1(:,2,3) * A1(:,3,1) 416 & + Atau1(:,2,4) * A1(:,4,1) 417 & + Atau1(:,2,5) * A1(:,5,1) 418c 419 A0(:,2,2)= 420 & Atau1(:,2,1) * A1(:,1,2) 421 & + Atau1(:,2,2) * A1(:,2,2) 422 & + Atau1(:,2,3) * A1(:,3,2) 423 & + Atau1(:,2,4) * A1(:,4,2) 424 & + Atau1(:,2,5) * A1(:,5,2) 425c 426 A0(:,2,3)= 427c & Atau1(:,2,1) * A1(:,1,3) 428c & + Atau1(:,2,2) * A1(:,2,3) 429 & Atau1(:,2,3) * A1(:,3,3) 430c & + Atau1(:,2,4) * A1(:,4,3) 431 & + Atau1(:,2,5) * A1(:,5,3) 432c 433 A0(:,2,4)= 434c & Atau1(:,2,1) * A1(:,1,4) 435c & + Atau1(:,2,2) * A1(:,2,4) 436c & + Atau1(:,2,3) * A1(:,3,4) 437 & Atau1(:,2,4) * A1(:,4,4) 438 & + Atau1(:,2,5) * A1(:,5,4) 439c 440 A0(:,2,5)= 441 & Atau1(:,2,1) * A1(:,1,5) 442 & + Atau1(:,2,2) * A1(:,2,5) 443 & + Atau1(:,2,3) * A1(:,3,5) 444 & + Atau1(:,2,4) * A1(:,4,5) 445 & + Atau1(:,2,5) * A1(:,5,5) 446c 447c row three 448c 449 A0(:,3,1)= 450 & Atau1(:,3,1) * A1(:,1,1) 451 & + Atau1(:,3,2) * A1(:,2,1) 452 & + Atau1(:,3,3) * A1(:,3,1) 453 & + Atau1(:,3,4) * A1(:,4,1) 454 & + Atau1(:,3,5) * A1(:,5,1) 455c 456 A0(:,3,2)= 457 & Atau1(:,3,1) * A1(:,1,2) 458 & + Atau1(:,3,2) * A1(:,2,2) 459 & + Atau1(:,3,3) * A1(:,3,2) 460 & + Atau1(:,3,4) * A1(:,4,2) 461 & + Atau1(:,3,5) * A1(:,5,2) 462c 463 A0(:,3,3)= 464c & Atau1(:,2,1) * A1(:,1,3) 465c & + Atau1(:,2,2) * A1(:,2,3) 466 & Atau1(:,3,3) * A1(:,3,3) 467c & + Atau1(:,2,4) * A1(:,4,3) 468 & + Atau1(:,3,5) * A1(:,5,3) 469c 470 A0(:,3,4)= 471c & Atau1(:,2,1) * A1(:,1,4) 472c & + Atau1(:,2,2) * A1(:,2,4) 473c & + Atau1(:,2,3) * A1(:,3,4) 474 & Atau1(:,3,4) * A1(:,4,4) 475 & + Atau1(:,3,5) * A1(:,5,4) 476c 477 A0(:,3,5)= 478 & Atau1(:,3,1) * A1(:,1,5) 479 & + Atau1(:,3,2) * A1(:,2,5) 480 & + Atau1(:,3,3) * A1(:,3,5) 481 & + Atau1(:,3,4) * A1(:,4,5) 482 & + Atau1(:,3,5) * A1(:,5,5) 483 484c 485c row four 486c 487 A0(:,4,1)= 488 & Atau1(:,4,1) * A1(:,1,1) 489 & + Atau1(:,4,2) * A1(:,2,1) 490 & + Atau1(:,4,3) * A1(:,3,1) 491 & + Atau1(:,4,4) * A1(:,4,1) 492 & + Atau1(:,4,5) * A1(:,5,1) 493c 494 A0(:,4,2)= 495 & Atau1(:,4,1) * A1(:,1,2) 496 & + Atau1(:,4,2) * A1(:,2,2) 497 & + Atau1(:,4,3) * A1(:,3,2) 498 & + Atau1(:,4,4) * A1(:,4,2) 499 & + Atau1(:,4,5) * A1(:,5,2) 500c 501 A0(:,4,3)= 502c & Atau1(:,2,1) * A1(:,1,3) 503c & + Atau1(:,2,2) * A1(:,2,3) 504 & Atau1(:,4,3) * A1(:,3,3) 505c & + Atau1(:,2,4) * A1(:,4,3) 506 & + Atau1(:,4,5) * A1(:,5,3) 507c 508 A0(:,4,4)= 509c & Atau1(:,2,1) * A1(:,1,4) 510c & + Atau1(:,2,2) * A1(:,2,4) 511c & + Atau1(:,2,3) * A1(:,3,4) 512 & Atau1(:,4,4) * A1(:,4,4) 513 & + Atau1(:,4,5) * A1(:,5,4) 514c 515 A0(:,4,5)= 516 & Atau1(:,4,1) * A1(:,1,5) 517 & + Atau1(:,4,2) * A1(:,2,5) 518 & + Atau1(:,4,3) * A1(:,3,5) 519 & + Atau1(:,4,4) * A1(:,4,5) 520 & + Atau1(:,4,5) * A1(:,5,5) 521c 522c row five 523c 524 525 A0(:,5,1)= 526 & Atau1(:,5,1) * A1(:,1,1) 527 & + Atau1(:,5,2) * A1(:,2,1) 528 & + Atau1(:,5,3) * A1(:,3,1) 529 & + Atau1(:,5,4) * A1(:,4,1) 530 & + Atau1(:,5,5) * A1(:,5,1) 531c 532 A0(:,5,2)= 533 & Atau1(:,5,1) * A1(:,1,2) 534 & + Atau1(:,5,2) * A1(:,2,2) 535 & + Atau1(:,5,3) * A1(:,3,2) 536 & + Atau1(:,5,4) * A1(:,4,2) 537 & + Atau1(:,5,5) * A1(:,5,2) 538c 539 A0(:,5,3)= 540c & Atau1(:,2,1) * A1(:,1,3) 541c & + Atau1(:,2,2) * A1(:,2,3) 542 & Atau1(:,5,3) * A1(:,3,3) 543c & + Atau1(:,2,4) * A1(:,4,3) 544 & + Atau1(:,3,5) * A1(:,5,3) 545c 546 A0(:,5,4)= 547c & Atau1(:,2,1) * A1(:,1,4) 548c & + Atau1(:,2,2) * A1(:,2,4) 549c & + Atau1(:,2,3) * A1(:,3,4) 550 & Atau1(:,5,4) * A1(:,4,4) 551 & + Atau1(:,5,5) * A1(:,5,4) 552c 553 A0(:,5,5)= 554 & Atau1(:,5,1) * A1(:,1,5) 555 & + Atau1(:,5,2) * A1(:,2,5) 556 & + Atau1(:,5,3) * A1(:,3,5) 557 & + Atau1(:,5,4) * A1(:,4,5) 558 & + Atau1(:,5,5) * A1(:,5,5) 559 560c 561 else 562 A0=zero 563 endif 564c 565 if (Navier .eq. 1) then 566 A0(:,2,2) = A0(:,2,2) + rlm2mu !rlm2mu 567 A0(:,3,3) = A0(:,3,3) + rmu !rmu 568 A0(:,4,4) = A0(:,4,4) + rmu !rmu 569 A0(:,5,2) = A0(:,5,2) + rlm2mu * u1 570 A0(:,5,3) = A0(:,5,3) + rmu * u2 571 A0(:,5,4) = A0(:,5,4) + rmu * u3 572 A0(:,5,5) = A0(:,5,5) + con !con 573 endif 574 do j = 1, nshl 575 tmp = WdetJ * shg(:,j,1) * shg(:,j,1) 576 do i=1,nflow 577 do k=1,nflow 578 BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k) 579 enddo 580 enddo 581 enddo 582c 583c start with 2 2 contribution 584c 585 if(ivart.ge.2) then 586 587c 588c row one 589c 590 A0(:,1,1) = 591 & Atau2(:,1,1) * A2(:,1,1) 592 & + Atau2(:,1,2) * A2(:,2,1) 593 & + Atau2(:,1,3) * A2(:,3,1) 594 & + Atau2(:,1,4) * A2(:,4,1) 595 & + Atau2(:,1,5) * A2(:,5,1) 596c 597 A0(:,1,2) = 598c & Atau2(:,1,1) * A2(:,1,2) 599 & Atau2(:,1,2) * A2(:,2,2) 600c & + Atau2(:,1,3) * A2(:,3,2) 601c & + Atau2(:,1,4) * A2(:,4,2) 602 & + Atau2(:,1,5) * A2(:,5,2) 603c 604 A0(:,1,3) = 605 & Atau2(:,1,1) * A2(:,1,3) 606 & + Atau2(:,1,2) * A2(:,2,3) 607 & + Atau2(:,1,3) * A2(:,3,3) 608 & + Atau2(:,1,4) * A2(:,4,3) 609 & + Atau2(:,1,5) * A2(:,5,3) 610c 611 A0(:,1,4)= 612c & Atau2(:,1,1) * A2(:,1,4) 613c & + Atau2(:,1,2) * A2(:,2,4) 614c & + Atau2(:,1,3) * A2(:,3,4) 615 & Atau2(:,1,4) * A2(:,4,4) 616 & + Atau2(:,1,5) * A2(:,5,4) 617c 618 A0(:,1,5)= 619 & Atau2(:,1,1) * A2(:,1,5) 620 & + Atau2(:,1,2) * A2(:,2,5) 621 & + Atau2(:,1,3) * A2(:,3,5) 622 & + Atau2(:,1,4) * A2(:,4,5) 623 & + Atau2(:,1,5) * A2(:,5,5) 624c 625c row two 626c 627 A0(:,2,1)= 628 & Atau2(:,2,1) * A2(:,1,1) 629 & + Atau2(:,2,2) * A2(:,2,1) 630 & + Atau2(:,2,3) * A2(:,3,1) 631 & + Atau2(:,2,4) * A2(:,4,1) 632 & + Atau2(:,2,5) * A2(:,5,1) 633c 634 A0(:,2,2)= 635c & Atau2(:,2,1) * A2(:,1,2) 636 & Atau2(:,2,2) * A2(:,2,2) 637c & + Atau2(:,2,3) * A2(:,3,2) 638c & + Atau2(:,2,4) * A2(:,4,2) 639 & + Atau2(:,2,5) * A2(:,5,2) 640c 641 A0(:,2,3)= 642 & Atau2(:,2,1) * A2(:,1,3) 643 & + Atau2(:,2,2) * A2(:,2,3) 644 & + Atau2(:,2,3) * A2(:,3,3) 645 & + Atau2(:,2,4) * A2(:,4,3) 646 & + Atau2(:,2,5) * A2(:,5,3) 647c 648 A0(:,2,4)= 649c & Atau2(:,2,1) * A2(:,1,4) 650c & + Atau2(:,2,2) * A2(:,2,4) 651c & + Atau2(:,2,3) * A2(:,3,4) 652 & Atau2(:,2,4) * A2(:,4,4) 653 & + Atau2(:,2,5) * A2(:,5,4) 654c 655 A0(:,2,5)= 656 & Atau2(:,2,1) * A2(:,1,5) 657 & + Atau2(:,2,2) * A2(:,2,5) 658 & + Atau2(:,2,3) * A2(:,3,5) 659 & + Atau2(:,2,4) * A2(:,4,5) 660 & + Atau2(:,2,5) * A2(:,5,5) 661c 662c row three 663c 664 A0(:,3,1)= 665 & Atau2(:,3,1) * A2(:,1,1) 666 & + Atau2(:,3,2) * A2(:,2,1) 667 & + Atau2(:,3,3) * A2(:,3,1) 668 & + Atau2(:,3,4) * A2(:,4,1) 669 & + Atau2(:,3,5) * A2(:,5,1) 670c 671 A0(:,3,2)= 672c & Atau2(:,3,1) * A2(:,1,2) 673 & Atau2(:,3,2) * A2(:,2,2) 674c & + Atau2(:,3,3) * A2(:,3,2) 675c & + Atau2(:,3,4) * A2(:,4,2) 676 & + Atau2(:,3,5) * A2(:,5,2) 677c 678 A0(:,3,3)= 679 & Atau2(:,2,1) * A2(:,1,3) 680 & + Atau2(:,2,2) * A2(:,2,3) 681 & + Atau2(:,3,3) * A2(:,3,3) 682 & + Atau2(:,2,4) * A2(:,4,3) 683 & + Atau2(:,3,5) * A2(:,5,3) 684c 685 A0(:,3,4)= 686c & Atau1(:,2,1) * A1(:,1,4) 687c & + Atau1(:,2,2) * A1(:,2,4) 688c & + Atau1(:,2,3) * A1(:,3,4) 689 & Atau2(:,3,4) * A2(:,4,4) 690 & + Atau2(:,3,5) * A2(:,5,4) 691c 692 A0(:,3,5)= 693 & Atau2(:,3,1) * A2(:,1,5) 694 & + Atau2(:,3,2) * A2(:,2,5) 695 & + Atau2(:,3,3) * A2(:,3,5) 696 & + Atau2(:,3,4) * A2(:,4,5) 697 & + Atau2(:,3,5) * A2(:,5,5) 698 699c 700c row four 701c 702 A0(:,4,1)= 703 & Atau2(:,4,1) * A2(:,1,1) 704 & + Atau2(:,4,2) * A2(:,2,1) 705 & + Atau2(:,4,3) * A2(:,3,1) 706 & + Atau2(:,4,4) * A2(:,4,1) 707 & + Atau2(:,4,5) * A2(:,5,1) 708c 709 A0(:,4,2)= 710c & Atau2(:,4,1) * A2(:,1,2) 711 & Atau2(:,4,2) * A2(:,2,2) 712c & + Atau2(:,4,3) * A2(:,3,2) 713c & + Atau2(:,4,4) * A2(:,4,2) 714 & + Atau2(:,4,5) * A2(:,5,2) 715c 716 A0(:,4,3)= 717 & Atau2(:,2,1) * A2(:,1,3) 718 & + Atau2(:,2,2) * A2(:,2,3) 719 & + Atau2(:,4,3) * A2(:,3,3) 720 & + Atau2(:,2,4) * A2(:,4,3) 721 & + Atau2(:,4,5) * A2(:,5,3) 722c 723 A0(:,4,4)= 724c & Atau1(:,2,1) * A1(:,1,4) 725c & + Atau1(:,2,2) * A1(:,2,4) 726c & + Atau1(:,2,3) * A1(:,3,4) 727 & Atau2(:,4,4) * A2(:,4,4) 728 & + Atau2(:,4,5) * A2(:,5,4) 729c 730 A0(:,4,5)= 731 & Atau2(:,4,1) * A2(:,1,5) 732 & + Atau2(:,4,2) * A2(:,2,5) 733 & + Atau2(:,4,3) * A2(:,3,5) 734 & + Atau2(:,4,4) * A2(:,4,5) 735 & + Atau2(:,4,5) * A2(:,5,5) 736c 737c row five 738c 739 740 A0(:,5,1)= 741 & Atau2(:,5,1) * A2(:,1,1) 742 & + Atau2(:,5,2) * A2(:,2,1) 743 & + Atau2(:,5,3) * A2(:,3,1) 744 & + Atau2(:,5,4) * A2(:,4,1) 745 & + Atau2(:,5,5) * A2(:,5,1) 746c 747 A0(:,5,2)= 748c & Atau2(:,5,1) * A2(:,1,2) 749 & Atau2(:,5,2) * A2(:,2,2) 750c & + Atau2(:,5,3) * A2(:,3,2) 751c & + Atau2(:,5,4) * A2(:,4,2) 752 & + Atau2(:,5,5) * A2(:,5,2) 753c 754 A0(:,5,3)= 755 & Atau2(:,2,1) * A2(:,1,3) 756 & + Atau2(:,2,2) * A2(:,2,3) 757 & + Atau2(:,5,3) * A2(:,3,3) 758 & + Atau2(:,2,4) * A2(:,4,3) 759 & + Atau2(:,3,5) * A2(:,5,3) 760c 761 A0(:,5,4)= 762c & Atau1(:,2,1) * A1(:,1,4) 763c & + Atau1(:,2,2) * A1(:,2,4) 764c & + Atau1(:,2,3) * A1(:,3,4) 765 & Atau2(:,5,4) * A2(:,4,4) 766 & + Atau2(:,5,5) * A2(:,5,4) 767c 768 A0(:,5,5)= 769 & Atau2(:,5,1) * A2(:,1,5) 770 & + Atau2(:,5,2) * A2(:,2,5) 771 & + Atau2(:,5,3) * A2(:,3,5) 772 & + Atau2(:,5,4) * A2(:,4,5) 773 & + Atau2(:,5,5) * A2(:,5,5) 774 775c 776 777 else 778 A0 = zero 779 endif 780c 781 if (Navier .eq. 1) then 782 A0(:,2,2) = A0(:,2,2) + rmu !rmu 783 A0(:,3,3) = A0(:,3,3) + rlm2mu !rlm2mu 784 A0(:,4,4) = A0(:,4,4) + rmu !rmu 785 A0(:,5,2) = A0(:,5,2) + rmu * u1 786 A0(:,5,3) = A0(:,5,3) + rlm2mu * u2 787 A0(:,5,4) = A0(:,5,4) + rmu * u3 788 A0(:,5,5) = A0(:,5,5) + con !con 789 endif 790c 791 do j = 1, nshl 792 tmp = WdetJ * shg(:,j,2) * shg(:,j,2) 793 do i=1,nflow 794 do k=1,nflow 795 BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k) 796 enddo 797 enddo 798 enddo 799c 800c start with 3 3 contribution 801c 802 if(ivart.ge.2) then 803 804c 805c row one 806c 807 A0(:,1,1) = 808 & Atau3(:,1,1) * A3(:,1,1) 809 & + Atau3(:,1,2) * A3(:,2,1) 810 & + Atau3(:,1,3) * A3(:,3,1) 811 & + Atau3(:,1,4) * A3(:,4,1) 812 & + Atau3(:,1,5) * A3(:,5,1) 813c 814 A0(:,1,2) = 815c & Atau2(:,1,1) * A2(:,1,2) 816 & Atau3(:,1,2) * A3(:,2,2) 817c & + Atau2(:,1,3) * A2(:,3,2) 818c & + Atau2(:,1,4) * A2(:,4,2) 819 & + Atau3(:,1,5) * A3(:,5,2) 820c 821 A0(:,1,3) = 822c & Atau2(:,1,1) * A2(:,1,3) 823c & + Atau2(:,1,2) * A2(:,2,3) 824 & Atau3(:,1,3) * A3(:,3,3) 825c & + Atau2(:,1,4) * A2(:,4,3) 826 & + Atau3(:,1,5) * A3(:,5,3) 827c 828 A0(:,1,4)= 829 & Atau3(:,1,1) * A3(:,1,4) 830 & + Atau3(:,1,2) * A3(:,2,4) 831 & + Atau3(:,1,3) * A3(:,3,4) 832 & + Atau3(:,1,4) * A3(:,4,4) 833 & + Atau3(:,1,5) * A3(:,5,4) 834c 835 A0(:,1,5)= 836 & Atau3(:,1,1) * A3(:,1,5) 837 & + Atau3(:,1,2) * A3(:,2,5) 838 & + Atau3(:,1,3) * A3(:,3,5) 839 & + Atau3(:,1,4) * A3(:,4,5) 840 & + Atau3(:,1,5) * A3(:,5,5) 841c 842c row two 843c 844 A0(:,2,1)= 845 & Atau3(:,2,1) * A3(:,1,1) 846 & + Atau3(:,2,2) * A3(:,2,1) 847 & + Atau3(:,2,3) * A3(:,3,1) 848 & + Atau3(:,2,4) * A3(:,4,1) 849 & + Atau3(:,2,5) * A3(:,5,1) 850c 851 A0(:,2,2)= 852c & Atau2(:,2,1) * A2(:,1,2) 853 & Atau3(:,2,2) * A3(:,2,2) 854c & + Atau2(:,2,3) * A2(:,3,2) 855c & + Atau2(:,2,4) * A2(:,4,2) 856 & + Atau3(:,2,5) * A3(:,5,2) 857c 858 A0(:,2,3)= 859c & Atau2(:,2,1) * A2(:,1,3) 860c & + Atau2(:,2,2) * A2(:,2,3) 861 & Atau3(:,2,3) * A3(:,3,3) 862c & + Atau2(:,2,4) * A2(:,4,3) 863 & + Atau3(:,2,5) * A3(:,5,3) 864c 865 A0(:,2,4)= 866 & Atau3(:,2,1) * A3(:,1,4) 867 & + Atau3(:,2,2) * A3(:,2,4) 868 & + Atau3(:,2,3) * A3(:,3,4) 869 & + Atau3(:,2,4) * A3(:,4,4) 870 & + Atau3(:,2,5) * A3(:,5,4) 871c 872 A0(:,2,5)= 873 & Atau3(:,2,1) * A3(:,1,5) 874 & + Atau3(:,2,2) * A3(:,2,5) 875 & + Atau3(:,2,3) * A3(:,3,5) 876 & + Atau3(:,2,4) * A3(:,4,5) 877 & + Atau3(:,2,5) * A3(:,5,5) 878c 879c row three 880c 881 A0(:,3,1)= 882 & Atau3(:,3,1) * A3(:,1,1) 883 & + Atau3(:,3,2) * A3(:,2,1) 884 & + Atau3(:,3,3) * A3(:,3,1) 885 & + Atau3(:,3,4) * A3(:,4,1) 886 & + Atau3(:,3,5) * A3(:,5,1) 887c 888 A0(:,3,2)= 889c & Atau2(:,3,1) * A2(:,1,2) 890 & Atau3(:,3,2) * A3(:,2,2) 891c & + Atau2(:,3,3) * A2(:,3,2) 892c & + Atau2(:,3,4) * A2(:,4,2) 893 & + Atau3(:,3,5) * A3(:,5,2) 894c 895 A0(:,3,3)= 896c & Atau2(:,2,1) * A2(:,1,3) 897c & + Atau2(:,2,2) * A2(:,2,3) 898 & Atau3(:,3,3) * A3(:,3,3) 899c & + Atau2(:,2,4) * A2(:,4,3) 900 & + Atau3(:,3,5) * A3(:,5,3) 901c 902 A0(:,3,4)= 903 & Atau3(:,2,1) * A3(:,1,4) 904 & + Atau3(:,2,2) * A3(:,2,4) 905 & + Atau3(:,2,3) * A3(:,3,4) 906 & + Atau3(:,3,4) * A3(:,4,4) 907 & + Atau3(:,3,5) * A3(:,5,4) 908c 909 A0(:,3,5)= 910 & Atau3(:,3,1) * A3(:,1,5) 911 & + Atau3(:,3,2) * A3(:,2,5) 912 & + Atau3(:,3,3) * A3(:,3,5) 913 & + Atau3(:,3,4) * A3(:,4,5) 914 & + Atau3(:,3,5) * A3(:,5,5) 915 916c 917c row four 918c 919 A0(:,4,1)= 920 & Atau3(:,4,1) * A3(:,1,1) 921 & + Atau3(:,4,2) * A3(:,2,1) 922 & + Atau3(:,4,3) * A3(:,3,1) 923 & + Atau3(:,4,4) * A3(:,4,1) 924 & + Atau3(:,4,5) * A3(:,5,1) 925c 926 A0(:,4,2)= 927c & Atau2(:,4,1) * A2(:,1,2) 928 & Atau3(:,4,2) * A3(:,2,2) 929c & + Atau2(:,4,3) * A2(:,3,2) 930c & + Atau2(:,4,4) * A2(:,4,2) 931 & + Atau3(:,4,5) * A3(:,5,2) 932c 933 A0(:,4,3)= 934c & Atau2(:,2,1) * A2(:,1,3) 935c & + Atau2(:,2,2) * A2(:,2,3) 936 & Atau3(:,4,3) * A3(:,3,3) 937c & + Atau2(:,2,4) * A2(:,4,3) 938 & + Atau3(:,4,5) * A3(:,5,3) 939c 940 A0(:,4,4)= 941 & Atau3(:,2,1) * A3(:,1,4) 942 & + Atau3(:,2,2) * A3(:,2,4) 943 & + Atau3(:,2,3) * A3(:,3,4) 944 & + Atau3(:,4,4) * A3(:,4,4) 945 & + Atau3(:,4,5) * A3(:,5,4) 946c 947 A0(:,4,5)= 948 & Atau3(:,4,1) * A3(:,1,5) 949 & + Atau3(:,4,2) * A3(:,2,5) 950 & + Atau3(:,4,3) * A3(:,3,5) 951 & + Atau3(:,4,4) * A3(:,4,5) 952 & + Atau3(:,4,5) * A3(:,5,5) 953c 954c row five 955c 956 957 A0(:,5,1)= 958 & Atau3(:,5,1) * A3(:,1,1) 959 & + Atau3(:,5,2) * A3(:,2,1) 960 & + Atau3(:,5,3) * A3(:,3,1) 961 & + Atau3(:,5,4) * A3(:,4,1) 962 & + Atau3(:,5,5) * A3(:,5,1) 963c 964 A0(:,5,2)= 965c & Atau2(:,5,1) * A2(:,1,2) 966 & Atau3(:,5,2) * A3(:,2,2) 967c & + Atau2(:,5,3) * A2(:,3,2) 968c & + Atau2(:,5,4) * A2(:,4,2) 969 & + Atau3(:,5,5) * A3(:,5,2) 970c 971 A0(:,5,3)= 972c & Atau2(:,2,1) * A2(:,1,3) 973c & + Atau2(:,2,2) * A2(:,2,3) 974 & Atau3(:,5,3) * A3(:,3,3) 975c & + Atau2(:,2,4) * A2(:,4,3) 976 & + Atau3(:,3,5) * A3(:,5,3) 977c 978 A0(:,5,4)= 979 & Atau3(:,2,1) * A3(:,1,4) 980 & + Atau3(:,2,2) * A3(:,2,4) 981 & + Atau3(:,2,3) * A3(:,3,4) 982 & + Atau3(:,5,4) * A3(:,4,4) 983 & + Atau3(:,5,5) * A3(:,5,4) 984c 985 A0(:,5,5)= 986 & Atau3(:,5,1) * A3(:,1,5) 987 & + Atau3(:,5,2) * A3(:,2,5) 988 & + Atau3(:,5,3) * A3(:,3,5) 989 & + Atau3(:,5,4) * A3(:,4,5) 990 & + Atau3(:,5,5) * A3(:,5,5) 991 992c 993 994 else 995 A0 = zero 996 endif 997c 998 if (Navier .eq. 1) then 999 A0(:,2,2) = A0(:,2,2) + rmu !rmu 1000 A0(:,3,3) = A0(:,3,3) + rmu !rmu 1001 A0(:,4,4) = A0(:,4,4) + rlm2mu !rlm2mu 1002 A0(:,5,2) = A0(:,5,2) + rmu * u1 1003 A0(:,5,3) = A0(:,5,3) + rmu * u2 1004 A0(:,5,4) = A0(:,5,4) + rlm2mu * u3 1005 A0(:,5,5) = A0(:,5,5) + con !con 1006 endif 1007c 1008 do j = 1, nshl 1009 tmp = WdetJ * shg(:,j,3) * shg(:,j,3) 1010 do i=1,nflow 1011 do k=1,nflow 1012 BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k) 1013 enddo 1014 enddo 1015 enddo 1016c 1017c now for the 1 2 plus 2 1 1018 if(ivart.ge.2) then 1019 1020 A0(:,1,1) = 1021 & Atau1(:,1,1) * A2(:,1,1) 1022 & + Atau1(:,1,2) * A2(:,2,1) 1023 & + Atau1(:,1,3) * A2(:,3,1) 1024 & + Atau1(:,1,4) * A2(:,4,1) 1025 & + Atau1(:,1,5) * A2(:,5,1) 1026 & + Atau2(:,1,1) * A1(:,1,1) 1027 & + Atau2(:,1,2) * A1(:,2,1) 1028 & + Atau2(:,1,3) * A1(:,3,1) 1029 & + Atau2(:,1,4) * A1(:,4,1) 1030 & + Atau2(:,1,5) * A1(:,5,1) 1031c 1032 A0(:,1,2) = 1033c & Atau1(:,1,1) * A2(:,1,2) 1034 & Atau1(:,1,2) * A2(:,2,2) 1035c & + Atau1(:,1,3) * A2(:,3,2) 1036c & + Atau1(:,1,4) * A2(:,4,2) 1037 & + Atau1(:,1,5) * A2(:,5,2) 1038 & + Atau2(:,1,1) * A1(:,1,2) 1039 & + Atau2(:,1,2) * A1(:,2,2) 1040 & + Atau2(:,1,3) * A1(:,3,2) 1041 & + Atau2(:,1,4) * A1(:,4,2) 1042 & + Atau2(:,1,5) * A1(:,5,2) 1043 1044c 1045 A0(:,1,3) = 1046 & Atau1(:,1,1) * A2(:,1,3) 1047 & + Atau1(:,1,2) * A2(:,2,3) 1048 & + Atau1(:,1,3) * A2(:,3,3) 1049 & + Atau1(:,1,4) * A2(:,4,3) 1050 & + Atau1(:,1,5) * A2(:,5,3) 1051c & + Atau2(:,1,1) * A1(:,1,3) 1052c & + Atau2(:,1,2) * A1(:,2,3) 1053 & + Atau2(:,1,3) * A1(:,3,3) 1054c & + Atau2(:,1,4) * A1(:,4,3) 1055 & + Atau2(:,1,5) * A1(:,5,3) 1056c 1057 A0(:,1,4)= 1058c & Atau1(:,1,1) * A2(:,1,4) 1059c & + Atau1(:,1,2) * A2(:,2,4) 1060c & + Atau1(:,1,3) * A2(:,3,4) 1061 & Atau1(:,1,4) * A2(:,4,4) 1062 & + Atau1(:,1,5) * A2(:,5,4) 1063c & + Atau2(:,1,1) * A1(:,1,4) 1064c & + Atau2(:,1,2) * A1(:,2,4) 1065c & + Atau2(:,1,3) * A1(:,3,4) 1066 & + Atau2(:,1,4) * A1(:,4,4) 1067 & + Atau2(:,1,5) * A1(:,5,4) 1068c 1069 A0(:,1,5)= 1070 & Atau1(:,1,1) * A2(:,1,5) 1071 & + Atau1(:,1,2) * A2(:,2,5) 1072 & + Atau1(:,1,3) * A2(:,3,5) 1073 & + Atau1(:,1,4) * A2(:,4,5) 1074 & + Atau1(:,1,5) * A2(:,5,5) 1075 & + Atau2(:,1,1) * A1(:,1,5) 1076 & + Atau2(:,1,2) * A1(:,2,5) 1077 & + Atau2(:,1,3) * A1(:,3,5) 1078 & + Atau2(:,1,4) * A1(:,4,5) 1079 & + Atau2(:,1,5) * A1(:,5,5) 1080c 1081c row two 1082c 1083 A0(:,2,1)= 1084 & Atau1(:,2,1) * A2(:,1,1) 1085 & + Atau1(:,2,2) * A2(:,2,1) 1086 & + Atau1(:,2,3) * A2(:,3,1) 1087 & + Atau1(:,2,4) * A2(:,4,1) 1088 & + Atau1(:,2,5) * A2(:,5,1) 1089 & + Atau2(:,2,1) * A1(:,1,1) 1090 & + Atau2(:,2,2) * A1(:,2,1) 1091 & + Atau2(:,2,3) * A1(:,3,1) 1092 & + Atau2(:,2,4) * A1(:,4,1) 1093 & + Atau2(:,2,5) * A1(:,5,1) 1094c 1095 A0(:,2,2)= 1096c & Atau1(:,2,1) * A2(:,1,2) 1097 & Atau1(:,2,2) * A2(:,2,2) 1098c & + Atau1(:,2,3) * A2(:,3,2) 1099c & + Atau1(:,2,4) * A2(:,4,2) 1100 & + Atau1(:,2,5) * A2(:,5,2) 1101 & + Atau2(:,2,1) * A1(:,1,2) 1102 & + Atau2(:,2,2) * A1(:,2,2) 1103 & + Atau2(:,2,3) * A1(:,3,2) 1104 & + Atau2(:,2,4) * A1(:,4,2) 1105 & + Atau2(:,2,5) * A1(:,5,2) 1106c 1107 A0(:,2,3)= 1108 & Atau1(:,2,1) * A2(:,1,3) 1109 & + Atau1(:,2,2) * A2(:,2,3) 1110 & + Atau1(:,2,3) * A2(:,3,3) 1111 & + Atau1(:,2,4) * A2(:,4,3) 1112 & + Atau1(:,2,5) * A2(:,5,3) 1113c & + Atau2(:,2,1) * A1(:,1,3) 1114c & + Atau2(:,2,2) * A1(:,2,3) 1115 & + Atau2(:,2,3) * A1(:,3,3) 1116c & + Atau2(:,2,4) * A1(:,4,3) 1117 & + Atau2(:,2,5) * A1(:,5,3) 1118c 1119 A0(:,2,4)= 1120c & Atau1(:,2,1) * A2(:,1,4) 1121c & + Atau1(:,2,2) * A2(:,2,4) 1122c & + Atau1(:,2,3) * A2(:,3,4) 1123 & Atau1(:,2,4) * A2(:,4,4) 1124 & + Atau1(:,2,5) * A2(:,5,4) 1125c & + Atau2(:,2,1) * A1(:,1,4) 1126c & + Atau2(:,2,2) * A1(:,2,4) 1127c & + Atau2(:,2,3) * A1(:,3,4) 1128 & + Atau2(:,2,4) * A1(:,4,4) 1129 & + Atau2(:,2,5) * A1(:,5,4) 1130 1131c 1132 A0(:,2,5)= 1133 & Atau1(:,2,1) * A2(:,1,5) 1134 & + Atau1(:,2,2) * A2(:,2,5) 1135 & + Atau1(:,2,3) * A2(:,3,5) 1136 & + Atau1(:,2,4) * A2(:,4,5) 1137 & + Atau1(:,2,5) * A2(:,5,5) 1138 & + Atau2(:,2,1) * A1(:,1,5) 1139 & + Atau2(:,2,2) * A1(:,2,5) 1140 & + Atau2(:,2,3) * A1(:,3,5) 1141 & + Atau2(:,2,4) * A1(:,4,5) 1142 & + Atau2(:,2,5) * A1(:,5,5) 1143 1144c 1145c row three 1146c 1147 A0(:,3,1)= 1148 & Atau1(:,3,1) * A2(:,1,1) 1149 & + Atau1(:,3,2) * A2(:,2,1) 1150 & + Atau1(:,3,3) * A2(:,3,1) 1151 & + Atau1(:,3,4) * A2(:,4,1) 1152 & + Atau1(:,3,5) * A2(:,5,1) 1153 & + Atau2(:,3,1) * A1(:,1,1) 1154 & + Atau2(:,3,2) * A1(:,2,1) 1155 & + Atau2(:,3,3) * A1(:,3,1) 1156 & + Atau2(:,3,4) * A1(:,4,1) 1157 & + Atau2(:,3,5) * A1(:,5,1) 1158 1159c 1160 A0(:,3,2)= 1161c & Atau1(:,3,1) * A2(:,1,2) 1162 & Atau1(:,3,2) * A2(:,2,2) 1163c & + Atau1(:,3,3) * A2(:,3,2) 1164c & + Atau1(:,3,4) * A2(:,4,2) 1165 & + Atau1(:,3,5) * A2(:,5,2) 1166 & + Atau2(:,3,1) * A1(:,1,2) 1167 & + Atau2(:,3,2) * A1(:,2,2) 1168 & + Atau2(:,3,3) * A1(:,3,2) 1169 & + Atau2(:,3,4) * A1(:,4,2) 1170 & + Atau2(:,3,5) * A1(:,5,2) 1171c 1172 A0(:,3,3)= 1173 & Atau1(:,3,1) * A2(:,1,3) 1174 & + Atau1(:,3,2) * A2(:,2,3) 1175 & + Atau1(:,3,3) * A2(:,3,3) 1176 & + Atau1(:,3,4) * A2(:,4,3) 1177 & + Atau1(:,3,5) * A2(:,5,3) 1178c & + Atau2(:,3,1) * A1(:,1,3) 1179c & + Atau2(:,3,2) * A1(:,2,3) 1180 & + Atau2(:,3,3) * A1(:,3,3) 1181c & + Atau2(:,3,4) * A1(:,4,3) 1182 & + Atau2(:,3,5) * A1(:,5,3) 1183c 1184 A0(:,3,4)= 1185c & Atau1(:,3,1) * A2(:,1,4) 1186c & + Atau1(:,3,2) * A2(:,2,4) 1187c & + Atau1(:,3,3) * A2(:,3,4) 1188 & Atau1(:,3,4) * A2(:,4,4) 1189 & + Atau1(:,3,5) * A2(:,5,4) 1190c & + Atau2(:,3,1) * A1(:,1,4) 1191c & + Atau2(:,3,2) * A1(:,2,4) 1192c & + Atau2(:,3,3) * A1(:,3,4) 1193 & + Atau2(:,3,4) * A1(:,4,4) 1194 & + Atau2(:,3,5) * A1(:,5,4) 1195c 1196 A0(:,3,5)= 1197 & Atau1(:,3,1) * A2(:,1,5) 1198 & + Atau1(:,3,2) * A2(:,2,5) 1199 & + Atau1(:,3,3) * A2(:,3,5) 1200 & + Atau1(:,3,4) * A2(:,4,5) 1201 & + Atau1(:,3,5) * A2(:,5,5) 1202 & + Atau2(:,3,1) * A1(:,1,5) 1203 & + Atau2(:,3,2) * A1(:,2,5) 1204 & + Atau2(:,3,3) * A1(:,3,5) 1205 & + Atau2(:,3,4) * A1(:,4,5) 1206 & + Atau2(:,3,5) * A1(:,5,5) 1207 1208 1209c 1210c row four 1211c 1212 A0(:,4,1)= 1213 & Atau1(:,4,1) * A2(:,1,1) 1214 & + Atau1(:,4,2) * A2(:,2,1) 1215 & + Atau1(:,4,3) * A2(:,3,1) 1216 & + Atau1(:,4,4) * A2(:,4,1) 1217 & + Atau1(:,4,5) * A2(:,5,1) 1218 & + Atau2(:,4,1) * A1(:,1,1) 1219 & + Atau2(:,4,2) * A1(:,2,1) 1220 & + Atau2(:,4,3) * A1(:,3,1) 1221 & + Atau2(:,4,4) * A1(:,4,1) 1222 & + Atau2(:,4,5) * A1(:,5,1) 1223c 1224 A0(:,4,2)= 1225c & Atau1(:,4,1) * A2(:,1,2) 1226 & Atau1(:,4,2) * A2(:,2,2) 1227c & + Atau1(:,4,3) * A2(:,3,2) 1228c & + Atau1(:,4,4) * A2(:,4,2) 1229 & + Atau1(:,4,5) * A2(:,5,2) 1230 & + Atau2(:,4,1) * A1(:,1,2) 1231 & + Atau2(:,4,2) * A1(:,2,2) 1232 & + Atau2(:,4,3) * A1(:,3,2) 1233 & + Atau2(:,4,4) * A1(:,4,2) 1234 & + Atau2(:,4,5) * A1(:,5,2) 1235c 1236 A0(:,4,3)= 1237 & Atau1(:,4,1) * A2(:,1,3) 1238 & + Atau1(:,4,2) * A2(:,2,3) 1239 & + Atau1(:,4,3) * A2(:,3,3) 1240 & + Atau1(:,4,4) * A2(:,4,3) 1241 & + Atau1(:,4,5) * A2(:,5,3) 1242c & + Atau2(:,4,1) * A1(:,1,3) 1243c & + Atau2(:,4,2) * A1(:,2,3) 1244 & + Atau2(:,4,3) * A1(:,3,3) 1245c & + Atau2(:,4,4) * A1(:,4,3) 1246 & + Atau2(:,4,5) * A1(:,5,3) 1247c 1248 A0(:,4,4)= 1249c & Atau1(:,4,1) * A2(:,1,4) 1250c & + Atau1(:,4,2) * A2(:,2,4) 1251c & + Atau1(:,4,3) * A2(:,3,4) 1252 & Atau1(:,4,4) * A2(:,4,4) 1253 & + Atau1(:,4,5) * A2(:,5,4) 1254c & + Atau2(:,4,1) * A1(:,1,4) 1255c & + Atau2(:,4,2) * A1(:,2,4) 1256c & + Atau2(:,4,3) * A1(:,3,4) 1257 & + Atau2(:,4,4) * A1(:,4,4) 1258 & + Atau2(:,4,5) * A1(:,5,4) 1259c 1260 A0(:,4,5)= 1261 & Atau1(:,4,1) * A2(:,1,5) 1262 & + Atau1(:,4,2) * A2(:,2,5) 1263 & + Atau1(:,4,3) * A2(:,3,5) 1264 & + Atau1(:,4,4) * A2(:,4,5) 1265 & + Atau1(:,4,5) * A2(:,5,5) 1266 & + Atau2(:,4,1) * A1(:,1,5) 1267 & + Atau2(:,4,2) * A1(:,2,5) 1268 & + Atau2(:,4,3) * A1(:,3,5) 1269 & + Atau2(:,4,4) * A1(:,4,5) 1270 & + Atau2(:,4,5) * A1(:,5,5) 1271c 1272c row five 1273c 1274 1275 A0(:,5,1)= 1276 & Atau1(:,5,1) * A2(:,1,1) 1277 & + Atau1(:,5,2) * A2(:,2,1) 1278 & + Atau1(:,5,3) * A2(:,3,1) 1279 & + Atau1(:,5,4) * A2(:,4,1) 1280 & + Atau1(:,5,5) * A2(:,5,1) 1281 & + Atau2(:,5,1) * A1(:,1,1) 1282 & + Atau2(:,5,2) * A1(:,2,1) 1283 & + Atau2(:,5,3) * A1(:,3,1) 1284 & + Atau2(:,5,4) * A1(:,4,1) 1285 & + Atau2(:,5,5) * A1(:,5,1) 1286c 1287 A0(:,5,2)= 1288c & Atau1(:,5,1) * A2(:,1,2) 1289 & Atau1(:,5,2) * A2(:,2,2) 1290c & + Atau1(:,5,3) * A2(:,3,2) 1291c & + Atau1(:,5,4) * A2(:,4,2) 1292 & + Atau1(:,5,5) * A2(:,5,2) 1293 & + Atau2(:,5,1) * A1(:,1,2) 1294 & + Atau2(:,5,2) * A1(:,2,2) 1295 & + Atau2(:,5,3) * A1(:,3,2) 1296 & + Atau2(:,5,4) * A1(:,4,2) 1297 & + Atau2(:,5,5) * A1(:,5,2) 1298c 1299 A0(:,5,3)= 1300 & Atau1(:,5,1) * A2(:,1,3) 1301 & + Atau1(:,5,2) * A2(:,2,3) 1302 & + Atau1(:,5,3) * A2(:,3,3) 1303 & + Atau1(:,5,4) * A2(:,4,3) 1304 & + Atau1(:,5,5) * A2(:,5,3) 1305c & + Atau2(:,5,1) * A1(:,1,3) 1306c & + Atau2(:,5,2) * A1(:,2,3) 1307 & + Atau2(:,5,3) * A1(:,3,3) 1308c & + Atau2(:,5,4) * A1(:,4,3) 1309 & + Atau2(:,5,5) * A1(:,5,3) 1310c 1311 A0(:,5,4)= 1312c & Atau1(:,5,1) * A2(:,1,4) 1313c & + Atau1(:,5,2) * A2(:,2,4) 1314c & + Atau1(:,5,3) * A2(:,3,4) 1315 & Atau1(:,5,4) * A2(:,4,4) 1316 & + Atau1(:,5,5) * A2(:,5,4) 1317c & + Atau2(:,5,1) * A1(:,1,4) 1318c & + Atau2(:,5,2) * A1(:,2,4) 1319c & + Atau2(:,5,3) * A1(:,3,4) 1320 & + Atau2(:,5,4) * A1(:,4,4) 1321 & + Atau2(:,5,5) * A1(:,5,4) 1322 1323c 1324 A0(:,5,5)= 1325 & Atau1(:,5,1) * A2(:,1,5) 1326 & + Atau1(:,5,2) * A2(:,2,5) 1327 & + Atau1(:,5,3) * A2(:,3,5) 1328 & + Atau1(:,5,4) * A2(:,4,5) 1329 & + Atau1(:,5,5) * A2(:,5,5) 1330 & + Atau2(:,5,1) * A1(:,1,5) 1331 & + Atau2(:,5,2) * A1(:,2,5) 1332 & + Atau2(:,5,3) * A1(:,3,5) 1333 & + Atau2(:,5,4) * A1(:,4,5) 1334 & + Atau2(:,5,5) * A1(:,5,5) 1335 1336c 1337 else 1338 A0 = zero 1339 endif 1340c 1341 if (Navier .eq. 1) then 1342 A0(:,2,3) = A0(:,2,3) + rlm2mu - rmu 1343 A0(:,3,2) = A0(:,3,2) + rlm2mu - rmu 1344 A0(:,5,2) = A0(:,5,2) + (rlm2mu - rmu) * u2 1345 A0(:,5,3) = A0(:,5,3) + (rlm2mu- rmu) * u1 1346 endif 1347c 1348 do j = 1, nshl 1349 tmp = WdetJ * shg(:,j,1) * shg(:,j,2) 1350 do i=1,nflow 1351 do k=1,nflow 1352 BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k) 1353 enddo 1354 enddo 1355 enddo 1356c 1357c now for the 1 3 plus 3 1 1358c 1359 if(ivart.ge.2) then 1360c 1361c row one 1362c 1363 A0(:,1,1)= 1364 & Atau1(:,1,1) * A3(:,1,1) 1365 & + Atau1(:,1,2) * A3(:,2,1) 1366 & + Atau1(:,1,3) * A3(:,3,1) 1367 & + Atau1(:,1,4) * A3(:,4,1) 1368 & + Atau1(:,1,5) * A3(:,5,1) 1369 & + Atau3(:,1,1) * A1(:,1,1) 1370 & + Atau3(:,1,2) * A1(:,2,1) 1371 & + Atau3(:,1,3) * A1(:,3,1) 1372 & + Atau3(:,1,4) * A1(:,4,1) 1373 & + Atau3(:,1,5) * A1(:,5,1) 1374c 1375 A0(:,1,2) = 1376c & Atau1(:,1,1) * A3(:,1,2) 1377 & Atau1(:,1,2) * A3(:,2,2) 1378c & + Atau1(:,1,3) * A3(:,3,2) 1379c & + Atau1(:,1,4) * A3(:,4,2) 1380 & + Atau1(:,1,5) * A3(:,5,2) 1381 & + Atau3(:,1,1) * A1(:,1,2) 1382 & + Atau3(:,1,2) * A1(:,2,2) 1383 & + Atau3(:,1,3) * A1(:,3,2) 1384 & + Atau3(:,1,4) * A1(:,4,2) 1385 & + Atau3(:,1,5) * A1(:,5,2) 1386 1387c 1388 A0(:,1,3) = 1389c & Atau1(:,1,1) * A3(:,1,3) 1390c & + Atau1(:,1,2) * A3(:,2,3) 1391 & Atau1(:,1,3) * A3(:,3,3) 1392c & + Atau1(:,1,4) * A3(:,4,3) 1393 & + Atau1(:,1,5) * A3(:,5,3) 1394c & + Atau3(:,1,1) * A1(:,1,3) 1395c & + Atau3(:,1,2) * A1(:,2,3) 1396 & + Atau3(:,1,3) * A1(:,3,3) 1397c & + Atau3(:,1,4) * A1(:,4,3) 1398 & + Atau3(:,1,5) * A1(:,5,3) 1399c 1400 A0(:,1,4)= 1401 & Atau1(:,1,1) * A3(:,1,4) 1402 & + Atau1(:,1,2) * A3(:,2,4) 1403 & + Atau1(:,1,3) * A3(:,3,4) 1404 & + Atau1(:,1,4) * A3(:,4,4) 1405 & + Atau1(:,1,5) * A3(:,5,4) 1406c & + Atau3(:,1,1) * A1(:,1,4) 1407c & + Atau3(:,1,2) * A1(:,2,4) 1408c & + Atau3(:,1,3) * A1(:,3,4) 1409 & + Atau3(:,1,4) * A1(:,4,4) 1410 & + Atau3(:,1,5) * A1(:,5,4) 1411c 1412 A0(:,1,5)= 1413 & Atau1(:,1,1) * A3(:,1,5) 1414 & + Atau1(:,1,2) * A3(:,2,5) 1415 & + Atau1(:,1,3) * A3(:,3,5) 1416 & + Atau1(:,1,4) * A3(:,4,5) 1417 & + Atau1(:,1,5) * A3(:,5,5) 1418 & + Atau3(:,1,1) * A1(:,1,5) 1419 & + Atau3(:,1,2) * A1(:,2,5) 1420 & + Atau3(:,1,3) * A1(:,3,5) 1421 & + Atau3(:,1,4) * A1(:,4,5) 1422 & + Atau3(:,1,5) * A1(:,5,5) 1423c 1424c row two 1425c 1426 A0(:,2,1)= 1427 & Atau1(:,2,1) * A3(:,1,1) 1428 & + Atau1(:,2,2) * A3(:,2,1) 1429 & + Atau1(:,2,3) * A3(:,3,1) 1430 & + Atau1(:,2,4) * A3(:,4,1) 1431 & + Atau1(:,2,5) * A3(:,5,1) 1432 & + Atau3(:,2,1) * A1(:,1,1) 1433 & + Atau3(:,2,2) * A1(:,2,1) 1434 & + Atau3(:,2,3) * A1(:,3,1) 1435 & + Atau3(:,2,4) * A1(:,4,1) 1436 & + Atau3(:,2,5) * A1(:,5,1) 1437c 1438 A0(:,2,2)= 1439c & Atau1(:,2,1) * A3(:,1,2) 1440 & Atau1(:,2,2) * A3(:,2,2) 1441c & + Atau1(:,2,3) * A3(:,3,2) 1442c & + Atau1(:,2,4) * A3(:,4,2) 1443 & + Atau1(:,2,5) * A3(:,5,2) 1444 & + Atau3(:,2,1) * A1(:,1,2) 1445 & + Atau3(:,2,2) * A1(:,2,2) 1446 & + Atau3(:,2,3) * A1(:,3,2) 1447 & + Atau3(:,2,4) * A1(:,4,2) 1448 & + Atau3(:,2,5) * A1(:,5,2) 1449c 1450 A0(:,2,3)= 1451c & Atau1(:,2,1) * A3(:,1,3) 1452c & + Atau1(:,2,2) * A3(:,2,3) 1453 & Atau1(:,2,3) * A3(:,3,3) 1454c & + Atau1(:,2,4) * A3(:,4,3) 1455 & + Atau1(:,2,5) * A3(:,5,3) 1456c & + Atau3(:,2,1) * A1(:,1,3) 1457c & + Atau3(:,2,2) * A1(:,2,3) 1458 & + Atau3(:,2,3) * A1(:,3,3) 1459c & + Atau3(:,2,4) * A1(:,4,3) 1460 & + Atau3(:,2,5) * A1(:,5,3) 1461c 1462 A0(:,2,4)= 1463 & Atau1(:,2,1) * A3(:,1,4) 1464 & + Atau1(:,2,2) * A3(:,2,4) 1465 & + Atau1(:,2,3) * A3(:,3,4) 1466 & + Atau1(:,2,4) * A3(:,4,4) 1467 & + Atau1(:,2,5) * A3(:,5,4) 1468c & + Atau3(:,2,1) * A1(:,1,4) 1469c & + Atau3(:,2,2) * A1(:,2,4) 1470c & + Atau3(:,2,3) * A1(:,3,4) 1471 & + Atau3(:,2,4) * A1(:,4,4) 1472 & + Atau3(:,2,5) * A1(:,5,4) 1473 1474c 1475 A0(:,2,5)= 1476 & Atau1(:,2,1) * A3(:,1,5) 1477 & + Atau1(:,2,2) * A3(:,2,5) 1478 & + Atau1(:,2,3) * A3(:,3,5) 1479 & + Atau1(:,2,4) * A3(:,4,5) 1480 & + Atau1(:,2,5) * A3(:,5,5) 1481 & + Atau3(:,2,1) * A1(:,1,5) 1482 & + Atau3(:,2,2) * A1(:,2,5) 1483 & + Atau3(:,2,3) * A1(:,3,5) 1484 & + Atau3(:,2,4) * A1(:,4,5) 1485 & + Atau3(:,2,5) * A1(:,5,5) 1486 1487c 1488c row three 1489c 1490 A0(:,3,1)= 1491 & Atau1(:,3,1) * A3(:,1,1) 1492 & + Atau1(:,3,2) * A3(:,2,1) 1493 & + Atau1(:,3,3) * A3(:,3,1) 1494 & + Atau1(:,3,4) * A3(:,4,1) 1495 & + Atau1(:,3,5) * A3(:,5,1) 1496 & + Atau3(:,3,1) * A1(:,1,1) 1497 & + Atau3(:,3,2) * A1(:,2,1) 1498 & + Atau3(:,3,3) * A1(:,3,1) 1499 & + Atau3(:,3,4) * A1(:,4,1) 1500 & + Atau3(:,3,5) * A1(:,5,1) 1501 1502c 1503 A0(:,3,2)= 1504c & Atau1(:,3,1) * A3(:,1,2) 1505 & Atau1(:,3,2) * A3(:,2,2) 1506c & + Atau1(:,3,3) * A3(:,3,2) 1507c & + Atau1(:,3,4) * A3(:,4,2) 1508 & + Atau1(:,3,5) * A3(:,5,2) 1509 & + Atau3(:,3,1) * A1(:,1,2) 1510 & + Atau3(:,3,2) * A1(:,2,2) 1511 & + Atau3(:,3,3) * A1(:,3,2) 1512 & + Atau3(:,3,4) * A1(:,4,2) 1513 & + Atau3(:,3,5) * A1(:,5,2) 1514c 1515 A0(:,3,3)= 1516c & Atau1(:,3,1) * A3(:,1,3) 1517c & + Atau1(:,3,2) * A3(:,2,3) 1518 & Atau1(:,3,3) * A3(:,3,3) 1519c & + Atau1(:,3,4) * A3(:,4,3) 1520 & + Atau1(:,3,5) * A3(:,5,3) 1521c & + Atau3(:,3,1) * A1(:,1,3) 1522c & + Atau3(:,3,2) * A1(:,2,3) 1523 & + Atau3(:,3,3) * A1(:,3,3) 1524c & + Atau3(:,3,4) * A1(:,4,3) 1525 & + Atau3(:,3,5) * A1(:,5,3) 1526c 1527 A0(:,3,4)= 1528 & Atau1(:,3,1) * A3(:,1,4) 1529 & + Atau1(:,3,2) * A3(:,2,4) 1530 & + Atau1(:,3,3) * A3(:,3,4) 1531 & + Atau1(:,3,4) * A3(:,4,4) 1532 & + Atau1(:,3,5) * A3(:,5,4) 1533c & + Atau3(:,3,1) * A1(:,1,4) 1534c & + Atau3(:,3,2) * A1(:,2,4) 1535c & + Atau3(:,3,3) * A1(:,3,4) 1536 & + Atau3(:,3,4) * A1(:,4,4) 1537 & + Atau3(:,3,5) * A1(:,5,4) 1538c 1539 A0(:,3,5)= 1540 & Atau1(:,3,1) * A3(:,1,5) 1541 & + Atau1(:,3,2) * A3(:,2,5) 1542 & + Atau1(:,3,3) * A3(:,3,5) 1543 & + Atau1(:,3,4) * A3(:,4,5) 1544 & + Atau1(:,3,5) * A3(:,5,5) 1545 & + Atau3(:,3,1) * A1(:,1,5) 1546 & + Atau3(:,3,2) * A1(:,2,5) 1547 & + Atau3(:,3,3) * A1(:,3,5) 1548 & + Atau3(:,3,4) * A1(:,4,5) 1549 & + Atau3(:,3,5) * A1(:,5,5) 1550 1551 1552c 1553c row four 1554c 1555 A0(:,4,1)= 1556 & Atau1(:,4,1) * A3(:,1,1) 1557 & + Atau1(:,4,2) * A3(:,2,1) 1558 & + Atau1(:,4,3) * A3(:,3,1) 1559 & + Atau1(:,4,4) * A3(:,4,1) 1560 & + Atau1(:,4,5) * A3(:,5,1) 1561 & + Atau3(:,4,1) * A1(:,1,1) 1562 & + Atau3(:,4,2) * A1(:,2,1) 1563 & + Atau3(:,4,3) * A1(:,3,1) 1564 & + Atau3(:,4,4) * A1(:,4,1) 1565 & + Atau3(:,4,5) * A1(:,5,1) 1566c 1567 A0(:,4,2)= 1568c & Atau1(:,4,1) * A3(:,1,2) 1569 & Atau1(:,4,2) * A3(:,2,2) 1570c & + Atau1(:,4,3) * A3(:,3,2) 1571c & + Atau1(:,4,4) * A3(:,4,2) 1572 & + Atau1(:,4,5) * A3(:,5,2) 1573 & + Atau3(:,4,1) * A1(:,1,2) 1574 & + Atau3(:,4,2) * A1(:,2,2) 1575 & + Atau3(:,4,3) * A1(:,3,2) 1576 & + Atau3(:,4,4) * A1(:,4,2) 1577 & + Atau3(:,4,5) * A1(:,5,2) 1578c 1579 A0(:,4,3)= 1580c & Atau1(:,4,1) * A3(:,1,3) 1581c & + Atau1(:,4,2) * A3(:,2,3) 1582 & Atau1(:,4,3) * A3(:,3,3) 1583c & + Atau1(:,4,4) * A3(:,4,3) 1584 & + Atau1(:,4,5) * A3(:,5,3) 1585c & + Atau3(:,4,1) * A1(:,1,3) 1586c & + Atau3(:,4,2) * A1(:,2,3) 1587 & + Atau3(:,4,3) * A1(:,3,3) 1588c & + Atau3(:,4,4) * A1(:,4,3) 1589 & + Atau3(:,4,5) * A1(:,5,3) 1590c 1591 A0(:,4,4)= 1592 & Atau1(:,4,1) * A3(:,1,4) 1593 & + Atau1(:,4,2) * A3(:,2,4) 1594 & + Atau1(:,4,3) * A3(:,3,4) 1595 & + Atau1(:,4,4) * A3(:,4,4) 1596 & + Atau1(:,4,5) * A3(:,5,4) 1597c & + Atau3(:,4,1) * A1(:,1,4) 1598c & + Atau3(:,4,2) * A1(:,2,4) 1599c & + Atau3(:,4,3) * A1(:,3,4) 1600 & + Atau3(:,4,4) * A1(:,4,4) 1601 & + Atau3(:,4,5) * A1(:,5,4) 1602c 1603 A0(:,4,5)= 1604 & Atau1(:,4,1) * A3(:,1,5) 1605 & + Atau1(:,4,2) * A3(:,2,5) 1606 & + Atau1(:,4,3) * A3(:,3,5) 1607 & + Atau1(:,4,4) * A3(:,4,5) 1608 & + Atau1(:,4,5) * A3(:,5,5) 1609 & + Atau3(:,4,1) * A1(:,1,5) 1610 & + Atau3(:,4,2) * A1(:,2,5) 1611 & + Atau3(:,4,3) * A1(:,3,5) 1612 & + Atau3(:,4,4) * A1(:,4,5) 1613 & + Atau3(:,4,5) * A1(:,5,5) 1614c 1615c row five 1616c 1617 1618 A0(:,5,1)= 1619 & Atau1(:,5,1) * A3(:,1,1) 1620 & + Atau1(:,5,2) * A3(:,2,1) 1621 & + Atau1(:,5,3) * A3(:,3,1) 1622 & + Atau1(:,5,4) * A3(:,4,1) 1623 & + Atau1(:,5,5) * A3(:,5,1) 1624 & + Atau3(:,5,1) * A1(:,1,1) 1625 & + Atau3(:,5,2) * A1(:,2,1) 1626 & + Atau3(:,5,3) * A1(:,3,1) 1627 & + Atau3(:,5,4) * A1(:,4,1) 1628 & + Atau3(:,5,5) * A1(:,5,1) 1629c 1630 A0(:,5,2)= 1631c & Atau1(:,5,1) * A3(:,1,2) 1632 & Atau1(:,5,2) * A3(:,2,2) 1633c & + Atau1(:,5,3) * A3(:,3,2) 1634c & + Atau1(:,5,4) * A3(:,4,2) 1635 & + Atau1(:,5,5) * A3(:,5,2) 1636 & + Atau3(:,5,1) * A1(:,1,2) 1637 & + Atau3(:,5,2) * A1(:,2,2) 1638 & + Atau3(:,5,3) * A1(:,3,2) 1639 & + Atau3(:,5,4) * A1(:,4,2) 1640 & + Atau3(:,5,5) * A1(:,5,2) 1641c 1642 A0(:,5,3)= 1643c & Atau1(:,5,1) * A3(:,1,3) 1644c & + Atau1(:,5,2) * A3(:,2,3) 1645 & Atau1(:,5,3) * A3(:,3,3) 1646c & + Atau1(:,5,4) * A3(:,4,3) 1647 & + Atau1(:,5,5) * A3(:,5,3) 1648c & + Atau3(:,5,1) * A1(:,1,3) 1649c & + Atau3(:,5,2) * A1(:,2,3) 1650 & + Atau3(:,5,3) * A1(:,3,3) 1651c & + Atau3(:,5,4) * A1(:,4,3) 1652 & + Atau3(:,5,5) * A1(:,5,3) 1653c 1654 A0(:,5,4)= 1655 & Atau1(:,5,1) * A3(:,1,4) 1656 & + Atau1(:,5,2) * A3(:,2,4) 1657 & + Atau1(:,5,3) * A3(:,3,4) 1658 & + Atau1(:,5,4) * A3(:,4,4) 1659 & + Atau1(:,5,5) * A3(:,5,4) 1660c & + Atau3(:,5,1) * A1(:,1,4) 1661c & + Atau3(:,5,2) * A1(:,2,4) 1662c & + Atau3(:,5,3) * A1(:,3,4) 1663 & + Atau3(:,5,4) * A1(:,4,4) 1664 & + Atau3(:,5,5) * A1(:,5,4) 1665 1666c 1667 A0(:,5,5)= 1668 & Atau1(:,5,1) * A3(:,1,5) 1669 & + Atau1(:,5,2) * A3(:,2,5) 1670 & + Atau1(:,5,3) * A3(:,3,5) 1671 & + Atau1(:,5,4) * A3(:,4,5) 1672 & + Atau1(:,5,5) * A3(:,5,5) 1673 & + Atau3(:,5,1) * A1(:,1,5) 1674 & + Atau3(:,5,2) * A1(:,2,5) 1675 & + Atau3(:,5,3) * A1(:,3,5) 1676 & + Atau3(:,5,4) * A1(:,4,5) 1677 & + Atau3(:,5,5) * A1(:,5,5) 1678 1679c 1680 else 1681 A0 = zero 1682 endif 1683c 1684 if (Navier .eq. 1) then 1685 A0(:,2,4) = A0(:,2,4) + rlm2mu - rmu 1686 A0(:,4,2) = A0(:,4,2) + rlm2mu - rmu 1687 A0(:,5,2) = A0(:,5,2) + (rlm2mu-rmu) * u3 1688 A0(:,5,4) = A0(:,5,4) + (rlm2mu-rmu) * u1 1689 endif 1690c 1691 do j = 1, nshl 1692 tmp = WdetJ * shg(:,j,1) * shg(:,j,3) 1693 do i=1,nflow 1694 do k=1,nflow 1695 BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k) 1696 enddo 1697 enddo 1698 enddo 1699c 1700c now for the 2 3 plus 3 2 1701c 1702 if(ivart.ge.2) then 1703c 1704c row one 1705c 1706 A0(:,1,1)= 1707 & Atau2(:,1,1) * A3(:,1,1) 1708 & + Atau2(:,1,2) * A3(:,2,1) 1709 & + Atau2(:,1,3) * A3(:,3,1) 1710 & + Atau2(:,1,4) * A3(:,4,1) 1711 & + Atau2(:,1,5) * A3(:,5,1) 1712 & + Atau3(:,1,1) * A2(:,1,1) 1713 & + Atau3(:,1,2) * A2(:,2,1) 1714 & + Atau3(:,1,3) * A2(:,3,1) 1715 & + Atau3(:,1,4) * A2(:,4,1) 1716 & + Atau3(:,1,5) * A2(:,5,1) 1717c 1718 A0(:,1,2) = 1719c & Atau2(:,1,1) * A3(:,1,2) 1720 & Atau2(:,1,2) * A3(:,2,2) 1721c & + Atau2(:,1,3) * A3(:,3,2) 1722c & + Atau2(:,1,4) * A3(:,4,2) 1723 & + Atau2(:,1,5) * A3(:,5,2) 1724c & + Atau3(:,1,1) * A2(:,1,2) 1725 & + Atau3(:,1,2) * A2(:,2,2) 1726c & + Atau3(:,1,3) * A2(:,3,2) 1727c & + Atau3(:,1,4) * A2(:,4,2) 1728 & + Atau3(:,1,5) * A2(:,5,2) 1729 1730c 1731 A0(:,1,3) = 1732c & Atau2(:,1,1) * A3(:,1,3) 1733c & + Atau2(:,1,2) * A3(:,2,3) 1734 & Atau2(:,1,3) * A3(:,3,3) 1735c & + Atau2(:,1,4) * A3(:,4,3) 1736 & + Atau2(:,1,5) * A3(:,5,3) 1737 & + Atau3(:,1,1) * A2(:,1,3) 1738 & + Atau3(:,1,2) * A2(:,2,3) 1739 & + Atau3(:,1,3) * A2(:,3,3) 1740 & + Atau3(:,1,4) * A2(:,4,3) 1741 & + Atau3(:,1,5) * A2(:,5,3) 1742c 1743 A0(:,1,4)= 1744 & Atau2(:,1,1) * A3(:,1,4) 1745 & + Atau2(:,1,2) * A3(:,2,4) 1746 & + Atau2(:,1,3) * A3(:,3,4) 1747 & + Atau2(:,1,4) * A3(:,4,4) 1748 & + Atau2(:,1,5) * A3(:,5,4) 1749c & + Atau3(:,1,1) * A2(:,1,4) 1750c & + Atau3(:,1,2) * A2(:,2,4) 1751c & + Atau3(:,1,3) * A2(:,3,4) 1752 & + Atau3(:,1,4) * A2(:,4,4) 1753 & + Atau3(:,1,5) * A2(:,5,4) 1754c 1755 A0(:,1,5)= 1756 & Atau2(:,1,1) * A3(:,1,5) 1757 & + Atau2(:,1,2) * A3(:,2,5) 1758 & + Atau2(:,1,3) * A3(:,3,5) 1759 & + Atau2(:,1,4) * A3(:,4,5) 1760 & + Atau2(:,1,5) * A3(:,5,5) 1761 & + Atau3(:,1,1) * A2(:,1,5) 1762 & + Atau3(:,1,2) * A2(:,2,5) 1763 & + Atau3(:,1,3) * A2(:,3,5) 1764 & + Atau3(:,1,4) * A2(:,4,5) 1765 & + Atau3(:,1,5) * A2(:,5,5) 1766c 1767c row two 1768c 1769 A0(:,2,1)= 1770 & Atau2(:,2,1) * A3(:,1,1) 1771 & + Atau2(:,2,2) * A3(:,2,1) 1772 & + Atau2(:,2,3) * A3(:,3,1) 1773 & + Atau2(:,2,4) * A3(:,4,1) 1774 & + Atau2(:,2,5) * A3(:,5,1) 1775 & + Atau3(:,2,1) * A2(:,1,1) 1776 & + Atau3(:,2,2) * A2(:,2,1) 1777 & + Atau3(:,2,3) * A2(:,3,1) 1778 & + Atau3(:,2,4) * A2(:,4,1) 1779 & + Atau3(:,2,5) * A2(:,5,1) 1780c 1781 A0(:,2,2)= 1782c & Atau2(:,2,1) * A3(:,1,2) 1783 & Atau2(:,2,2) * A3(:,2,2) 1784c & + Atau2(:,2,3) * A3(:,3,2) 1785c & + Atau2(:,2,4) * A3(:,4,2) 1786 & + Atau2(:,2,5) * A3(:,5,2) 1787c & + Atau3(:,2,1) * A2(:,1,2) 1788 & + Atau3(:,2,2) * A2(:,2,2) 1789c & + Atau3(:,2,3) * A2(:,3,2) 1790c & + Atau3(:,2,4) * A2(:,4,2) 1791 & + Atau3(:,2,5) * A2(:,5,2) 1792c 1793 A0(:,2,3)= 1794c & Atau2(:,2,1) * A3(:,1,3) 1795c & + Atau2(:,2,2) * A3(:,2,3) 1796 & Atau2(:,2,3) * A3(:,3,3) 1797c & + Atau2(:,2,4) * A3(:,4,3) 1798 & + Atau2(:,2,5) * A3(:,5,3) 1799 & + Atau3(:,2,1) * A2(:,1,3) 1800 & + Atau3(:,2,2) * A2(:,2,3) 1801 & + Atau3(:,2,3) * A2(:,3,3) 1802 & + Atau3(:,2,4) * A2(:,4,3) 1803 & + Atau3(:,2,5) * A2(:,5,3) 1804c 1805 A0(:,2,4)= 1806 & Atau2(:,2,1) * A3(:,1,4) 1807 & + Atau2(:,2,2) * A3(:,2,4) 1808 & + Atau2(:,2,3) * A3(:,3,4) 1809 & + Atau2(:,2,4) * A3(:,4,4) 1810 & + Atau2(:,2,5) * A3(:,5,4) 1811c & + Atau3(:,2,1) * A2(:,1,4) 1812c & + Atau3(:,2,2) * A2(:,2,4) 1813c & + Atau3(:,2,3) * A2(:,3,4) 1814 & + Atau3(:,2,4) * A2(:,4,4) 1815 & + Atau3(:,2,5) * A2(:,5,4) 1816 1817c 1818 A0(:,2,5)= 1819 & Atau2(:,2,1) * A3(:,1,5) 1820 & + Atau2(:,2,2) * A3(:,2,5) 1821 & + Atau2(:,2,3) * A3(:,3,5) 1822 & + Atau2(:,2,4) * A3(:,4,5) 1823 & + Atau2(:,2,5) * A3(:,5,5) 1824 & + Atau3(:,2,1) * A2(:,1,5) 1825 & + Atau3(:,2,2) * A2(:,2,5) 1826 & + Atau3(:,2,3) * A2(:,3,5) 1827 & + Atau3(:,2,4) * A2(:,4,5) 1828 & + Atau3(:,2,5) * A2(:,5,5) 1829 1830c 1831c row three 1832c 1833 A0(:,3,1)= 1834 & Atau2(:,3,1) * A3(:,1,1) 1835 & + Atau2(:,3,2) * A3(:,2,1) 1836 & + Atau2(:,3,3) * A3(:,3,1) 1837 & + Atau2(:,3,4) * A3(:,4,1) 1838 & + Atau2(:,3,5) * A3(:,5,1) 1839 & + Atau3(:,3,1) * A2(:,1,1) 1840 & + Atau3(:,3,2) * A2(:,2,1) 1841 & + Atau3(:,3,3) * A2(:,3,1) 1842 & + Atau3(:,3,4) * A2(:,4,1) 1843 & + Atau3(:,3,5) * A2(:,5,1) 1844 1845c 1846 A0(:,3,2)= 1847c & Atau2(:,3,1) * A3(:,1,2) 1848 & Atau2(:,3,2) * A3(:,2,2) 1849c & + Atau2(:,3,3) * A3(:,3,2) 1850c & + Atau2(:,3,4) * A3(:,4,2) 1851 & + Atau2(:,3,5) * A3(:,5,2) 1852c & + Atau3(:,3,1) * A2(:,1,2) 1853 & + Atau3(:,3,2) * A2(:,2,2) 1854c & + Atau3(:,3,3) * A2(:,3,2) 1855c & + Atau3(:,3,4) * A2(:,4,2) 1856 & + Atau3(:,3,5) * A2(:,5,2) 1857c 1858 A0(:,3,3)= 1859c & Atau2(:,3,1) * A3(:,1,3) 1860c & + Atau2(:,3,2) * A3(:,2,3) 1861 & Atau2(:,3,3) * A3(:,3,3) 1862c & + Atau2(:,3,4) * A3(:,4,3) 1863 & + Atau2(:,3,5) * A3(:,5,3) 1864 & + Atau3(:,3,1) * A2(:,1,3) 1865 & + Atau3(:,3,2) * A2(:,2,3) 1866 & + Atau3(:,3,3) * A2(:,3,3) 1867 & + Atau3(:,3,4) * A2(:,4,3) 1868 & + Atau3(:,3,5) * A2(:,5,3) 1869c 1870 A0(:,3,4)= 1871 & Atau2(:,3,1) * A3(:,1,4) 1872 & + Atau2(:,3,2) * A3(:,2,4) 1873 & + Atau2(:,3,3) * A3(:,3,4) 1874 & + Atau2(:,3,4) * A3(:,4,4) 1875 & + Atau2(:,3,5) * A3(:,5,4) 1876c & + Atau3(:,3,1) * A2(:,1,4) 1877c & + Atau3(:,3,2) * A2(:,2,4) 1878c & + Atau3(:,3,3) * A2(:,3,4) 1879 & + Atau3(:,3,4) * A2(:,4,4) 1880 & + Atau3(:,3,5) * A2(:,5,4) 1881c 1882 A0(:,3,5)= 1883 & Atau2(:,3,1) * A3(:,1,5) 1884 & + Atau2(:,3,2) * A3(:,2,5) 1885 & + Atau2(:,3,3) * A3(:,3,5) 1886 & + Atau2(:,3,4) * A3(:,4,5) 1887 & + Atau2(:,3,5) * A3(:,5,5) 1888 & + Atau3(:,3,1) * A2(:,1,5) 1889 & + Atau3(:,3,2) * A2(:,2,5) 1890 & + Atau3(:,3,3) * A2(:,3,5) 1891 & + Atau3(:,3,4) * A2(:,4,5) 1892 & + Atau3(:,3,5) * A2(:,5,5) 1893 1894 1895c 1896c row four 1897c 1898 A0(:,4,1)= 1899 & Atau2(:,4,1) * A3(:,1,1) 1900 & + Atau2(:,4,2) * A3(:,2,1) 1901 & + Atau2(:,4,3) * A3(:,3,1) 1902 & + Atau2(:,4,4) * A3(:,4,1) 1903 & + Atau2(:,4,5) * A3(:,5,1) 1904 & + Atau3(:,4,1) * A2(:,1,1) 1905 & + Atau3(:,4,2) * A2(:,2,1) 1906 & + Atau3(:,4,3) * A2(:,3,1) 1907 & + Atau3(:,4,4) * A2(:,4,1) 1908 & + Atau3(:,4,5) * A2(:,5,1) 1909c 1910 A0(:,4,2)= 1911c & Atau2(:,4,1) * A3(:,1,2) 1912 & Atau2(:,4,2) * A3(:,2,2) 1913c & + Atau2(:,4,3) * A3(:,3,2) 1914c & + Atau2(:,4,4) * A3(:,4,2) 1915 & + Atau2(:,4,5) * A3(:,5,2) 1916c & + Atau3(:,4,1) * A2(:,1,2) 1917 & + Atau3(:,4,2) * A2(:,2,2) 1918c & + Atau3(:,4,3) * A2(:,3,2) 1919c & + Atau3(:,4,4) * A2(:,4,2) 1920 & + Atau3(:,4,5) * A2(:,5,2) 1921c 1922 A0(:,4,3)= 1923c & Atau2(:,4,1) * A3(:,1,3) 1924c & + Atau2(:,4,2) * A3(:,2,3) 1925 & Atau2(:,4,3) * A3(:,3,3) 1926c & + Atau2(:,4,4) * A3(:,4,3) 1927 & + Atau2(:,4,5) * A3(:,5,3) 1928 & + Atau3(:,4,1) * A2(:,1,3) 1929 & + Atau3(:,4,2) * A2(:,2,3) 1930 & + Atau3(:,4,3) * A2(:,3,3) 1931 & + Atau3(:,4,4) * A2(:,4,3) 1932 & + Atau3(:,4,5) * A2(:,5,3) 1933c 1934 A0(:,4,4)= 1935 & Atau2(:,4,1) * A3(:,1,4) 1936 & + Atau2(:,4,2) * A3(:,2,4) 1937 & + Atau2(:,4,3) * A3(:,3,4) 1938 & + Atau2(:,4,4) * A3(:,4,4) 1939 & + Atau2(:,4,5) * A3(:,5,4) 1940c & + Atau3(:,4,1) * A2(:,1,4) 1941c & + Atau3(:,4,2) * A2(:,2,4) 1942c & + Atau3(:,4,3) * A2(:,3,4) 1943 & + Atau3(:,4,4) * A2(:,4,4) 1944 & + Atau3(:,4,5) * A2(:,5,4) 1945c 1946 A0(:,4,5)= 1947 & Atau2(:,4,1) * A3(:,1,5) 1948 & + Atau2(:,4,2) * A3(:,2,5) 1949 & + Atau2(:,4,3) * A3(:,3,5) 1950 & + Atau2(:,4,4) * A3(:,4,5) 1951 & + Atau2(:,4,5) * A3(:,5,5) 1952 & + Atau3(:,4,1) * A2(:,1,5) 1953 & + Atau3(:,4,2) * A2(:,2,5) 1954 & + Atau3(:,4,3) * A2(:,3,5) 1955 & + Atau3(:,4,4) * A2(:,4,5) 1956 & + Atau3(:,4,5) * A2(:,5,5) 1957c 1958c row five 1959c 1960 1961 A0(:,5,1)= 1962 & Atau2(:,5,1) * A3(:,1,1) 1963 & + Atau2(:,5,2) * A3(:,2,1) 1964 & + Atau2(:,5,3) * A3(:,3,1) 1965 & + Atau2(:,5,4) * A3(:,4,1) 1966 & + Atau2(:,5,5) * A3(:,5,1) 1967 & + Atau3(:,5,1) * A2(:,1,1) 1968 & + Atau3(:,5,2) * A2(:,2,1) 1969 & + Atau3(:,5,3) * A2(:,3,1) 1970 & + Atau3(:,5,4) * A2(:,4,1) 1971 & + Atau3(:,5,5) * A2(:,5,1) 1972c 1973 A0(:,5,2)= 1974c & Atau2(:,5,1) * A3(:,1,2) 1975 & Atau2(:,5,2) * A3(:,2,2) 1976c & + Atau2(:,5,3) * A3(:,3,2) 1977c & + Atau2(:,5,4) * A3(:,4,2) 1978 & + Atau2(:,5,5) * A3(:,5,2) 1979c & + Atau3(:,5,1) * A2(:,1,2) 1980 & + Atau3(:,5,2) * A2(:,2,2) 1981c & + Atau3(:,5,3) * A2(:,3,2) 1982c & + Atau3(:,5,4) * A2(:,4,2) 1983 & + Atau3(:,5,5) * A2(:,5,2) 1984c 1985 A0(:,5,3)= 1986c & Atau2(:,5,1) * A3(:,1,3) 1987c & + Atau2(:,5,2) * A3(:,2,3) 1988 & Atau2(:,5,3) * A3(:,3,3) 1989c & + Atau2(:,5,4) * A3(:,4,3) 1990 & + Atau2(:,5,5) * A3(:,5,3) 1991 & + Atau3(:,5,1) * A2(:,1,3) 1992 & + Atau3(:,5,2) * A2(:,2,3) 1993 & + Atau3(:,5,3) * A2(:,3,3) 1994 & + Atau3(:,5,4) * A2(:,4,3) 1995 & + Atau3(:,5,5) * A2(:,5,3) 1996c 1997 A0(:,5,4)= 1998 & Atau2(:,5,1) * A3(:,1,4) 1999 & + Atau2(:,5,2) * A3(:,2,4) 2000 & + Atau2(:,5,3) * A3(:,3,4) 2001 & + Atau2(:,5,4) * A3(:,4,4) 2002 & + Atau2(:,5,5) * A3(:,5,4) 2003c & + Atau3(:,5,1) * A2(:,1,4) 2004c & + Atau3(:,5,2) * A2(:,2,4) 2005c & + Atau3(:,5,3) * A2(:,3,4) 2006 & + Atau3(:,5,4) * A2(:,4,4) 2007 & + Atau3(:,5,5) * A2(:,5,4) 2008 2009c 2010 A0(:,5,5)= 2011 & Atau2(:,5,1) * A3(:,1,5) 2012 & + Atau2(:,5,2) * A3(:,2,5) 2013 & + Atau2(:,5,3) * A3(:,3,5) 2014 & + Atau2(:,5,4) * A3(:,4,5) 2015 & + Atau2(:,5,5) * A3(:,5,5) 2016 & + Atau3(:,5,1) * A2(:,1,5) 2017 & + Atau3(:,5,2) * A2(:,2,5) 2018 & + Atau3(:,5,3) * A2(:,3,5) 2019 & + Atau3(:,5,4) * A2(:,4,5) 2020 & + Atau3(:,5,5) * A2(:,5,5) 2021 2022c 2023 else 2024 A0 = zero 2025 endif 2026c 2027 if (Navier .eq. 1) then 2028 A0(:,3,4) = A0(:,3,4) + rlm2mu - rmu 2029 A0(:,4,3) = A0(:,4,3) + rlm2mu - rmu 2030 A0(:,5,3) = A0(:,5,3) + (rlm2mu-rmu) * u3 2031 A0(:,5,4) = A0(:,5,4) + (rlm2mu-rmu) * u2 2032 endif 2033 2034 do j = 1, nshl 2035 tmp = WdetJ * shg(:,j,2) * shg(:,j,3) 2036 do i=1,nflow 2037 do k=1,nflow 2038 BDiagl(:,j,i,k) = BDiagl(:,j,i,k) + tmp * A0(:,i,k) 2039 enddo 2040 enddo 2041 enddo 2042 2043 ttim(30) = ttim(30) + tmr() 2044 2045 return 2046 end 2047 2048 2049