Changeset 1942
- Timestamp:
- Jun 14, 2016 12:18:18 PM (8 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
advec_ws.f90 (modified) (3 diffs)
-
init_grid.f90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r1877 r1942 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Initialization of flags for ws-scheme moved from init_grid. 22 22 ! 23 23 ! Former revisions: … … 167 167 PUBLIC advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc, & 168 168 advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc, & 169 ws_init, ws_ statistics169 ws_init, ws_init_flags, ws_statistics 170 170 171 171 INTERFACE ws_init 172 172 MODULE PROCEDURE ws_init 173 173 END INTERFACE ws_init 174 175 INTERFACE ws_init_flags 176 MODULE PROCEDURE ws_init_flags 177 END INTERFACE ws_init_flags 174 178 175 179 INTERFACE ws_statistics … … 364 368 365 369 END SUBROUTINE ws_init 370 371 !------------------------------------------------------------------------------! 372 ! Description: 373 ! ------------ 374 !> Initialization of flags for WS-scheme used to degrade the order of the scheme 375 !> near walls. 376 !------------------------------------------------------------------------------! 377 SUBROUTINE ws_init_flags 378 379 USE control_parameters, & 380 ONLY: inflow_l, inflow_n, inflow_r, inflow_s, momentum_advec, & 381 nest_bound_l, nest_bound_n, nest_bound_r, nest_bound_s, & 382 outflow_l, outflow_n, outflow_r, outflow_s, scalar_advec 383 384 USE indices, & 385 ONLY: nbgp, nxl, nxlu, nxr, nyn, nys, nysv, nzb, nzb_s_inner, & 386 nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt, wall_flags_0, & 387 wall_flags_00 388 389 USE kinds 390 391 IMPLICIT NONE 392 393 INTEGER(iwp) :: i !< index variable along x 394 INTEGER(iwp) :: j !< index variable along y 395 INTEGER(iwp) :: k !< index variable along z 396 397 LOGICAL :: flag_set !< steering variable for advection flags 398 399 400 IF ( scalar_advec == 'ws-scheme' .OR. & 401 scalar_advec == 'ws-scheme-mono' ) THEN 402 ! 403 !-- Set flags to steer the degradation of the advection scheme in advec_ws 404 !-- near topography, inflow- and outflow boundaries as well as bottom and 405 !-- top of model domain. wall_flags_0 remains zero for all non-prognostic 406 !-- grid points. 407 DO i = nxl, nxr 408 DO j = nys, nyn 409 DO k = nzb_s_inner(j,i)+1, nzt 410 ! 411 !-- scalar - x-direction 412 !-- WS1 (0), WS3 (1), WS5 (2) 413 IF ( ( k <= nzb_s_inner(j,i+1) .OR. k <= nzb_s_inner(j,i+2) & 414 .OR. k <= nzb_s_inner(j,i-1) ) & 415 .OR. ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 416 .AND. i == nxl ) .OR. & 417 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 418 .AND. i == nxr ) ) & 419 THEN 420 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 ) 421 ELSEIF ( ( k <= nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i+1)& 422 .AND. k > nzb_s_inner(j,i+2)& 423 .AND. k > nzb_s_inner(j,i-1)& 424 ) .OR. & 425 ( k <= nzb_s_inner(j,i-2) .AND. k > nzb_s_inner(j,i+1)& 426 .AND. k > nzb_s_inner(j,i+2)& 427 .AND. k > nzb_s_inner(j,i-1)& 428 ) & 429 .OR. & 430 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 431 .AND. i == nxr-1 ) .OR. & 432 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 433 .AND. i == nxlu ) ) & 434 THEN 435 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 ) 436 ELSEIF ( k > nzb_s_inner(j,i+1) .AND. k > nzb_s_inner(j,i+2)& 437 .AND. k > nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i-1)& 438 .AND. k > nzb_s_inner(j,i-2) ) & 439 THEN 440 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 ) 441 ENDIF 442 ! 443 !-- scalar - y-direction 444 !-- WS1 (3), WS3 (4), WS5 (5) 445 IF ( ( k <= nzb_s_inner(j+1,i) .OR. k <= nzb_s_inner(j+2,i) & 446 .OR. k <= nzb_s_inner(j-1,i) ) & 447 .OR. & 448 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 449 .AND. j == nys ) .OR. & 450 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 451 .AND. j == nyn ) ) & 452 THEN 453 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 ) 454 ! 455 !-- WS3 456 ELSEIF ( ( k <= nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j+1,i)& 457 .AND. k > nzb_s_inner(j+2,i)& 458 .AND. k > nzb_s_inner(j-1,i)& 459 ) .OR. & 460 ( k <= nzb_s_inner(j-2,i) .AND. k > nzb_s_inner(j+1,i)& 461 .AND. k > nzb_s_inner(j+2,i)& 462 .AND. k > nzb_s_inner(j-1,i)& 463 ) & 464 .OR. & 465 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 466 .AND. j == nysv ) .OR. & 467 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 468 .AND. j == nyn-1 ) ) & 469 THEN 470 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 ) 471 ! 472 !-- WS5 473 ELSEIF ( k > nzb_s_inner(j+1,i) .AND. k > nzb_s_inner(j+2,i)& 474 .AND. k > nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j-1,i)& 475 .AND. k > nzb_s_inner(j-2,i) ) & 476 THEN 477 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 ) 478 ENDIF 479 ! 480 !-- scalar - z-direction 481 !-- WS1 (6), WS3 (7), WS5 (8) 482 flag_set = .FALSE. 483 IF ( k == nzb_s_inner(j,i) + 1 .OR. k == nzt ) THEN 484 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 ) 485 flag_set = .TRUE. 486 ELSEIF ( k == nzb_s_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 487 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 7 ) 488 flag_set = .TRUE. 489 ELSEIF ( k > nzb_s_inner(j,i) .AND. .NOT. flag_set ) THEN 490 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 ) 491 ENDIF 492 493 ENDDO 494 ENDDO 495 ENDDO 496 ENDIF 497 498 IF ( momentum_advec == 'ws-scheme' ) THEN 499 ! 500 !-- Set wall_flags_0 to steer the degradation of the advection scheme in advec_ws 501 !-- near topography, inflow- and outflow boundaries as well as bottom and 502 !-- top of model domain. wall_flags_0 remains zero for all non-prognostic 503 !-- grid points. 504 DO i = nxl, nxr 505 DO j = nys, nyn 506 DO k = nzb+1, nzt 507 ! 508 !-- At first, set flags to WS1. 509 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 510 !-- in order to handle the left/south flux. 511 !-- near vertical walls. 512 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 ) 513 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 ) 514 ! 515 !-- u component - x-direction 516 !-- WS1 (9), WS3 (10), WS5 (11) 517 IF ( k <= nzb_u_inner(j,i+1) .OR. & 518 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 519 .AND. i <= nxlu ) .OR. & 520 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 521 .AND. i == nxr ) ) & 522 THEN 523 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 ) 524 ELSEIF ( ( k <= nzb_u_inner(j,i+2) .AND. & 525 k > nzb_u_inner(j,i+1) ) .OR. & 526 k <= nzb_u_inner(j,i-1) & 527 .OR. & 528 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 529 .AND. i == nxr-1 ) .OR. & 530 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 531 .AND. i == nxlu+1) ) & 532 THEN 533 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 ) 534 ! 535 !-- Clear flag for WS1 536 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 ) 537 ELSEIF ( k > nzb_u_inner(j,i+1) .AND. k > nzb_u_inner(j,i+2) & 538 .AND. k > nzb_u_inner(j,i-1) ) & 539 THEN 540 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 ) 541 ! 542 !-- Clear flag for WS1 543 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 ) 544 ENDIF 545 ! 546 !-- u component - y-direction 547 !-- WS1 (12), WS3 (13), WS5 (14) 548 IF ( k <= nzb_u_inner(j+1,i) .OR. & 549 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 550 .AND. j == nys ) .OR. & 551 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 552 .AND. j == nyn ) ) & 553 THEN 554 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 ) 555 ELSEIF ( ( k <= nzb_u_inner(j+2,i) .AND. & 556 k > nzb_u_inner(j+1,i) ) .OR. & 557 k <= nzb_u_inner(j-1,i) & 558 .OR. & 559 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 560 .AND. j == nysv ) .OR. & 561 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 562 .AND. j == nyn-1 ) ) & 563 THEN 564 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 ) 565 ! 566 !-- Clear flag for WS1 567 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 ) 568 ELSEIF ( k > nzb_u_inner(j+1,i) .AND. k > nzb_u_inner(j+2,i) & 569 .AND. k > nzb_u_inner(j-1,i) ) & 570 THEN 571 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 ) 572 ! 573 !-- Clear flag for WS1 574 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 ) 575 ENDIF 576 ! 577 !-- u component - z-direction 578 !-- WS1 (15), WS3 (16), WS5 (17) 579 flag_set = .FALSE. 580 IF ( k == nzb_u_inner(j,i) + 1 .OR. k == nzt ) THEN 581 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 ) 582 flag_set = .TRUE. 583 ELSEIF ( k == nzb_u_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 584 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 ) 585 flag_set = .TRUE. 586 ELSEIF ( k > nzb_u_inner(j,i) .AND. .NOT. flag_set ) THEN 587 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 ) 588 ENDIF 589 590 ENDDO 591 ENDDO 592 ENDDO 593 594 DO i = nxl, nxr 595 DO j = nys, nyn 596 DO k = nzb+1, nzt 597 ! 598 !-- At first, set flags to WS1. 599 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 600 !-- in order to handle the left/south flux. 601 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 ) 602 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 ) 603 ! 604 !-- v component - x-direction 605 !-- WS1 (18), WS3 (19), WS5 (20) 606 IF ( k <= nzb_v_inner(j,i+1) .OR. & 607 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 608 .AND. i == nxl ) .OR. & 609 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 610 .AND. i == nxr ) ) & 611 THEN 612 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 ) 613 ! 614 !-- WS3 615 ELSEIF ( ( k <= nzb_v_inner(j,i+2) .AND. & 616 k > nzb_v_inner(j,i+1) ) .OR. & 617 k <= nzb_v_inner(j,i-1) & 618 .OR. & 619 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 620 .AND. i == nxr-1 ) .OR. & 621 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 622 .AND. i == nxlu ) ) & 623 THEN 624 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 ) 625 ! 626 !-- Clear flag for WS1 627 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 ) 628 ELSEIF ( k > nzb_v_inner(j,i+1) .AND. k > nzb_v_inner(j,i+2) & 629 .AND. k > nzb_v_inner(j,i-1) ) & 630 THEN 631 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 ) 632 ! 633 !-- Clear flag for WS1 634 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 ) 635 ENDIF 636 ! 637 !-- v component - y-direction 638 !-- WS1 (21), WS3 (22), WS5 (23) 639 IF ( k <= nzb_v_inner(j+1,i) .OR. & 640 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 641 .AND. j <= nysv ) .OR. & 642 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 643 .AND. j == nyn ) ) & 644 THEN 645 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 ) 646 ELSEIF ( ( k <= nzb_v_inner(j+2,i) .AND. & 647 k > nzb_v_inner(j+1,i) ) .OR. & 648 k <= nzb_v_inner(j-1,i) & 649 .OR. & 650 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 651 .AND. j == nysv+1) .OR. & 652 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 653 .AND. j == nyn-1 ) ) & 654 THEN 655 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 ) 656 ! 657 !-- Clear flag for WS1 658 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 ) 659 ELSEIF ( k > nzb_v_inner(j+1,i) .AND. k > nzb_v_inner(j+2,i) & 660 .AND. k > nzb_v_inner(j-1,i) ) & 661 THEN 662 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 ) 663 ! 664 !-- Clear flag for WS1 665 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 ) 666 ENDIF 667 ! 668 !-- v component - z-direction 669 !-- WS1 (24), WS3 (25), WS5 (26) 670 flag_set = .FALSE. 671 IF ( k == nzb_v_inner(j,i) + 1 .OR. k == nzt ) THEN 672 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 ) 673 flag_set = .TRUE. 674 ELSEIF ( k == nzb_v_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 675 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 ) 676 flag_set = .TRUE. 677 ELSEIF ( k > nzb_v_inner(j,i) .AND. .NOT. flag_set ) THEN 678 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 ) 679 ENDIF 680 681 ENDDO 682 ENDDO 683 ENDDO 684 DO i = nxl, nxr 685 DO j = nys, nyn 686 DO k = nzb+1, nzt 687 ! 688 !-- At first, set flags to WS1. 689 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 690 !-- in order to handle the left/south flux. 691 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 ) 692 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 ) 693 ! 694 !-- w component - x-direction 695 !-- WS1 (27), WS3 (28), WS5 (29) 696 IF ( k <= nzb_w_inner(j,i+1) .OR. & 697 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 698 .AND. i == nxl ) .OR. & 699 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 700 .AND. i == nxr ) ) & 701 THEN 702 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 ) 703 ELSEIF ( ( k <= nzb_w_inner(j,i+2) .AND. & 704 k > nzb_w_inner(j,i+1) ) .OR. & 705 k <= nzb_w_inner(j,i-1) & 706 .OR. & 707 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 708 .AND. i == nxr-1 ) .OR. & 709 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 710 .AND. i == nxlu ) ) & 711 THEN 712 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 ) 713 ! 714 !-- Clear flag for WS1 715 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 ) 716 ELSEIF ( k > nzb_w_inner(j,i+1) .AND. k > nzb_w_inner(j,i+2) & 717 .AND. k > nzb_w_inner(j,i-1) ) & 718 THEN 719 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i),29 ) 720 ! 721 !-- Clear flag for WS1 722 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 ) 723 ENDIF 724 ! 725 !-- w component - y-direction 726 !-- WS1 (30), WS3 (31), WS5 (32) 727 IF ( k <= nzb_w_inner(j+1,i) .OR. & 728 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 729 .AND. j == nys ) .OR. & 730 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 731 .AND. j == nyn ) ) & 732 THEN 733 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 ) 734 ELSEIF ( ( k <= nzb_w_inner(j+2,i) .AND. & 735 k > nzb_w_inner(j+1,i) ) .OR. & 736 k <= nzb_w_inner(j-1,i) & 737 .OR. & 738 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 739 .AND. j == nysv ) .OR. & 740 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 741 .AND. j == nyn-1 ) ) & 742 THEN 743 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 ) 744 ! 745 !-- Clear flag for WS1 746 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 ) 747 ELSEIF ( k > nzb_w_inner(j+1,i) .AND. k > nzb_w_inner(j+2,i) & 748 .AND. k > nzb_w_inner(j-1,i) ) & 749 THEN 750 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 0 ) 751 ! 752 !-- Clear flag for WS1 753 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 ) 754 ENDIF 755 ! 756 !-- w component - z-direction 757 !-- WS1 (33), WS3 (34), WS5 (35) 758 flag_set = .FALSE. 759 IF ( k == nzb_w_inner(j,i) .OR. k == nzb_w_inner(j,i) + 1 & 760 .OR. k == nzt ) THEN 761 ! 762 !-- Please note, at k == nzb_w_inner(j,i) a flag is explictely 763 !-- set, although this is not a prognostic level. However, 764 !-- contrary to the advection of u,v and s this is necessary 765 !-- because flux_t(nzb_w_inner(j,i)) is used for the tendency 766 !-- at k == nzb_w_inner(j,i)+1. 767 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 1 ) 768 flag_set = .TRUE. 769 ELSEIF ( k == nzb_w_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 770 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 2 ) 771 flag_set = .TRUE. 772 ELSEIF ( k > nzb_w_inner(j,i) .AND. .NOT. flag_set ) THEN 773 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 3 ) 774 ENDIF 775 776 ENDDO 777 ENDDO 778 ENDDO 779 780 ENDIF 781 782 783 ! 784 !-- Exchange 3D integer wall_flags. 785 IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' & 786 .OR. scalar_advec == 'ws-scheme-mono' ) THEN 787 ! 788 !-- Exchange ghost points for advection flags 789 CALL exchange_horiz_int( wall_flags_0, nbgp ) 790 CALL exchange_horiz_int( wall_flags_00, nbgp ) 791 ! 792 !-- Set boundary flags at inflow and outflow boundary in case of 793 !-- non-cyclic boundary conditions. 794 IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN 795 wall_flags_0(:,:,nxl-1) = wall_flags_0(:,:,nxl) 796 wall_flags_00(:,:,nxl-1) = wall_flags_00(:,:,nxl) 797 ENDIF 798 799 IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN 800 wall_flags_0(:,:,nxr+1) = wall_flags_0(:,:,nxr) 801 wall_flags_00(:,:,nxr+1) = wall_flags_00(:,:,nxr) 802 ENDIF 803 804 IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN 805 wall_flags_0(:,nyn+1,:) = wall_flags_0(:,nyn,:) 806 wall_flags_00(:,nyn+1,:) = wall_flags_00(:,nyn,:) 807 ENDIF 808 809 IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN 810 wall_flags_0(:,nys-1,:) = wall_flags_0(:,nys,:) 811 wall_flags_00(:,nys-1,:) = wall_flags_00(:,nys,:) 812 ENDIF 813 814 ENDIF 815 816 817 END SUBROUTINE ws_init_flags 366 818 367 819 -
palm/trunk/SOURCE/init_grid.f90
r1932 r1942 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Topography filter implemented to fill holes resolved by only one grid point. 22 ! Initialization of flags for ws-scheme moved to advec_ws. 22 23 ! 23 24 ! Former revisions: … … 171 172 SUBROUTINE init_grid 172 173 174 USE advec_ws, & 175 ONLY: ws_init_flags 173 176 174 177 USE arrays_3d, & … … 185 188 io_group, inflow_l, inflow_n, inflow_r, inflow_s, & 186 189 masking_method, maximum_grid_level, message_string, & 187 momentum_advec, nest_domain, nest_bound_l, nest_bound_n, & 188 nest_bound_r, nest_bound_s, ocean, outflow_l, outflow_n, & 190 momentum_advec, nest_domain, ocean, outflow_l, outflow_n, & 189 191 outflow_r, outflow_s, psolver, scalar_advec, topography, & 190 192 topography_grid_convention, use_surface_fluxes, use_top_fluxes, & … … 197 199 198 200 USE indices, & 199 ONLY: flags, nbgp, nx, nxl, nxlg, nxl u, nxl_mg, nxr, nxrg, nxr_mg,&200 ny, nyn, nyng, nyn_mg, nys, nys v, nys_mg, nysg, nz, nzb,&201 ONLY: flags, nbgp, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, & 202 ny, nyn, nyng, nyn_mg, nys, nys_mg, nysg, nz, nzb, & 201 203 nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u, & 202 204 nzb_diff_v, nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner, & … … 214 216 IMPLICIT NONE 215 217 216 INTEGER(iwp) :: bh !< temporary vertical index of building height 217 INTEGER(iwp) :: blx !< grid point number of building size along x 218 INTEGER(iwp) :: bly !< grid point number of building size along y 219 INTEGER(iwp) :: bxl !< index for left building wall 220 INTEGER(iwp) :: bxr !< index for right building wall 221 INTEGER(iwp) :: byn !< index for north building wall 222 INTEGER(iwp) :: bys !< index for south building wall 223 INTEGER(iwp) :: ch !< temporary vertical index for canyon height 224 INTEGER(iwp) :: cwx !< grid point number of canyon size along x 225 INTEGER(iwp) :: cwy !< grid point number of canyon size along y 226 INTEGER(iwp) :: cxl !< index for left canyon wall 227 INTEGER(iwp) :: cxr !< index for right canyon wall 228 INTEGER(iwp) :: cyn !< index for north canyon wall 229 INTEGER(iwp) :: cys !< index for south canyon wall 230 INTEGER(iwp) :: gls !< number of lateral ghost points at total model domain boundaries required for multigrid solver 231 INTEGER(iwp) :: i !< index variable along x 232 INTEGER(iwp) :: ii !< loop variable for reading topography file 233 INTEGER(iwp) :: inc !< incremental parameter for coarsening grid level 234 INTEGER(iwp) :: j !< index variable along y 235 INTEGER(iwp) :: k !< index variable along z 236 INTEGER(iwp) :: l !< loop variable 237 INTEGER(iwp) :: nxl_l !< index of left PE boundary for multigrid level 238 INTEGER(iwp) :: nxr_l !< index of right PE boundary for multigrid level 239 INTEGER(iwp) :: nyn_l !< index of north PE boundary for multigrid level 240 INTEGER(iwp) :: nys_l !< index of south PE boundary for multigrid level 241 INTEGER(iwp) :: nzb_si !< dummy index for local nzb_s_inner 242 INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level 243 INTEGER(iwp) :: vi !< dummy for vertical influence 218 INTEGER(iwp) :: bh !< temporary vertical index of building height 219 INTEGER(iwp) :: blx !< grid point number of building size along x 220 INTEGER(iwp) :: bly !< grid point number of building size along y 221 INTEGER(iwp) :: bxl !< index for left building wall 222 INTEGER(iwp) :: bxr !< index for right building wall 223 INTEGER(iwp) :: byn !< index for north building wall 224 INTEGER(iwp) :: bys !< index for south building wall 225 INTEGER(iwp) :: ch !< temporary vertical index for canyon height 226 INTEGER(iwp) :: cwx !< grid point number of canyon size along x 227 INTEGER(iwp) :: cwy !< grid point number of canyon size along y 228 INTEGER(iwp) :: cxl !< index for left canyon wall 229 INTEGER(iwp) :: cxr !< index for right canyon wall 230 INTEGER(iwp) :: cyn !< index for north canyon wall 231 INTEGER(iwp) :: cys !< index for south canyon wall 232 INTEGER(iwp) :: gls !< number of lateral ghost points at total model domain boundaries required for multigrid solver 233 INTEGER(iwp) :: i !< index variable along x 234 INTEGER(iwp) :: ii !< loop variable for reading topography file 235 INTEGER(iwp) :: inc !< incremental parameter for coarsening grid level 236 INTEGER(iwp) :: j !< index variable along y 237 INTEGER(iwp) :: k !< index variable along z 238 INTEGER(iwp) :: l !< loop variable 239 INTEGER(iwp) :: nxl_l !< index of left PE boundary for multigrid level 240 INTEGER(iwp) :: nxr_l !< index of right PE boundary for multigrid level 241 INTEGER(iwp) :: nyn_l !< index of north PE boundary for multigrid level 242 INTEGER(iwp) :: nys_l !< index of south PE boundary for multigrid level 243 INTEGER(iwp) :: nzb_si !< dummy index for local nzb_s_inner 244 INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level 245 INTEGER(iwp) :: num_hole !< number of holes (in topography) resolved by only one grid point 246 INTEGER(iwp) :: num_wall !< number of surrounding vertical walls for a single grid point 247 INTEGER(iwp) :: vi !< dummy for vertical influence 244 248 245 249 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: & … … 259 263 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp !< dummy to calculate topography indices on u- and v-grid 260 264 261 LOGICAL :: flag_set = .FALSE. !< steering variable for advection flags265 LOGICAL :: hole = .FALSE. !< flag to check if any holes resolved by only 1 grid point were filled 262 266 263 267 REAL(wp) :: dx_l !< grid spacing along x on different multigrid level … … 744 748 745 749 DEALLOCATE ( topo_height ) 746 750 ! 751 !-- Filter topography, i.e. fill holes resolved by only one grid point. 752 !-- Such holes are suspected to lead to velocity blow-ups as continuity 753 !-- equation on discrete grid cannot be fulfilled in such case. 754 !-- For now, check only for holes and fill them to the lowest height level 755 !-- of the directly adjoining grid points along x- and y- direction. 756 !-- Before checking for holes, set lateral boundary conditions for 757 !-- topography. After hole-filling, boundary conditions must be set again! 758 IF ( bc_ns_cyc ) THEN 759 nzb_local(-1,:) = nzb_local(ny,:) 760 nzb_local(ny+1,:) = nzb_local(0,:) 761 ELSE 762 nzb_local(-1,:) = nzb_local(0,:) 763 nzb_local(ny+1,:) = nzb_local(ny,:) 764 ENDIF 765 766 IF ( bc_lr_cyc ) THEN 767 nzb_local(:,-1) = nzb_local(:,nx) 768 nzb_local(:,nx+1) = nzb_local(:,0) 769 ELSE 770 nzb_local(:,-1) = nzb_local(:,0) 771 nzb_local(:,nx+1) = nzb_local(:,nx) 772 ENDIF 773 774 num_hole = 0 775 DO i = 0, nx 776 DO j = 0, ny 777 778 num_wall = 0 779 780 IF ( nzb_local(j-1,i) > nzb_local(j,i) ) & 781 num_wall = num_wall + 1 782 IF ( nzb_local(j+1,i) > nzb_local(j,i) ) & 783 num_wall = num_wall + 1 784 IF ( nzb_local(j,i-1) > nzb_local(j,i) ) & 785 num_wall = num_wall + 1 786 IF ( nzb_local(j,i+1) > nzb_local(j,i) ) & 787 num_wall = num_wall + 1 788 789 IF ( num_wall == 4 ) THEN 790 hole = .TRUE. 791 nzb_local(j,i) = MIN( nzb_local(j-1,i), nzb_local(j+1,i), & 792 nzb_local(j,i-1), nzb_local(j,i+1) ) 793 num_hole = num_hole + 1 794 ENDIF 795 ENDDO 796 ENDDO 797 ! 798 !-- Create an informative message if any hole was removed. 799 IF ( hole ) THEN 800 WRITE( message_string, * ) num_hole, 'hole(s) resolved by only '//& 801 'one grid point were filled' 802 CALL message( 'init_grid', 'PA0430', 0, 0, 0, 6, 0 ) 803 ENDIF 747 804 ! 748 805 !-- Add cyclic or Neumann boundary conditions (additional layers are for … … 1237 1294 ENDIF 1238 1295 1239 !1240 !-- Test output of flag arrays1241 ! i = nxl_l1242 ! WRITE (9,*) ' '1243 ! WRITE (9,*) '*** mg level ', l, ' ***', mg_switch_to_pe0_level1244 ! WRITE (9,*) ' inc=', inc, ' i =', nxl_l1245 ! WRITE (9,*) ' nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l1246 ! DO k = nzt_l+1, nzb, -11247 ! WRITE (9,'(194(1X,I2))') ( flags(k,j,i), j = nys_l-1, nyn_l+1 )1248 ! ENDDO1249 1250 1296 inc = inc * 2 1251 1297 … … 1254 1300 ENDIF 1255 1301 ! 1256 !-- Allocate flags needed for masking walls. 1302 !-- Allocate flags needed for masking walls. Even though these flags are only 1303 !-- required in the ws-scheme, the arrays need to be allocated as they are 1304 !-- used in OpenACC directives. 1257 1305 ALLOCATE( wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 1258 1306 wall_flags_00(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1259 1307 wall_flags_0 = 0 1260 1308 wall_flags_00 = 0 1261 1262 IF ( scalar_advec == 'ws-scheme' .OR. & 1263 scalar_advec == 'ws-scheme-mono' ) THEN 1264 ! 1265 !-- Set flags to steer the degradation of the advection scheme in advec_ws 1266 !-- near topography, inflow- and outflow boundaries as well as bottom and 1267 !-- top of model domain. wall_flags_0 remains zero for all non-prognostic 1268 !-- grid points. 1269 DO i = nxl, nxr 1270 DO j = nys, nyn 1271 DO k = nzb_s_inner(j,i)+1, nzt 1272 ! 1273 !-- scalar - x-direction 1274 !-- WS1 (0), WS3 (1), WS5 (2) 1275 IF ( ( k <= nzb_s_inner(j,i+1) .OR. k <= nzb_s_inner(j,i+2) & 1276 .OR. k <= nzb_s_inner(j,i-1) ) & 1277 .OR. ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1278 .AND. i == nxl ) .OR. & 1279 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1280 .AND. i == nxr ) ) & 1281 THEN 1282 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 ) 1283 ELSEIF ( ( k <= nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i+1)& 1284 .AND. k > nzb_s_inner(j,i+2)& 1285 .AND. k > nzb_s_inner(j,i-1)& 1286 ) .OR. & 1287 ( k <= nzb_s_inner(j,i-2) .AND. k > nzb_s_inner(j,i+1)& 1288 .AND. k > nzb_s_inner(j,i+2)& 1289 .AND. k > nzb_s_inner(j,i-1)& 1290 ) & 1291 .OR. & 1292 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1293 .AND. i == nxr-1 ) .OR. & 1294 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1295 .AND. i == nxlu ) ) & 1296 THEN 1297 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 ) 1298 ELSEIF ( k > nzb_s_inner(j,i+1) .AND. k > nzb_s_inner(j,i+2) & 1299 .AND. k > nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i-1) & 1300 .AND. k > nzb_s_inner(j,i-2) ) & 1301 THEN 1302 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 ) 1303 ENDIF 1304 ! 1305 !-- scalar - y-direction 1306 !-- WS1 (3), WS3 (4), WS5 (5) 1307 IF ( ( k <= nzb_s_inner(j+1,i) .OR. k <= nzb_s_inner(j+2,i) & 1308 .OR. k <= nzb_s_inner(j-1,i) ) & 1309 .OR. & 1310 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1311 .AND. j == nys ) .OR. & 1312 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1313 .AND. j == nyn ) ) & 1314 THEN 1315 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 ) 1316 ! 1317 !-- WS3 1318 ELSEIF ( ( k <= nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j+1,i)& 1319 .AND. k > nzb_s_inner(j+2,i)& 1320 .AND. k > nzb_s_inner(j-1,i)& 1321 ) .OR. & 1322 ( k <= nzb_s_inner(j-2,i) .AND. k > nzb_s_inner(j+1,i)& 1323 .AND. k > nzb_s_inner(j+2,i)& 1324 .AND. k > nzb_s_inner(j-1,i)& 1325 ) & 1326 .OR. & 1327 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1328 .AND. j == nysv ) .OR. & 1329 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1330 .AND. j == nyn-1 ) ) & 1331 THEN 1332 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 ) 1333 ! 1334 !-- WS5 1335 ELSEIF ( k > nzb_s_inner(j+1,i) .AND. k > nzb_s_inner(j+2,i) & 1336 .AND. k > nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j-1,i) & 1337 .AND. k > nzb_s_inner(j-2,i) ) & 1338 THEN 1339 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 ) 1340 ENDIF 1341 ! 1342 !-- scalar - z-direction 1343 !-- WS1 (6), WS3 (7), WS5 (8) 1344 flag_set = .FALSE. 1345 IF ( k == nzb_s_inner(j,i) + 1 .OR. k == nzt ) THEN 1346 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 ) 1347 flag_set = .TRUE. 1348 ELSEIF ( k == nzb_s_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 1349 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 7 ) 1350 flag_set = .TRUE. 1351 ELSEIF ( k > nzb_s_inner(j,i) .AND. .NOT. flag_set ) THEN 1352 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 ) 1353 ENDIF 1354 1355 ENDDO 1356 ENDDO 1357 ENDDO 1358 ENDIF 1359 1360 IF ( momentum_advec == 'ws-scheme' ) THEN 1361 ! 1362 !-- Set wall_flags_0 to steer the degradation of the advection scheme in advec_ws 1363 !-- near topography, inflow- and outflow boundaries as well as bottom and 1364 !-- top of model domain. wall_flags_0 remains zero for all non-prognostic 1365 !-- grid points. 1366 DO i = nxl, nxr 1367 DO j = nys, nyn 1368 DO k = nzb+1, nzt 1369 ! 1370 !-- At first, set flags to WS1. 1371 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 1372 !-- in order to handle the left/south flux. 1373 !-- near vertical walls. 1374 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 ) 1375 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 ) 1376 ! 1377 !-- u component - x-direction 1378 !-- WS1 (9), WS3 (10), WS5 (11) 1379 IF ( k <= nzb_u_inner(j,i+1) .OR. & 1380 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1381 .AND. i <= nxlu ) .OR. & 1382 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1383 .AND. i == nxr ) ) & 1384 THEN 1385 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 ) 1386 ELSEIF ( ( k <= nzb_u_inner(j,i+2) .AND. & 1387 k > nzb_u_inner(j,i+1) ) .OR. & 1388 k <= nzb_u_inner(j,i-1) & 1389 .OR. & 1390 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1391 .AND. i == nxr-1 ) .OR. & 1392 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1393 .AND. i == nxlu+1) ) & 1394 THEN 1395 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 ) 1396 ! 1397 !-- Clear flag for WS1 1398 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 ) 1399 ELSEIF ( k > nzb_u_inner(j,i+1) .AND. k > nzb_u_inner(j,i+2) & 1400 .AND. k > nzb_u_inner(j,i-1) ) & 1401 THEN 1402 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 ) 1403 ! 1404 !-- Clear flag for WS1 1405 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 ) 1406 ENDIF 1407 ! 1408 !-- u component - y-direction 1409 !-- WS1 (12), WS3 (13), WS5 (14) 1410 IF ( k <= nzb_u_inner(j+1,i) .OR. & 1411 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1412 .AND. j == nys ) .OR. & 1413 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1414 .AND. j == nyn ) ) & 1415 THEN 1416 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 ) 1417 ELSEIF ( ( k <= nzb_u_inner(j+2,i) .AND. & 1418 k > nzb_u_inner(j+1,i) ) .OR. & 1419 k <= nzb_u_inner(j-1,i) & 1420 .OR. & 1421 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1422 .AND. j == nysv ) .OR. & 1423 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1424 .AND. j == nyn-1 ) ) & 1425 THEN 1426 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 ) 1427 ! 1428 !-- Clear flag for WS1 1429 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 ) 1430 ELSEIF ( k > nzb_u_inner(j+1,i) .AND. k > nzb_u_inner(j+2,i) & 1431 .AND. k > nzb_u_inner(j-1,i) ) & 1432 THEN 1433 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 ) 1434 ! 1435 !-- Clear flag for WS1 1436 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 ) 1437 ENDIF 1438 ! 1439 !-- u component - z-direction 1440 !-- WS1 (15), WS3 (16), WS5 (17) 1441 flag_set = .FALSE. 1442 IF ( k == nzb_u_inner(j,i) + 1 .OR. k == nzt ) THEN 1443 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 ) 1444 flag_set = .TRUE. 1445 ELSEIF ( k == nzb_u_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 1446 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 ) 1447 flag_set = .TRUE. 1448 ELSEIF ( k > nzb_u_inner(j,i) .AND. .NOT. flag_set ) THEN 1449 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 ) 1450 ENDIF 1451 1452 ENDDO 1453 ENDDO 1454 ENDDO 1455 1456 DO i = nxl, nxr 1457 DO j = nys, nyn 1458 DO k = nzb+1, nzt 1459 ! 1460 !-- At first, set flags to WS1. 1461 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 1462 !-- in order to handle the left/south flux. 1463 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 ) 1464 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 ) 1465 ! 1466 !-- v component - x-direction 1467 !-- WS1 (18), WS3 (19), WS5 (20) 1468 IF ( k <= nzb_v_inner(j,i+1) .OR. & 1469 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1470 .AND. i == nxl ) .OR. & 1471 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1472 .AND. i == nxr ) ) & 1473 THEN 1474 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 ) 1475 ! 1476 !-- WS3 1477 ELSEIF ( ( k <= nzb_v_inner(j,i+2) .AND. & 1478 k > nzb_v_inner(j,i+1) ) .OR. & 1479 k <= nzb_v_inner(j,i-1) & 1480 .OR. & 1481 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1482 .AND. i == nxr-1 ) .OR. & 1483 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1484 .AND. i == nxlu ) ) & 1485 THEN 1486 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 ) 1487 ! 1488 !-- Clear flag for WS1 1489 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 ) 1490 ELSEIF ( k > nzb_v_inner(j,i+1) .AND. k > nzb_v_inner(j,i+2) & 1491 .AND. k > nzb_v_inner(j,i-1) ) & 1492 THEN 1493 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 ) 1494 ! 1495 !-- Clear flag for WS1 1496 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 ) 1497 ENDIF 1498 ! 1499 !-- v component - y-direction 1500 !-- WS1 (21), WS3 (22), WS5 (23) 1501 IF ( k <= nzb_v_inner(j+1,i) .OR. & 1502 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1503 .AND. j <= nysv ) .OR. & 1504 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1505 .AND. j == nyn ) ) & 1506 THEN 1507 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 ) 1508 ELSEIF ( ( k <= nzb_v_inner(j+2,i) .AND. & 1509 k > nzb_v_inner(j+1,i) ) .OR. & 1510 k <= nzb_v_inner(j-1,i) & 1511 .OR. & 1512 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1513 .AND. j == nysv+1) .OR. & 1514 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1515 .AND. j == nyn-1 ) ) & 1516 THEN 1517 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 ) 1518 ! 1519 !-- Clear flag for WS1 1520 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 ) 1521 ELSEIF ( k > nzb_v_inner(j+1,i) .AND. k > nzb_v_inner(j+2,i) & 1522 .AND. k > nzb_v_inner(j-1,i) ) & 1523 THEN 1524 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 ) 1525 ! 1526 !-- Clear flag for WS1 1527 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 ) 1528 ENDIF 1529 ! 1530 !-- v component - z-direction 1531 !-- WS1 (24), WS3 (25), WS5 (26) 1532 flag_set = .FALSE. 1533 IF ( k == nzb_v_inner(j,i) + 1 .OR. k == nzt ) THEN 1534 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 ) 1535 flag_set = .TRUE. 1536 ELSEIF ( k == nzb_v_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 1537 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 ) 1538 flag_set = .TRUE. 1539 ELSEIF ( k > nzb_v_inner(j,i) .AND. .NOT. flag_set ) THEN 1540 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 ) 1541 ENDIF 1542 1543 ENDDO 1544 ENDDO 1545 ENDDO 1546 DO i = nxl, nxr 1547 DO j = nys, nyn 1548 DO k = nzb+1, nzt 1549 ! 1550 !-- At first, set flags to WS1. 1551 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 1552 !-- in order to handle the left/south flux. 1553 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 ) 1554 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 ) 1555 ! 1556 !-- w component - x-direction 1557 !-- WS1 (27), WS3 (28), WS5 (29) 1558 IF ( k <= nzb_w_inner(j,i+1) .OR. & 1559 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1560 .AND. i == nxl ) .OR. & 1561 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1562 .AND. i == nxr ) ) & 1563 THEN 1564 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 ) 1565 ELSEIF ( ( k <= nzb_w_inner(j,i+2) .AND. & 1566 k > nzb_w_inner(j,i+1) ) .OR. & 1567 k <= nzb_w_inner(j,i-1) & 1568 .OR. & 1569 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1570 .AND. i == nxr-1 ) .OR. & 1571 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1572 .AND. i == nxlu ) ) & 1573 THEN 1574 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 ) 1575 ! 1576 !-- Clear flag for WS1 1577 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 ) 1578 ELSEIF ( k > nzb_w_inner(j,i+1) .AND. k > nzb_w_inner(j,i+2) & 1579 .AND. k > nzb_w_inner(j,i-1) ) & 1580 THEN 1581 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i),29 ) 1582 ! 1583 !-- Clear flag for WS1 1584 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 ) 1585 ENDIF 1586 ! 1587 !-- w component - y-direction 1588 !-- WS1 (30), WS3 (31), WS5 (32) 1589 IF ( k <= nzb_w_inner(j+1,i) .OR. & 1590 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1591 .AND. j == nys ) .OR. & 1592 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1593 .AND. j == nyn ) ) & 1594 THEN 1595 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 ) 1596 ELSEIF ( ( k <= nzb_w_inner(j+2,i) .AND. & 1597 k > nzb_w_inner(j+1,i) ) .OR. & 1598 k <= nzb_w_inner(j-1,i) & 1599 .OR. & 1600 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1601 .AND. j == nysv ) .OR. & 1602 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1603 .AND. j == nyn-1 ) ) & 1604 THEN 1605 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 ) 1606 ! 1607 !-- Clear flag for WS1 1608 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 ) 1609 ELSEIF ( k > nzb_w_inner(j+1,i) .AND. k > nzb_w_inner(j+2,i) & 1610 .AND. k > nzb_w_inner(j-1,i) ) & 1611 THEN 1612 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 0 ) 1613 ! 1614 !-- Clear flag for WS1 1615 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 ) 1616 ENDIF 1617 ! 1618 !-- w component - z-direction 1619 !-- WS1 (33), WS3 (34), WS5 (35) 1620 flag_set = .FALSE. 1621 IF ( k == nzb_w_inner(j,i) .OR. k == nzb_w_inner(j,i) + 1 & 1622 .OR. k == nzt ) THEN 1623 ! 1624 !-- Please note, at k == nzb_w_inner(j,i) a flag is explictely 1625 !-- set, although this is not a prognostic level. However, 1626 !-- contrary to the advection of u,v and s this is necessary 1627 !-- because flux_t(nzb_w_inner(j,i)) is used for the tendency 1628 !-- at k == nzb_w_inner(j,i)+1. 1629 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 1 ) 1630 flag_set = .TRUE. 1631 ELSEIF ( k == nzb_w_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 1632 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 2 ) 1633 flag_set = .TRUE. 1634 ELSEIF ( k > nzb_w_inner(j,i) .AND. .NOT. flag_set ) THEN 1635 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 3 ) 1636 ENDIF 1637 1638 ENDDO 1639 ENDDO 1640 ENDDO 1641 1642 ENDIF 1643 1644 ! 1645 !-- Exchange 3D integer wall_flags. 1646 IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' & 1647 .OR. scalar_advec == 'ws-scheme-mono' ) THEN 1648 ! 1649 !-- Exchange ghost points for advection flags 1650 CALL exchange_horiz_int( wall_flags_0, nbgp ) 1651 CALL exchange_horiz_int( wall_flags_00, nbgp ) 1652 ! 1653 !-- Set boundary flags at inflow and outflow boundary in case of 1654 !-- non-cyclic boundary conditions. 1655 IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN 1656 wall_flags_0(:,:,nxl-1) = wall_flags_0(:,:,nxl) 1657 wall_flags_00(:,:,nxl-1) = wall_flags_00(:,:,nxl) 1658 ENDIF 1659 1660 IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN 1661 wall_flags_0(:,:,nxr+1) = wall_flags_0(:,:,nxr) 1662 wall_flags_00(:,:,nxr+1) = wall_flags_00(:,:,nxr) 1663 ENDIF 1664 1665 IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN 1666 wall_flags_0(:,nyn+1,:) = wall_flags_0(:,nyn,:) 1667 wall_flags_00(:,nyn+1,:) = wall_flags_00(:,nyn,:) 1668 ENDIF 1669 1670 IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN 1671 wall_flags_0(:,nys-1,:) = wall_flags_0(:,nys,:) 1672 wall_flags_00(:,nys-1,:) = wall_flags_00(:,nys,:) 1673 ENDIF 1674 1309 ! 1310 !-- Init flags for ws-scheme to degrade order near walls 1311 IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' .OR.& 1312 scalar_advec == 'ws-scheme-mono' ) THEN 1313 CALL ws_init_flags 1675 1314 ENDIF 1676 1315
Note: See TracChangeset
for help on using the changeset viewer.
