1 subroutine e3stab (rho, u1, u2, 2 & u3, dxidx, rLui, 3 & rmu, tauC, tauM, 4 & tauBar, uBar ) 5c 6c---------------------------------------------------------------------- 7c 8c This routine computes the diagonal Tau for least-squares operator. 9c Diagonal tau proposed by Shakib. 10c 11c input: 12c u1 (npro) : x1-velocity component 13c u2 (npro) : x2-velocity component 14c u3 (npro) : x3-velocity component 15c dxidx (npro,nsd,nsd) : inverse of deformation gradient 16c rLui (npro,nsd) : least-squares residual vector 17c 18c output: 19c tauC (npro) : continuity tau 20c tauM (npro) : momentum tau 21c tauBar (npro) : additional tau 22c uBar (npro,nsd) : modified velocity 23c 24c Zdenek Johan, Summer 1990. (Modified from e2tau.f) 25c Zdenek Johan, Winter 1991. (Fortran 90) 26c---------------------------------------------------------------------- 27c 28 include "common.h" 29c 30 dimension rho(npro), u1(npro), 31 & u2(npro), u3(npro), 32 & dxidx(npro,nsd,nsd), 33 & rLui(npro,nsd), 34 & tauC(npro), tauM(npro), tauBar(npro), 35 & rmu(npro), uBar(npro,3), unorm(npro) 36 37c 38 dimension gijd(npro,6), fact(npro), rnu(npro), 39 & rhoinv(npro) 40c 41c 42c.... get the metric tensor 43c 44 call e3gijd( dxidx, gijd ) 45c 46c... higher order element diffusive correction 47c 48 if (ipord == 1) then 49 fff = 36.0d0 50 else if (ipord == 2) then 51 fff = 60.0d0 52c fff = 36.0d0 53 else if (ipord == 3) then 54 fff = 128.0d0 55c fff = 144.0d0 56 endif 57 58 omegasq=zero 59 if(matflg(6,1).eq.1) omegasq = datmat(1,6,1)**2 60 . +datmat(2,6,1)**2 61 . +datmat(3,6,1)**2 62 rhoinv=one/rho 63 rnu=rmu*rhoinv 64 65 if(itau.eq.0) then ! original tau 66c 67c... momentum tau 68c 69 70!MR CHANGE 71 dts= Dtgl*dtsfct ! Dtgl = (time step)^-1 72! dts= min(Dtgl,28800.0d0)*dtsfct ! Dtgl = (time step)^-1 !28800 = 1600*180 / 10 73! dts= min(Dtgl,2880.0d0)*dtsfct ! Dtgl = (time step)^-1 !2880 = 1600*180 / 100 74! dts= min(Dtgl,288.0d0)*dtsfct ! Dtgl = (time step)^-1 !288 = 1600*180 / 1000 75!MR CHANGE 76 77 tauM = ( (two*dts)**2 78 3 + ( u1 * ( gijd(:,1) * u1 79 4 + gijd(:,4) * u2 80 5 + gijd(:,6) * u3 ) 81 6 + u2 * ( gijd(:,4) * u1 82 7 + gijd(:,2) * u2 83 8 + gijd(:,5) * u3 ) 84 9 + u3 * ( gijd(:,6) * u1 85 a + gijd(:,5) * u2 86 1 + gijd(:,3) * u3 ) ) ) 87 2 + fff * rnu** 2 88 3 * ( gijd(:,1) ** 2 89 4 + gijd(:,2) ** 2 90 5 + gijd(:,3) ** 2 91 6 + 2. 92 7 * ( gijd(:,4) ** 2 93 8 + gijd(:,5) ** 2 94 9 + gijd(:,6) ** 2 ) 95 b +omegasq) 96 97 fact = sqrt(tauM) 98 dtsi=one/dts 99 ff=taucfct/dtsfct 100 tauC =rho* pt125*fact/(gijd(:,1)+gijd(:,2)+gijd(:,3))*ff 101 tauM = one/fact 102 else if(itau.eq.1) then ! new tau 103 104c 105c determinant of gijd 106c 107 fact = gijd(:,1) * gijd(:,2) * gijd(:,3) 108 & - gijd(:,2) * gijd(:,6) * gijd(:,6) 109 & - gijd(:,1) * gijd(:,5) * gijd(:,5) 110 & - gijd(:,3) * gijd(:,4) * gijd(:,4) 111 & + gijd(:,6) * gijd(:,4) * gijd(:,5) * two 112 113c 114c put 1/2u*h1 = sqrt(u_i g^{ij} u_j) into tau_M note inverse is calculated 115c on the fly here from cofactors over the determinent dotted from left and 116c right with u 117c 118 119 tauM = 120 1 u1 * ( (gijd(:,2)*gijd(:,3)-gijd(:,5)*gijd(:,5)) * u1 121 2 + two * (gijd(:,5)*gijd(:,6)-gijd(:,4)*gijd(:,3)) * u2 122 3 + two * (gijd(:,4)*gijd(:,5)-gijd(:,6)*gijd(:,2)) * u3) 123 1 + u2 * ( (gijd(:,1)*gijd(:,3)-gijd(:,6)*gijd(:,6)) * u2 124 3 + two * (gijd(:,4)*gijd(:,6)-gijd(:,1)*gijd(:,5)) * u3) 125 1 + u3 * ( (gijd(:,1)*gijd(:,2)-gijd(:,4)*gijd(:,4)) * u3) 126 tauM=fact/taum ! here we have (u_i g^{ij} u^j)^{-1} approx 4/u^2h^2 127c 128c we can calculate tauC more efficiently now 129c 130 tauC=tauM*(one+tauM*rmu*rmu) 131 tauC=one/tauC 132 tauC=taucfct*sqrt(tauC) 133c 134c 135c... momentum tau 136c 137c 138c this tau needs a u/h instead of a u*h so we contract with g_{ij} as 139c follows (i.e. u_i g_{ij} u_j approx u^2/(h^2)/4) 140c 141 fact = 142 3 u1 * ( gijd(:,1) * u1 143 4 + gijd(:,4) * u2 144 5 + gijd(:,6) * u3 ) 145 6 + u2 * ( gijd(:,4) * u1 146 7 + gijd(:,2) * u2 147 8 + gijd(:,5) * u3 ) 148 9 + u3 * ( gijd(:,6) * u1 149 a + gijd(:,5) * u2 150 1 + gijd(:,3) * u3 ) 151c 152c first limit dt effect on tau from causing trouble if user drops CFL below 153c .05 (this could cause loss of spatial stability) 154c 155 velsq=vel*vel 156 unorm = (u1*u1+u2*u2+u3*u3)/velsq 157 dtsfsq=dtsfct*dtsfct 158 dt=one/Dtgl 159 taubar= dtsfsq/( dt*dt + .01*unorm/fact) ! never gets above (C_1 20*u_inf/h)^2 160c 161c this means tau will never get below h/(20*C_1*u) no matter what time step 162c you choose. The 0.01 constant comes from minCFL=.05=> .05*.05*4 (where the 163c 4 comes from the bi-unit mapping). If you want to limit sooner the formula 164c would be ".01-factor"=minCFL^2*4 165c 166 167 tauM = rho ** 2 168 1 * ( four*taubar + fact 169 2 + fff * rmu** 2 170 3 * ( gijd(:,1) ** 2 171 4 + gijd(:,2) ** 2 172 5 + gijd(:,3) ** 2 173 6 + 2. 174 7 * ( gijd(:,4) ** 2 175 8 + gijd(:,5) ** 2 176 9 + gijd(:,6) ** 2 ) ) 177 b +omegasq) 178 fact=sqrt(tauM) 179cdebugcheck tauBar = pt125*fact/(gijd(:,1)+gijd(:,2)+gijd(:,3)) !*dtsi 180 181 tauM=one/fact ! turn it right side up. 182 else if(itau.eq.2) then ! new tau different continuity h 183 184 unorm = (u1*u1+u2*u2+u3*u3) 185 186 tauM=(gijd(:,1)+gijd(:,2)+gijd(:,3))/unorm ! here we have 4/u^2h^2 187c 188c we can calculate tauC more efficiently now 189c 190 tauC=tauM*(one+tauM*rmu*rmu) 191 tauC=one/tauC 192 tauC=sqrt(tauC)*taucfct 193c 194c 195c... momentum tau 196c 197c 198c this tau needs a u/h instead of a u*h so we contract with g_{ij} as 199c follows (i.e. u_i g_{ij} u_j approx u^2/(h^2)/4) 200c 201 fact = 202 3 u1 * ( gijd(:,1) * u1 203 4 + gijd(:,4) * u2 204 5 + gijd(:,6) * u3 ) 205 6 + u2 * ( gijd(:,4) * u1 206 7 + gijd(:,2) * u2 207 8 + gijd(:,5) * u3 ) 208 9 + u3 * ( gijd(:,6) * u1 209 a + gijd(:,5) * u2 210 1 + gijd(:,3) * u3 ) 211c 212c first limit dt effect on tau from causing trouble if user drops CFL below 213c .05 (this could cause loss of spatial stability) 214c 215 velsq=vel*vel 216 dtsfsq=dtsfct*dtsfct 217 dt=one/Dtgl 218 unorm=unorm/velsq 219 taubar= dtsfsq/( dt*dt + .01*unorm/fact) ! never gets above (C_1 20*u_inf/h)^2 220c 221c this means tau will never get below h/(20*C_1*u) no matter what time step 222c you choose. The 0.01 constant comes from minCFL=.05=> .05*.05*4 (where the 223c 4 comes from the bi-unit mapping). If you want to limit sooner the formula 224c would be ".01-factor"=minCFL^2*4 225c 226 227 tauM = rho ** 2 228 1 * ( four*taubar + fact 229 2 + fff * rmu** 2 230 3 * ( gijd(:,1) ** 2 231 4 + gijd(:,2) ** 2 232 5 + gijd(:,3) ** 2 233 6 + 2. 234 7 * ( gijd(:,4) ** 2 235 8 + gijd(:,5) ** 2 236 9 + gijd(:,6) ** 2 ) ) 237 b +omegasq) 238 fact=sqrt(tauM) 239c tauBar = pt125*fact/(gijd(:,1)+gijd(:,2)+gijd(:,3)) !*dtsi 240 241 tauM=one/fact ! turn it right side up. 242 else if(itau.eq.3) then ! compressible tau 243 244c 245c determinant of gijd 246c 247 fact = gijd(:,1) * gijd(:,2) * gijd(:,3) 248 & - gijd(:,2) * gijd(:,6) * gijd(:,6) 249 & - gijd(:,1) * gijd(:,5) * gijd(:,5) 250 & - gijd(:,3) * gijd(:,4) * gijd(:,4) 251 & + gijd(:,6) * gijd(:,4) * gijd(:,5) * two 252 253c 254c put 1/2u*h1 = sqrt(u_i g^{ij} u_j) into tau_M note inverse is calculated 255c on the fly here from cofactors over the determinent dotted from left and 256c right with u 257c 258 259 tauM = 260 1 u1 * ( (gijd(:,2)*gijd(:,3)-gijd(:,5)*gijd(:,5)) * u1 261 2 + two * (gijd(:,5)*gijd(:,6)-gijd(:,4)*gijd(:,3)) * u2 262 3 + two * (gijd(:,4)*gijd(:,5)-gijd(:,6)*gijd(:,2)) * u3) 263 1 + u2 * ( (gijd(:,1)*gijd(:,3)-gijd(:,6)*gijd(:,6)) * u2 264 3 + two * (gijd(:,4)*gijd(:,6)-gijd(:,1)*gijd(:,5)) * u3) 265 1 + u3 * ( (gijd(:,1)*gijd(:,2)-gijd(:,4)*gijd(:,4)) * u3) 266c 267c we can calculate tauC more efficiently now 268c 269 tauM=sqrt(tauM/fact)*two 270 tauC=pt5*tauM*min(one,pt5*tauM/rmu)*taucfct 271c 272c 273c... momentum tau 274c 275c 276c this tau needs a u/h instead of a u*h so we contract with g_{ij} as 277c follows (i.e. u_i g_{ij} u_j approx u^2/(h^2)/4) 278c 279 fact = 280 3 u1 * ( gijd(:,1) * u1 281 4 + gijd(:,4) * u2 282 5 + gijd(:,6) * u3 ) 283 6 + u2 * ( gijd(:,4) * u1 284 7 + gijd(:,2) * u2 285 8 + gijd(:,5) * u3 ) 286 9 + u3 * ( gijd(:,6) * u1 287 a + gijd(:,5) * u2 288 1 + gijd(:,3) * u3 ) 289 fact=one/sqrt(fact) 290 291 unorm = (u1*u1+u2*u2+u3*u3) 292 293 dts= one/( Dtgl*dtsfct) 294 tauM =min(dts,min(fact,fact*fact*unorm*pt33/rmu)) 295 endif 296c 297c.... calculate tauBar 298c 299 tauBar = rLui(:,1) * ( gijd(:,1) * rLui(:,1) 300 & + gijd(:,4) * rLui(:,2) 301 & + gijd(:,6) * rLui(:,3) ) 302 & + rLui(:,2) * ( gijd(:,4) * rLui(:,1) 303 & + gijd(:,2) * rLui(:,2) 304 & + gijd(:,5) * rLui(:,3) ) 305 & + rLui(:,3) * ( gijd(:,6) * rLui(:,1) 306 & + gijd(:,5) * rLui(:,2) 307 & + gijd(:,3) * rLui(:,3) ) 308 where ( tauBar .ne. 0.0 ) 309 tauBar = tauM / sqrt(tauBar) 310 endwhere 311 312c 313c.... compute the modified velocity, uBar 314c 315 uBar(:,1) = u1 - tauM * rLui(:,1)*rhoinv 316 uBar(:,2) = u2 - tauM * rLui(:,2)*rhoinv 317 uBar(:,3) = u3 - tauM * rLui(:,3)*rhoinv 318c 319c.... return 320c 321 return 322 end 323 324c----------------------------------------------------------------------- 325c 326c Momentum tau 327c 328c----------------------------------------------------------------------- 329 subroutine e3uBar (rho, ui, dxidx, 330 & rLui, rmu, uBar ) 331 332 include "common.h" 333 334 real*8 rho(npro), ui(npro,nsd), 335 & dxidx(npro,nsd,nsd), rLui(npro,nsd), 336 & rmu(npro), uBar(npro,nsd) 337 338 real*8 gijd(npro,6), tauM(npro) 339 340c 341c.... get the metric tensor 342c 343 call e3gijd( dxidx, gijd ) 344c 345c.... higher order element diffusive correction 346c 347 if (ipord == 1) then 348 fff = 36.0d0 349 else if (ipord == 2) then 350 fff = 60.0d0 351 else if (ipord == 3) then 352 fff = 128.0d0 353 endif 354 355!MR CHANGE 356 dts = (Dtgl*dtsfct) 357! dts= min(Dtgl,28800.0d0)*dtsfct !28800 = 1600*180 / 10 358! dts= min(Dtgl,2880.0d0)*dtsfct !2880 = 1600*180 / 100 359! dts= min(Dtgl,288.0d0)*dtsfct !288 = 1600*180 / 1000 360!MR CHANGE 361 362 tauM = rho ** 2 363 1 * ( (two*dts)**2 364 3 + ( ui(:,1) * ( gijd(:,1) * ui(:,1) 365 4 + gijd(:,4) * ui(:,2) 366 5 + gijd(:,6) * ui(:,3) ) 367 6 + ui(:,2) * ( gijd(:,4) * ui(:,1) 368 7 + gijd(:,2) * ui(:,2) 369 8 + gijd(:,5) * ui(:,3) ) 370 9 + ui(:,3) * ( gijd(:,6) * ui(:,1) 371 a + gijd(:,5) * ui(:,2) 372 1 + gijd(:,3) * ui(:,3) ) ) ) 373 2 + fff * rmu** 2 374 3 * ( gijd(:,1) ** 2 375 4 + gijd(:,2) ** 2 376 5 + gijd(:,3) ** 2 377 6 + 2. 378 7 * ( gijd(:,4) ** 2 379 8 + gijd(:,5) ** 2 380 9 + gijd(:,6) ** 2 ) ) 381 382 tauM = one/sqrt(tauM) 383c 384c.... compute the modified velocity, uBar 385c 386 uBar(:,1) = ui(:,1) - tauM * rLui(:,1) 387 uBar(:,2) = ui(:,2) - tauM * rLui(:,2) 388 uBar(:,3) = ui(:,3) - tauM * rLui(:,3) 389 390 return 391 end 392 393c----------------------------------------------------------------------- 394c get the metric tensor g_{ij}=xi_{k,i} xi_{k,j}. 395c----------------------------------------------------------------------- 396 subroutine e3gijd( dxidx, gijd ) 397 398 include "common.h" 399 400 real*8 dxidx(npro,nsd,nsd), gijd(npro,6), 401 & tmp1(npro), tmp2(npro), 402 & tmp3(npro) 403c 404c form metric tensor g_{ij}=xi_{k,i} xi_{k,j}. It is a symmetric 405c tensor so we only form 6 components and use symmetric matrix numbering. 406c 407 if (lcsyst .ge. 2) then ! note this makes wedges like hexs..should 408c be corrected later 409 410 gijd(:,1) = dxidx(:,1,1) * dxidx(:,1,1) 411 & + dxidx(:,2,1) * dxidx(:,2,1) 412 & + dxidx(:,3,1) * dxidx(:,3,1) 413c 414 gijd(:,4) = dxidx(:,1,1) * dxidx(:,1,2) 415 & + dxidx(:,2,1) * dxidx(:,2,2) 416 & + dxidx(:,3,1) * dxidx(:,3,2) 417c 418 gijd(:,2) = dxidx(:,1,2) * dxidx(:,1,2) 419 & + dxidx(:,2,2) * dxidx(:,2,2) 420 & + dxidx(:,3,2) * dxidx(:,3,2) 421c 422 gijd(:,5) = dxidx(:,1,2) * dxidx(:,1,3) 423 & + dxidx(:,2,2) * dxidx(:,2,3) 424 & + dxidx(:,3,2) * dxidx(:,3,3) 425c 426 gijd(:,6) = dxidx(:,1,1) * dxidx(:,1,3) 427 & + dxidx(:,2,1) * dxidx(:,2,3) 428 & + dxidx(:,3,1) * dxidx(:,3,3) 429c 430 gijd(:,3) = dxidx(:,1,3) * dxidx(:,1,3) 431 & + dxidx(:,2,3) * dxidx(:,2,3) 432 & + dxidx(:,3,3) * dxidx(:,3,3) 433c 434 else if (lcsyst .eq. 1) then 435c 436c There is an invariance problem with tets 437c It is fixed by the following modifications to gijd 438c 439 440 c1 = 1.259921049894873D+00 441 c2 = 6.299605249474365D-01 442c 443 tmp1(:) = c1 * dxidx(:,1,1)+c2 *(dxidx(:,2,1)+dxidx(:,3,1)) 444 tmp2(:) = c1 * dxidx(:,2,1)+c2 *(dxidx(:,1,1)+dxidx(:,3,1)) 445 tmp3(:) = c1 * dxidx(:,3,1)+c2 *(dxidx(:,1,1)+dxidx(:,2,1)) 446 gijd(:,1) = dxidx(:,1,1) * tmp1 447 1 + dxidx(:,2,1) * tmp2 448 2 + dxidx(:,3,1) * tmp3 449c 450 tmp1(:) = c1 * dxidx(:,1,2)+c2 *(dxidx(:,2,2)+dxidx(:,3,2)) 451 tmp2(:) = c1 * dxidx(:,2,2)+c2 *(dxidx(:,1,2)+dxidx(:,3,2)) 452 tmp3(:) = c1 * dxidx(:,3,2)+c2 *(dxidx(:,1,2)+dxidx(:,2,2)) 453 gijd(:,2) = dxidx(:,1,2) * tmp1 454 1 + dxidx(:,2,2) * tmp2 455 2 + dxidx(:,3,2) * tmp3 456c 457 gijd(:,4) = dxidx(:,1,1) * tmp1 458 1 + dxidx(:,2,1) * tmp2 459 2 + dxidx(:,3,1) * tmp3 460c 461 tmp1(:) = c1 * dxidx(:,1,3)+c2 *(dxidx(:,2,3)+dxidx(:,3,3)) 462 tmp2(:) = c1 * dxidx(:,2,3)+c2 *(dxidx(:,1,3)+dxidx(:,3,3)) 463 tmp3(:) = c1 * dxidx(:,3,3)+c2 *(dxidx(:,1,3)+dxidx(:,2,3)) 464 gijd(:,3) = dxidx(:,1,3) * tmp1 465 1 + dxidx(:,2,3) * tmp2 466 2 + dxidx(:,3,3) * tmp3 467c 468 gijd(:,5) = dxidx(:,1,2) * tmp1 469 1 + dxidx(:,2,2) * tmp2 470 2 + dxidx(:,3,2) * tmp3 471c 472 gijd(:,6) = dxidx(:,1,1) * tmp1 473 1 + dxidx(:,2,1) * tmp2 474 2 + dxidx(:,3,1) * tmp3 475c 476 else 477 write(*,*) 'lcsyst eq',lcsyst,'not supported' 478 stop 479 endif 480 481 return 482 end 483 484c------------------------------------------------------------------------ 485c 486c calculate the stabilization for the advection-diffusion equation 487c 488c------------------------------------------------------------------------ 489 subroutine e3StabSclr (uMod, dxidx, tauT, diffus, srcR, giju, 490 & srcRat ) 491c 492c 493 include "common.h" 494c 495 real*8 rho(npro), uMod(npro,nsd), 496 & dxidx(npro,nsd,nsd), diffus(npro), 497 & tauT(npro), srcR(npro) 498 499c 500 real*8 gijd(npro,6), giju(npro,6), 501 & tmp1(npro), tmp2(npro), 502 & tmp3(npro), fact(npro), 503 & srcRat(npro) 504 505 real*8 fff 506 if(ivart.eq.1) then 507 tauT=zero 508 return 509 endif 510c 511c.... get the metric tensor 512c 513 call e3gijd( dxidx, gijd ) 514c 515c... momentum tau 516c 517c 518c... higher order element diffusive correction 519c 520 if (ipord == 1) then 521 fff = 9.0d0 522 else if (ipord == 2) then 523 fff = 36.0d0 524 else if (ipord == 3) then 525 fff = 64.0d0 526 endif 527 528!MR CHANGE 529 dts = (Dtgl*dtsfct) 530! dts= min(Dtgl,28800.0d0)*dtsfct !28800 = 1600*180 / 10 531! dts= min(Dtgl,2880.0d0)*dtsfct !2880 = 1600*180 / 100 532! dts= min(Dtgl,288.0d0)*dtsfct !288 = 1600*180 / 1000 533!MR CHANGE 534 535c if(iRANS.ne.-2) srcRat=srcR 536 tauT = 537 1 (two*dts)**2 538 2 + srcRat ** 2 539 3 + uMod(:,1) * ( gijd(:,1) * uMod(:,1) 540 4 + gijd(:,4) * uMod(:,2) 541 5 + gijd(:,6) * uMod(:,3) ) 542 6 + uMod(:,2) * ( gijd(:,4) * uMod(:,1) 543 7 + gijd(:,2) * uMod(:,2) 544 8 + gijd(:,5) * uMod(:,3) ) 545 9 + uMod(:,3) * ( gijd(:,6) * uMod(:,1) 546 a + gijd(:,5) * uMod(:,2) 547 1 + gijd(:,3) * uMod(:,3) ) 548 2 + fff * diffus(:)** 2 549 3 * ( gijd(:,1) ** 2 550 4 + gijd(:,2) ** 2 551 5 + gijd(:,3) ** 2 552 6 + 2. 553 7 * ( gijd(:,4) ** 2 554 8 + gijd(:,5) ** 2 555 9 + gijd(:,6) ** 2 ) ) 556 557 tauT = one/sqrt(tauT) 558c 559 if(idcsclr(1) .ne. 0) then 560 if ((idcsclr(2).eq.1 .and. isclr.eq.1) .or. 561 & (idcsclr(2).eq.2 .and. isclr.eq.2)) then ! scalar with dc 562c 563c determinant of gijd 564c 565 fact = one/(gijd(:,1) * gijd(:,2) * gijd(:,3) 566 & - gijd(:,2) * gijd(:,6) * gijd(:,6) 567 & - gijd(:,1) * gijd(:,5) * gijd(:,5) 568 & - gijd(:,3) * gijd(:,4) * gijd(:,4) 569 & + gijd(:,6) * gijd(:,4) * gijd(:,5) * two) 570c 571c ... note between compressible and incompressible 5 and 6 of giju 572c are switched 573c 574 giju(:,1) = fact * (gijd(:,2)*gijd(:,3) 575 & - gijd(:,5)**2) 576 giju(:,2) = fact * (gijd(:,1)*gijd(:,3) 577 & - gijd(:,6)**2) 578 giju(:,3) = fact * (gijd(:,1)*gijd(:,2) 579 & - gijd(:,4)**2) 580 giju(:,4) = fact * (gijd(:,5)*gijd(:,6) 581 & - gijd(:,4)*gijd(:,3) ) 582 giju(:,5) = fact * (gijd(:,4)*gijd(:,6) 583 & - gijd(:,1)*gijd(:,5) ) 584 giju(:,6) = fact * (gijd(:,4)*gijd(:,5) 585 & - gijd(:,6)*gijd(:,2) ) 586 587c 588 endif 589 endif ! end of idcsclr.ne.0 590c 591c.... return 592c 593 return 594 end 595 596