1c 2c file mud3pn.f 3c 4c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5c . . 6c . copyright (c) 1999 by UCAR . 7c . . 8c . UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH . 9c . . 10c . all rights reserved . 11c . . 12c . . 13c . MUDPACK version 5.0 . 14c . . 15c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16c 17c ... author and specialist 18c 19c John C. Adams (National Center for Atmospheric Research) 20c email: johnad@ucar.edu, phone: 303-497-1213 21 22c ... For MUDPACK 5.0 information, visit the website: 23c (http://www.scd.ucar.edu/css/software/mudpack) 24c 25c ... purpose 26c 27c mud3pn.f contains subroutines for planar relaxation in the x or 28c y or z direction. This file must be loaded with any of the real 29c 3-d mudpack solvers except mud3sp. 30c 31 subroutine planxy(wk) 32c 33c planar relaxation on (x,y) planes in z direction 34c 35 implicit none 36 real wk(*) 37 integer iparm(16),mgo(4) 38 real fparm(6) 39 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 40 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 41 +kcycle,iprer,ipost,intpol 42 real xa,xb,yc,yd,ze,zf,tolmax,relmax 43 integer kpbgn,kcbgn,ktxbgn,ktybgn, 44 +ktzbgn,nxk,nyk,nzk,ngrid,klevel,kcur,kps 45 integer nx,ny,nz,k,kb,kset,kl,ixy,icxy,ip,ic,isx,jsy,kst,kfn 46 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 47 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 48 +kcycle,iprer,ipost,intpol 49 common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax 50 common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 51 +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps 52c 53c set integer parameters for xy plane 54c 55 iparm(1) = 1 56 iparm(2) = nxa 57 iparm(3) = nxb 58 iparm(4) = nyc 59 iparm(5) = nyd 60 iparm(6) = ixp 61 iparm(7) = jyq 62 iparm(8) = max0(klevel+iex-ngrid,1) 63 iparm(9) = max0(klevel+jey-ngrid,1) 64 iparm(10) = nxk(klevel) 65 iparm(11) = nyk(klevel) 66 iparm(12) = 1 67 iparm(13) = 1 68c 69c set relaxation method for 2-d 70c 71 iparm(14) = meth2 72 73 fparm(1) = xa 74 fparm(2) = xb 75 fparm(3) = yc 76 fparm(4) = yd 77 fparm(5) = 0.0 78c 79c set 2-d multigrid options 80c 81 mgo(1) = kcycle 82 mgo(2) = iprer 83 mgo(3) = ipost 84 mgo(4) = intpol 85 isx = 0 86 if (meth2.eq.1.or.meth2.eq.3) then 87 isx = 3 88 if (nxa.eq.0) isx = 5 89 end if 90 jsy = 0 91 if (meth2.eq.2.or.meth2.eq.3) then 92 jsy = 3 93 if (nyc.eq.0) jsy = 5 94 end if 95 nz = nzk(klevel) 96 kst = 2 97 if (nze.ne.1) kst = 1 98 kfn = nz-1 99 if (nzf.ne.1) kfn = nz 100 do k=kst,kfn 101 kset = k 102 ixy = kps 103 ip = kpbgn(klevel) 104 ic = kcbgn(klevel) 105 do kb=1,klevel 106 kl = klevel-kb+1 107 nx = nxk(kl) 108 ny = nyk(kl) 109c 110c set 2-d coefficient pointer 111c 112 icxy = ixy + (nx+2)*(ny+2) 113c 114c transfer coefficients from 3-d to 2-d for current kset 115c 116 call trscxy(kset,nx,ny,nz,wk(ic),wk(icxy),wk(ip)) 117 ixy = ixy+(6+isx+jsy)*nx*ny + (nx+2)*(ny+2) 118 ip = ip + (nx+2)*(ny+2)*(nz+2) 119 ic = ic + 8*nx*ny*nz 120 end do 121c 122c set 2-d solution for current kset for passage to mup2 123c 124 nx = nxk(klevel) 125 ny = nyk(klevel) 126 ip = kpbgn(klevel) 127 ixy = 0 128 call setpxy(kset,nx,ny,nz,wk(ip),wk(kps),ixy) 129c 130c solve on current z plane with full 2-d multigrid cycling 131c 132 call mup2(iparm,fparm,wk(kps),mgo) 133c 134c reset approx from mup2 in wk(ip) 135c 136 ixy = 1 137 call setpxy(kset,nx,ny,nz,wk(ip),wk(kps),ixy) 138 end do 139c 140c set periodic virtual boundaries if necessary 141c 142 if (nxa*nyc*nze.eq.0) then 143 nx = nxk(klevel) 144 ny = nyk(klevel) 145 ip = kpbgn(klevel) 146 call per3vb(nx,ny,nz,wk(ip),nxa,nyc,nze) 147 end if 148 return 149 end 150 151 subroutine trscxy(kset,nx,ny,nz,cof,cofxy,phi) 152c 153c transfer coefficients from 3-d to 2-d 154c to allow 2-d relaxation 155c 156 implicit none 157 integer kset,nx,ny,nz,i,j,k 158 real cof(nx,ny,nz,8),cofxy(nx,ny,6) 159 real phi(0:nx+1,0:ny+1,0:nz+1) 160 k = kset 161 do j=1,ny 162 do i=1,nx 163 cofxy(i,j,1) = cof(i,j,k,1) 164 cofxy(i,j,2) = cof(i,j,k,2) 165 cofxy(i,j,3) = cof(i,j,k,3) 166 cofxy(i,j,4) = cof(i,j,k,4) 167 cofxy(i,j,5) = cof(i,j,k,7) 168 cofxy(i,j,6) = cof(i,j,k,8)-(cof(i,j,k,5)*phi(i,j,k-1)+ 169 + cof(i,j,k,6)*phi(i,j,k+1)) 170 end do 171 end do 172 return 173 end 174 175 subroutine setpxy(kset,nx,ny,nz,phi,phxy,ixy) 176 implicit none 177 integer kset,nx,ny,nz,i,j,ixy 178 real phi(0:nx+1,0:ny+1,0:nz+1),phxy(0:nx+1,0:ny+1) 179 if (ixy.eq.0) then 180 do j=0,ny+1 181 do i=0,nx+1 182 phxy(i,j) = phi(i,j,kset) 183 end do 184 end do 185 else 186 do j=0,ny+1 187 do i=0,nx+1 188 phi(i,j,kset) = phxy(i,j) 189 end do 190 end do 191 end if 192 return 193 end 194 195 subroutine planxz(wk) 196c 197c planar relaxation on (x,z) planes in y direction 198c 199 implicit none 200 real wk(*) 201 integer iparm(16),mgo(4) 202 real fparm(6) 203 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 204 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 205 +kcycle,iprer,ipost,intpol 206 real xa,xb,yc,yd,ze,zf,tolmax,relmax 207 integer kpbgn,kcbgn,ktxbgn,ktybgn, 208 +ktzbgn,nxk,nyk,nzk,ngrid,klevel,kcur,kps 209 integer nx,ny,nz,j,kb,jset,kl,ixz,icxz,ip,ic,isx,ksz,jst,jfn 210 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 211 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 212 +kcycle,iprer,ipost,intpol 213 common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax 214 common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 215 +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps 216c 217c set integer parameters for xz plane 218c 219 iparm(1) = 1 220 iparm(2) = nxa 221 iparm(3) = nxb 222 iparm(4) = nze 223 iparm(5) = nzf 224 iparm(6) = ixp 225 iparm(7) = kzr 226 iparm(8) = max0(klevel+iex-ngrid,1) 227 iparm(9) = max0(klevel+kez-ngrid,1) 228 iparm(10) = nxk(klevel) 229 iparm(11) = nzk(klevel) 230 iparm(12) = 1 231 iparm(13) = 1 232c 233c set relaxation method for 2-d 234c 235 iparm(14) = meth2 236 237 fparm(1) = xa 238 fparm(2) = xb 239 fparm(3) = ze 240 fparm(4) = zf 241 fparm(5) = 0.0 242c 243c set 2-d multigrid options 244c 245 mgo(1) = kcycle 246 mgo(2) = iprer 247 mgo(3) = ipost 248 mgo(4) = intpol 249 isx = 0 250 if (meth2.eq.1.or.meth2.eq.3) then 251 isx = 3 252 if (nxa.eq.0) isx = 5 253 end if 254 ksz = 0 255 if (meth2.eq.2.or.meth2.eq.3) then 256 ksz = 3 257 if (nze.eq.0) ksz = 5 258 end if 259 ny = nyk(klevel) 260 jst = 2 261 if (nyc.ne.1) jst = 1 262 jfn = ny-1 263 if (nyd.ne.1) jfn = ny 264 do j=jst,jfn 265 jset = j 266 ixz = kps 267 ip = kpbgn(klevel) 268 ic = kcbgn(klevel) 269 do kb=1,klevel 270 kl = klevel-kb+1 271 nx = nxk(kl) 272 nz = nzk(kl) 273c 274c set 2-d coefficient pointer 275c 276 icxz = ixz + (nx+2)*(nz+2) 277c 278c transfer coefficients from 3-d to 2-d for current jset 279c 280 call trscxz(jset,nx,ny,nz,wk(ic),wk(icxz),wk(ip)) 281 ixz = ixz+(6+isx+ksz)*nx*nz + (nx+2)*(nz+2) 282 ip = ip + (nx+2)*(ny+2)*(nz+2) 283 ic = ic + 8*nx*ny*nz 284 end do 285c 286c set 2-d solution for current jset for passage to mup2 287c 288 nx = nxk(klevel) 289 nz = nzk(klevel) 290 ip = kpbgn(klevel) 291 ixz = 0 292 call setpxz(jset,nx,ny,nz,wk(ip),wk(kps),ixz) 293c 294c solve on current y plane with full 2-d multigrid cycling 295c 296 call mup2(iparm,fparm,wk(kps),mgo) 297 ixz = 1 298 call setpxz(jset,nx,ny,nz,wk(ip),wk(kps),ixz) 299 end do 300c 301c set periodic virtual boundaries if necessary 302c 303 if (nxa*nyc*nze.eq.0) then 304 nx = nxk(klevel) 305 nz = nzk(klevel) 306 ip = kpbgn(klevel) 307 call per3vb(nx,ny,nz,wk(ip),nxa,nyc,nze) 308 end if 309 return 310 end 311 312 subroutine trscxz(jset,nx,ny,nz,cof,cofxz,phi) 313c 314c transfer coefficients from 3-d to 2-d 315c to allow 2-d relaxation 316c 317 implicit none 318 integer jset,nx,ny,nz,i,j,k 319 real cof(nx,ny,nz,8),cofxz(nx,nz,6) 320 real phi(0:nx+1,0:ny+1,0:nz+1) 321 j = jset 322 do k=1,nz 323 do i=1,nx 324 cofxz(i,k,1) = cof(i,j,k,1) 325 cofxz(i,k,2) = cof(i,j,k,2) 326 cofxz(i,k,3) = cof(i,j,k,5) 327 cofxz(i,k,4) = cof(i,j,k,6) 328 cofxz(i,k,5) = cof(i,j,k,7) 329 cofxz(i,k,6) = cof(i,j,k,8)-(cof(i,j,k,3)*phi(i,j-1,k)+ 330 + cof(i,j,k,4)*phi(i,j+1,k)) 331 end do 332 end do 333 return 334 end 335 336 subroutine setpxz(jset,nx,ny,nz,phi,phxz,ixz) 337 implicit none 338 integer jset,nx,ny,nz,i,k,ixz 339 real phi(0:nx+1,0:ny+1,0:nz+1),phxz(0:nx+1,0:nz+1) 340 if (ixz.eq.0) then 341 do k=0,nz+1 342 do i=0,nx+1 343 phxz(i,k) = phi(i,jset,k) 344 end do 345 end do 346 else 347 do k=0,nz+1 348 do i=0,nx+1 349 phi(i,jset,k) = phxz(i,k) 350 end do 351 end do 352 end if 353 return 354 end 355 356 subroutine planyz(wk) 357c 358c planar relaxation on (y,z) planes in x direction 359c 360 implicit none 361 real wk(*) 362 integer iparm(16),mgo(4) 363 real fparm(6) 364 integer intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 365 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 366 +kcycle,iprer,ipost,intpol 367 real xa,xb,yc,yd,ze,zf,tolmax,relmax 368 integer kpbgn,kcbgn,ktxbgn,ktybgn, 369 +ktzbgn,nxk,nyk,nzk,ngrid,klevel,kcur,kps 370 integer nx,ny,nz,i,kb,iset,kl,iyz,icyz,ip,ic,jsy,ksz,ist,ifn 371 common/imud3/intl,nxa,nxb,nyc,nyd,nze,nzf,ixp,jyq,kzr,iex,jey,kez, 372 +nfx,nfy,nfz,iguess,maxcy,method,meth2,nwork,lwork,itero, 373 +kcycle,iprer,ipost,intpol 374 common/fmud3/xa,xb,yc,yd,ze,zf,tolmax,relmax 375 common/mud3c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 376 +ktzbgn(50),nxk(50),nyk(50),nzk(50),ngrid,klevel,kcur,kps 377c 378c set integer parameters for yz plane 379c 380 iparm(1) = 1 381 iparm(2) = nyc 382 iparm(3) = nyd 383 iparm(4) = nze 384 iparm(5) = nzf 385 iparm(6) = jyq 386 iparm(7) = kzr 387 iparm(8) = max0(klevel+jey-ngrid,1) 388 iparm(9) = max0(klevel+kez-ngrid,1) 389 iparm(10) = nyk(klevel) 390 iparm(11) = nzk(klevel) 391 iparm(12) = 1 392 iparm(13) = 1 393c 394c set relaxation method for 2-d 395c 396 iparm(14) = meth2 397 398 fparm(1) = yc 399 fparm(2) = yd 400 fparm(3) = ze 401 fparm(4) = zf 402 fparm(5) = 0.0 403c 404c set 2-d multigrid options 405c 406 mgo(1) = kcycle 407 mgo(2) = iprer 408 mgo(3) = ipost 409 mgo(4) = intpol 410 jsy = 0 411 if (meth2.eq.1.or.meth2.eq.3) then 412 jsy = 3 413 if (nyc.eq.0) jsy = 5 414 end if 415 ksz = 0 416 if (meth2.eq.2.or.meth2.eq.3) then 417 ksz = 3 418 if (nze.eq.0) ksz = 5 419 end if 420 nx = nxk(klevel) 421 ist = 2 422 if (nxa.ne.1) ist = 1 423 ifn = nx-1 424 if (nxb.ne.1) ifn = nx 425 do i=ist,ifn 426 iset = i 427 iyz = kps 428 ip = kpbgn(klevel) 429 ic = kcbgn(klevel) 430 do kb=1,klevel 431 kl = klevel-kb+1 432 ny = nyk(kl) 433 nz = nzk(kl) 434c 435c set 2-d coefficient pointer 436c 437 icyz = iyz + (ny+2)*(nz+2) 438c 439c transfer coefficients from 3-d to 2-d for current iset 440c 441 call trscyz(iset,nx,ny,nz,wk(ic),wk(icyz),wk(ip)) 442 iyz = iyz+(6+jsy+ksz)*ny*nz + (ny+2)*(nz+2) 443 ip = ip + (nx+2)*(ny+2)*(nz+2) 444 ic = ic + 8*nx*ny*nz 445 end do 446c 447c set 2-d solution for current iset for passage to mup2 448c 449 ny = nyk(klevel) 450 nz = nzk(klevel) 451 ip = kpbgn(klevel) 452 iyz = 0 453 call setpyz(iset,nx,ny,nz,wk(ip),wk(kps),iyz) 454c 455c solve on current y plane with full 2-d multigrid cycling 456c 457 call mup2(iparm,fparm,wk(kps),mgo) 458 iyz = 1 459 call setpyz(iset,nx,ny,nz,wk(ip),wk(kps),iyz) 460 end do 461c 462c set periodic virtual boundaries if necessary 463c 464 if (nxa*nyc*nze.eq.0) then 465 ny = nyk(klevel) 466 nz = nzk(klevel) 467 ip = kpbgn(klevel) 468 call per3vb(nx,ny,nz,wk(ip),nxa,nyc,nze) 469 end if 470 return 471 end 472 473 subroutine trscyz(iset,nx,ny,nz,cof,cofyz,phi) 474c 475c transfer coefficients from 3-d to 2-d 476c to allow 2-d relaxation 477c 478 implicit none 479 integer iset,nx,ny,nz,i,j,k 480 real cof(nx,ny,nz,8),cofyz(ny,nz,6) 481 real phi(0:nx+1,0:ny+1,0:nz+1) 482 i = iset 483 do k=1,nz 484 do j=1,ny 485 cofyz(j,k,1) = cof(i,j,k,3) 486 cofyz(j,k,2) = cof(i,j,k,4) 487 cofyz(j,k,3) = cof(i,j,k,5) 488 cofyz(j,k,4) = cof(i,j,k,6) 489 cofyz(j,k,5) = cof(i,j,k,7) 490 cofyz(j,k,6) = cof(i,j,k,8)-(cof(i,j,k,1)*phi(i-1,j,k)+ 491 + cof(i,j,k,2)*phi(i+1,j,k)) 492 end do 493 end do 494 return 495 end 496 497 subroutine setpyz(iset,nx,ny,nz,phi,phyz,iyz) 498 implicit none 499 integer iset,nx,ny,nz,iyz,j,k 500 real phi(0:nx+1,0:ny+1,0:nz+1),phyz(0:ny+1,0:nz+1) 501 if (iyz.eq.0) then 502 do k=0,nz+1 503 do j=0,ny+1 504 phyz(j,k) = phi(iset,j,k) 505 end do 506 end do 507 else 508 do k=0,nz+1 509 do j=0,ny+1 510 phi(iset,j,k) = phyz(j,k) 511 end do 512 end do 513 end if 514 return 515 end 516 517 subroutine mup2(iparm,fparm,work,mgopt) 518c 519c modification of mud2 for planar with mud3 520c whenever mup2 is called coefficients from discretization 521c have already been set but matrices for line (if flagged 522c have not been set) 523c 524 implicit none 525 integer iparm,mgopt 526 integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 527 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 528 + kcycle,iprer,ipost,intpol,kps 529 real fparm,xa,xb,yc,yd,tolmax,relmax 530 integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy 531 integer iw,k,kb,nx,ny,ic,itx,ity 532 dimension iparm(16),fparm(6),mgopt(4) 533 real work(*) 534 common/imup2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 535 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 536 + kcycle,iprer,ipost,intpol,kps 537 common/fmup2/xa,xb,yc,yd,tolmax,relmax 538 common/mup2c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 539 +nxk(50),nyk(50),isx,jsy 540 intl = 1 541 nxa = iparm(2) 542 nxb = iparm(3) 543 nyc = iparm(4) 544 nyd = iparm(5) 545 ixp = iparm(6) 546 jyq = iparm(7) 547 iex = iparm(8) 548 jey = iparm(9) 549 ngrid = max0(iex,jey) 550 nfx = iparm(10) 551 nfy = iparm(11) 552 iguess = iparm(12) 553 maxcy = iparm(13) 554 method = iparm(14) 555 nwork = iparm(15) 556 kcycle = mgopt(1) 557 iprer = mgopt(2) 558 ipost = mgopt(3) 559 intpol = mgopt(4) 560 xa = fparm(1) 561 xb = fparm(2) 562 yc = fparm(3) 563 yd = fparm(4) 564 tolmax = fparm(5) 565 isx = 0 566 jsy = 0 567 if ((method-1)*(method-3).eq.0) then 568 isx = 3 569 if (nxa.eq.0) isx = 5 570 end if 571 if ((method-2)*(method-3).eq.0) then 572 jsy = 3 573 if (nyc.eq.0) jsy = 5 574 end if 575 kps = 1 576 do k=1,ngrid 577c set subgrid sizes 578 nxk(k) = ixp*2**(max0(k+iex-ngrid,1)-1)+1 579 nyk(k) = jyq*2**(max0(k+jey-ngrid,1)-1)+1 580 nx = nxk(k) 581 ny = nyk(k) 582 kps = kps+(nx+2)*(ny+2)+nx*ny*(6+isx+jsy) 583 end do 584c 585c set work space pointers and discretize pde at each grid level 586c 587 iw = 1 588 do kb=1,ngrid 589 k = ngrid-kb+1 590 nx = nxk(k) 591 ny = nyk(k) 592 kpbgn(k) = iw 593 kcbgn(k) = kpbgn(k)+(nx+2)*(ny+2) 594 ktxbgn(k) = kcbgn(k)+6*nx*ny 595 ktybgn(k) = ktxbgn(k)+isx*nx*ny 596 iw = ktybgn(k)+jsy*nx*ny 597 ic = kcbgn(k) 598 itx = ktxbgn(k) 599 ity = ktybgn(k) 600 klevel = k 601 call dismp2(nx,ny,work(ic),work(itx),work(ity),work(kps)) 602 end do 603 call mup21(work) 604 return 605 end 606 607 subroutine mup21(wk) 608 implicit none 609 real wk(*) 610 integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy 611 integer iter 612 integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 613 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 614 + kcycle,iprer,ipost,intpol,kps 615 common/imup2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 616 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 617 + kcycle,iprer,ipost,intpol,kps 618 common/mup2c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 619 +nxk(50),nyk(50),isx,jsy 620c 621c cycle from ngrid level 622c 623 kcur = ngrid 624 do iter=1,maxcy 625 itero = iter 626 call kcymp2(wk) 627 end do 628 return 629 end 630 631 subroutine kcymp2(wk) 632c 633c execute multigrid k cycle from kcur grid level 634c kcycle=1 for v cycles, kcycle=2 for w cycles 635c 636 implicit none 637 real wk(*) 638 integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 639 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 640 + kcycle,iprer,ipost,intpol,kps 641 integer nx,ny,ip,ic,ipc,irc,itx,ity,ncx,ncy,l,nrel 642 real xa,xb,yc,yd,tolmax,relmax 643 integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy 644 common/imup2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 645 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 646 + kcycle,iprer,ipost,intpol,kps 647 common/fmup2/xa,xb,yc,yd,tolmax,relmax 648 common/mup2c/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), 649 +nxk(50),nyk(50),isx,jsy 650 integer kount(50) 651 klevel = kcur 652 nx = nxk(klevel) 653 ny = nyk(klevel) 654 ip = kpbgn(klevel) 655 ic = kcbgn(klevel) 656 itx = ktxbgn(klevel) 657 ity = ktybgn(klevel) 658c 659c prerelax at current finest grid level 660c 661 do l=1,iprer 662 call relmp2(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps)) 663 end do 664 if (kcur .eq. 1) go to 5 665c 666c restrict residual to kcur-1 level 667c 668 ipc = kpbgn(klevel-1) 669 ncx = nxk(klevel-1) 670 ncy = nyk(klevel-1) 671 irc = kcbgn(klevel-1)+5*ncx*ncy 672 call resmp2(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic),wk(kps)) 673c 674c set counter for grid levels to zero 675c 676 do l = 1,kcur 677 kount(l) = 0 678 end do 679c 680c set new grid level and continue k-cycling 681c 682 klevel = kcur-1 683 nrel = iprer 684c 685c kcycle control point 686c 687 10 continue 688c 689c post relax when kcur revisited 690c 691 if (klevel .eq. kcur) go to 5 692c 693c count hit at current level 694c 695 kount(klevel) = kount(klevel)+1 696c 697c relax at current level 698c 699 nx = nxk(klevel) 700 ny = nyk(klevel) 701 ip = kpbgn(klevel) 702 ic = kcbgn(klevel) 703 itx = ktxbgn(klevel) 704 ity = ktybgn(klevel) 705 do l=1,nrel 706 call relmp2(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps)) 707 end do 708 if (kount(klevel) .eq. kcycle+1) then 709c 710c kcycle complete at klevel 711c 712 ipc = ip 713 ip = kpbgn(klevel+1) 714 ncx = nxk(klevel) 715 ncy = nyk(klevel) 716 nx = nxk(klevel+1) 717 ny = nyk(klevel+1) 718c 719c inject correction to finer grid 720c 721 call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd, 722 + intpol,wk(kps)) 723c 724c reset counter to zero 725c 726 kount(klevel) = 0 727c 728c ascend to next higher level and set to postrelax there 729c 730 klevel = klevel+1 731 nrel = ipost 732 go to 10 733 else 734 if (klevel .gt. 1) then 735c 736c kcycle not complete so descend unless at coarsest grid 737c 738 ipc = kpbgn(klevel-1) 739 ncx = nxk(klevel-1) 740 ncy = nyk(klevel-1) 741 irc = kcbgn(klevel-1)+5*ncx*ncy 742 call resmp2(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic), 743 + wk(kps)) 744c 745c prerelax at next coarser level 746c 747 klevel = klevel-1 748 nrel = iprer 749 go to 10 750 else 751c 752c postrelax at coarsest level 753c 754 do l=1,ipost 755 call relmp2(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps)) 756 end do 757 ipc = ip 758 ip = kpbgn(2) 759 ncx = nxk(1) 760 ncy = nyk(1) 761 nx = nxk(2) 762 ny = nyk(2) 763c 764c inject correction to level 2 765c 766 call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd, 767 + intpol,wk(kps)) 768c 769c set to postrelax at level 2 770c 771 nrel = ipost 772 klevel = 2 773 go to 10 774 end if 775 end if 776 5 continue 777c 778c post relax at current finest grid level 779c 780 nx = nxk(kcur) 781 ny = nyk(kcur) 782 ip = kpbgn(kcur) 783 ic = kcbgn(kcur) 784 itx = ktxbgn(kcur) 785 ity = ktybgn(kcur) 786 do l=1,ipost 787 call relmp2(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps)) 788 end do 789 return 790 end 791 792 subroutine dismp2(nx,ny,cof,tx,ty,sum) 793c 794c set tridiagonal matrices for line relaxation if necessary 795c 796 implicit none 797 integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 798 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 799 + kcycle,iprer,ipost,intpol,kps 800 integer nx,ny,i,j,im1,jm1 801 real tx(nx,ny,*),ty(ny,nx,*),cof(nx,ny,6),sum(*) 802 common/imup2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 803 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 804 + kcycle,iprer,ipost,intpol,kps 805 if (method.eq.1.or.method.eq.3) then 806 if (nxa.ne.0) then 807c 808c nonperiodic x line relaxation 809c 810 do i=1,nx 811 im1 = max0(i-1,1) 812 do j=1,ny 813 tx(im1,j,1) = cof(i,j,1) 814 tx(i,j,2) = cof(i,j,5) 815 tx(i,j,3) = cof(i,j,2) 816 end do 817 end do 818 call factri(ny,nx,tx(1,1,1),tx(1,1,2),tx(1,1,3)) 819 else 820c 821c periodic x line relaxation 822c 823 if (nx .gt. 3) then 824c 825c set and factor iff nx > 3 826c 827 do i=1,nx-1 828 do j=1,ny 829 tx(i,j,1) = cof(i,j,1) 830 tx(i,j,2) = cof(i,j,5) 831 tx(i,j,3) = cof(i,j,2) 832 end do 833 end do 834 call factrp(ny,nx,tx,tx(1,1,2),tx(1,1,3),tx(1,1,4), 835 + tx(1,1,5),sum) 836 end if 837 end if 838 end if 839 840 if (method.eq.2.or.method.eq.3) then 841 if (nyc.ne.0) then 842c 843c nonperiodic y line relaxation 844c 845 do j=1,ny 846 jm1 = max0(j-1,1) 847 do i=1,nx 848 ty(jm1,i,1) = cof(i,j,3) 849 ty(j,i,2) = cof(i,j,5) 850 ty(j,i,3) = cof(i,j,4) 851 end do 852 end do 853 call factri(nx,ny,ty(1,1,1),ty(1,1,2),ty(1,1,3)) 854 else 855c 856c periodic y line relaxation 857c 858 if (ny .gt. 3) then 859c 860c set and factor iff ny > 3 861c 862 do j=1,ny-1 863 do i=1,nx 864 ty(j,i,1) = cof(i,j,3) 865 ty(j,i,2) = cof(i,j,5) 866 ty(j,i,3) = cof(i,j,4) 867 end do 868 end do 869 call factrp(nx,ny,ty,ty(1,1,2),ty(1,1,3),ty(1,1,4), 870 + ty(1,1,5),sum) 871 end if 872 end if 873 end if 874 return 875 end 876 877 subroutine resmp2(nx,ny,phi,ncx,ncy,phic,rhsc,cof,resf) 878c 879c restrict residual from fine to coarse mesh using fully weighted 880c residual restriction 881c 882 implicit none 883 integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 884 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 885 + kcycle,iprer,ipost,intpol,kps 886 integer nx,ny,ncx,ncy,i,j,ic,jc 887 common/imup2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 888 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 889 + kcycle,iprer,ipost,intpol,kps 890 real rhsc(ncx,ncy),resf(nx,ny) 891 real phi(0:nx+1,0:ny+1),phic(0:ncx+1,0:ncy+1) 892 real cof(nx,ny,6) 893c 894c set phic zero 895c 896 do jc=0,ncy+1 897 do ic=0,ncx+1 898 phic(ic,jc) = 0.0 899 end do 900 end do 901c 902c compute residual on fine mesh in resf 903c 904!$OMP PARALLEL DO SHARED(resf,cof,phi,nx,ny) PRIVATE(i,j) 905 do j=1,ny 906 do i=1,nx 907 resf(i,j) = cof(i,j,6)-( 908 + cof(i,j,1)*phi(i-1,j)+ 909 + cof(i,j,2)*phi(i+1,j)+ 910 + cof(i,j,3)*phi(i,j-1)+ 911 + cof(i,j,4)*phi(i,j+1)+ 912 + cof(i,j,5)*phi(i,j)) 913 end do 914 end do 915!$OMP END PARALLEL DO 916c 917c restrict resf to coarse mesh in rhsc 918c 919 call res2(nx,ny,resf,ncx,ncy,rhsc,nxa,nxb,nyc,nyd) 920 return 921 end 922 923 subroutine relmp2(nx,ny,phi,cof,tx,ty,sum) 924c 925c relaxation for mud2 926c 927 implicit none 928 integer nx,ny 929 real phi(*),cof(*),tx(*),ty(*),sum(*) 930 integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 931 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 932 + kcycle,iprer,ipost,intpol,kps 933 common/imup2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 934 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 935 + kcycle,iprer,ipost,intpol,kps 936 if (method.eq.0) then ! point relaxation 937 call relmp2p(nx,ny,phi,cof) 938 else if (method.eq.1) then ! line x relaxation 939 call slxmp2(nx,ny,phi,cof,tx,sum) 940 else if (method.eq.2) then ! line y relaxation 941 call slymp2(nx,ny,phi,cof,ty,sum) 942 else if (method.eq.3) then ! line x&y relaxation 943 call slxmp2(nx,ny,phi,cof,tx,sum) 944 call slymp2(nx,ny,phi,cof,ty,sum) 945 end if 946 return 947 end 948 949 subroutine relmp2p(nx,ny,phi,cof) 950c 951c Gauss-Seidel red/black point relaxation 952c 953 implicit none 954 integer nx,ny,i,j 955 integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 956 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 957 + kcycle,iprer,ipost,intpol,kps 958 common/imup2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 959 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 960 + kcycle,iprer,ipost,intpol,kps 961 real phi(0:nx+1,0:ny+1),cof(nx,ny,6) 962c 963c periodic adjustment bypass block 964c 965 if (nxa*nyc.ne.0) then 966c 967c relax on red grid points 968c 969!$OMP PARALLEL DO SHARED(cof,phi,nx,ny) PRIVATE(i,j) 970 do i=1,nx,2 971 do j=1,ny,2 972 phi(i,j) = (cof(i,j,6) - 973 + (cof(i,j,1)*phi(i-1,j)+cof(i,j,2)*phi(i+1,j) + 974 + cof(i,j,3)*phi(i,j-1)+cof(i,j,4)*phi(i,j+1)))/ 975 + cof(i,j,5) 976 end do 977 end do 978!$OMP END PARALLEL DO 979c 980!$OMP PARALLEL DO SHARED(cof,phi,nx,ny) PRIVATE(i,j) 981 do i=2,nx,2 982 do j=2,ny,2 983 phi(i,j) = (cof(i,j,6) - 984 + (cof(i,j,1)*phi(i-1,j)+cof(i,j,2)*phi(i+1,j) + 985 + cof(i,j,3)*phi(i,j-1)+cof(i,j,4)*phi(i,j+1)))/ 986 + cof(i,j,5) 987 end do 988 end do 989!$OMP END PARALLEL DO 990c 991c relax on black grid points 992c 993c 994!$OMP PARALLEL DO SHARED(cof,phi,nx,ny) PRIVATE(i,j) 995 do i=1,nx,2 996 do j=2,ny,2 997 phi(i,j) = (cof(i,j,6) - 998 + (cof(i,j,1)*phi(i-1,j)+cof(i,j,2)*phi(i+1,j) + 999 + cof(i,j,3)*phi(i,j-1)+cof(i,j,4)*phi(i,j+1)))/ 1000 + cof(i,j,5) 1001 end do 1002 end do 1003!$OMP END PARALLEL DO 1004c 1005!$OMP PARALLEL DO SHARED(cof,phi,nx,ny) PRIVATE(i,j) 1006 do i=2,nx,2 1007 do j=1,ny,2 1008 phi(i,j) = (cof(i,j,6) - 1009 + (cof(i,j,1)*phi(i-1,j)+cof(i,j,2)*phi(i+1,j) + 1010 + cof(i,j,3)*phi(i,j-1)+cof(i,j,4)*phi(i,j+1)))/ 1011 + cof(i,j,5) 1012 end do 1013 end do 1014!$OMP END PARALLEL DO 1015 return 1016 end if 1017c 1018c set periodic virtual boundaries 1019c 1020 if (nxa.eq.0) then 1021 do j=1,ny 1022 phi(0,j) = phi(nx-1,j) 1023 phi(nx+1,j) = phi(2,j) 1024 end do 1025 end if 1026 if (nyc.eq.0) then 1027 do i=1,nx 1028 phi(i,0) = phi(i,ny-1) 1029 phi(i,ny+1) = phi(i,2) 1030 end do 1031 end if 1032c 1033c relax on red grid points 1034c 1035!$OMP PARALLEL DO SHARED(cof,phi,nx,ny) PRIVATE(i,j) 1036 do i=1,nx,2 1037 do j=1,ny,2 1038 phi(i,j) = (cof(i,j,6) - 1039 + (cof(i,j,1)*phi(i-1,j)+cof(i,j,2)*phi(i+1,j) + 1040 + cof(i,j,3)*phi(i,j-1)+cof(i,j,4)*phi(i,j+1)))/ 1041 + cof(i,j,5) 1042 end do 1043 end do 1044!$OMP END PARALLEL DO 1045!$OMP PARALLEL DO SHARED(cof,phi,nx,ny) PRIVATE(i,j) 1046 do i=2,nx,2 1047 do j=2,ny,2 1048 phi(i,j) = (cof(i,j,6) - 1049 + (cof(i,j,1)*phi(i-1,j)+cof(i,j,2)*phi(i+1,j) + 1050 + cof(i,j,3)*phi(i,j-1)+cof(i,j,4)*phi(i,j+1)))/ 1051 + cof(i,j,5) 1052 end do 1053 end do 1054!$OMP END PARALLEL DO 1055c 1056c ensure periodic virtual boundary red points are set 1057c 1058 if (nxa.eq.0) then 1059 do j=1,ny 1060 phi(0,j) = phi(nx-1,j) 1061 phi(nx+1,j) = phi(2,j) 1062 end do 1063 end if 1064 if (nyc.eq.0) then 1065 do i=1,nx 1066 phi(i,0) = phi(i,ny-1) 1067 phi(i,ny+1) = phi(i,2) 1068 end do 1069 end if 1070c 1071c relax on black grid points 1072c 1073!$OMP PARALLEL DO SHARED(cof,phi,nx,ny) PRIVATE(i,j) 1074 do i=1,nx,2 1075 do j=2,ny,2 1076 phi(i,j) = (cof(i,j,6) - 1077 + (cof(i,j,1)*phi(i-1,j)+cof(i,j,2)*phi(i+1,j) + 1078 + cof(i,j,3)*phi(i,j-1)+cof(i,j,4)*phi(i,j+1)))/ 1079 + cof(i,j,5) 1080 end do 1081 end do 1082!$OMP END PARALLEL DO 1083!$OMP PARALLEL DO SHARED(cof,phi,nx,ny) PRIVATE(i,j) 1084 do i=2,nx,2 1085 do j=1,ny,2 1086 phi(i,j) = (cof(i,j,6) - 1087 + (cof(i,j,1)*phi(i-1,j)+cof(i,j,2)*phi(i+1,j) + 1088 + cof(i,j,3)*phi(i,j-1)+cof(i,j,4)*phi(i,j+1)))/ 1089 + cof(i,j,5) 1090 end do 1091 end do 1092!$OMP END PARALLEL DO 1093c 1094c final set of periodic virtual boundaries 1095c 1096 if (nxa.eq.0) then 1097 do j=1,ny 1098 phi(0,j) = phi(nx-1,j) 1099 phi(nx+1,j) = phi(2,j) 1100 end do 1101 end if 1102 if (nyc.eq.0) then 1103 do i=1,nx 1104 phi(i,0) = phi(i,ny-1) 1105 phi(i,ny+1) = phi(i,2) 1106 end do 1107 end if 1108 return 1109 end 1110 1111 subroutine slxmp2(nx,ny,phi,cof,tx,sum) 1112c 1113c line relaxation in the x direction (periodic or nonperiodic) 1114c 1115 implicit none 1116 integer nx,ny,i,ib,j 1117 integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 1118 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 1119 + kcycle,iprer,ipost,intpol,kps 1120 common/imup2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 1121 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 1122 + kcycle,iprer,ipost,intpol,kps 1123 real phi(0:nx+1,0:ny+1),cof(nx,ny,6),tx(nx,ny,*),sum(ny) 1124c 1125c replace line x with point gauss-seidel if 1126c x direction is periodic and nx = 3 (coarsest) 1127c 1128 if (nxa .eq. 0 .and. nx .eq. 3) then 1129 call relmp2p(nx,ny,phi,cof) 1130 return 1131 end if 1132c 1133c set periodic y virtual boundary if necessary 1134c 1135 if (nyc.eq.0) then 1136 do i=1,nx 1137 phi(i,0) = phi(i,ny-1) 1138 phi(i,ny+1) = phi(i,2) 1139 end do 1140 end if 1141 1142 if (nxa.ne.0) then 1143c 1144c x direction not periodic 1145c 1146c sweep odd j lines 1147c 1148!$OMP PARALLEL DO SHARED(cof,phi,tx,nx,ny) PRIVATE(i,ib,j) 1149 do j=1,ny,2 1150 do i=1,nx 1151 phi(i,j)=cof(i,j,6)-cof(i,j,3)*phi(i,j-1)-cof(i,j,4)* 1152 + phi(i,j+1) 1153 end do 1154c 1155c forward sweep 1156c 1157 do i=2,nx 1158 phi(i,j) = phi(i,j)-tx(i-1,j,1)*phi(i-1,j) 1159 end do 1160c 1161c backward sweep 1162c 1163 phi(nx,j) = phi(nx,j)/tx(nx,j,2) 1164 do ib=2,nx 1165 i = nx-ib+1 1166 phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j))/tx(i,j,2) 1167 end do 1168 end do 1169!$OMP END PARALLEL DO 1170c 1171c sweep even j lines forward and back 1172c 1173!$OMP PARALLEL DO SHARED(cof,phi,tx,nx,ny) PRIVATE(i,j,ib) 1174 do j=2,ny,2 1175 do i=1,nx 1176 phi(i,j)=cof(i,j,6)-cof(i,j,3)*phi(i,j-1)-cof(i,j,4)* 1177 + phi(i,j+1) 1178 end do 1179 do i=2,nx 1180 phi(i,j) = phi(i,j)-tx(i-1,j,1)*phi(i-1,j) 1181 end do 1182 phi(nx,j) = phi(nx,j)/tx(nx,j,2) 1183 do ib=2,nx 1184 i = nx-ib+1 1185 phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j))/tx(i,j,2) 1186 end do 1187 end do 1188!$OMP END PARALLEL DO 1189 else 1190c 1191c x direction periodic 1192c 1193 do j=1,ny 1194 sum(j) = 0.0 1195 phi(0,j) = phi(nx-1,j) 1196 phi(nx+1,j) = phi(2,j) 1197 end do 1198c 1199c sweep odd lines forward and back 1200c 1201!$OMP PARALLEL DO SHARED(sum,cof,phi,tx,nx,ny) PRIVATE(i,j,ib) 1202 do j=1,ny,2 1203 do i=1,nx-1 1204 phi(i,j)=cof(i,j,6)-cof(i,j,3)*phi(i,j-1)-cof(i,j,4)* 1205 + phi(i,j+1) 1206 end do 1207c 1208c forward sweep 1209 do i=2,nx-2 1210 phi(i,j) = phi(i,j)-tx(i,j,1)*phi(i-1,j) 1211 end do 1212 do i=1,nx-2 1213 sum(j) = sum(j)+tx(i,j,5)*phi(i,j) 1214 end do 1215 phi(nx-1,j) = phi(nx-1,j)-sum(j) 1216c 1217c backward sweep 1218c 1219 phi(nx-1,j) = phi(nx-1,j)/tx(nx-1,j,2) 1220 phi(nx-2,j) = (phi(nx-2,j)-tx(nx-2,j,4)*phi(nx-1,j))/ 1221 + tx(nx-2,j,2) 1222 do ib=4,nx 1223 i = nx-ib+1 1224 phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j)-tx(i,j,4)* 1225 + phi(nx-1,j))/tx(i,j,2) 1226 end do 1227 end do 1228!$OMP END PARALLEL DO 1229c 1230c set periodic and virtual points for j odd 1231c 1232 do j=1,ny,2 1233 phi(nx,j) = phi(1,j) 1234 phi(0,j) = phi(nx-1,j) 1235 phi(nx+1,j) = phi(2,j) 1236 end do 1237c 1238c sweep even j lines 1239c 1240!$OMP PARALLEL DO SHARED(sum,cof,phi,tx,nx,ny) PRIVATE(i,j,ib) 1241 do j=2,ny,2 1242 do i=1,nx-1 1243 phi(i,j)=cof(i,j,6)-cof(i,j,3)*phi(i,j-1)-cof(i,j,4)* 1244 + phi(i,j+1) 1245 end do 1246c 1247c forward sweep 1248c 1249 do i=2,nx-2 1250 phi(i,j) = phi(i,j)-tx(i,j,1)*phi(i-1,j) 1251 end do 1252 do i=1,nx-2 1253 sum(j) = sum(j)+tx(i,j,5)*phi(i,j) 1254 end do 1255 phi(nx-1,j) = phi(nx-1,j)-sum(j) 1256c 1257c backward sweep 1258c 1259 phi(nx-1,j) = phi(nx-1,j)/tx(nx-1,j,2) 1260 phi(nx-2,j) = (phi(nx-2,j)-tx(nx-2,j,4)*phi(nx-1,j))/ 1261 + tx(nx-2,j,2) 1262 do ib=4,nx 1263 i = nx-ib+1 1264 phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j)-tx(i,j,4)* 1265 + phi(nx-1,j))/tx(i,j,2) 1266 end do 1267 end do 1268!$OMP END PARALLEL DO 1269c 1270c set periodic and virtual points for j even 1271c 1272 do j=2,ny,2 1273 phi(nx,j) = phi(1,j) 1274 phi(0,j) = phi(nx-1,j) 1275 phi(nx+1,j) = phi(2,j) 1276 end do 1277 end if 1278c 1279c set periodic y virtual boundaries if necessary 1280c 1281 if (nyc.eq.0) then 1282 do i=1,nx 1283 phi(i,0) = phi(i,ny-1) 1284 phi(i,ny+1) = phi(i,2) 1285 end do 1286 end if 1287 return 1288 end 1289 1290 subroutine slymp2(nx,ny,phi,cof,ty,sum) 1291 implicit none 1292 integer nx,ny,i,j,jb 1293 integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 1294 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 1295 + kcycle,iprer,ipost,intpol,kps 1296 common/imup2/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, 1297 + maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, 1298 + kcycle,iprer,ipost,intpol,kps 1299 real phi(0:nx+1,0:ny+1),cof(nx,ny,6),ty(ny,nx,*),sum(nx) 1300c 1301c replace line y with point gauss-seidel if 1302c y direction is periodic and ny = 3 1303c 1304 if (nyc .eq. 0 .and. ny .eq. 3) then 1305 call relmp2p(nx,ny,phi,cof) 1306 return 1307 end if 1308c 1309c set periodic and virtual x boundaries if necessary 1310c 1311 if (nxa.eq.0) then 1312 do j=1,ny 1313 phi(0,j) = phi(nx-1,j) 1314 phi(nx,j) = phi(1,j) 1315 phi(nx+1,j) = phi(2,j) 1316 end do 1317 end if 1318 1319 if (nyc.ne.0) then 1320c 1321c y direction not periodic 1322c 1323c 1324c sweep odd x lines 1325c 1326!$OMP PARALLEL DO SHARED(cof,phi,ty,nx,ny) PRIVATE(i,j,jb) 1327 do i=1,nx,2 1328 do j=1,ny 1329 phi(i,j)=cof(i,j,6)-cof(i,j,1)*phi(i-1,j)-cof(i,j,2)* 1330 + phi(i+1,j) 1331 end do 1332c 1333c forward sweep thru odd x lines 1334c 1335 do j=2,ny 1336 phi(i,j) = phi(i,j)-ty(j-1,i,1)*phi(i,j-1) 1337 end do 1338c 1339c backward sweep 1340c 1341 phi(i,ny) = phi(i,ny)/ty(ny,i,2) 1342 do jb=2,ny 1343 j = ny-jb+1 1344 phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1))/ty(j,i,2) 1345 end do 1346 end do 1347!$OMP END PARALLEL DO 1348c 1349c forward sweep even x lines 1350c 1351!$OMP PARALLEL DO SHARED(cof,phi,ty,nx,ny) PRIVATE(i,j,jb) 1352 do i=2,nx,2 1353 do j=1,ny 1354 phi(i,j)=cof(i,j,6)-cof(i,j,1)*phi(i-1,j)-cof(i,j,2)* 1355 + phi(i+1,j) 1356 end do 1357 do j=2,ny 1358 phi(i,j) = phi(i,j)-ty(j-1,i,1)*phi(i,j-1) 1359 end do 1360c 1361c backward sweep 1362c 1363 phi(i,ny) = phi(i,ny)/ty(ny,i,2) 1364 do jb=2,ny 1365 j = ny-jb+1 1366 phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1))/ty(j,i,2) 1367 end do 1368 end do 1369!$OMP END PARALLEL DO 1370 else 1371c 1372c y direction periodic 1373c 1374 do i=1,nx 1375 sum(i) = 0.0 1376 phi(i,0) = phi(i,ny-1) 1377 phi(i,ny) = phi(i,1) 1378 phi(i,ny+1) = phi(i,2) 1379 end do 1380c 1381!$OMP PARALLEL DO SHARED(cof,phi,ty,sum,nx,ny) PRIVATE(i,j,jb) 1382 do i=1,nx,2 1383 do j=1,ny-1 1384 phi(i,j)=cof(i,j,6)-cof(i,j,1)*phi(i-1,j)-cof(i,j,2)* 1385 + phi(i+1,j) 1386 end do 1387 do j=2,ny-2 1388 phi(i,j) = phi(i,j)-ty(j,i,1)*phi(i,j-1) 1389 end do 1390 do j=1,ny-2 1391 sum(i) = sum(i)+ty(j,i,5)*phi(i,j) 1392 end do 1393 phi(i,ny-1) = phi(i,ny-1)-sum(i) 1394c 1395c backward sweep 1396c 1397 phi(i,ny-1) = phi(i,ny-1)/ty(ny-1,i,2) 1398 phi(i,ny-2) = (phi(i,ny-2)-ty(ny-2,i,4)*phi(i,ny-1))/ 1399 + ty(ny-2,i,2) 1400 do jb=4,ny 1401 j = ny-jb+1 1402 phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1)-ty(j,i,4)* 1403 + phi(i,ny-1))/ty(j,i,2) 1404 end do 1405 end do 1406!$OMP END PARALLEL DO 1407c 1408c set odd periodic and virtual y boundaries 1409c 1410 do i=1,nx,2 1411 phi(i,0) = phi(i,ny-1) 1412 phi(i,ny) = phi(i,1) 1413 phi(i,ny+1) = phi(i,2) 1414 end do 1415c 1416c forward sweep even x lines 1417c 1418c 1419!$OMP PARALLEL DO SHARED(sum,cof,phi,ty,nx,ny) PRIVATE(i,j,jb) 1420 do i=2,nx,2 1421 do j=1,ny-1 1422 phi(i,j)=cof(i,j,6)-cof(i,j,1)*phi(i-1,j)-cof(i,j,2)* 1423 + phi(i+1,j) 1424 1425 end do 1426 do j=2,ny-2 1427 phi(i,j) = phi(i,j)-ty(j,i,1)*phi(i,j-1) 1428 end do 1429 do j=1,ny-2 1430 sum(i) = sum(i)+ty(j,i,5)*phi(i,j) 1431 end do 1432 phi(i,ny-1) = phi(i,ny-1)-sum(i) 1433c 1434c backward sweep 1435c 1436 phi(i,ny-1) = phi(i,ny-1)/ty(ny-1,i,2) 1437 phi(i,ny-2) = (phi(i,ny-2)-ty(ny-2,i,4)*phi(i,ny-1))/ 1438 + ty(ny-2,i,2) 1439 do jb=4,ny 1440 j = ny-jb+1 1441 phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1)-ty(j,i,4)* 1442 + phi(i,ny-1))/ty(j,i,2) 1443 end do 1444 end do 1445!$OMP END PARALLEL DO 1446c 1447c set even periodic and virtual y boundaries 1448c 1449 do i=2,nx,2 1450 phi(i,0) = phi(i,ny-1) 1451 phi(i,ny) = phi(i,1) 1452 phi(i,ny+1) = phi(i,2) 1453 end do 1454 end if 1455c 1456c set periodic and virtual x boundaries if necessary 1457c 1458 if (nxa.eq.0) then 1459 do j=1,ny 1460 phi(0,j) = phi(nx-1,j) 1461 phi(nx+1,j) = phi(2,j) 1462 end do 1463 end if 1464 return 1465 end 1466 1467