Changeset 114 for palm/trunk/SOURCE/poismg.f90
- Timestamp:
- Oct 10, 2007 12:03:15 AM (17 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/poismg.f90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/poismg.f90
r77 r114 8 8 ! Actual revisions: 9 9 ! ----------------- 10 ! 10 ! Boundary conditions at walls are implicitly set using flag arrays. Only 11 ! Neumann BC is allowed. Upper walls are still not realized. 12 ! Bottom and top BCs for array f_mg in restrict removed because boundary 13 ! values are not needed (right hand side of SOR iteration) 11 14 ! 12 15 ! Former revisions: … … 150 153 l = grid_level 151 154 155 ! 156 !-- Choose flag array of this level 157 SELECT CASE ( l ) 158 CASE ( 1 ) 159 flags => wall_flags_1 160 CASE ( 2 ) 161 flags => wall_flags_2 162 CASE ( 3 ) 163 flags => wall_flags_3 164 CASE ( 4 ) 165 flags => wall_flags_4 166 CASE ( 5 ) 167 flags => wall_flags_5 168 CASE ( 6 ) 169 flags => wall_flags_6 170 CASE ( 7 ) 171 flags => wall_flags_7 172 CASE ( 8 ) 173 flags => wall_flags_8 174 CASE ( 9 ) 175 flags => wall_flags_9 176 CASE ( 10 ) 177 flags => wall_flags_10 178 END SELECT 179 152 180 !$OMP PARALLEL PRIVATE (i,j,k) 153 181 !$OMP DO … … 155 183 DO j = nys_mg(l), nyn_mg(l) 156 184 DO k = nzb+1, nzt_mg(l) 157 r(k,j,i) = f_mg(k,j,i) & 158 - ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 159 - ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 160 - f2_mg(k,l) * p_mg(k+1,j,i) & 161 - f3_mg(k,l) * p_mg(k-1,j,i) & 185 r(k,j,i) = f_mg(k,j,i) & 186 - ddx2_mg(l) * & 187 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 188 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 189 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 190 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 191 - ddy2_mg(l) * & 192 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 193 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 194 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 195 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 196 - f2_mg(k,l) * p_mg(k+1,j,i) & 197 - f3_mg(k,l) * & 198 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 199 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 162 200 + f1_mg(k,l) * p_mg(k,j,i) 201 ! 202 !-- Residual within topography should be zero 203 r(k,j,i) = r(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) ) 163 204 ENDDO 164 205 ENDDO … … 181 222 182 223 ! 183 !-- Bottom and top boundary conditions 184 IF ( ibc_p_b == 1 ) THEN 185 r(nzb,:,: ) = r(nzb+1,:,:) 186 ELSE 187 r(nzb,:,: ) = 0.0 188 ENDIF 189 224 !-- Top boundary condition 225 !-- A Neumann boundary condition for r is implicitly set in routine restrict 190 226 IF ( ibc_p_t == 1 ) THEN 191 227 r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:) … … 217 253 INTEGER :: i, ic, j, jc, k, kc, l 218 254 255 REAL :: rkjim, rkjip, rkjmi, rkjmim, rkjmip, rkjpi, rkjpim, rkjpip, & 256 rkmji, rkmjim, rkmjip, rkmjmi, rkmjmim, rkmjmip, rkmjpi, rkmjpim, & 257 rkmjpip 258 219 259 REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, & 220 260 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & … … 229 269 l = grid_level 230 270 271 ! 272 !-- Choose flag array of the upper level 273 SELECT CASE ( l ) 274 CASE ( 1 ) 275 flags => wall_flags_1 276 CASE ( 2 ) 277 flags => wall_flags_2 278 CASE ( 3 ) 279 flags => wall_flags_3 280 CASE ( 4 ) 281 flags => wall_flags_4 282 CASE ( 5 ) 283 flags => wall_flags_5 284 CASE ( 6 ) 285 flags => wall_flags_6 286 CASE ( 7 ) 287 flags => wall_flags_7 288 CASE ( 8 ) 289 flags => wall_flags_8 290 CASE ( 9 ) 291 flags => wall_flags_9 292 CASE ( 10 ) 293 flags => wall_flags_10 294 END SELECT 295 231 296 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc) 232 297 !$OMP DO … … 237 302 DO kc = nzb+1, nzt_mg(l) 238 303 k = 2*kc-1 304 ! 305 !-- Use implicit Neumann BCs if the respective gridpoint is inside 306 !-- the building 307 rkjim = r(k,j,i-1) + IBITS( flags(k,j,i-1), 6, 1 ) * & 308 ( r(k,j,i) - r(k,j,i-1) ) 309 rkjip = r(k,j,i+1) + IBITS( flags(k,j,i+1), 6, 1 ) * & 310 ( r(k,j,i) - r(k,j,i+1) ) 311 rkjpi = r(k,j+1,i) + IBITS( flags(k,j+1,i), 6, 1 ) * & 312 ( r(k,j,i) - r(k,j+1,i) ) 313 rkjmi = r(k,j-1,i) + IBITS( flags(k,j-1,i), 6, 1 ) * & 314 ( r(k,j,i) - r(k,j-1,i) ) 315 rkjmim = r(k,j-1,i-1) + IBITS( flags(k,j-1,i-1), 6, 1 ) * & 316 ( r(k,j,i) - r(k,j-1,i-1) ) 317 rkjpim = r(k,j+1,i-1) + IBITS( flags(k,j+1,i-1), 6, 1 ) * & 318 ( r(k,j,i) - r(k,j+1,i-1) ) 319 rkjmip = r(k,j-1,i+1) + IBITS( flags(k,j-1,i+1), 6, 1 ) * & 320 ( r(k,j,i) - r(k,j-1,i+1) ) 321 rkjpip = r(k,j+1,i+1) + IBITS( flags(k,j+1,i+1), 6, 1 ) * & 322 ( r(k,j,i) - r(k,j+1,i+1) ) 323 rkmji = r(k-1,j,i) + IBITS( flags(k-1,j,i), 6, 1 ) * & 324 ( r(k,j,i) - r(k-1,j,i) ) 325 rkmjim = r(k-1,j,i-1) + IBITS( flags(k-1,j,i-1), 6, 1 ) * & 326 ( r(k,j,i) - r(k-1,j,i-1) ) 327 rkmjip = r(k-1,j,i+1) + IBITS( flags(k-1,j,i+1), 6, 1 ) * & 328 ( r(k,j,i) - r(k-1,j,i+1) ) 329 rkmjpi = r(k-1,j+1,i) + IBITS( flags(k-1,j+1,i), 6, 1 ) * & 330 ( r(k,j,i) - r(k-1,j+1,i) ) 331 rkmjmi = r(k-1,j-1,i) + IBITS( flags(k-1,j-1,i), 6, 1 ) * & 332 ( r(k,j,i) - r(k-1,j-1,i) ) 333 rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) * & 334 ( r(k,j,i) - r(k-1,j-1,i-1) ) 335 rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) * & 336 ( r(k,j,i) - r(k-1,j+1,i-1) ) 337 rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) * & 338 ( r(k,j,i) - r(k-1,j-1,i+1) ) 339 rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) * & 340 ( r(k,j,i) - r(k-1,j+1,i+1) ) 341 239 342 f_mg(kc,jc,ic) = 1.0 / 64.0 * ( & 240 343 8.0 * r(k,j,i) & 241 + 4.0 * ( r (k,j,i-1) + r(k,j,i+1) +&242 r (k,j+1,i) + r(k,j-1,i) )&243 + 2.0 * ( r (k,j-1,i-1) + r(k,j+1,i-1) +&244 r (k,j-1,i+1) + r(k,j+1,i+1) )&245 + 4.0 * r (k-1,j,i)&246 + 2.0 * ( r (k-1,j,i-1) + r(k-1,j,i+1) +&247 r (k-1,j+1,i) + r(k-1,j-1,i) )&248 + ( r (k-1,j-1,i-1) + r(k-1,j+1,i-1) +&249 r (k-1,j-1,i+1) + r(k-1,j+1,i+1) )&344 + 4.0 * ( rkjim + rkjip + & 345 rkjpi + rkjmi ) & 346 + 2.0 * ( rkjmim + rkjpim + & 347 rkjmip + rkjpip ) & 348 + 4.0 * rkmji & 349 + 2.0 * ( rkmjim + rkmjim + & 350 rkmjpi + rkmjmi ) & 351 + ( rkmjmim + rkmjpim + & 352 rkmjmip + rkmjpip ) & 250 353 + 4.0 * r(k+1,j,i) & 251 354 + 2.0 * ( r(k+1,j,i-1) + r(k+1,j,i+1) + & … … 254 357 r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) & 255 358 ) 359 360 ! f_mg(kc,jc,ic) = 1.0 / 64.0 * ( & 361 ! 8.0 * r(k,j,i) & 362 ! + 4.0 * ( r(k,j,i-1) + r(k,j,i+1) + & 363 ! r(k,j+1,i) + r(k,j-1,i) ) & 364 ! + 2.0 * ( r(k,j-1,i-1) + r(k,j+1,i-1) + & 365 ! r(k,j-1,i+1) + r(k,j+1,i+1) ) & 366 ! + 4.0 * r(k-1,j,i) & 367 ! + 2.0 * ( r(k-1,j,i-1) + r(k-1,j,i+1) + & 368 ! r(k-1,j+1,i) + r(k-1,j-1,i) ) & 369 ! + ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + & 370 ! r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) & 371 ! + 4.0 * r(k+1,j,i) & 372 ! + 2.0 * ( r(k+1,j,i-1) + r(k+1,j,i+1) + & 373 ! r(k+1,j+1,i) + r(k+1,j-1,i) ) & 374 ! + ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + & 375 ! r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) & 376 ! ) 256 377 ENDDO 257 378 ENDDO … … 275 396 ! 276 397 !-- Bottom and top boundary conditions 277 IF ( ibc_p_b == 1 ) THEN278 f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)279 ELSE280 f_mg(nzb,:,: ) = 0.0281 ENDIF282 283 IF ( ibc_p_t == 1 ) THEN284 f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)285 ELSE286 f_mg(nzt_mg(l)+1,:,: ) = 0.0287 ENDIF398 ! IF ( ibc_p_b == 1 ) THEN 399 ! f_mg(nzb,:,: ) = f_mg(nzb+1,:,:) 400 ! ELSE 401 ! f_mg(nzb,:,: ) = 0.0 402 ! ENDIF 403 ! 404 ! IF ( ibc_p_t == 1 ) THEN 405 ! f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:) 406 ! ELSE 407 ! f_mg(nzt_mg(l)+1,:,: ) = 0.0 408 ! ENDIF 288 409 289 410 … … 412 533 LOGICAL :: unroll 413 534 535 REAL :: wall_left, wall_north, wall_right, wall_south, wall_total, wall_top 536 414 537 REAL, DIMENSION(nzb:nzt_mg(grid_level)+1, & 415 538 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & … … 418 541 419 542 l = grid_level 543 544 ! 545 !-- Choose flag array of this level 546 SELECT CASE ( l ) 547 CASE ( 1 ) 548 flags => wall_flags_1 549 CASE ( 2 ) 550 flags => wall_flags_2 551 CASE ( 3 ) 552 flags => wall_flags_3 553 CASE ( 4 ) 554 flags => wall_flags_4 555 CASE ( 5 ) 556 flags => wall_flags_5 557 CASE ( 6 ) 558 flags => wall_flags_6 559 CASE ( 7 ) 560 flags => wall_flags_7 561 CASE ( 8 ) 562 flags => wall_flags_8 563 CASE ( 9 ) 564 flags => wall_flags_9 565 CASE ( 10 ) 566 flags => wall_flags_10 567 END SELECT 420 568 421 569 unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. & … … 434 582 DO j = nys_mg(l) + 2 - colour, nyn_mg(l), 2 435 583 DO k = nzb+1, nzt_mg(l), 2 584 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 585 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 586 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 587 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 588 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 589 ! ) 590 436 591 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 437 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 438 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 439 + f2_mg(k,l) * p_mg(k+1,j,i) & 440 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 441 ) 592 ddx2_mg(l) * & 593 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 594 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 595 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 596 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 597 + ddy2_mg(l) * & 598 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 599 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 600 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 601 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 602 + f2_mg(k,l) * p_mg(k+1,j,i) & 603 + f3_mg(k,l) * & 604 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 605 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 606 - f_mg(k,j,i) ) 442 607 ENDDO 443 608 ENDDO … … 448 613 DO k = nzb+1, nzt_mg(l), 2 449 614 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 450 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 451 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 452 + f2_mg(k,l) * p_mg(k+1,j,i) & 453 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 454 ) 615 ddx2_mg(l) * & 616 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 617 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 618 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 619 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 620 + ddy2_mg(l) * & 621 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 622 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 623 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 624 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 625 + f2_mg(k,l) * p_mg(k+1,j,i) & 626 + f3_mg(k,l) * & 627 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 628 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 629 - f_mg(k,j,i) ) 455 630 ENDDO 456 631 ENDDO … … 461 636 DO k = nzb+2, nzt_mg(l), 2 462 637 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 463 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 464 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 465 + f2_mg(k,l) * p_mg(k+1,j,i) & 466 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 467 ) 638 ddx2_mg(l) * & 639 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 640 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 641 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 642 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 643 + ddy2_mg(l) * & 644 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 645 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 646 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 647 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 648 + f2_mg(k,l) * p_mg(k+1,j,i) & 649 + f3_mg(k,l) * & 650 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 651 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 652 - f_mg(k,j,i) ) 468 653 ENDDO 469 654 ENDDO … … 474 659 DO k = nzb+2, nzt_mg(l), 2 475 660 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 476 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 477 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 478 + f2_mg(k,l) * p_mg(k+1,j,i) & 479 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 480 ) 661 ddx2_mg(l) * & 662 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 663 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 664 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 665 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 666 + ddy2_mg(l) * & 667 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 668 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 669 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 670 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 671 + f2_mg(k,l) * p_mg(k+1,j,i) & 672 + f3_mg(k,l) * & 673 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 674 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 675 - f_mg(k,j,i) ) 481 676 ENDDO 482 677 ENDDO … … 496 691 j = jj 497 692 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 498 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 499 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 500 + f2_mg(k,l) * p_mg(k+1,j,i) & 501 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 502 ) 693 ddx2_mg(l) * & 694 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 695 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 696 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 697 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 698 + ddy2_mg(l) * & 699 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 700 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 701 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 702 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 703 + f2_mg(k,l) * p_mg(k+1,j,i) & 704 + f3_mg(k,l) * & 705 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 706 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 707 - f_mg(k,j,i) ) 503 708 j = jj+2 504 709 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 505 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 506 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 507 + f2_mg(k,l) * p_mg(k+1,j,i) & 508 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 509 ) 510 ! j = jj+4 511 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 512 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 513 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 514 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 515 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 516 ! ) 517 ! j = jj+6 518 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 519 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 520 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 521 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 522 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 523 ! ) 710 ddx2_mg(l) * & 711 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 712 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 713 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 714 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 715 + ddy2_mg(l) * & 716 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 717 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 718 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 719 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 720 + f2_mg(k,l) * p_mg(k+1,j,i) & 721 + f3_mg(k,l) * & 722 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 723 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 724 - f_mg(k,j,i) ) 524 725 ENDDO 525 726 … … 529 730 j =jj 530 731 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 531 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 532 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 533 + f2_mg(k,l) * p_mg(k+1,j,i) & 534 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 535 ) 732 ddx2_mg(l) * & 733 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 734 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 735 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 736 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 737 + ddy2_mg(l) * & 738 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 739 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 740 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 741 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 742 + f2_mg(k,l) * p_mg(k+1,j,i) & 743 + f3_mg(k,l) * & 744 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 745 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 746 - f_mg(k,j,i) ) 536 747 j = jj+2 537 748 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 538 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 539 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 540 + f2_mg(k,l) * p_mg(k+1,j,i) & 541 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 542 ) 543 ! j = jj+4 544 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 545 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 546 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 547 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 548 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 549 ! ) 550 ! j = jj+6 551 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 552 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 553 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 554 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 555 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 556 ! ) 749 ddx2_mg(l) * & 750 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 751 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 752 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 753 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 754 + ddy2_mg(l) * & 755 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 756 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 757 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 758 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 759 + f2_mg(k,l) * p_mg(k+1,j,i) & 760 + f3_mg(k,l) * & 761 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 762 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 763 - f_mg(k,j,i) ) 557 764 ENDDO 558 765 … … 562 769 j =jj 563 770 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 564 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 565 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 566 + f2_mg(k,l) * p_mg(k+1,j,i) & 567 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 568 ) 771 ddx2_mg(l) * & 772 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 773 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 774 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 775 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 776 + ddy2_mg(l) * & 777 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 778 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 779 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 780 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 781 + f2_mg(k,l) * p_mg(k+1,j,i) & 782 + f3_mg(k,l) * & 783 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 784 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 785 - f_mg(k,j,i) ) 569 786 j = jj+2 570 787 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 571 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 572 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 573 + f2_mg(k,l) * p_mg(k+1,j,i) & 574 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 575 ) 576 ! j = jj+4 577 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 578 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 579 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 580 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 581 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 582 ! ) 583 ! j = jj+6 584 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 585 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 586 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 587 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 588 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 589 ! ) 788 ddx2_mg(l) * & 789 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 790 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 791 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 792 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 793 + ddy2_mg(l) * & 794 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 795 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 796 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 797 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 798 + f2_mg(k,l) * p_mg(k+1,j,i) & 799 + f3_mg(k,l) * & 800 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 801 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 802 - f_mg(k,j,i) ) 590 803 ENDDO 591 804 … … 595 808 j =jj 596 809 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 597 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 598 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 599 + f2_mg(k,l) * p_mg(k+1,j,i) & 600 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 601 ) 810 ddx2_mg(l) * & 811 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 812 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 813 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 814 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 815 + ddy2_mg(l) * & 816 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 817 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 818 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 819 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 820 + f2_mg(k,l) * p_mg(k+1,j,i) & 821 + f3_mg(k,l) * & 822 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 823 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 824 - f_mg(k,j,i) ) 602 825 j = jj+2 603 826 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 604 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 605 + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 606 + f2_mg(k,l) * p_mg(k+1,j,i) & 607 + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 608 ) 609 ! j = jj+4 610 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 611 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 612 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 613 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 614 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 615 ! ) 616 ! j = jj+6 617 ! p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 618 ! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 619 ! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 620 ! + f2_mg(k,l) * p_mg(k+1,j,i) & 621 ! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) & 622 ! ) 827 ddx2_mg(l) * & 828 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 829 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 830 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 831 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 832 + ddy2_mg(l) * & 833 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 834 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 835 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 836 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 837 + f2_mg(k,l) * p_mg(k+1,j,i) & 838 + f3_mg(k,l) * & 839 ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 840 ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) & 841 - f_mg(k,j,i) ) 623 842 ENDDO 624 843 … … 669 888 ENDDO 670 889 890 ! 891 !-- Set pressure within topography and at the topography surfaces 892 !$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total) 893 !$OMP DO 894 DO i = nxl_mg(l), nxr_mg(l) 895 DO j = nys_mg(l), nyn_mg(l) 896 DO k = nzb, nzt_mg(l) 897 ! 898 !-- First, set pressure inside topography to zero 899 p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) ) 900 ! 901 !-- Second, determine if the gridpoint inside topography is adjacent 902 !-- to a wall and set its value to a value given by the average of 903 !-- those values obtained from Neumann boundary condition 904 wall_left = IBITS( flags(k,j,i-1), 5, 1 ) 905 wall_right = IBITS( flags(k,j,i+1), 4, 1 ) 906 wall_south = IBITS( flags(k,j-1,i), 3, 1 ) 907 wall_north = IBITS( flags(k,j+1,i), 2, 1 ) 908 wall_top = IBITS( flags(k+1,j,i), 0, 1 ) 909 wall_total = wall_left + wall_right + wall_south + wall_north + & 910 wall_top 911 912 IF ( wall_total > 0.0 ) THEN 913 p_mg(k,j,i) = 1.0 / wall_total * & 914 ( wall_left * p_mg(k,j,i-1) + & 915 wall_right * p_mg(k,j,i+1) + & 916 wall_south * p_mg(k,j-1,i) + & 917 wall_north * p_mg(k,j+1,i) + & 918 wall_top * p_mg(k+1,j,i) ) 919 ENDIF 920 ENDDO 921 ENDDO 922 ENDDO 923 !$OMP END PARALLEL 924 925 ! 926 !-- One more time horizontal boundary conditions 927 CALL exchange_horiz( p_mg ) 671 928 672 929 END SUBROUTINE redblack
Note: See TracChangeset
for help on using the changeset viewer.
