- Timestamp:
- May 3, 2016 11:27:17 AM (9 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/poismg_fast_mod.f90 (modified) (34 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/poismg_fast_mod.f90
r1851 r1898 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Bugfix: bottom and top boundary condition in resid_fast 22 ! Bugfix: restriction at nzb+1 23 ! formatting adjustments, variable descriptions added in some declaration blocks 22 24 ! 23 25 ! Former revisions: … … 247 249 USE control_parameters, & 248 250 ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,& 249 inflow_n, inflow_r, inflow_s, masking_method, nest_bound_l,&250 nest_bound_ n, nest_bound_r, nest_bound_s, outflow_l,&251 outflow_ n, outflow_r, outflow_s, topography251 inflow_n, inflow_r, inflow_s, nest_bound_l, nest_bound_n, & 252 nest_bound_r, nest_bound_s, outflow_l, outflow_n, outflow_r, & 253 outflow_s 252 254 253 255 USE grid_variables, & … … 255 257 256 258 USE indices, & 257 ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, & 258 wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7, & 259 wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, & 260 nys_mg, nyn_mg, nzb, nzt_mg 259 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 261 260 262 261 IMPLICIT NONE 263 262 264 INTEGER(iwp) :: i 265 INTEGER(iwp) :: j 266 INTEGER(iwp) :: k 267 INTEGER(iwp) :: l 268 INTEGER(iwp) :: km1 !<269 INTEGER(iwp) :: kp1 !<263 INTEGER(iwp) :: i !< index variable along x 264 INTEGER(iwp) :: j !< index variable along y 265 INTEGER(iwp) :: k !< index variable along z 266 INTEGER(iwp) :: l !< index indicating grid level 267 INTEGER(iwp) :: km1 !< index variable along z dimension (k-1) 268 INTEGER(iwp) :: kp1 !< index variable along z dimension (k+1) 270 269 271 270 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 272 271 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 273 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< 272 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< velocity divergence 274 273 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 275 274 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 276 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< 275 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< perturbation pressure 277 276 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 278 277 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 279 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !< 278 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !< residuum of perturbation pressure 280 279 281 280 ! … … 284 283 285 284 CALL cpu_log( log_point_s(53), 'resid', 'start' ) 286 IF ( topography == 'flat' .OR. masking_method ) THEN 287 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 288 !$OMP DO 289 DO i = nxl_mg(l), nxr_mg(l) 290 DO j = nys_mg(l), nyn_mg(l) 285 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 286 !$OMP DO 287 DO i = nxl_mg(l), nxr_mg(l) 288 DO j = nys_mg(l), nyn_mg(l) 291 289 !DIR$ IVDEP 292 293 294 295 r(k,j,i) = f_mg(k,j,i)&296 - ddx2_mg(l) * &297 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &298 - ddy2_mg(l) * &299 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &300 - f2_mg_b(k,l) * p_mg(kp1,j,i) &301 - f3_mg_b(k,l) * p_mg(km1,j,i) &290 DO k = ind_even_odd+1, nzt_mg(l) 291 km1 = k-ind_even_odd-1 292 kp1 = k-ind_even_odd 293 r(k,j,i) = f_mg(k,j,i) & 294 - ddx2_mg(l) * & 295 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 296 - ddy2_mg(l) * & 297 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 298 - f2_mg_b(k,l) * p_mg(kp1,j,i) & 299 - f3_mg_b(k,l) * p_mg(km1,j,i) & 302 300 + f1_mg_b(k,l) * p_mg(k,j,i) 303 304 305 306 307 308 r(k,j,i) = f_mg(k,j,i)&309 - ddx2_mg(l) * &310 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &311 - ddy2_mg(l) * &312 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &313 - f2_mg_b(k,l) * p_mg(kp1,j,i) &314 - f3_mg_b(k,l) * p_mg(km1,j,i) &301 ENDDO 302 !DIR$ IVDEP 303 DO k = nzb+1, ind_even_odd 304 km1 = k+ind_even_odd 305 kp1 = k+ind_even_odd+1 306 r(k,j,i) = f_mg(k,j,i) & 307 - ddx2_mg(l) * & 308 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 309 - ddy2_mg(l) * & 310 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 311 - f2_mg_b(k,l) * p_mg(kp1,j,i) & 312 - f3_mg_b(k,l) * p_mg(km1,j,i) & 315 313 + f1_mg_b(k,l) * p_mg(k,j,i) 316 ENDDO 317 ENDDO 318 ENDDO 319 !$OMP END PARALLEL 320 ELSE 321 ! 322 !-- Choose flag array of this level 323 SELECT CASE ( l ) 324 CASE ( 1 ) 325 flags => wall_flags_1 326 CASE ( 2 ) 327 flags => wall_flags_2 328 CASE ( 3 ) 329 flags => wall_flags_3 330 CASE ( 4 ) 331 flags => wall_flags_4 332 CASE ( 5 ) 333 flags => wall_flags_5 334 CASE ( 6 ) 335 flags => wall_flags_6 336 CASE ( 7 ) 337 flags => wall_flags_7 338 CASE ( 8 ) 339 flags => wall_flags_8 340 CASE ( 9 ) 341 flags => wall_flags_9 342 CASE ( 10 ) 343 flags => wall_flags_10 344 END SELECT 345 346 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 347 !$OMP DO 348 DO i = nxl_mg(l), nxr_mg(l) 349 DO j = nys_mg(l), nyn_mg(l) 350 351 !DIR$ IVDEP 352 DO k = ind_even_odd+1, nzt_mg(l) 353 km1 = k-ind_even_odd-1 354 kp1 = k-ind_even_odd 355 356 r(k,j,i) = f_mg(k,j,i) & 357 - ddx2_mg(l) * & 358 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 )) & 359 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) & 360 - ddy2_mg(l) & 361 * ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) & 362 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) & 363 - f2_mg(k,l) * p_mg(kp1,j,i) & 364 - f3_mg(k,l) * & 365 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 )) & 366 + f1_mg(k,l) * p_mg(k,j,i) 367 ENDDO 368 369 !DIR$ IVDEP 370 DO k = nzb+1, ind_even_odd 371 km1 = k+ind_even_odd 372 kp1 = k+ind_even_odd+1 373 374 r(k,j,i) = f_mg(k,j,i) & 375 - ddx2_mg(l) * & 376 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 )) & 377 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) & 378 - ddy2_mg(l) & 379 * ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) & 380 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) & 381 - f2_mg(k,l) * p_mg(kp1,j,i) & 382 - f3_mg(k,l) * & 383 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 )) & 384 + f1_mg(k,l) * p_mg(k,j,i) 385 ENDDO 386 ! 387 !-- The residual within topography should be zero 388 WHERE( BTEST(flags(nzb+1:nzt_mg(l),j,i), 6) ) 389 r(nzb+1:nzt_mg(l),j,i) = 0.0 390 END WHERE 391 392 ENDDO 393 ENDDO 394 !$OMP END PARALLEL 395 396 ENDIF 397 314 ENDDO 315 ENDDO 316 ENDDO 317 !$OMP END PARALLEL 398 318 ! 399 319 !-- Horizontal boundary conditions … … 419 339 420 340 ! 421 !-- Boundary conditions at bottom and top of the domain.outflow_l, outflow_n, 422 !-- These points are not handled by the above loop. Points may be within 341 !-- Boundary conditions at bottom and top of the domain. Points may be within 423 342 !-- buildings, but that doesn't matter. 424 343 IF ( ibc_p_b == 1 ) THEN 425 r(nzb,:,: ) = r(nzb+1,:,:) 344 ! 345 !-- equivalent to r(nzb,:,: ) = r(nzb+1,:,:) 346 r(nzb,:,: ) = r(ind_even_odd+1,:,:) 426 347 ELSE 427 348 r(nzb,:,: ) = 0.0_wp … … 429 350 430 351 IF ( ibc_p_t == 1 ) THEN 431 r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:) 352 ! 353 !-- equivalent to r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:) 354 r(nzt_mg(l)+1,:,: ) = r(ind_even_odd,:,:) 432 355 ELSE 433 356 r(nzt_mg(l)+1,:,: ) = 0.0_wp … … 450 373 USE control_parameters, & 451 374 ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,& 452 inflow_n, inflow_r, inflow_s, masking_method, nest_bound_l,&453 nest_bound_ n, nest_bound_r, nest_bound_s, outflow_l,&454 outflow_ n, outflow_r, outflow_s, topography375 inflow_n, inflow_r, inflow_s, nest_bound_l, nest_bound_n, & 376 nest_bound_r, nest_bound_s, outflow_l, outflow_n, outflow_r, & 377 outflow_s 455 378 456 379 USE indices, & 457 ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, & 458 wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7, & 459 wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, & 460 nys_mg, nyn_mg, nzb, nzt_mg 380 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 461 381 462 382 IMPLICIT NONE 463 383 464 INTEGER(iwp) :: i !< 465 INTEGER(iwp) :: ic !< 466 INTEGER(iwp) :: j !< 467 INTEGER(iwp) :: jc !< 468 INTEGER(iwp) :: k !< 469 INTEGER(iwp) :: kc !< 470 INTEGER(iwp) :: l !< 471 INTEGER(iwp) :: km1 !< 472 INTEGER(iwp) :: kp1 !< 473 474 REAL(wp) :: rkjim !< 475 REAL(wp) :: rkjip !< 476 REAL(wp) :: rkjmi !< 477 REAL(wp) :: rkjmim !< 478 REAL(wp) :: rkjmip !< 479 REAL(wp) :: rkjpi !< 480 REAL(wp) :: rkjpim !< 481 REAL(wp) :: rkjpip !< 482 REAL(wp) :: rkmji !< 483 REAL(wp) :: rkmjim !< 484 REAL(wp) :: rkmjip !< 485 REAL(wp) :: rkmjmi !< 486 REAL(wp) :: rkmjmim !< 487 REAL(wp) :: rkmjmip !< 488 REAL(wp) :: rkmjpi !< 489 REAL(wp) :: rkmjpim !< 490 REAL(wp) :: rkmjpip !< 384 INTEGER(iwp) :: i !< index variable along x on finer grid 385 INTEGER(iwp) :: ic !< index variable along x on coarser grid 386 INTEGER(iwp) :: j !< index variable along y on finer grid 387 INTEGER(iwp) :: jc !< index variable along y on coarser grid 388 INTEGER(iwp) :: k !< index variable along z on finer grid 389 INTEGER(iwp) :: kc !< index variable along z on coarser grid 390 INTEGER(iwp) :: l !< index indicating finer grid level 391 INTEGER(iwp) :: km1 !< index variable along z dimension (k-1 on finer level) 392 INTEGER(iwp) :: kp1 !< index variable along z dimension (k+1 on finer level) 393 491 394 492 395 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 493 396 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 494 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< 397 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 398 f_mg !< Residual on coarser grid level 495 399 496 400 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1, & 497 401 nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, & 498 nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: r !< 402 nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: & 403 r !< Residual on finer grid level 499 404 500 405 ! … … 503 408 504 409 CALL cpu_log( log_point_s(54), 'restrict', 'start' ) 505 506 IF ( topography == 'flat' .OR. masking_method ) THEN 507 ! 508 !-- No wall treatment 509 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1) 510 DO ic = nxl_mg(l), nxr_mg(l) 511 i = 2*ic 512 !$OMP DO SCHEDULE( STATIC ) 513 DO jc = nys_mg(l), nyn_mg(l) 514 ! 515 !-- Calculation for the first point along k 516 j = 2*jc 517 k = ind_even_odd+1 410 ! 411 !-- No wall treatment 412 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1) 413 DO ic = nxl_mg(l), nxr_mg(l) 414 i = 2*ic 415 !$OMP DO SCHEDULE( STATIC ) 416 DO jc = nys_mg(l), nyn_mg(l) 417 ! 418 !-- Calculation for the first point along k 419 j = 2*jc 420 ! 421 !-- Calculation for the other points along k 422 !DIR$ IVDEP 423 DO k = ind_even_odd+1, nzt_mg(l+1) ! Fine grid at this point 424 km1 = k-ind_even_odd-1 518 425 kp1 = k-ind_even_odd 519 kc = k-ind_even_odd 520 f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( & 521 8.0_wp * r(k,j,i) & 522 + 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) + & 523 r(k,j+1,i) + r(k,j-1,i) ) & 524 + 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) + & 525 r(k,j-1,i+1) + r(k,j+1,i+1) ) & 526 + 16.0_wp * r(k,j,i) & 527 + 4.0_wp * r(kp1,j,i) & 528 + 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) + & 529 r(kp1,j+1,i) + r(kp1,j-1,i) ) & 530 + ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + & 531 r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) & 532 ) 533 ! 534 !-- Calculation for the other points along k 535 !DIR$ IVDEP 536 DO k = ind_even_odd+2, nzt_mg(l+1) ! Fine grid at this point 537 km1 = k-ind_even_odd-1 538 kp1 = k-ind_even_odd 539 kc = k-ind_even_odd ! Coarse grid index 540 541 f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( & 542 8.0_wp * r(k,j,i) & 543 + 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) + & 544 r(k,j+1,i) + r(k,j-1,i) ) & 545 + 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) + & 546 r(k,j-1,i+1) + r(k,j+1,i+1) ) & 547 + 4.0_wp * r(km1,j,i) & 548 + 2.0_wp * ( r(km1,j,i-1) + r(km1,j,i+1) + & 549 r(km1,j+1,i) + r(km1,j-1,i) ) & 550 + ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + & 551 r(km1,j-1,i+1) + r(km1,j+1,i+1) ) & 552 + 4.0_wp * r(kp1,j,i) & 553 + 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) + & 554 r(kp1,j+1,i) + r(kp1,j-1,i) ) & 555 + ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + & 556 r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) & 557 ) 558 ENDDO 559 ENDDO 560 !$OMP ENDDO nowait 561 ENDDO 562 !$OMP END PARALLEL 563 564 ELSE 565 ! 566 !-- Choose flag array of the upper level 567 SELECT CASE ( l+1 ) 568 CASE ( 1 ) 569 flags => wall_flags_1 570 CASE ( 2 ) 571 flags => wall_flags_2 572 CASE ( 3 ) 573 flags => wall_flags_3 574 CASE ( 4 ) 575 flags => wall_flags_4 576 CASE ( 5 ) 577 flags => wall_flags_5 578 CASE ( 6 ) 579 flags => wall_flags_6 580 CASE ( 7 ) 581 flags => wall_flags_7 582 CASE ( 8 ) 583 flags => wall_flags_8 584 CASE ( 9 ) 585 flags => wall_flags_9 586 CASE ( 10 ) 587 flags => wall_flags_10 588 END SELECT 589 590 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1,rkjim,rkjip,rkjpi, & 591 !$OMP rkjmim,rkjpim,rkjmip,rkjpip,rkmji,rkmjim, & 592 !$OMP rkmjip,rkmjpi,rkmjmi,rkmjmim,rkmjpim, & 593 !$OMP rkmjmip,rkmjpip) 594 !$OMP DO 595 DO ic = nxl_mg(l), nxr_mg(l) 596 i = 2*ic 597 DO jc = nys_mg(l), nyn_mg(l) 598 j = 2*jc 599 !DIR$ IVDEP 600 DO k = ind_even_odd+1, nzt_mg(l+1) !fine grid at this point 601 602 km1 = k-ind_even_odd-1 603 kp1 = k-ind_even_odd 604 kc = k-ind_even_odd !Coarse grid index 605 ! 606 !-- Use implicit Neumann BCs if the respective gridpoint is 607 !-- inside the building 608 rkjim = MERGE(r(k,j,i),r(k,j,i-1),BTEST( flags(k,j,i-1), 6)) 609 rkjip = MERGE(r(k,j,i),r(k,j,i+1),BTEST( flags(k,j,i+1), 6)) 610 rkjpi = MERGE(r(k,j,i),r(k,j+1,i),BTEST( flags(k,j+1,i), 6)) 611 rkjmi = MERGE(r(k,j,i),r(k,j-1,i),BTEST( flags(k,j-1,i), 6)) 612 613 rkjmim = MERGE(r(k,j,i),r(k,j-1,i-1),BTEST( flags(k,j-1,i-1), 6)) 614 rkjpim = MERGE(r(k,j,i),r(k,j+1,i-1),BTEST( flags(k,j+1,i-1), 6)) 615 rkjmip = MERGE(r(k,j,i),r(k,j-1,i+1),BTEST( flags(k,j-1,i+1), 6)) 616 rkjpip = MERGE(r(k,j,i),r(k,j+1,i+1),BTEST( flags(k,j+1,i+1), 6)) 617 618 rkmji = MERGE(r(k,j,i),r(km1,j,i) ,BTEST( flags(km1,j,i) , 6)) 619 rkmjim = MERGE(r(k,j,i),r(km1,j,i-1),BTEST( flags(km1,j,i-1), 6)) 620 rkmjip = MERGE(r(k,j,i),r(km1,j,i+1),BTEST( flags(km1,j,i+1), 6)) 621 rkmjpi = MERGE(r(k,j,i),r(km1,j+1,i),BTEST( flags(km1,j+1,i), 6)) 622 rkmjmi = MERGE(r(k,j,i),r(km1,j-1,i),BTEST( flags(km1,j-1,i), 6)) 623 624 rkmjmim = MERGE(r(k,j,i),r(km1,j-1,i-1),BTEST( flags(km1,j-1,i-1), 6)) 625 rkmjpim = MERGE(r(k,j,i),r(km1,j+1,i-1),BTEST( flags(km1,j+1,i-1), 6)) 626 rkmjmip = MERGE(r(k,j,i),r(km1,j-1,i+1),BTEST( flags(km1,j-1,i+1), 6)) 627 rkmjpip = MERGE(r(k,j,i),r(km1,j+1,i+1),BTEST( flags(km1,j+1,i+1), 6)) 628 629 f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( & 630 8.0_wp * r(k,j,i) & 631 + 4.0_wp * ( rkjim + rkjip + rkjpi + rkjmi ) & 632 + 2.0_wp * ( rkjmim + rkjpim + rkjmip + rkjpip ) & 633 + 4.0_wp * rkmji & 634 + 2.0_wp * ( rkmjim + rkmjip + rkmjpi + rkmjmi ) & 635 + ( rkmjmim + rkmjpim + rkmjmip + rkmjpip ) & 636 + 4.0_wp * r(kp1,j,i) & 637 + 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) + & 638 r(kp1,j+1,i) + r(kp1,j-1,i) ) & 639 + ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + & 640 r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) & 641 ) 642 ENDDO 643 ENDDO 644 ENDDO 645 !$OMP END PARALLEL 646 647 ENDIF 648 426 kc = k-ind_even_odd ! Coarse grid index 427 428 f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( & 429 8.0_wp * r(k,j,i) & 430 + 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) + & 431 r(k,j+1,i) + r(k,j-1,i) ) & 432 + 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) + & 433 r(k,j-1,i+1) + r(k,j+1,i+1) ) & 434 + 4.0_wp * r(km1,j,i) & 435 + 2.0_wp * ( r(km1,j,i-1) + r(km1,j,i+1) + & 436 r(km1,j+1,i) + r(km1,j-1,i) ) & 437 + ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + & 438 r(km1,j-1,i+1) + r(km1,j+1,i+1) ) & 439 + 4.0_wp * r(kp1,j,i) & 440 + 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) + & 441 r(kp1,j+1,i) + r(kp1,j-1,i) ) & 442 + ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + & 443 r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) & 444 ) 445 ENDDO 446 ENDDO 447 !$OMP ENDDO nowait 448 ENDDO 449 !$OMP END PARALLEL 450 451 ! 452 !-- Ghost point exchange 453 CALL exchange_horiz( f_mg, 1) 649 454 ! 650 455 !-- Horizontal boundary conditions 651 CALL exchange_horiz( f_mg, 1)652 653 456 IF ( .NOT. bc_lr_cyc ) THEN 654 457 IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN … … 672 475 !-- Boundary conditions at bottom and top of the domain. 673 476 !-- These points are not handled by the above loop. Points may be within 674 !-- buildings, but that doesn't matter. 477 !-- buildings, but that doesn't matter. Remark: f_mg is ordered sequentielly 478 !-- after interpolation on coarse grid (is ordered in odd-even blocks further 479 !-- below). 675 480 IF ( ibc_p_b == 1 ) THEN 676 481 f_mg(nzb,:,: ) = f_mg(nzb+1,:,:) … … 687 492 CALL cpu_log( log_point_s(54), 'restrict', 'stop' ) 688 493 ! 689 !-- Why do we need a sorting here? 494 !-- Since residual is in sequential order after interpolation, an additional 495 !-- sorting in odd-even blocks along z dimension is required at this point. 690 496 CALL sort_k_to_even_odd_blocks( f_mg , l) 691 497 … … 713 519 IMPLICIT NONE 714 520 715 INTEGER(iwp) :: i !< 716 INTEGER(iwp) :: j !< 717 INTEGER(iwp) :: k !< 718 INTEGER(iwp) :: l !< 719 INTEGER(iwp) :: kp1 !< 720 INTEGER(iwp) :: ke !< Index for prolog even721 INTEGER(iwp) :: ko !< Index for prolog odd521 INTEGER(iwp) :: i !< index variable along x on coarser grid level 522 INTEGER(iwp) :: j !< index variable along y on coarser grid level 523 INTEGER(iwp) :: k !< index variable along z on coarser grid level 524 INTEGER(iwp) :: l !< index indicating finer grid level 525 INTEGER(iwp) :: kp1 !< index variable along z 526 INTEGER(iwp) :: ke !< index for prolog even 527 INTEGER(iwp) :: ko !< index for prolog odd 722 528 723 529 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, & 724 530 nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, & 725 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: p !< 531 nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: & 532 p !< perturbation pressure on coarser grid level 726 533 727 534 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 728 535 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 729 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: temp !< 536 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 537 temp !< perturbation pressure on finer grid level 730 538 731 539 … … 772 580 p(kp1,j,i) + p(kp1,j,i+1) + & 773 581 p(kp1,j+1,i) + p(kp1,j+1,i+1) ) 582 774 583 ENDDO 775 584 … … 804 613 p(kp1,j,i) + p(kp1,j,i+1) + & 805 614 p(kp1,j+1,i) + p(kp1,j+1,i+1) ) 615 806 616 ENDDO 807 617 … … 836 646 !-- Bottom and top boundary conditions 837 647 IF ( ibc_p_b == 1 ) THEN 648 ! 649 !-- equivalent to temp(nzb,:,: ) = temp(nzb+1,:,:) 838 650 temp(nzb,:,: ) = temp(ind_even_odd+1,:,:) 839 651 ELSE … … 842 654 843 655 IF ( ibc_p_t == 1 ) THEN 656 ! 657 !-- equivalent to temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:) 844 658 temp(nzt_mg(l)+1,:,: ) = temp(ind_even_odd,:,:) 845 659 ELSE … … 866 680 USE control_parameters, & 867 681 ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,& 868 inflow_n, inflow_r, inflow_s, masking_method, nest_bound_l,&869 nest_bound_ n, nest_bound_r, nest_bound_s, ngsrb,&870 outflow_ l, outflow_n, outflow_r, outflow_s, topography682 inflow_n, inflow_r, inflow_s, nest_bound_l, nest_bound_n, & 683 nest_bound_r, nest_bound_s, ngsrb, outflow_l, outflow_n, & 684 outflow_r, outflow_s 871 685 872 686 USE grid_variables, & … … 874 688 875 689 USE indices, & 876 ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, & 877 wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7, & 878 wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, & 879 nys_mg, nyn_mg, nzb, nzt_mg 690 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 880 691 881 692 IMPLICIT NONE 882 693 883 INTEGER(iwp) :: color !< 884 INTEGER(iwp) :: i !< 885 INTEGER(iwp) :: ic !< 886 INTEGER(iwp) :: j !< 887 INTEGER(iwp) :: jc !< 888 INTEGER(iwp) :: jj !< 889 INTEGER(iwp) :: k !< 890 INTEGER(iwp) :: l !< 891 INTEGER(iwp) :: n !< 892 INTEGER(iwp) :: km1 !< 893 INTEGER(iwp) :: kp1 !< 894 895 LOGICAL :: unroll !< 896 897 REAL(wp) :: wall_left !< 898 REAL(wp) :: wall_north !< 899 REAL(wp) :: wall_right !< 900 REAL(wp) :: wall_south !< 901 REAL(wp) :: wall_total !< 902 REAL(wp) :: wall_top !< 694 INTEGER(iwp) :: color !< grid point color, either red or black 695 INTEGER(iwp) :: i !< index variable along x 696 INTEGER(iwp) :: ic !< index variable along x 697 INTEGER(iwp) :: j !< index variable along y 698 INTEGER(iwp) :: jc !< index variable along y 699 INTEGER(iwp) :: jj !< index variable along y 700 INTEGER(iwp) :: k !< index variable along z 701 INTEGER(iwp) :: l !< grid level 702 INTEGER(iwp) :: n !< loop variable GauÃ-Seidel iterations 703 INTEGER(iwp) :: km1 !< index variable (k-1) 704 INTEGER(iwp) :: kp1 !< index variable (k+1) 705 706 LOGICAL :: unroll !< flag indicating whether loop unrolling is possible 903 707 904 708 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 905 709 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 906 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !< 710 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 711 f_mg !< residual of perturbation pressure 907 712 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 908 713 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 909 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< 714 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 715 p_mg !< perturbation pressure 910 716 911 717 l = grid_level 912 718 913 ! 914 !-- Choose flag array of this level 915 IF ( topography /= 'flat' .AND. .NOT. masking_method ) THEN 916 917 SELECT CASE ( l ) 918 CASE ( 1 ) 919 flags => wall_flags_1 920 CASE ( 2 ) 921 flags => wall_flags_2 922 CASE ( 3 ) 923 flags => wall_flags_3 924 CASE ( 4 ) 925 flags => wall_flags_4 926 CASE ( 5 ) 927 flags => wall_flags_5 928 CASE ( 6 ) 929 flags => wall_flags_6 930 CASE ( 7 ) 931 flags => wall_flags_7 932 CASE ( 8 ) 933 flags => wall_flags_8 934 CASE ( 9 ) 935 flags => wall_flags_9 936 CASE ( 10 ) 937 flags => wall_flags_10 938 END SELECT 939 940 ENDIF 941 942 unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. & 719 unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. & 943 720 MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 ) 944 721 … … 946 723 947 724 DO color = 1, 2 948 ! 949 !-- No wall treatment in case of flat topography 950 IF ( topography == 'flat' .OR. masking_method ) THEN 951 952 IF ( .NOT. unroll ) THEN 953 954 CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' ) 955 ! 956 !-- Without unrolling of loops, no cache optimization 957 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 958 !$OMP DO 959 DO i = nxl_mg(l), nxr_mg(l), 2 960 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 961 !DIR$ IVDEP 962 DO k = ind_even_odd+1, nzt_mg(l) 963 km1 = k-ind_even_odd-1 964 kp1 = k-ind_even_odd 965 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 966 ddx2_mg(l) * & 967 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 968 + ddy2_mg(l) * & 969 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 970 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 971 + f3_mg_b(k,l) * p_mg(km1,j,i) & 972 - f_mg(k,j,i) ) 973 ENDDO 974 ENDDO 975 ENDDO 976 977 !$OMP DO 978 DO i = nxl_mg(l)+1, nxr_mg(l), 2 979 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 980 !DIR$ IVDEP 981 DO k = ind_even_odd+1, nzt_mg(l) 982 km1 = k-ind_even_odd-1 983 kp1 = k-ind_even_odd 984 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 985 ddx2_mg(l) * & 986 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 987 + ddy2_mg(l) * & 988 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 989 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 990 + f3_mg_b(k,l) * p_mg(km1,j,i) & 991 - f_mg(k,j,i) ) 992 ENDDO 993 ENDDO 994 ENDDO 995 996 !$OMP DO 997 DO i = nxl_mg(l), nxr_mg(l), 2 998 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 999 !DIR$ IVDEP 1000 DO k = nzb+1, ind_even_odd 1001 km1 = k+ind_even_odd 1002 kp1 = k+ind_even_odd+1 1003 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1004 ddx2_mg(l) * & 725 726 IF ( .NOT. unroll ) THEN 727 728 CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' ) 729 ! 730 !-- Without unrolling of loops, no cache optimization 731 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 732 !$OMP DO 733 DO i = nxl_mg(l), nxr_mg(l), 2 734 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 735 !DIR$ IVDEP 736 DO k = ind_even_odd+1, nzt_mg(l) 737 km1 = k-ind_even_odd-1 738 kp1 = k-ind_even_odd 739 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 740 ddx2_mg(l) * & 1005 741 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1006 742 + ddy2_mg(l) * & … … 1008 744 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1009 745 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1010 - f_mg(k,j,i) ) 1011 ENDDO 746 - f_mg(k,j,i) ) 1012 747 ENDDO 1013 748 ENDDO 1014 1015 !$OMP DO 1016 DO i = nxl_mg(l)+1, nxr_mg(l), 2 1017 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 1018 !DIR$ IVDEP 1019 DO k = nzb+1, ind_even_odd 1020 km1 = k+ind_even_odd 1021 kp1 = k+ind_even_odd+1 1022 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1023 ddx2_mg(l) * & 749 ENDDO 750 751 !$OMP DO 752 DO i = nxl_mg(l)+1, nxr_mg(l), 2 753 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 754 !DIR$ IVDEP 755 DO k = ind_even_odd+1, nzt_mg(l) 756 km1 = k-ind_even_odd-1 757 kp1 = k-ind_even_odd 758 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 759 ddx2_mg(l) * & 1024 760 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1025 761 + ddy2_mg(l) * & … … 1027 763 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1028 764 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1029 - f_mg(k,j,i) ) 1030 ENDDO 765 - f_mg(k,j,i) ) 1031 766 ENDDO 1032 767 ENDDO 1033 !$OMP END PARALLEL 1034 1035 CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' ) 1036 1037 ELSE 1038 ! 1039 !-- Loop unrolling along y, only one i loop for better cache use 1040 CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' ) 1041 1042 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj) 1043 !$OMP DO 1044 DO ic = nxl_mg(l), nxr_mg(l), 2 1045 DO jc = nys_mg(l), nyn_mg(l), 4 1046 i = ic 1047 jj = jc+2-color 1048 !DIR$ IVDEP 1049 DO k = ind_even_odd+1, nzt_mg(l) 1050 km1 = k-ind_even_odd-1 1051 kp1 = k-ind_even_odd 1052 j = jj 1053 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1054 ddx2_mg(l) * & 768 ENDDO 769 770 !$OMP DO 771 DO i = nxl_mg(l), nxr_mg(l), 2 772 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 773 !DIR$ IVDEP 774 DO k = nzb+1, ind_even_odd 775 km1 = k+ind_even_odd 776 kp1 = k+ind_even_odd+1 777 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 778 ddx2_mg(l) * & 1055 779 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1056 780 + ddy2_mg(l) * & … … 1058 782 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1059 783 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1060 - f_mg(k,j,i) ) 1061 j = jj+2 1062 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1063 ddx2_mg(l) * & 784 - f_mg(k,j,i) ) 785 ENDDO 786 ENDDO 787 ENDDO 788 789 !$OMP DO 790 DO i = nxl_mg(l)+1, nxr_mg(l), 2 791 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 792 !DIR$ IVDEP 793 DO k = nzb+1, ind_even_odd 794 km1 = k+ind_even_odd 795 kp1 = k+ind_even_odd+1 796 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 797 ddx2_mg(l) * & 1064 798 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1065 799 + ddy2_mg(l) * & … … 1067 801 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1068 802 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1069 - f_mg(k,j,i) ) 1070 ENDDO 1071 1072 i = ic+1 1073 jj = jc+color-1 1074 !DIR$ IVDEP 1075 DO k = ind_even_odd+1, nzt_mg(l) 1076 km1 = k-ind_even_odd-1 1077 kp1 = k-ind_even_odd 1078 j = jj 1079 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1080 ddx2_mg(l) * & 1081 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1082 + ddy2_mg(l) * & 1083 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 1084 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1085 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1086 - f_mg(k,j,i) ) 1087 j = jj+2 1088 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1089 ddx2_mg(l) * & 1090 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1091 + ddy2_mg(l) * & 1092 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 1093 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1094 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1095 - f_mg(k,j,i) ) 1096 ENDDO 1097 1098 i = ic 1099 jj = jc+color-1 1100 !DIR$ IVDEP 1101 DO k = nzb+1, ind_even_odd 1102 km1 = k+ind_even_odd 1103 kp1 = k+ind_even_odd+1 1104 j = jj 1105 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1106 ddx2_mg(l) * & 1107 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1108 + ddy2_mg(l) * & 1109 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 1110 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1111 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1112 - f_mg(k,j,i) ) 1113 j = jj+2 1114 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1115 ddx2_mg(l) * & 803 - f_mg(k,j,i) ) 804 ENDDO 805 ENDDO 806 ENDDO 807 !$OMP END PARALLEL 808 809 CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' ) 810 811 ELSE 812 ! 813 !-- Loop unrolling along y, only one i loop for better cache use 814 CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' ) 815 816 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj) 817 !$OMP DO 818 DO ic = nxl_mg(l), nxr_mg(l), 2 819 DO jc = nys_mg(l), nyn_mg(l), 4 820 i = ic 821 jj = jc+2-color 822 !DIR$ IVDEP 823 DO k = ind_even_odd+1, nzt_mg(l) 824 km1 = k-ind_even_odd-1 825 kp1 = k-ind_even_odd 826 j = jj 827 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 828 ddx2_mg(l) * & 1116 829 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1117 830 + ddy2_mg(l) * & … … 1119 832 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1120 833 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1121 - f_mg(k,j,i) ) 1122 ENDDO 1123 1124 i = ic+1 1125 jj = jc+2-color 1126 !DIR$ IVDEP 1127 DO k = nzb+1, ind_even_odd 1128 km1 = k+ind_even_odd 1129 kp1 = k+ind_even_odd+1 1130 j = jj 1131 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1132 ddx2_mg(l) * & 1133 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 834 - f_mg(k,j,i) ) 835 j = jj+2 836 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 837 ddx2_mg(l) * & 838 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1134 839 + ddy2_mg(l) * & 1135 840 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 1136 841 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1137 842 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1138 - f_mg(k,j,i) ) 1139 j = jj+2 1140 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 1141 ddx2_mg(l) * & 843 - f_mg(k,j,i) ) 844 ENDDO 845 846 i = ic+1 847 jj = jc+color-1 848 !DIR$ IVDEP 849 DO k = ind_even_odd+1, nzt_mg(l) 850 km1 = k-ind_even_odd-1 851 kp1 = k-ind_even_odd 852 j = jj 853 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 854 ddx2_mg(l) * & 1142 855 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1143 856 + ddy2_mg(l) * & 1144 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )&857 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 1145 858 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 1146 859 + f3_mg_b(k,l) * p_mg(km1,j,i) & 1147 - f_mg(k,j,i) ) 1148 ENDDO 1149 860 - f_mg(k,j,i) ) 861 j = jj+2 862 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 863 ddx2_mg(l) * & 864 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 865 + ddy2_mg(l) * & 866 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 867 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 868 + f3_mg_b(k,l) * p_mg(km1,j,i) & 869 - f_mg(k,j,i) ) 1150 870 ENDDO 1151 ENDDO 1152 !$OMP END PARALLEL 1153 1154 CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' ) 1155 1156 ENDIF 1157 1158 ELSE 1159 1160 IF ( .NOT. unroll ) THEN 1161 1162 CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' ) 1163 ! 1164 !-- Without unrolling of loops, no cache optimization / 1165 !-- vectorization. Therefore, the non-vectorized IBITS function 1166 !-- is still used. 1167 !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1) 1168 !$OMP DO 1169 DO i = nxl_mg(l), nxr_mg(l), 2 1170 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 1171 !DIR$ IVDEP 1172 DO k = ind_even_odd+1, nzt_mg(l) 1173 km1 = k-ind_even_odd-1 1174 kp1 = k-ind_even_odd 1175 1176 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 1177 ddx2_mg(l) * & 1178 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 1179 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 1180 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 1181 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 871 872 i = ic 873 jj = jc+color-1 874 !DIR$ IVDEP 875 DO k = nzb+1, ind_even_odd 876 km1 = k+ind_even_odd 877 kp1 = k+ind_even_odd+1 878 j = jj 879 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 880 ddx2_mg(l) * & 881 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1182 882 + ddy2_mg(l) * & 1183 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 1184 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 1185 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 1186 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 1187 + f2_mg(k,l) * p_mg(kp1,j,i) & 1188 + f3_mg(k,l) * & 1189 ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 1190 ( p_mg(k,j,i) - p_mg(km1,j,i) ) ) & 1191 - f_mg(k,j,i) ) 1192 ENDDO 883 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 884 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 885 + f3_mg_b(k,l) * p_mg(km1,j,i) & 886 - f_mg(k,j,i) ) 887 j = jj+2 888 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 889 ddx2_mg(l) * & 890 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 891 + ddy2_mg(l) * & 892 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 893 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 894 + f3_mg_b(k,l) * p_mg(km1,j,i) & 895 - f_mg(k,j,i) ) 1193 896 ENDDO 1194 ENDDO 1195 1196 !$OMP DO 1197 DO i = nxl_mg(l)+1, nxr_mg(l), 2 1198 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 1199 !DIR$ IVDEP 1200 DO k = ind_even_odd+1, nzt_mg(l) 1201 km1 = k-ind_even_odd-1 1202 kp1 = k-ind_even_odd 1203 1204 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 1205 ddx2_mg(l) * & 1206 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 1207 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 1208 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 1209 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 897 898 i = ic+1 899 jj = jc+2-color 900 !DIR$ IVDEP 901 DO k = nzb+1, ind_even_odd 902 km1 = k+ind_even_odd 903 kp1 = k+ind_even_odd+1 904 j = jj 905 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 906 ddx2_mg(l) * & 907 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 1210 908 + ddy2_mg(l) * & 1211 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 1212 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 1213 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 1214 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 1215 + f2_mg(k,l) * p_mg(kp1,j,i) & 1216 + f3_mg(k,l) * & 1217 ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 1218 ( p_mg(k,j,i) - p_mg(km1,j,i) ) ) & 1219 - f_mg(k,j,i) ) 1220 ENDDO 909 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 910 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 911 + f3_mg_b(k,l) * p_mg(km1,j,i) & 912 - f_mg(k,j,i) ) 913 j = jj+2 914 p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( & 915 ddx2_mg(l) * & 916 ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) & 917 + ddy2_mg(l) * & 918 ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) & 919 + f2_mg_b(k,l) * p_mg(kp1,j,i) & 920 + f3_mg_b(k,l) * p_mg(km1,j,i) & 921 - f_mg(k,j,i) ) 1221 922 ENDDO 1222 ENDDO 1223 1224 !$OMP DO 1225 DO i = nxl_mg(l), nxr_mg(l), 2 1226 DO j = nys_mg(l) + (color-1), nyn_mg(l), 2 1227 !DIR$ IVDEP 1228 DO k = nzb+1, ind_even_odd 1229 km1 = k+ind_even_odd 1230 kp1 = k+ind_even_odd+1 1231 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 1232 ddx2_mg(l) * & 1233 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 1234 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 1235 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 1236 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 1237 + ddy2_mg(l) * & 1238 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 1239 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 1240 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 1241 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 1242 + f2_mg(k,l) * p_mg(kp1,j,i) & 1243 + f3_mg(k,l) * & 1244 ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 1245 ( p_mg(k,j,i) - p_mg(km1,j,i) ) ) & 1246 - f_mg(k,j,i) ) 1247 ENDDO 1248 ENDDO 1249 ENDDO 1250 1251 !$OMP DO 1252 DO i = nxl_mg(l)+1, nxr_mg(l), 2 1253 DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2 1254 !DIR$ IVDEP 1255 DO k = nzb+1, ind_even_odd 1256 km1 = k+ind_even_odd 1257 kp1 = k+ind_even_odd+1 1258 1259 p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( & 1260 ddx2_mg(l) * & 1261 ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * & 1262 ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + & 1263 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * & 1264 ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) & 1265 + ddy2_mg(l) * & 1266 ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * & 1267 ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + & 1268 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * & 1269 ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) & 1270 + f2_mg(k,l) * p_mg(kp1,j,i) & 1271 + f3_mg(k,l) * & 1272 ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * & 1273 ( p_mg(k,j,i) - p_mg(km1,j,i) ) ) & 1274 - f_mg(k,j,i) ) 1275 ENDDO 1276 ENDDO 1277 ENDDO 1278 !$OMP END PARALLEL 1279 1280 CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' ) 1281 1282 ELSE 1283 ! 1284 !-- Loop unrolling along y, only one i loop for better cache use 1285 CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' ) 1286 1287 !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj) 1288 !$OMP DO 1289 DO ic = nxl_mg(l), nxr_mg(l), 2 1290 DO jc = nys_mg(l), nyn_mg(l), 4 1291 i = ic 1292 jj = jc+2-color 1293 !DIR$ IVDEP 1294 DO k = ind_even_odd+1, nzt_mg(l) 1295 1296 km1 = k-ind_even_odd-1 1297 kp1 = k-ind_even_odd 1298 1299 j = jj 1300 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 1301 ddx2_mg(l) * & 1302 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) & 1303 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) & 1304 + ddy2_mg(l) * & 1305 ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) & 1306 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) & 1307 + f2_mg(k,l) * p_mg(kp1,j,i) & 1308 + f3_mg(k,l) * & 1309 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) & 1310 - f_mg(k,j,i) ) 1311 j = jj+2 1312 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 1313 ddx2_mg(l) * & 1314 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) & 1315 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) & 1316 + ddy2_mg(l) * & 1317 ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) & 1318 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) & 1319 + f2_mg(k,l) * p_mg(kp1,j,i) & 1320 + f3_mg(k,l) * & 1321 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) & 1322 - f_mg(k,j,i) ) 1323 1324 ENDDO 1325 1326 i = ic+1 1327 jj = jc+color-1 1328 !DIR$ IVDEP 1329 DO k = ind_even_odd+1, nzt_mg(l) 1330 1331 km1 = k-ind_even_odd-1 1332 kp1 = k-ind_even_odd 1333 1334 j =jj 1335 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 1336 ddx2_mg(l) * & 1337 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) & 1338 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) & 1339 + ddy2_mg(l) * & 1340 ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) & 1341 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) & 1342 + f2_mg(k,l) * p_mg(kp1,j,i) & 1343 + f3_mg(k,l) * & 1344 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) & 1345 - f_mg(k,j,i) ) 1346 1347 j = jj+2 1348 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 1349 ddx2_mg(l) * & 1350 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) & 1351 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) & 1352 + ddy2_mg(l) * & 1353 ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) & 1354 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) & 1355 + f2_mg(k,l) * p_mg(kp1,j,i) & 1356 + f3_mg(k,l) * & 1357 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) & 1358 - f_mg(k,j,i) ) 1359 1360 ENDDO 1361 1362 i = ic 1363 jj = jc+color-1 1364 !DIR$ IVDEP 1365 DO k = nzb+1, ind_even_odd 1366 1367 km1 = k+ind_even_odd 1368 kp1 = k+ind_even_odd+1 1369 1370 j =jj 1371 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 1372 ddx2_mg(l) * & 1373 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) & 1374 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) & 1375 + ddy2_mg(l) * & 1376 ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) & 1377 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) & 1378 + f2_mg(k,l) * p_mg(kp1,j,i) & 1379 + f3_mg(k,l) * & 1380 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) & 1381 - f_mg(k,j,i) ) 1382 1383 j = jj+2 1384 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 1385 ddx2_mg(l) * & 1386 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) & 1387 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) & 1388 + ddy2_mg(l) * & 1389 ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) & 1390 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) & 1391 + f2_mg(k,l) * p_mg(kp1,j,i) & 1392 + f3_mg(k,l) * & 1393 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) & 1394 - f_mg(k,j,i) ) 1395 1396 ENDDO 1397 1398 i = ic+1 1399 jj = jc+2-color 1400 !DIR$ IVDEP 1401 DO k = nzb+1, ind_even_odd 1402 1403 km1 = k+ind_even_odd 1404 kp1 = k+ind_even_odd+1 1405 1406 j =jj 1407 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 1408 ddx2_mg(l) * & 1409 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) & 1410 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) & 1411 + ddy2_mg(l) * & 1412 ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) & 1413 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) & 1414 + f2_mg(k,l) * p_mg(kp1,j,i) & 1415 + f3_mg(k,l) * & 1416 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) & 1417 - f_mg(k,j,i) ) 1418 1419 j = jj+2 1420 p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( & 1421 ddx2_mg(l) * & 1422 ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) & 1423 + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) & 1424 + ddy2_mg(l) * & 1425 ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) & 1426 + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) & 1427 + f2_mg(k,l) * p_mg(kp1,j,i) & 1428 + f3_mg(k,l) * & 1429 MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) & 1430 - f_mg(k,j,i) ) 1431 1432 ENDDO 1433 1434 ENDDO 1435 ENDDO 1436 !$OMP END PARALLEL 1437 1438 CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' ) 1439 1440 ENDIF 923 924 ENDDO 925 ENDDO 926 !$OMP END PARALLEL 927 928 CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' ) 1441 929 1442 930 ENDIF … … 1467 955 !-- Bottom and top boundary conditions 1468 956 IF ( ibc_p_b == 1 ) THEN 957 ! 958 !-- equivalent to p_mg(nzb,:,: ) = p_mg(nzb+1,:,:) 1469 959 p_mg(nzb,:,: ) = p_mg(ind_even_odd+1,:,:) 1470 960 ELSE … … 1473 963 1474 964 IF ( ibc_p_t == 1 ) THEN 965 ! 966 !-- equivalent to p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:) 1475 967 p_mg(nzt_mg(l)+1,:,: ) = p_mg(ind_even_odd,:,:) 1476 968 ELSE … … 1479 971 1480 972 ENDDO 973 1481 974 ENDDO 1482 !1483 !-- Set pressure within topography and at the topography surfaces1484 IF ( topography /= 'flat' .AND. .NOT. masking_method ) THEN1485 1486 !$OMP PARALLEL PRIVATE (i,j,k,kp1,wall_left,wall_north,wall_right, &1487 !$OMP wall_south,wall_top,wall_total)1488 !$OMP DO1489 DO i = nxl_mg(l), nxr_mg(l)1490 DO j = nys_mg(l), nyn_mg(l)1491 !DIR$ IVDEP1492 DO k = ind_even_odd+1, nzt_mg(l)1493 km1 = k-ind_even_odd-11494 kp1 = k-ind_even_odd1495 !1496 !-- First, set pressure inside topography to zero1497 p_mg(k,j,i) = p_mg(k,j,i) * &1498 ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )1499 !1500 !-- Second, determine if the gridpoint inside topography is1501 !-- adjacent to a wall and set its value to a value given by the1502 !-- average of those values obtained from Neumann boundary1503 !-- condition1504 wall_left = IBITS( flags(k,j,i-1), 5, 1 )1505 wall_right = IBITS( flags(k,j,i+1), 4, 1 )1506 wall_south = IBITS( flags(k,j-1,i), 3, 1 )1507 wall_north = IBITS( flags(k,j+1,i), 2, 1 )1508 wall_top = IBITS( flags(kp1,j,i), 0, 1 )1509 wall_total = wall_left + wall_right + wall_south + &1510 wall_north + wall_top1511 1512 IF ( wall_total > 0.0_wp ) THEN1513 p_mg(k,j,i) = 1.0_wp / wall_total * &1514 ( wall_left * p_mg(k,j,i-1) + &1515 wall_right * p_mg(k,j,i+1) + &1516 wall_south * p_mg(k,j-1,i) + &1517 wall_north * p_mg(k,j+1,i) + &1518 wall_top * p_mg(kp1,j,i) )1519 ENDIF1520 ENDDO1521 1522 !DIR$ IVDEP1523 DO k = nzb+1, ind_even_odd1524 km1 = k+ind_even_odd1525 kp1 = k+ind_even_odd+11526 !1527 !-- First, set pressure inside topography to zero1528 p_mg(k,j,i) = p_mg(k,j,i) * &1529 ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )1530 !1531 !-- Second, determine if the gridpoint inside topography is1532 !-- adjacent to a wall and set its value to a value given by the1533 !-- average of those values obtained from Neumann boundary1534 !-- condition1535 wall_left = IBITS( flags(k,j,i-1), 5, 1 )1536 wall_right = IBITS( flags(k,j,i+1), 4, 1 )1537 wall_south = IBITS( flags(k,j-1,i), 3, 1 )1538 wall_north = IBITS( flags(k,j+1,i), 2, 1 )1539 wall_top = IBITS( flags(kp1,j,i), 0, 1 )1540 wall_total = wall_left + wall_right + wall_south + &1541 wall_north + wall_top1542 1543 IF ( wall_total > 0.0_wp ) THEN1544 p_mg(k,j,i) = 1.0_wp / wall_total * &1545 ( wall_left * p_mg(k,j,i-1) + &1546 wall_right * p_mg(k,j,i+1) + &1547 wall_south * p_mg(k,j-1,i) + &1548 wall_north * p_mg(k,j+1,i) + &1549 wall_top * p_mg(kp1,j,i) )1550 ENDIF1551 ENDDO1552 ENDDO1553 ENDDO1554 !$OMP END PARALLEL1555 1556 !1557 !-- One more time horizontal boundary conditions1558 CALL exchange_horiz( p_mg, 1)1559 1560 ENDIF1561 975 1562 976 END SUBROUTINE redblack_fast … … 1581 995 IMPLICIT NONE 1582 996 1583 INTEGER(iwp), INTENT(IN) :: glevel 997 INTEGER(iwp), INTENT(IN) :: glevel !< grid level 1584 998 1585 999 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1, & 1586 1000 nys_mg(glevel)-1:nyn_mg(glevel)+1, & 1587 nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: p_mg !< 1001 nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: & 1002 p_mg !< array to be sorted 1588 1003 ! 1589 1004 !-- Local variables 1590 INTEGER(iwp) :: i !< 1591 INTEGER(iwp) :: j !< 1592 INTEGER(iwp) :: k !< 1593 INTEGER(iwp) :: l !< 1594 INTEGER(iwp) :: ind !< 1595 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< 1005 INTEGER(iwp) :: i !< index variable along x 1006 INTEGER(iwp) :: j !< index variable along y 1007 INTEGER(iwp) :: k !< index variable along z 1008 INTEGER(iwp) :: l !< grid level 1009 INTEGER(iwp) :: ind !< index variable along z 1010 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< odd-even sorted temporary array 1596 1011 1597 1012 … … 1646 1061 IMPLICIT NONE 1647 1062 1648 INTEGER(iwp), INTENT(IN) :: glevel 1649 1650 REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) :: f_mg 1651 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: f_mg_b 1063 INTEGER(iwp), INTENT(IN) :: glevel !< grid level 1064 1065 REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) :: f_mg !< 1D input array 1066 REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: f_mg_b !< 1D output array 1652 1067 1653 1068 ! 1654 1069 !-- Local variables 1655 INTEGER(iwp) :: ind !< 1656 INTEGER(iwp) :: k !< 1070 INTEGER(iwp) :: ind !< index variable along z 1071 INTEGER(iwp) :: k !< index variable along z 1657 1072 1658 1073 … … 1696 1111 IMPLICIT NONE 1697 1112 1698 INTEGER(iwp), INTENT(IN) :: glevel 1113 INTEGER(iwp), INTENT(IN) :: glevel !< grid level 1699 1114 1700 1115 INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1, & 1701 1116 nys_mg(glevel)-1:nyn_mg(glevel)+1, & 1702 nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: i_mg !< 1117 nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: & 1118 i_mg !< array to be sorted 1703 1119 ! 1704 1120 !-- Local variables 1705 INTEGER(iwp) :: i !< 1706 INTEGER(iwp) :: j !< 1707 INTEGER(iwp) :: k !< 1708 INTEGER(iwp) :: l !< 1709 INTEGER(iwp) :: ind !< 1710 INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp 1121 INTEGER(iwp) :: i !< index variabel along x 1122 INTEGER(iwp) :: j !< index variable along y 1123 INTEGER(iwp) :: k !< index variable along z 1124 INTEGER(iwp) :: l !< grid level 1125 INTEGER(iwp) :: ind !< index variable along z 1126 INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !< temporary odd-even sorted array 1711 1127 1712 1128 … … 1770 1186 REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, & 1771 1187 nys_mg(grid_level)-1:nyn_mg(grid_level)+1, & 1772 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !< 1188 nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: & 1189 p_mg !< array to be sorted 1773 1190 ! 1774 1191 !-- Local variables 1775 INTEGER(iwp) :: i !< 1776 INTEGER(iwp) :: j !< 1777 INTEGER(iwp) :: k !< 1778 INTEGER(iwp) :: l !< 1779 INTEGER(iwp) :: ind !< 1192 INTEGER(iwp) :: i !< index variable along x 1193 INTEGER(iwp) :: j !< index variable along y 1194 INTEGER(iwp) :: k !< index variable along z 1195 INTEGER(iwp) :: l !< grid level 1196 INTEGER(iwp) :: ind !< index variable along z 1780 1197 1781 1198 REAL(wp),DIMENSION(nzb:nzt_mg(grid_level)+1) :: tmp … … 1958 1375 IMPLICIT NONE 1959 1376 1960 INTEGER(iwp) :: i !< 1961 INTEGER(iwp) :: j !< 1962 INTEGER(iwp) :: k !< 1377 INTEGER(iwp) :: i !< index variable along x 1378 INTEGER(iwp) :: j !< index variable along y 1379 INTEGER(iwp) :: k !< index variable along z 1963 1380 INTEGER(iwp) :: nxl_mg_save !< 1964 1381 INTEGER(iwp) :: nxr_mg_save !< … … 2308 1725 2309 1726 USE control_parameters, & 2310 ONLY: grid_level, ma sking_method, maximum_grid_level, topography1727 ONLY: grid_level, maximum_grid_level 2311 1728 2312 1729 USE indices, & … … 2314 1731 2315 1732 USE indices, & 2316 ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, & 2317 wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7, & 2318 wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, & 2319 nys_mg, nyn_mg, nzb, nzt_mg 1733 ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg 2320 1734 2321 1735 IMPLICIT NONE 2322 1736 ! 2323 1737 !-- Local variables 2324 INTEGER(iwp) :: i !< 1738 INTEGER(iwp) :: i !< 2325 1739 INTEGER(iwp) :: l !< 2326 1740 … … 2357 1771 ENDDO 2358 1772 2359 !2360 !-- Sort flags array for all levels to block (even/odd) structure2361 IF ( topography /= 'flat' .AND. .NOT. masking_method ) THEN2362 2363 DO l = maximum_grid_level, 1 , -12364 !2365 !-- Assign the flag level to be calculated2366 SELECT CASE ( l )2367 CASE ( 1 )2368 flags => wall_flags_12369 CASE ( 2 )2370 flags => wall_flags_22371 CASE ( 3 )2372 flags => wall_flags_32373 CASE ( 4 )2374 flags => wall_flags_42375 CASE ( 5 )2376 flags => wall_flags_52377 CASE ( 6 )2378 flags => wall_flags_62379 CASE ( 7 )2380 flags => wall_flags_72381 CASE ( 8 )2382 flags => wall_flags_82383 CASE ( 9 )2384 flags => wall_flags_92385 CASE ( 10 )2386 flags => wall_flags_102387 END SELECT2388 2389 CALL sort_k_to_even_odd_blocks( flags, l )2390 2391 ENDDO2392 2393 ENDIF2394 2395 1773 lfirst = .FALSE. 2396 1774
Note: See TracChangeset
for help on using the changeset viewer.
