463 | | ! |
464 | | !-- u component - x-direction |
465 | | !-- WS1 (0), WS3 (1), WS5 (2) |
466 | | IF ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),1) .OR. & |
467 | | ( ( bc_dirichlet_l .OR. bc_radiation_l ) & |
468 | | .AND. i <= nxlu ) .OR. & |
469 | | ( ( bc_dirichlet_r .OR. bc_radiation_r ) & |
470 | | .AND. i == nxr ) ) & |
471 | | THEN |
472 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 ) |
473 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),1) .AND. & |
474 | | BTEST(wall_flags_total_0(k,j,i+1),1) .OR. & |
475 | | .NOT. BTEST(wall_flags_total_0(k,j,i-1),1) ) & |
476 | | .OR. & |
477 | | ( ( bc_dirichlet_r .OR. bc_radiation_r ) & |
478 | | .AND. i == nxr-1 ) .OR. & |
479 | | ( ( bc_dirichlet_l .OR. bc_radiation_l ) & |
480 | | .AND. i == nxlu+1) ) & |
481 | | THEN |
482 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 ) |
483 | | ! |
484 | | !-- Clear flag for WS1 |
485 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 ) |
486 | | ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),1) .AND. & |
487 | | BTEST(wall_flags_total_0(k,j,i+2),1) .AND. & |
488 | | BTEST(wall_flags_total_0(k,j,i-1),1) ) & |
489 | | THEN |
490 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 ) |
491 | | ! |
492 | | !-- Clear flag for WS1 |
493 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 ) |
494 | | ENDIF |
495 | | ! |
496 | | !-- u component - y-direction |
497 | | !-- WS1 (3), WS3 (4), WS5 (5) |
498 | | IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),1) .OR. & |
499 | | ( ( bc_dirichlet_s .OR. bc_radiation_s ) & |
500 | | .AND. j == nys ) .OR. & |
501 | | ( ( bc_dirichlet_n .OR. bc_radiation_n ) & |
502 | | .AND. j == nyn ) ) & |
503 | | THEN |
504 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 ) |
505 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),1) .AND. & |
506 | | BTEST(wall_flags_total_0(k,j+1,i),1) .OR. & |
507 | | .NOT. BTEST(wall_flags_total_0(k,j-1,i),1) ) & |
508 | | .OR. & |
509 | | ( ( bc_dirichlet_s .OR. bc_radiation_s ) & |
510 | | .AND. j == nysv ) .OR. & |
511 | | ( ( bc_dirichlet_n .OR. bc_radiation_n ) & |
512 | | .AND. j == nyn-1 ) ) & |
513 | | THEN |
514 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 ) |
515 | | ! |
516 | | !-- Clear flag for WS1 |
517 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) |
518 | | ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),1) .AND. & |
519 | | BTEST(wall_flags_total_0(k,j+2,i),1) .AND. & |
520 | | BTEST(wall_flags_total_0(k,j-1,i),1) ) & |
521 | | THEN |
522 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 ) |
523 | | ! |
524 | | !-- Clear flag for WS1 |
525 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) |
526 | | ENDIF |
527 | | ! |
528 | | !-- u component - z-direction. Fluxes are calculated on w-grid |
529 | | !-- level. Boundary u-values at/within walls aren't used. |
530 | | !-- WS1 (6), WS3 (7), WS5 (8) |
531 | | IF ( k == nzb+1 ) THEN |
532 | | k_mm = nzb |
533 | | ELSE |
534 | | k_mm = k - 2 |
535 | | ENDIF |
536 | | IF ( k > nzt-1 ) THEN |
537 | | k_pp = nzt+1 |
538 | | ELSE |
539 | | k_pp = k + 2 |
540 | | ENDIF |
541 | | IF ( k > nzt-2 ) THEN |
542 | | k_ppp = nzt+1 |
543 | | ELSE |
544 | | k_ppp = k + 3 |
545 | | ENDIF |
546 | | |
547 | | flag_set = .FALSE. |
548 | | IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),1) .AND. & |
549 | | BTEST(wall_flags_total_0(k,j,i),1) .AND. & |
550 | | BTEST(wall_flags_total_0(k+1,j,i),1) ) .OR. & |
551 | | ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),1) .AND. & |
552 | | BTEST(wall_flags_total_0(k+1,j,i),1) .AND. & |
553 | | BTEST(wall_flags_total_0(k,j,i),1) ) .OR. & |
554 | | ( k == nzt .AND. symmetry_flag == 0 ) ) & |
555 | | THEN |
556 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 ) |
557 | | flag_set = .TRUE. |
558 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),1) .OR. & |
559 | | .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),1) ) .AND.& |
560 | | BTEST(wall_flags_total_0(k-1,j,i),1) .AND.& |
561 | | BTEST(wall_flags_total_0(k,j,i),1) .AND.& |
562 | | BTEST(wall_flags_total_0(k+1,j,i),1) .AND.& |
563 | | BTEST(wall_flags_total_0(k_pp,j,i),1) .AND.& |
564 | | .NOT. flag_set .OR.& |
565 | | ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) & |
566 | | THEN |
567 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 ) |
568 | | flag_set = .TRUE. |
569 | | ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),1) .AND.& |
570 | | BTEST(wall_flags_total_0(k-1,j,i),1) .AND.& |
571 | | BTEST(wall_flags_total_0(k,j,i),1) .AND.& |
572 | | BTEST(wall_flags_total_0(k+1,j,i),1) .AND.& |
573 | | BTEST(wall_flags_total_0(k_pp,j,i),1) .AND.& |
574 | | BTEST(wall_flags_total_0(k_ppp,j,i),1) .AND.& |
575 | | .NOT. flag_set ) & |
576 | | THEN |
577 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 8 ) |
578 | | ENDIF |
579 | | |
580 | | ENDDO |
| 488 | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),1) .AND. & |
| 489 | BTEST(wall_flags_total_0(k,j+1,i),1) .OR. & |
| 490 | .NOT. BTEST(wall_flags_total_0(k,j-1,i),1) ) & |
| 491 | .OR. & |
| 492 | ( ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j == nysv ) .OR. & |
| 493 | ( ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j == nyn-1 ) ) & |
| 494 | THEN |
| 495 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 ) |
| 496 | ! |
| 497 | !-- Clear flag for WS1 |
| 498 | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) |
| 499 | ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),1) .AND. & |
| 500 | BTEST(wall_flags_total_0(k,j+2,i),1) .AND. & |
| 501 | BTEST(wall_flags_total_0(k,j-1,i),1) ) & |
| 502 | THEN |
| 503 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 ) |
| 504 | ! |
| 505 | !-- Clear flag for WS1 |
| 506 | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) |
| 507 | ENDIF |
| 508 | ! |
| 509 | !-- u component - z-direction. Fluxes are calculated on w-grid level. Boundary u-values |
| 510 | !-- at/within walls aren't used. |
| 511 | !-- WS1 (6), WS3 (7), WS5 (8) |
| 512 | IF ( k == nzb+1 ) THEN |
| 513 | k_mm = nzb |
| 514 | ELSE |
| 515 | k_mm = k - 2 |
| 516 | ENDIF |
| 517 | IF ( k > nzt-1 ) THEN |
| 518 | k_pp = nzt+1 |
| 519 | ELSE |
| 520 | k_pp = k + 2 |
| 521 | ENDIF |
| 522 | IF ( k > nzt-2 ) THEN |
| 523 | k_ppp = nzt+1 |
| 524 | ELSE |
| 525 | k_ppp = k + 3 |
| 526 | ENDIF |
| 527 | |
| 528 | flag_set = .FALSE. |
| 529 | IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),1) .AND. & |
| 530 | BTEST(wall_flags_total_0(k,j,i),1) .AND. & |
| 531 | BTEST(wall_flags_total_0(k+1,j,i),1) ) .OR. & |
| 532 | ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),1) .AND. & |
| 533 | BTEST(wall_flags_total_0(k+1,j,i),1) .AND. & |
| 534 | BTEST(wall_flags_total_0(k,j,i),1) ) .OR. & |
| 535 | ( k == nzt .AND. symmetry_flag == 0 ) ) & |
| 536 | THEN |
| 537 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 ) |
| 538 | flag_set = .TRUE. |
| 539 | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),1) .OR. & |
| 540 | .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),1) ) .AND. & |
| 541 | BTEST(wall_flags_total_0(k-1,j,i),1) .AND. & |
| 542 | BTEST(wall_flags_total_0(k,j,i),1) .AND. & |
| 543 | BTEST(wall_flags_total_0(k+1,j,i),1) .AND. & |
| 544 | BTEST(wall_flags_total_0(k_pp,j,i),1) .AND. & |
| 545 | .NOT. flag_set .OR. & |
| 546 | ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) & |
| 547 | THEN |
| 548 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 ) |
| 549 | flag_set = .TRUE. |
| 550 | ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),1) .AND. & |
| 551 | BTEST(wall_flags_total_0(k-1,j,i),1) .AND. & |
| 552 | BTEST(wall_flags_total_0(k,j,i),1) .AND. & |
| 553 | BTEST(wall_flags_total_0(k+1,j,i),1) .AND. & |
| 554 | BTEST(wall_flags_total_0(k_pp,j,i),1) .AND. & |
| 555 | BTEST(wall_flags_total_0(k_ppp,j,i),1) .AND. & |
| 556 | .NOT. flag_set ) & |
| 557 | THEN |
| 558 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 8 ) |
| 559 | ENDIF |
| 560 | |
594 | | ! |
595 | | !-- v component - x-direction |
596 | | !-- WS1 (9), WS3 (10), WS5 (11) |
597 | | IF ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),2) .OR. & |
598 | | ( ( bc_dirichlet_l .OR. bc_radiation_l ) & |
599 | | .AND. i == nxl ) .OR. & |
600 | | ( ( bc_dirichlet_r .OR. bc_radiation_r ) & |
601 | | .AND. i == nxr ) ) & |
602 | | THEN |
603 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 ) |
604 | | ! |
605 | | !-- WS3 |
606 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),2) .AND. & |
607 | | BTEST(wall_flags_total_0(k,j,i+1),2) ) .OR. & |
608 | | .NOT. BTEST(wall_flags_total_0(k,j,i-1),2) & |
609 | | .OR. & |
610 | | ( ( bc_dirichlet_r .OR. bc_radiation_r ) & |
611 | | .AND. i == nxr-1 ) .OR. & |
612 | | ( ( bc_dirichlet_l .OR. bc_radiation_l ) & |
613 | | .AND. i == nxlu ) ) & |
614 | | THEN |
615 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 ) |
616 | | ! |
617 | | !-- Clear flag for WS1 |
618 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 ) |
619 | | ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),2) .AND. & |
620 | | BTEST(wall_flags_total_0(k,j,i+2),2) .AND. & |
621 | | BTEST(wall_flags_total_0(k,j,i-1),2) ) & |
622 | | THEN |
623 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 ) |
624 | | ! |
625 | | !-- Clear flag for WS1 |
626 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 ) |
627 | | ENDIF |
628 | | ! |
629 | | !-- v component - y-direction |
630 | | !-- WS1 (12), WS3 (13), WS5 (14) |
631 | | IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),2) .OR. & |
632 | | ( ( bc_dirichlet_s .OR. bc_radiation_s ) & |
633 | | .AND. j <= nysv ) .OR. & |
634 | | ( ( bc_dirichlet_n .OR. bc_radiation_n ) & |
635 | | .AND. j == nyn ) ) & |
636 | | THEN |
637 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 ) |
638 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),2) .AND. & |
639 | | BTEST(wall_flags_total_0(k,j+1,i),2) .OR. & |
640 | | .NOT. BTEST(wall_flags_total_0(k,j-1,i),2) ) & |
641 | | .OR. & |
642 | | ( ( bc_dirichlet_s .OR. bc_radiation_s ) & |
643 | | .AND. j == nysv+1) .OR. & |
644 | | ( ( bc_dirichlet_n .OR. bc_radiation_n ) & |
645 | | .AND. j == nyn-1 ) ) & |
646 | | THEN |
647 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 ) |
648 | | ! |
649 | | !-- Clear flag for WS1 |
650 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 ) |
651 | | ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),2) .AND. & |
652 | | BTEST(wall_flags_total_0(k,j+2,i),2) .AND. & |
653 | | BTEST(wall_flags_total_0(k,j-1,i),2) ) & |
654 | | THEN |
655 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 14 ) |
656 | | ! |
657 | | !-- Clear flag for WS1 |
658 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 ) |
659 | | ENDIF |
660 | | ! |
661 | | !-- v component - z-direction. Fluxes are calculated on w-grid |
662 | | !-- level. Boundary v-values at/within walls aren't used. |
663 | | !-- WS1 (15), WS3 (16), WS5 (17) |
664 | | IF ( k == nzb+1 ) THEN |
665 | | k_mm = nzb |
666 | | ELSE |
667 | | k_mm = k - 2 |
668 | | ENDIF |
669 | | IF ( k > nzt-1 ) THEN |
670 | | k_pp = nzt+1 |
671 | | ELSE |
672 | | k_pp = k + 2 |
673 | | ENDIF |
674 | | IF ( k > nzt-2 ) THEN |
675 | | k_ppp = nzt+1 |
676 | | ELSE |
677 | | k_ppp = k + 3 |
678 | | ENDIF |
679 | | |
680 | | flag_set = .FALSE. |
681 | | IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),2) .AND.& |
682 | | BTEST(wall_flags_total_0(k,j,i),2) .AND.& |
683 | | BTEST(wall_flags_total_0(k+1,j,i),2) ) .OR. & |
684 | | ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),2) .AND.& |
685 | | BTEST(wall_flags_total_0(k+1,j,i),2) .AND.& |
686 | | BTEST(wall_flags_total_0(k,j,i),2) ) .OR. & |
687 | | ( k == nzt .AND. symmetry_flag == 0 ) ) & |
688 | | THEN |
689 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 ) |
690 | | flag_set = .TRUE. |
691 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),2) .OR. & |
692 | | .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),2) ) .AND.& |
693 | | BTEST(wall_flags_total_0(k-1,j,i),2) .AND.& |
694 | | BTEST(wall_flags_total_0(k,j,i),2) .AND.& |
695 | | BTEST(wall_flags_total_0(k+1,j,i),2) .AND.& |
696 | | BTEST(wall_flags_total_0(k_pp,j,i),2) .AND.& |
697 | | .NOT. flag_set .OR.& |
698 | | ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) & |
699 | | THEN |
700 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 ) |
701 | | flag_set = .TRUE. |
702 | | ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),2) .AND.& |
703 | | BTEST(wall_flags_total_0(k-1,j,i),2) .AND.& |
704 | | BTEST(wall_flags_total_0(k,j,i),2) .AND.& |
705 | | BTEST(wall_flags_total_0(k+1,j,i),2) .AND.& |
706 | | BTEST(wall_flags_total_0(k_pp,j,i),2) .AND.& |
707 | | BTEST(wall_flags_total_0(k_ppp,j,i),2) .AND.& |
708 | | .NOT. flag_set ) & |
709 | | THEN |
710 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 17 ) |
711 | | ENDIF |
712 | | |
713 | | ENDDO |
| 613 | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),2) .AND. & |
| 614 | BTEST(wall_flags_total_0(k,j+1,i),2) .OR. & |
| 615 | .NOT. BTEST(wall_flags_total_0(k,j-1,i),2) ) & |
| 616 | .OR. & |
| 617 | ( ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j == nysv+1) .OR. & |
| 618 | ( ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j == nyn-1 ) ) & |
| 619 | THEN |
| 620 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 ) |
| 621 | ! |
| 622 | !-- Clear flag for WS1 |
| 623 | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 ) |
| 624 | ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),2) .AND. & |
| 625 | BTEST(wall_flags_total_0(k,j+2,i),2) .AND. & |
| 626 | BTEST(wall_flags_total_0(k,j-1,i),2) ) & |
| 627 | THEN |
| 628 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 14 ) |
| 629 | ! |
| 630 | !-- Clear flag for WS1 |
| 631 | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 ) |
| 632 | ENDIF |
| 633 | ! |
| 634 | !-- v component - z-direction. Fluxes are calculated on w-grid level. Boundary v-values |
| 635 | !-- at/within walls aren't used. |
| 636 | !-- WS1 (15), WS3 (16), WS5 (17) |
| 637 | IF ( k == nzb+1 ) THEN |
| 638 | k_mm = nzb |
| 639 | ELSE |
| 640 | k_mm = k - 2 |
| 641 | ENDIF |
| 642 | IF ( k > nzt-1 ) THEN |
| 643 | k_pp = nzt+1 |
| 644 | ELSE |
| 645 | k_pp = k + 2 |
| 646 | ENDIF |
| 647 | IF ( k > nzt-2 ) THEN |
| 648 | k_ppp = nzt+1 |
| 649 | ELSE |
| 650 | k_ppp = k + 3 |
| 651 | ENDIF |
| 652 | |
| 653 | flag_set = .FALSE. |
| 654 | IF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),2) .AND. & |
| 655 | BTEST(wall_flags_total_0(k,j,i),2) .AND. & |
| 656 | BTEST(wall_flags_total_0(k+1,j,i),2) ) .OR. & |
| 657 | ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),2) .AND. & |
| 658 | BTEST(wall_flags_total_0(k+1,j,i),2) .AND. & |
| 659 | BTEST(wall_flags_total_0(k,j,i),2) ) .OR. & |
| 660 | ( k == nzt .AND. symmetry_flag == 0 ) ) & |
| 661 | THEN |
| 662 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 ) |
| 663 | flag_set = .TRUE. |
| 664 | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k_mm,j,i),2) .OR. & |
| 665 | .NOT. BTEST(wall_flags_total_0(k_ppp,j,i),2) ) .AND. & |
| 666 | BTEST(wall_flags_total_0(k-1,j,i),2) .AND. & |
| 667 | BTEST(wall_flags_total_0(k,j,i),2) .AND. & |
| 668 | BTEST(wall_flags_total_0(k+1,j,i),2) .AND. & |
| 669 | BTEST(wall_flags_total_0(k_pp,j,i),2) .AND. & |
| 670 | .NOT. flag_set .OR. & |
| 671 | ( k == nzt - 1 .AND. symmetry_flag == 0 ) ) & |
| 672 | THEN |
| 673 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 ) |
| 674 | flag_set = .TRUE. |
| 675 | ELSEIF ( BTEST(wall_flags_total_0(k_mm,j,i),2) .AND. & |
| 676 | BTEST(wall_flags_total_0(k-1,j,i),2) .AND. & |
| 677 | BTEST(wall_flags_total_0(k,j,i),2) .AND. & |
| 678 | BTEST(wall_flags_total_0(k+1,j,i),2) .AND. & |
| 679 | BTEST(wall_flags_total_0(k_pp,j,i),2) .AND. & |
| 680 | BTEST(wall_flags_total_0(k_ppp,j,i),2) .AND. & |
| 681 | .NOT. flag_set ) & |
| 682 | THEN |
| 683 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 17 ) |
| 684 | ENDIF |
| 685 | |
727 | | ! |
728 | | !-- w component - x-direction |
729 | | !-- WS1 (18), WS3 (19), WS5 (20) |
730 | | IF ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),3) .OR. & |
731 | | ( ( bc_dirichlet_l .OR. bc_radiation_l ) & |
732 | | .AND. i == nxl ) .OR. & |
733 | | ( ( bc_dirichlet_r .OR. bc_radiation_r ) & |
734 | | .AND. i == nxr ) ) & |
735 | | THEN |
736 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 ) |
737 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+2),3) .AND. & |
738 | | BTEST(wall_flags_total_0(k,j,i+1),3) .OR. & |
739 | | .NOT. BTEST(wall_flags_total_0(k,j,i-1),3) ) & |
740 | | .OR. & |
741 | | ( ( bc_dirichlet_r .OR. bc_radiation_r ) & |
742 | | .AND. i == nxr-1 ) .OR. & |
743 | | ( ( bc_dirichlet_l .OR. bc_radiation_l ) & |
744 | | .AND. i == nxlu ) ) & |
745 | | THEN |
746 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 ) |
747 | | ! |
748 | | !-- Clear flag for WS1 |
749 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 ) |
750 | | ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),3) .AND. & |
751 | | BTEST(wall_flags_total_0(k,j,i+2),3) .AND. & |
752 | | BTEST(wall_flags_total_0(k,j,i-1),3) ) & |
753 | | THEN |
754 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 ) |
755 | | ! |
756 | | !-- Clear flag for WS1 |
757 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 ) |
758 | | ENDIF |
759 | | ! |
760 | | !-- w component - y-direction |
761 | | !-- WS1 (21), WS3 (22), WS5 (23) |
762 | | IF ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),3) .OR. & |
763 | | ( ( bc_dirichlet_s .OR. bc_radiation_s ) & |
764 | | .AND. j == nys ) .OR. & |
765 | | ( ( bc_dirichlet_n .OR. bc_radiation_n ) & |
766 | | .AND. j == nyn ) ) & |
767 | | THEN |
768 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 ) |
769 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),3) .AND. & |
770 | | BTEST(wall_flags_total_0(k,j+1,i),3) .OR. & |
771 | | .NOT. BTEST(wall_flags_total_0(k,j-1,i),3) ) & |
772 | | .OR. & |
773 | | ( ( bc_dirichlet_s .OR. bc_radiation_s ) & |
774 | | .AND. j == nysv ) .OR. & |
775 | | ( ( bc_dirichlet_n .OR. bc_radiation_n ) & |
776 | | .AND. j == nyn-1 ) ) & |
777 | | THEN |
778 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 ) |
779 | | ! |
780 | | !-- Clear flag for WS1 |
781 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 ) |
782 | | ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),3) .AND. & |
783 | | BTEST(wall_flags_total_0(k,j+2,i),3) .AND. & |
784 | | BTEST(wall_flags_total_0(k,j-1,i),3) ) & |
785 | | THEN |
786 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 23 ) |
787 | | ! |
788 | | !-- Clear flag for WS1 |
789 | | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 ) |
790 | | ENDIF |
791 | | ! |
792 | | !-- w component - z-direction. Fluxes are calculated on scalar grid |
793 | | !-- level. Boundary w-values at walls are used. Flux at k=i is |
794 | | !-- defined at scalar position k=i+1 with i being an integer. |
795 | | !-- WS1 (24), WS3 (25), WS5 (26) |
796 | | IF ( k == nzb+1 ) THEN |
797 | | k_mm = nzb |
798 | | ELSE |
799 | | k_mm = k - 2 |
800 | | ENDIF |
801 | | IF ( k > nzt-1 ) THEN |
802 | | k_pp = nzt+1 |
803 | | ELSE |
804 | | k_pp = k + 2 |
805 | | ENDIF |
806 | | IF ( k > nzt-2 ) THEN |
807 | | k_ppp = nzt+1 |
808 | | ELSE |
809 | | k_ppp = k + 3 |
810 | | ENDIF |
811 | | |
812 | | flag_set = .FALSE. |
813 | | IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i),3) .AND. & |
814 | | BTEST(wall_flags_total_0(k+1,j,i),3) ) .OR. & |
815 | | ( .NOT. BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & |
816 | | BTEST(wall_flags_total_0(k,j,i),3) ) .OR. & |
817 | | k == nzt -1 ) & |
818 | | THEN |
819 | | ! |
820 | | !-- Please note, at k == nzb_w_inner(j,i) a flag is explicitly |
821 | | !-- set, although this is not a prognostic level. However, |
822 | | !-- contrary to the advection of u,v and s this is necessary |
823 | | !-- because flux_t(nzb_w_inner(j,i)) is used for the tendency |
824 | | !-- at k == nzb_w_inner(j,i)+1. |
825 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 ) |
826 | | flag_set = .TRUE. |
827 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),3) .AND. & |
828 | | BTEST(wall_flags_total_0(k,j,i),3) .AND. & |
829 | | BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & |
830 | | BTEST(wall_flags_total_0(k_pp,j,i),3) ) .OR. & |
831 | | ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),3) .AND. & |
832 | | BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & |
833 | | BTEST(wall_flags_total_0(k,j,i),3) .AND. & |
834 | | BTEST(wall_flags_total_0(k-1,j,i),3) ) .AND. & |
835 | | .NOT. flag_set .OR. & |
836 | | k == nzt - 2 ) & |
837 | | THEN |
838 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 ) |
839 | | flag_set = .TRUE. |
840 | | ELSEIF ( BTEST(wall_flags_total_0(k-1,j,i),3) .AND. & |
841 | | BTEST(wall_flags_total_0(k,j,i),3) .AND. & |
842 | | BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & |
843 | | BTEST(wall_flags_total_0(k_pp,j,i),3) .AND. & |
844 | | .NOT. flag_set ) & |
845 | | THEN |
846 | | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 26 ) |
847 | | ENDIF |
848 | | |
849 | | ENDDO |
| 736 | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+2,i),3) .AND. & |
| 737 | BTEST(wall_flags_total_0(k,j+1,i),3) .OR. & |
| 738 | .NOT. BTEST(wall_flags_total_0(k,j-1,i),3) ) & |
| 739 | .OR. & |
| 740 | ( ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j == nysv ) .OR. & |
| 741 | ( ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j == nyn-1 ) ) & |
| 742 | THEN |
| 743 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 ) |
| 744 | ! |
| 745 | !-- Clear flag for WS1 |
| 746 | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 ) |
| 747 | ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),3) .AND. & |
| 748 | BTEST(wall_flags_total_0(k,j+2,i),3) .AND. & |
| 749 | BTEST(wall_flags_total_0(k,j-1,i),3) ) & |
| 750 | THEN |
| 751 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 23 ) |
| 752 | ! |
| 753 | !-- Clear flag for WS1 |
| 754 | advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 ) |
| 755 | ENDIF |
| 756 | ! |
| 757 | !-- w component - z-direction. Fluxes are calculated on scalar grid level. Boundary |
| 758 | !-- w-values at walls are used. Flux at k=i is defined at scalar position k=i+1 with i |
| 759 | !-- being an integer. |
| 760 | !-- WS1 (24), WS3 (25), WS5 (26) |
| 761 | IF ( k == nzb+1 ) THEN |
| 762 | k_mm = nzb |
| 763 | ELSE |
| 764 | k_mm = k - 2 |
| 765 | ENDIF |
| 766 | IF ( k > nzt-1 ) THEN |
| 767 | k_pp = nzt+1 |
| 768 | ELSE |
| 769 | k_pp = k + 2 |
| 770 | ENDIF |
| 771 | IF ( k > nzt-2 ) THEN |
| 772 | k_ppp = nzt+1 |
| 773 | ELSE |
| 774 | k_ppp = k + 3 |
| 775 | ENDIF |
| 776 | |
| 777 | flag_set = .FALSE. |
| 778 | IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i),3) .AND. & |
| 779 | BTEST(wall_flags_total_0(k+1,j,i),3) ) .OR. & |
| 780 | ( .NOT. BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & |
| 781 | BTEST(wall_flags_total_0(k,j,i),3) ) .OR. & |
| 782 | k == nzt -1 ) & |
| 783 | THEN |
| 784 | ! |
| 785 | !-- Please note, at k == nzb_w_inner(j,i) a flag is explicitly set, although this is not |
| 786 | !-- a prognostic level. However, contrary to the advection of u,v and s this is |
| 787 | !-- necessary because flux_t(nzb_w_inner(j,i)) is used for the tendency at k == |
| 788 | !-- 0nzb_w_inner(j,i)+1. |
| 789 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 ) |
| 790 | flag_set = .TRUE. |
| 791 | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k-1,j,i),3) .AND. & |
| 792 | BTEST(wall_flags_total_0(k,j,i),3) .AND. & |
| 793 | BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & |
| 794 | BTEST(wall_flags_total_0(k_pp,j,i),3) ) .OR. & |
| 795 | ( .NOT. BTEST(wall_flags_total_0(k_pp,j,i),3) .AND. & |
| 796 | BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & |
| 797 | BTEST(wall_flags_total_0(k,j,i),3) .AND. & |
| 798 | BTEST(wall_flags_total_0(k-1,j,i),3) ) .AND. & |
| 799 | .NOT. flag_set .OR. & |
| 800 | k == nzt - 2 ) & |
| 801 | THEN |
| 802 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 ) |
| 803 | flag_set = .TRUE. |
| 804 | ELSEIF ( BTEST(wall_flags_total_0(k-1,j,i),3) .AND. & |
| 805 | BTEST(wall_flags_total_0(k,j,i),3) .AND. & |
| 806 | BTEST(wall_flags_total_0(k+1,j,i),3) .AND. & |
| 807 | BTEST(wall_flags_total_0(k_pp,j,i),3) .AND. & |
| 808 | .NOT. flag_set ) & |
| 809 | THEN |
| 810 | advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 26 ) |
| 811 | ENDIF |
| 812 | |
880 | | !> Initialization of flags to control the order of the advection scheme near |
881 | | !> solid walls and non-cyclic inflow boundaries, where the order is sucessively |
882 | | !> degraded. |
883 | | !------------------------------------------------------------------------------! |
884 | | SUBROUTINE ws_init_flags_scalar( non_cyclic_l, non_cyclic_n, non_cyclic_r, & |
885 | | non_cyclic_s, advc_flag, extensive_degrad ) |
886 | | |
887 | | |
888 | | INTEGER(iwp) :: i !< index variable along x |
889 | | INTEGER(iwp) :: j !< index variable along y |
890 | | INTEGER(iwp) :: k !< index variable along z |
891 | | INTEGER(iwp) :: k_mm !< dummy index along z |
892 | | INTEGER(iwp) :: k_pp !< dummy index along z |
893 | | INTEGER(iwp) :: k_ppp !< dummy index along z |
894 | | |
895 | | INTEGER(iwp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::& |
896 | | advc_flag !< flag array to control order of scalar advection |
897 | | |
898 | | LOGICAL :: flag_set !< steering variable for advection flags |
899 | | LOGICAL :: non_cyclic_l !< flag that indicates non-cyclic boundary on the left |
900 | | LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north |
901 | | LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right |
902 | | LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south |
903 | | |
904 | | LOGICAL, OPTIONAL :: extensive_degrad !< flag indicating that extensive degradation is required, e.g. for |
905 | | !< passive scalars nearby topography along the horizontal directions, |
906 | | !< as no monotonic limiter can be applied there |
907 | | ! |
908 | | !-- Set flags to steer the degradation of the advection scheme in advec_ws |
909 | | !-- near topography, inflow- and outflow boundaries as well as bottom and |
910 | | !-- top of model domain. advc_flags_m remains zero for all non-prognostic |
911 | | !-- grid points. |
912 | | DO i = nxl, nxr |
913 | | DO j = nys, nyn |
914 | | DO k = nzb+1, nzt |
915 | | IF ( .NOT. BTEST(wall_flags_total_0(k,j,i),0) ) CYCLE |
916 | | ! |
917 | | !-- scalar - x-direction |
918 | | !-- WS1 (0), WS3 (1), WS5 (2) |
919 | | IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),0) .OR. & |
920 | | .NOT. BTEST(wall_flags_total_0(k,j,i+2),0) .OR. & |
921 | | .NOT. BTEST(wall_flags_total_0(k,j,i-1),0) ) .OR. & |
922 | | ( non_cyclic_l .AND. i == 0 ) .OR. & |
923 | | ( non_cyclic_r .AND. i == nx ) ) THEN |
924 | | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) |
925 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+3),0) .AND. & |
926 | | BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & |
927 | | BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & |
928 | | BTEST(wall_flags_total_0(k,j,i-1),0) & |
929 | | ) .OR. & |
930 | | ( .NOT. BTEST(wall_flags_total_0(k,j,i-2),0) .AND. & |
931 | | BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & |
932 | | BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & |
933 | | BTEST(wall_flags_total_0(k,j,i-1),0) & |
934 | | ) & |
935 | | .OR. & |
936 | | ( non_cyclic_r .AND. i == nx-1 ) .OR. & |
937 | | ( non_cyclic_l .AND. i == 1 ) ) THEN |
938 | | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 ) |
939 | | ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & |
940 | | BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & |
941 | | BTEST(wall_flags_total_0(k,j,i+3),0) .AND. & |
942 | | BTEST(wall_flags_total_0(k,j,i-1),0) .AND. & |
943 | | BTEST(wall_flags_total_0(k,j,i-2),0) ) & |
944 | | THEN |
945 | | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 ) |
946 | | ENDIF |
947 | | ! |
948 | | !-- scalar - y-direction |
949 | | !-- WS1 (3), WS3 (4), WS5 (5) |
950 | | IF ( ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),0) .OR. & |
951 | | .NOT. BTEST(wall_flags_total_0(k,j+2,i),0) .OR. & |
952 | | .NOT. BTEST(wall_flags_total_0(k,j-1,i),0)) .OR. & |
953 | | ( non_cyclic_s .AND. j == 0 ) .OR. & |
954 | | ( non_cyclic_n .AND. j == ny ) ) THEN |
955 | | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) |
956 | | ! |
957 | | !-- WS3 |
958 | | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+3,i),0) .AND. & |
959 | | BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & |
960 | | BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & |
961 | | BTEST(wall_flags_total_0(k,j-1,i),0) & |
962 | | ) .OR. & |
963 | | ( .NOT. BTEST(wall_flags_total_0(k,j-2,i),0) .AND. & |
964 | | BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & |
965 | | BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & |
966 | | BTEST(wall_flags_total_0(k,j-1,i),0) & |
967 | | ) & |
968 | | .OR. & |
969 | | ( non_cyclic_s .AND. j == 1 ) .OR. & |
970 | | ( non_cyclic_n .AND. j == ny-1 ) ) THEN |
971 | | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 ) |
972 | | ! |
973 | | !-- WS5 |
974 | | ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & |
975 | | BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & |
976 | | BTEST(wall_flags_total_0(k,j+3,i),0) .AND. & |
977 | | BTEST(wall_flags_total_0(k,j-1,i),0) .AND. & |
978 | | BTEST(wall_flags_total_0(k,j-2,i),0) ) & |
979 | | THEN |
980 | | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 ) |
981 | | ENDIF |
982 | | ! |
983 | | !-- Near topography, set horizontal advection scheme to 1st order |
984 | | !-- for passive scalars, even if only one direction may be |
985 | | !-- blocked by topography. These locations will be identified |
986 | | !-- by wall_flags_total_0 bit 31. Note, since several modules define |
987 | | !-- advection flags but may apply different scalar boundary |
988 | | !-- conditions, bit 31 is temporarily stored on advc_flags. |
989 | | !-- Moreover, note that this extended degradtion for passive |
990 | | !-- scalars is not required for the vertical direction as there |
991 | | !-- the monotonic limiter can be applied. |
992 | | IF ( PRESENT( extensive_degrad ) ) THEN |
993 | | IF ( extensive_degrad ) THEN |
994 | | ! |
995 | | !-- At all grid points that are within a three-grid point |
996 | | !-- range to topography, set 1st-order scheme. |
997 | | IF( BTEST( advc_flag(k,j,i), 31 ) ) THEN |
998 | | ! |
999 | | !-- Clear flags that might indicate higher-order |
1000 | | !-- advection along x- and y-direction. |
| 844 | !> Initialization of flags to control the order of the advection scheme near solid walls and |
| 845 | !> non-cyclic inflow boundaries, where the order is sucessively degraded. |
| 846 | !--------------------------------------------------------------------------------------------------! |
| 847 | SUBROUTINE ws_init_flags_scalar( non_cyclic_l, non_cyclic_n, non_cyclic_r, non_cyclic_s, & |
| 848 | advc_flag, extensive_degrad ) |
| 849 | |
| 850 | |
| 851 | INTEGER(iwp) :: i !< index variable along x |
| 852 | INTEGER(iwp) :: j !< index variable along y |
| 853 | INTEGER(iwp) :: k !< index variable along z |
| 854 | INTEGER(iwp) :: k_mm !< dummy index along z |
| 855 | INTEGER(iwp) :: k_pp !< dummy index along z |
| 856 | INTEGER(iwp) :: k_ppp !< dummy index along z |
| 857 | |
| 858 | INTEGER(iwp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: & |
| 859 | advc_flag !< flag array to control order of scalar advection |
| 860 | |
| 861 | LOGICAL :: flag_set !< steering variable for advection flags |
| 862 | LOGICAL :: non_cyclic_l !< flag that indicates non-cyclic boundary on the left |
| 863 | LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north |
| 864 | LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right |
| 865 | LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south |
| 866 | |
| 867 | LOGICAL, OPTIONAL :: extensive_degrad !< flag indicating that extensive degradation is required, e.g. for |
| 868 | !< passive scalars nearby topography along the horizontal directions, |
| 869 | !< as no monotonic limiter can be applied there |
| 870 | ! |
| 871 | !-- Set flags to steer the degradation of the advection scheme in advec_ws near topography, inflow- |
| 872 | !-- and outflow boundaries as well as bottom and top of model domain. advc_flags_m remains zero for |
| 873 | !-- all non-prognostic grid points. |
| 874 | DO i = nxl, nxr |
| 875 | DO j = nys, nyn |
| 876 | DO k = nzb+1, nzt |
| 877 | IF ( .NOT. BTEST(wall_flags_total_0(k,j,i),0) ) CYCLE |
| 878 | ! |
| 879 | !-- scalar - x-direction |
| 880 | !-- WS1 (0), WS3 (1), WS5 (2) |
| 881 | IF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+1),0) .OR. & |
| 882 | .NOT. BTEST(wall_flags_total_0(k,j,i+2),0) .OR. & |
| 883 | .NOT. BTEST(wall_flags_total_0(k,j,i-1),0) ) .OR. & |
| 884 | ( non_cyclic_l .AND. i == 0 ) .OR. & |
| 885 | ( non_cyclic_r .AND. i == nx ) ) & |
| 886 | THEN |
| 887 | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) |
| 888 | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j,i+3),0) .AND. & |
| 889 | BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & |
| 890 | BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & |
| 891 | BTEST(wall_flags_total_0(k,j,i-1),0) & |
| 892 | ) .OR. & |
| 893 | ( .NOT. BTEST(wall_flags_total_0(k,j,i-2),0) .AND. & |
| 894 | BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & |
| 895 | BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & |
| 896 | BTEST(wall_flags_total_0(k,j,i-1),0) & |
| 897 | ) .OR. & |
| 898 | ( non_cyclic_r .AND. i == nx-1 ) .OR. & |
| 899 | ( non_cyclic_l .AND. i == 1 ) ) & |
| 900 | THEN |
| 901 | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 ) |
| 902 | ELSEIF ( BTEST(wall_flags_total_0(k,j,i+1),0) .AND. & |
| 903 | BTEST(wall_flags_total_0(k,j,i+2),0) .AND. & |
| 904 | BTEST(wall_flags_total_0(k,j,i+3),0) .AND. & |
| 905 | BTEST(wall_flags_total_0(k,j,i-1),0) .AND. & |
| 906 | BTEST(wall_flags_total_0(k,j,i-2),0) ) & |
| 907 | THEN |
| 908 | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 ) |
| 909 | ENDIF |
| 910 | ! |
| 911 | !-- scalar - y-direction |
| 912 | !-- WS1 (3), WS3 (4), WS5 (5) |
| 913 | IF ( ( .NOT. BTEST(wall_flags_total_0(k,j+1,i),0) .OR. & |
| 914 | .NOT. BTEST(wall_flags_total_0(k,j+2,i),0) .OR. & |
| 915 | .NOT. BTEST(wall_flags_total_0(k,j-1,i),0)) .OR. & |
| 916 | ( non_cyclic_s .AND. j == 0 ) .OR. & |
| 917 | ( non_cyclic_n .AND. j == ny ) ) & |
| 918 | THEN |
| 919 | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) |
| 920 | ! |
| 921 | !-- WS3 |
| 922 | ELSEIF ( ( .NOT. BTEST(wall_flags_total_0(k,j+3,i),0) .AND. & |
| 923 | BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & |
| 924 | BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & |
| 925 | BTEST(wall_flags_total_0(k,j-1,i),0) & |
| 926 | ) .OR. & |
| 927 | ( .NOT. BTEST(wall_flags_total_0(k,j-2,i),0) .AND. & |
| 928 | BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & |
| 929 | BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & |
| 930 | BTEST(wall_flags_total_0(k,j-1,i),0) & |
| 931 | ) .OR. & |
| 932 | ( non_cyclic_s .AND. j == 1 ) .OR. & |
| 933 | ( non_cyclic_n .AND. j == ny-1 ) ) & |
| 934 | THEN |
| 935 | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 ) |
| 936 | ! |
| 937 | !-- WS5 |
| 938 | ELSEIF ( BTEST(wall_flags_total_0(k,j+1,i),0) .AND. & |
| 939 | BTEST(wall_flags_total_0(k,j+2,i),0) .AND. & |
| 940 | BTEST(wall_flags_total_0(k,j+3,i),0) .AND. & |
| 941 | BTEST(wall_flags_total_0(k,j-1,i),0) .AND. & |
| 942 | BTEST(wall_flags_total_0(k,j-2,i),0) ) & |
| 943 | THEN |
| 944 | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 ) |
| 945 | ENDIF |
| 946 | ! |
| 947 | !-- Near topography, set horizontal advection scheme to 1st order for passive scalars, even |
| 948 | !-- if only one direction may be blocked by topography. These locations will be identified |
| 949 | !-- by wall_flags_total_0 bit 31. Note, since several modules define advection flags but |
| 950 | !-- may apply different scalar boundary conditions, bit 31 is temporarily stored on |
| 951 | !-- advc_flags. |
| 952 | !-- Moreover, note that this extended degradtion for passive scalars is not required for |
| 953 | !-- the vertical direction as there the monotonic limiter can be applied. |
| 954 | IF ( PRESENT( extensive_degrad ) ) THEN |
| 955 | IF ( extensive_degrad ) THEN |
| 956 | ! |
| 957 | !-- At all grid points that are within a three-grid point range to topography, set |
| 958 | !-- 1st-order scheme. |
| 959 | IF( BTEST( advc_flag(k,j,i), 31 ) ) THEN |
| 960 | ! |
| 961 | !-- Clear flags that might indicate higher-order advection along x- and |
| 962 | !-- y-direction. |
| 963 | advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) |
| 964 | advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) |
| 965 | advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) |
| 966 | advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) |
| 967 | ! |
| 968 | !-- Set flags that indicate 1st-order advection along x- and y-direction. |
| 969 | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) |
| 970 | advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) |
| 971 | ENDIF |
| 972 | ! |
| 973 | !-- Adjacent to this extended degradation zone, successively upgrade the order of the |
| 974 | !-- scheme if this grid point isn't flagged with bit 31 (indicating extended |
| 975 | !-- degradation zone). |
| 976 | IF ( .NOT. BTEST( advc_flag(k,j,i), 31 ) ) THEN |
| 977 | ! |
| 978 | !-- x-direction. First, clear all previous settings, then set flag for 3rd-order |
| 979 | !-- scheme. |
| 980 | IF ( BTEST( advc_flag(k,j,i-1), 31 ) .AND. & |
| 981 | BTEST( advc_flag(k,j,i+1), 31 ) ) THEN |
| 982 | advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 ) |
1235 | | !------------------------------------------------------------------------------! |
1236 | | SUBROUTINE advec_s_ws_ij( advc_flag, i, j, sk, sk_char, swap_flux_y_local, & |
1237 | | swap_diss_y_local, swap_flux_x_local, & |
1238 | | swap_diss_x_local, i_omp, tn, & |
1239 | | non_cyclic_l, non_cyclic_n, & |
1240 | | non_cyclic_r, non_cyclic_s, & |
1241 | | flux_limitation ) |
1242 | | |
1243 | | |
1244 | | CHARACTER (LEN = *), INTENT(IN) :: sk_char !< string identifier, used for assign fluxes to the |
1245 | | !<correct dimension in the analysis array |
1246 | | |
1247 | | INTEGER(iwp) :: i !< grid index along x-direction |
1248 | | INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread |
1249 | | INTEGER(iwp) :: j !< grid index along y-direction |
1250 | | INTEGER(iwp) :: k !< grid index along z-direction |
1251 | | INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults |
1252 | | INTEGER(iwp) :: k_mmm !< k-3 index in disretization, can be modified to avoid segmentation faults |
1253 | | INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults |
1254 | | INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults |
1255 | | INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms |
1256 | | INTEGER(iwp) :: tn !< number of OpenMP thread |
1257 | | |
1258 | | INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: & |
1259 | | advc_flag !< flag array to control order of scalar advection |
1260 | | |
1261 | | LOGICAL :: non_cyclic_l !< flag that indicates non-cyclic boundary on the left |
1262 | | LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north |
1263 | | LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right |
1264 | | LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south |
1265 | | LOGICAL, OPTIONAL :: flux_limitation !< flag indicating flux limitation of the vertical advection |
1266 | | LOGICAL :: limiter !< control flag indicating the application of flux limitation |
1267 | | |
1268 | | REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom |
1269 | | REAL(wp) :: div !< velocity diverence on scalar grid |
1270 | | REAL(wp) :: div_in !< vertical flux divergence of ingoing fluxes |
1271 | | REAL(wp) :: div_out !< vertical flux divergence of outgoing fluxes |
1272 | | REAL(wp) :: f_corr_t !< correction flux at grid-cell top, i.e. the difference between high and low-order flux |
1273 | | REAL(wp) :: f_corr_d !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux |
1274 | | REAL(wp) :: f_corr_t_in !< correction flux of ingoing flux part at grid-cell top |
1275 | | REAL(wp) :: f_corr_d_in !< correction flux of ingoing flux part at grid-cell bottom |
1276 | | REAL(wp) :: f_corr_t_out !< correction flux of outgoing flux part at grid-cell top |
1277 | | REAL(wp) :: f_corr_d_out !< correction flux of outgoing flux part at grid-cell bottom |
1278 | | REAL(wp) :: fac_correction!< factor to limit the in- and outgoing fluxes |
1279 | | REAL(wp) :: flux_d !< 6th-order flux at grid box bottom |
1280 | | REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction |
1281 | | REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction |
1282 | | REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction |
1283 | | REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction |
1284 | | REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction |
1285 | | REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction |
1286 | | REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction |
1287 | | REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction |
1288 | | REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction |
1289 | | REAL(wp) :: max_val !< maximum value of the quanitity along the numerical stencil (in vertical direction) |
1290 | | REAL(wp) :: min_val !< maximum value of the quanitity along the numerical stencil (in vertical direction) |
1291 | | REAL(wp) :: mon !< monotone solution of the advection equation using 1st-order fluxes |
1292 | | REAL(wp) :: u_comp !< advection velocity along x-direction |
1293 | | REAL(wp) :: v_comp !< advection velocity along y-direction |
1294 | | ! |
1295 | | !-- sk is an array from parameter list. It should not be a pointer, because |
1296 | | !-- in that case the compiler can not assume a stride 1 and cannot perform |
1297 | | !-- a strided one vector load. Adding the CONTIGUOUS keyword makes things |
1298 | | !-- even worse, because the compiler cannot assume strided one in the |
1299 | | !-- caller side. |
1300 | | REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar |
1301 | | |
1302 | | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side |
1303 | | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side |
1304 | | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top |
1305 | | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side |
1306 | | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side |
1307 | | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top |
1308 | | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t_1st !< discretized 1st-order flux at top |
1309 | | |
1310 | | REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !< discretized artificial dissipation at southward-side |
1311 | | REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local !< discretized 6th-order flux at northward-side |
1312 | | REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !< discretized artificial dissipation at leftward-side |
1313 | | REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !< discretized 6th-order flux at leftward-side |
1314 | | ! |
1315 | | !-- Used local modified copy of nzb_max (used to degrade order of |
1316 | | !-- discretization) at non-cyclic boundaries. Modify only at relevant points |
1317 | | !-- instead of the entire subdomain. This should lead to better |
1318 | | !-- load balance between boundary and non-boundary PEs. |
1319 | | IF( non_cyclic_l .AND. i <= nxl + 2 .OR. & |
1320 | | non_cyclic_r .AND. i >= nxr - 2 .OR. & |
1321 | | non_cyclic_s .AND. j <= nys + 2 .OR. & |
1322 | | non_cyclic_n .AND. j >= nyn - 2 ) THEN |
1323 | | nzb_max_l = nzt |
1324 | | ELSE |
1325 | | nzb_max_l = nzb_max |
1326 | | END IF |
1327 | | ! |
1328 | | !-- Set control flag for flux limiter |
1329 | | limiter = .FALSE. |
1330 | | IF ( PRESENT( flux_limitation) ) limiter = flux_limitation |
1331 | | ! |
1332 | | !-- Compute southside fluxes of the respective PE bounds. |
1333 | | IF ( j == nys ) THEN |
1334 | | ! |
1335 | | !-- Up to the top of the highest topography. |
1336 | | DO k = nzb+1, nzb_max_l |
1337 | | |
1338 | | ibit5 = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) |
1339 | | ibit4 = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) |
1340 | | ibit3 = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) |
1341 | | |
1342 | | v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) |
1343 | | swap_flux_y_local(k,tn) = v_comp * ( & |
1344 | | ( 37.0_wp * ibit5 * adv_sca_5 & |
1345 | | + 7.0_wp * ibit4 * adv_sca_3 & |
1346 | | + ibit3 * adv_sca_1 & |
1347 | | ) * & |
1348 | | ( sk(k,j,i) + sk(k,j-1,i) ) & |
1349 | | - ( 8.0_wp * ibit5 * adv_sca_5 & |
1350 | | + ibit4 * adv_sca_3 & |
1351 | | ) * & |
1352 | | ( sk(k,j+1,i) + sk(k,j-2,i) ) & |
1353 | | + ( ibit5 * adv_sca_5 & |
1354 | | ) * & |
1355 | | ( sk(k,j+2,i) + sk(k,j-3,i) ) & |
1356 | | ) |
1357 | | |
1358 | | swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( & |
1359 | | ( 10.0_wp * ibit5 * adv_sca_5 & |
1360 | | + 3.0_wp * ibit4 * adv_sca_3 & |
1361 | | + ibit3 * adv_sca_1 & |
1362 | | ) * & |
1363 | | ( sk(k,j,i) - sk(k,j-1,i) ) & |
1364 | | - ( 5.0_wp * ibit5 * adv_sca_5 & |
1365 | | + ibit4 * adv_sca_3 & |
1366 | | ) * & |
1367 | | ( sk(k,j+1,i) - sk(k,j-2,i) ) & |
1368 | | + ( ibit5 * adv_sca_5 & |
1369 | | ) * & |
1370 | | ( sk(k,j+2,i) - sk(k,j-3,i) ) & |
1371 | | ) |
1372 | | |
| 1193 | !--------------------------------------------------------------------------------------------------! |
| 1194 | SUBROUTINE advec_s_ws_ij( advc_flag, i, j, sk, sk_char, swap_flux_y_local, swap_diss_y_local, & |
| 1195 | swap_flux_x_local, swap_diss_x_local, i_omp, tn, non_cyclic_l, & |
| 1196 | non_cyclic_n, non_cyclic_r, non_cyclic_s, flux_limitation ) |
| 1197 | |
| 1198 | |
| 1199 | CHARACTER (LEN = *), INTENT(IN) :: sk_char !< string identifier, used for assign fluxes to the |
| 1200 | !<correct dimension in the analysis array |
| 1201 | |
| 1202 | INTEGER(iwp) :: i !< grid index along x-direction |
| 1203 | INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread |
| 1204 | INTEGER(iwp) :: j !< grid index along y-direction |
| 1205 | INTEGER(iwp) :: k !< grid index along z-direction |
| 1206 | INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults |
| 1207 | INTEGER(iwp) :: k_mmm !< k-3 index in disretization, can be modified to avoid segmentation faults |
| 1208 | INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults |
| 1209 | INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults |
| 1210 | INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms |
| 1211 | INTEGER(iwp) :: tn !< number of OpenMP thread |
| 1212 | |
| 1213 | INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: & |
| 1214 | advc_flag !< flag array to control order of scalar advection |
| 1215 | |
| 1216 | LOGICAL :: limiter !< control flag indicating the application of flux limitation |
| 1217 | LOGICAL :: non_cyclic_l !< flag that indicates non-cyclic boundary on the left |
| 1218 | LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north |
| 1219 | LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right |
| 1220 | LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south |
| 1221 | LOGICAL, OPTIONAL :: flux_limitation !< flag indicating flux limitation of the vertical advection |
| 1222 | |
| 1223 | REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom |
| 1224 | REAL(wp) :: div !< velocity diverence on scalar grid |
| 1225 | REAL(wp) :: div_in !< vertical flux divergence of ingoing fluxes |
| 1226 | REAL(wp) :: div_out !< vertical flux divergence of outgoing fluxes |
| 1227 | REAL(wp) :: f_corr_d !< correction flux at grid-cell bottom, i.e. the difference between high and low-order flux |
| 1228 | REAL(wp) :: f_corr_t !< correction flux at grid-cell top, i.e. the difference between high and low-order flux |
| 1229 | REAL(wp) :: f_corr_d_in !< correction flux of ingoing flux part at grid-cell bottom |
| 1230 | REAL(wp) :: f_corr_t_in !< correction flux of ingoing flux part at grid-cell top |
| 1231 | REAL(wp) :: f_corr_d_out !< correction flux of outgoing flux part at grid-cell bottom |
| 1232 | REAL(wp) :: f_corr_t_out !< correction flux of outgoing flux part at grid-cell top |
| 1233 | REAL(wp) :: fac_correction!< factor to limit the in- and outgoing fluxes |
| 1234 | REAL(wp) :: flux_d !< 6th-order flux at grid box bottom |
| 1235 | REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction |
| 1236 | REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction |
| 1237 | REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction |
| 1238 | REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction |
| 1239 | REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction |
| 1240 | REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction |
| 1241 | REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction |
| 1242 | REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction |
| 1243 | REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction |
| 1244 | REAL(wp) :: max_val !< maximum value of the quanitity along the numerical stencil (in vertical direction) |
| 1245 | REAL(wp) :: min_val !< maximum value of the quanitity along the numerical stencil (in vertical direction) |
| 1246 | REAL(wp) :: mon !< monotone solution of the advection equation using 1st-order fluxes |
| 1247 | REAL(wp) :: u_comp !< advection velocity along x-direction |
| 1248 | REAL(wp) :: v_comp !< advection velocity along y-direction |
| 1249 | |
| 1250 | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side |
| 1251 | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side |
| 1252 | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top |
| 1253 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side |
| 1254 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side |
| 1255 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top |
| 1256 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t_1st !< discretized 1st-order flux at top |
| 1257 | |
| 1258 | REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local !< discretized artificial dissipation at southward-side |
| 1259 | REAL(wp), DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local !< discretized 6th-order flux at northward-side |
| 1260 | REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_diss_x_local !< discretized artificial dissipation at leftward-side |
| 1261 | REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: swap_flux_x_local !< discretized 6th-order flux at leftward-side |
| 1262 | |
| 1263 | ! |
| 1264 | !-- sk is an array from parameter list. It should not be a pointer, because in that case the |
| 1265 | !-- compiler can not assume a stride 1 and cannot perform a strided one vector load. Adding the |
| 1266 | !-- CONTIGUOUS keyword makes things even worse, because the compiler cannot assume strided one in |
| 1267 | !-- the caller side. |
| 1268 | REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar |
| 1269 | |
| 1270 | ! |
| 1271 | !-- Used local modified copy of nzb_max (used to degrade order of discretization) at non-cyclic |
| 1272 | !-- boundaries. Modify only at relevant points instead of the entire subdomain. This should lead to |
| 1273 | !-- betterload balance between boundary and non-boundary PEs. |
| 1274 | IF( non_cyclic_l .AND. i <= nxl + 2 .OR. & |
| 1275 | non_cyclic_r .AND. i >= nxr - 2 .OR. & |
| 1276 | non_cyclic_s .AND. j <= nys + 2 .OR. & |
| 1277 | non_cyclic_n .AND. j >= nyn - 2 ) THEN |
| 1278 | nzb_max_l = nzt |
| 1279 | ELSE |
| 1280 | nzb_max_l = nzb_max |
| 1281 | END IF |
| 1282 | ! |
| 1283 | !-- Set control flag for flux limiter |
| 1284 | limiter = .FALSE. |
| 1285 | IF ( PRESENT( flux_limitation) ) limiter = flux_limitation |
| 1286 | ! |
| 1287 | !-- Compute southside fluxes of the respective PE bounds. |
| 1288 | IF ( j == nys ) THEN |
| 1289 | ! |
| 1290 | !-- Up to the top of the highest topography. |
| 1291 | DO k = nzb+1, nzb_max_l |
| 1292 | |
| 1293 | ibit5 = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) |
| 1294 | ibit4 = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) |
| 1295 | ibit3 = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) |
| 1296 | |
| 1297 | v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) |
| 1298 | swap_flux_y_local(k,tn) = v_comp * ( & |
| 1299 | ( 37.0_wp * ibit5 * adv_sca_5 & |
| 1300 | + 7.0_wp * ibit4 * adv_sca_3 & |
| 1301 | + ibit3 * adv_sca_1 & |
| 1302 | ) * ( sk(k,j,i) + sk(k,j-1,i) ) & |
| 1303 | - ( 8.0_wp * ibit5 * adv_sca_5 & |
| 1304 | + ibit4 * adv_sca_3 & |
| 1305 | ) * ( sk(k,j+1,i) + sk(k,j-2,i) ) & |
| 1306 | + ( ibit5 * adv_sca_5 ) & |
| 1307 | * ( sk(k,j+2,i) + sk(k,j-3,i) ) & |
| 1308 | ) |
| 1309 | |
| 1310 | swap_diss_y_local(k,tn) = - ABS( v_comp ) * ( & |
| 1311 | ( 10.0_wp * ibit5 * adv_sca_5 & |
| 1312 | + 3.0_wp * ibit4 * adv_sca_3 & |
| 1313 | + ibit3 * adv_sca_1 & |
| 1314 | ) * ( sk(k,j,i) - sk(k,j-1,i) ) & |
| 1315 | - ( 5.0_wp * ibit5 * adv_sca_5 & |
| 1316 | + ibit4 * adv_sca_3 & |
| 1317 | ) * ( sk(k,j+1,i) - sk(k,j-2,i) ) & |
| 1318 | + ( ibit5 * adv_sca_5 ) & |
| 1319 | * ( sk(k,j+2,i) - sk(k,j-3,i) ) & |
| 1320 | ) |
| 1321 | |
| 1322 | ENDDO |
| 1323 | ! |
| 1324 | !-- Above to the top of the highest topography. No degradation necessary. |
| 1325 | DO k = nzb_max_l+1, nzt |
| 1326 | |
| 1327 | v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) |
| 1328 | swap_flux_y_local(k,tn) = v_comp * ( 37.0_wp * ( sk(k,j,i) + sk(k,j-1,i) ) & |
| 1329 | - 8.0_wp * ( sk(k,j+1,i) + sk(k,j-2,i) ) & |
| 1330 | + ( sk(k,j+2,i) + sk(k,j-3,i) ) & |
| 1331 | ) * adv_sca_5 |
| 1332 | swap_diss_y_local(k,tn) = - ABS( v_comp ) * ( & |
| 1333 | 10.0_wp * ( sk(k,j,i) - sk(k,j-1,i) ) & |
| 1334 | - 5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) ) & |
| 1335 | + sk(k,j+2,i) - sk(k,j-3,i) & |
| 1336 | ) * adv_sca_5 |
| 1337 | |
| 1338 | ENDDO |
| 1339 | |
| 1340 | ENDIF |
| 1341 | ! |
| 1342 | !-- Compute leftside fluxes of the respective PE bounds. |
| 1343 | IF ( i == i_omp ) THEN |
| 1344 | |
| 1345 | DO k = nzb+1, nzb_max_l |
| 1346 | |
| 1347 | ibit2 = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) |
| 1348 | ibit1 = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) |
| 1349 | ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) |
| 1350 | |
| 1351 | u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) |
| 1352 | swap_flux_x_local(k,j,tn) = u_comp * ( & |
| 1353 | ( 37.0_wp * ibit2 * adv_sca_5 & |
| 1354 | + 7.0_wp * ibit1 * adv_sca_3 & |
| 1355 | + ibit0 * adv_sca_1 & |
| 1356 | ) * ( sk(k,j,i) + sk(k,j,i-1) ) & |
| 1357 | - ( 8.0_wp * ibit2 * adv_sca_5 & |
| 1358 | + ibit1 * adv_sca_3 & |
| 1359 | ) * ( sk(k,j,i+1) + sk(k,j,i-2) ) & |
| 1360 | + ( ibit2 * adv_sca_5 & |
| 1361 | ) * ( sk(k,j,i+2) + sk(k,j,i-3) ) & |
| 1362 | ) |
| 1363 | |
| 1364 | swap_diss_x_local(k,j,tn) = - ABS( u_comp ) * ( & |
| 1365 | ( 10.0_wp * ibit2 * adv_sca_5 & |
| 1366 | + 3.0_wp * ibit1 * adv_sca_3 & |
| 1367 | + ibit0 * adv_sca_1 & |
| 1368 | ) * ( sk(k,j,i) - sk(k,j,i-1) ) & |
| 1369 | - ( 5.0_wp * ibit2 * adv_sca_5 & |
| 1370 | + ibit1 * adv_sca_3 & |
| 1371 | ) * ( sk(k,j,i+1) - sk(k,j,i-2) ) & |
| 1372 | + ( ibit2 * adv_sca_5 & |
| 1373 | ) * ( sk(k,j,i+2) - sk(k,j,i-3) ) & |
| 1374 | ) |
| 1375 | |
| 1376 | ENDDO |
| 1377 | |
| 1378 | DO k = nzb_max_l+1, nzt |
| 1379 | |
| 1380 | u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) |
| 1381 | swap_flux_x_local(k,j,tn) = u_comp * ( & |
| 1382 | 37.0_wp * ( sk(k,j,i) + sk(k,j,i-1) ) & |
| 1383 | - 8.0_wp * ( sk(k,j,i+1) + sk(k,j,i-2) ) & |
| 1384 | + ( sk(k,j,i+2) + sk(k,j,i-3) ) & |
| 1385 | ) * adv_sca_5 |
| 1386 | |
| 1387 | swap_diss_x_local(k,j,tn) = - ABS( u_comp ) * ( & |
| 1388 | 10.0_wp * ( sk(k,j,i) - sk(k,j,i-1) ) & |
| 1389 | - 5.0_wp * ( sk(k,j,i+1) - sk(k,j,i-2) ) & |
| 1390 | + ( sk(k,j,i+2) - sk(k,j,i-3) ) & |
| 1391 | ) * adv_sca_5 |
| 1392 | |
| 1393 | ENDDO |
| 1394 | |
| 1395 | ENDIF |
| 1396 | ! |
| 1397 | !-- Now compute the fluxes for the horizontal termns up to the highest |
| 1398 | !-- topography. |
| 1399 | DO k = nzb+1, nzb_max_l |
| 1400 | |
| 1401 | ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp ) |
| 1402 | ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp ) |
| 1403 | ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp ) |
| 1404 | |
| 1405 | u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) |
| 1406 | flux_r(k) = u_comp * ( & |
| 1407 | ( 37.0_wp * ibit2 * adv_sca_5 & |
| 1408 | + 7.0_wp * ibit1 * adv_sca_3 & |
| 1409 | + ibit0 * adv_sca_1 ) * ( sk(k,j,i+1) + sk(k,j,i) ) & |
| 1410 | - ( 8.0_wp * ibit2 * adv_sca_5 & |
| 1411 | + ibit1 * adv_sca_3 ) * ( sk(k,j,i+2) + sk(k,j,i-1) ) & |
| 1412 | + ( ibit2 * adv_sca_5 ) * ( sk(k,j,i+3) + sk(k,j,i-2) ) & |
| 1413 | ) |
| 1414 | |
| 1415 | diss_r(k) = - ABS( u_comp ) * ( & |
| 1416 | ( 10.0_wp * ibit2 * adv_sca_5 & |
| 1417 | + 3.0_wp * ibit1 * adv_sca_3 & |
| 1418 | + ibit0 * adv_sca_1 ) * ( sk(k,j,i+1) - sk(k,j,i) ) & |
| 1419 | - ( 5.0_wp * ibit2 * adv_sca_5 & |
| 1420 | + ibit1 * adv_sca_3 ) * ( sk(k,j,i+2) - sk(k,j,i-1) ) & |
| 1421 | + ( ibit2 * adv_sca_5 ) * ( sk(k,j,i+3) - sk(k,j,i-2) ) & |
| 1422 | ) |
| 1423 | |
| 1424 | ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp ) |
| 1425 | ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp ) |
| 1426 | ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp ) |
| 1427 | |
| 1428 | v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) |
| 1429 | flux_n(k) = v_comp * ( & |
| 1430 | ( 37.0_wp * ibit5 * adv_sca_5 & |
| 1431 | + 7.0_wp * ibit4 * adv_sca_3 & |
| 1432 | + ibit3 * adv_sca_1 ) * ( sk(k,j+1,i) + sk(k,j,i) ) & |
| 1433 | - ( 8.0_wp * ibit5 * adv_sca_5 & |
| 1434 | + ibit4 * adv_sca_3 ) * ( sk(k,j+2,i) + sk(k,j-1,i) ) & |
| 1435 | + ( ibit5 * adv_sca_5 ) * ( sk(k,j+3,i) + sk(k,j-2,i) ) & |
| 1436 | ) |
| 1437 | |
| 1438 | diss_n(k) = - ABS( v_comp ) * ( & |
| 1439 | ( 10.0_wp * ibit5 * adv_sca_5 & |
| 1440 | + 3.0_wp * ibit4 * adv_sca_3 & |
| 1441 | + ibit3 * adv_sca_1 ) * ( sk(k,j+1,i) - sk(k,j,i) ) & |
| 1442 | - ( 5.0_wp * ibit5 * adv_sca_5 & |
| 1443 | + ibit4 * adv_sca_3 ) * ( sk(k,j+2,i) - sk(k,j-1,i) ) & |
| 1444 | + ( ibit5 * adv_sca_5 ) * ( sk(k,j+3,i) - sk(k,j-2,i) ) & |
| 1445 | ) |
| 1446 | ENDDO |
| 1447 | ! |
| 1448 | !-- Now compute the fluxes for the horizontal terms above the topography |
| 1449 | !-- where no degradation along the horizontal parts is necessary (except |
| 1450 | !-- for the non-cyclic lateral boundaries treated by nzb_max_l). |
| 1451 | DO k = nzb_max_l+1, nzt |
| 1452 | |
| 1453 | u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) |
| 1454 | flux_r(k) = u_comp * ( & |
| 1455 | 37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) & |
| 1456 | - 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) & |
| 1457 | + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5 |
| 1458 | diss_r(k) = - ABS( u_comp ) * ( & |
| 1459 | 10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) & |
| 1460 | - 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) & |
| 1461 | + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 |
| 1462 | |
| 1463 | v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) |
| 1464 | flux_n(k) = v_comp * ( & |
| 1465 | 37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) & |
| 1466 | - 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) & |
| 1467 | + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5 |
| 1468 | diss_n(k) = - ABS( v_comp ) * ( & |
| 1469 | 10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) & |
| 1470 | - 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) & |
| 1471 | + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 |
| 1472 | |
| 1473 | ENDDO |
| 1474 | ! |
| 1475 | !-- Now, compute vertical fluxes. Split loop into a part treating the lowest grid points with |
| 1476 | !-- indirect indexing, a main loop without indirect indexing, and a loop for the uppermost grid |
| 1477 | !-- points with indirect indexing. This allows better vectorization for the main loop. |
| 1478 | !-- First, compute the flux at model surface, which need has to be calculated explicetely for the |
| 1479 | !-- tendency at the first w-level. For topography wall this is done implicitely by advc_flag. |
| 1480 | flux_t(nzb) = 0.0_wp |
| 1481 | diss_t(nzb) = 0.0_wp |
| 1482 | |
| 1483 | DO k = nzb+1, nzb+1 |
| 1484 | ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) |
| 1485 | ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) |
| 1486 | ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) |
| 1487 | ! |
| 1488 | !-- k index has to be modified near bottom and top, else array subscripts will be exceeded. |
| 1489 | k_ppp = k + 3 * ibit8 |
| 1490 | k_pp = k + 2 * ( 1 - ibit6 ) |
| 1491 | k_mm = k - 2 * ibit8 |
| 1492 | |
| 1493 | flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & |
| 1494 | ( 37.0_wp * ibit8 * adv_sca_5 & |
| 1495 | + 7.0_wp * ibit7 * adv_sca_3 & |
| 1496 | + ibit6 * adv_sca_1 ) * ( sk(k+1,j,i) + sk(k,j,i) ) & |
| 1497 | - ( 8.0_wp * ibit8 * adv_sca_5 & |
| 1498 | + ibit7 * adv_sca_3 ) * ( sk(k_pp,j,i) + sk(k-1,j,i) ) & |
| 1499 | + ( ibit8 * adv_sca_5 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & |
| 1500 | ) |
| 1501 | |
| 1502 | diss_t(k) = - ABS( w(k,j,i) ) * rho_air_zw(k) * ( & |
| 1503 | ( 10.0_wp * ibit8 * adv_sca_5 & |
| 1504 | + 3.0_wp * ibit7 * adv_sca_3 & |
| 1505 | + ibit6 * adv_sca_1 ) * ( sk(k+1,j,i) - sk(k,j,i) ) & |
| 1506 | - ( 5.0_wp * ibit8 * adv_sca_5 & |
| 1507 | + ibit7 * adv_sca_3 ) * ( sk(k_pp,j,i) - sk(k-1,j,i) ) & |
| 1508 | + ( ibit8 * adv_sca_5 ) * ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & |
| 1509 | ) |
| 1510 | ENDDO |
| 1511 | |
| 1512 | DO k = nzb+2, nzt-2 |
| 1513 | ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) |
| 1514 | ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) |
| 1515 | ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) |
| 1516 | |
| 1517 | flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & |
| 1518 | ( 37.0_wp * ibit8 * adv_sca_5 & |
| 1519 | + 7.0_wp * ibit7 * adv_sca_3 & |
| 1520 | + ibit6 * adv_sca_1 ) * ( sk(k+1,j,i) + sk(k,j,i) ) & |
| 1521 | - ( 8.0_wp * ibit8 * adv_sca_5 & |
| 1522 | + ibit7 * adv_sca_3 ) * ( sk(k+2,j,i) + sk(k-1,j,i) ) & |
| 1523 | + ( ibit8 * adv_sca_5 ) * ( sk(k+3,j,i) + sk(k-2,j,i) ) & |
| 1524 | ) |
| 1525 | |
| 1526 | diss_t(k) = - ABS( w(k,j,i) ) * rho_air_zw(k) * ( & |
| 1527 | ( 10.0_wp * ibit8 * adv_sca_5 & |
| 1528 | + 3.0_wp * ibit7 * adv_sca_3 & |
| 1529 | + ibit6 * adv_sca_1 ) * ( sk(k+1,j,i) - sk(k,j,i) ) & |
| 1530 | - ( 5.0_wp * ibit8 * adv_sca_5 & |
| 1531 | + ibit7 * adv_sca_3 ) * ( sk(k+2,j,i) - sk(k-1,j,i) ) & |
| 1532 | + ( ibit8 * adv_sca_5 ) * ( sk(k+3,j,i) - sk(k-2,j,i) ) & |
| 1533 | ) |
| 1534 | ENDDO |
| 1535 | |
| 1536 | DO k = nzt-1, nzt-symmetry_flag |
| 1537 | ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) |
| 1538 | ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) |
| 1539 | ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) |
| 1540 | ! |
| 1541 | !-- k index has to be modified near bottom and top, else array subscripts will be exceeded. |
| 1542 | k_ppp = k + 3 * ibit8 |
| 1543 | k_pp = k + 2 * ( 1 - ibit6 ) |
| 1544 | k_mm = k - 2 * ibit8 |
| 1545 | |
| 1546 | |
| 1547 | flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & |
| 1548 | ( 37.0_wp * ibit8 * adv_sca_5 & |
| 1549 | + 7.0_wp * ibit7 * adv_sca_3 & |
| 1550 | + ibit6 * adv_sca_1 ) * ( sk(k+1,j,i) + sk(k,j,i) ) & |
| 1551 | - ( 8.0_wp * ibit8 * adv_sca_5 & |
| 1552 | + ibit7 * adv_sca_3 ) * ( sk(k_pp,j,i) + sk(k-1,j,i) ) & |
| 1553 | + ( ibit8 * adv_sca_5 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & |
| 1554 | ) |
| 1555 | |
| 1556 | diss_t(k) = - ABS( w(k,j,i) ) * rho_air_zw(k) * ( & |
| 1557 | ( 10.0_wp * ibit8 * adv_sca_5 & |
| 1558 | + 3.0_wp * ibit7 * adv_sca_3 & |
| 1559 | + ibit6 * adv_sca_1 ) * ( sk(k+1,j,i) - sk(k,j,i) ) & |
| 1560 | - ( 5.0_wp * ibit8 * adv_sca_5 & |
| 1561 | + ibit7 * adv_sca_3 ) * ( sk(k_pp,j,i) - sk(k-1,j,i) ) & |
| 1562 | + ( ibit8 * adv_sca_5 ) * ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & |
| 1563 | ) |
| 1564 | ENDDO |
| 1565 | |
| 1566 | ! |
| 1567 | !-- Set resolved/turbulent flux at model top to zero (w-level). In case that a symmetric behavior |
| 1568 | !-- between bottom and top shall be guaranteed (closed channel flow), the flux at nzt is also set to |
| 1569 | !-- zero. |
| 1570 | IF ( symmetry_flag == 1 ) THEN |
| 1571 | flux_t(nzt) = 0.0_wp |
| 1572 | diss_t(nzt) = 0.0_wp |
| 1573 | ENDIF |
| 1574 | flux_t(nzt+1) = 0.0_wp |
| 1575 | diss_t(nzt+1) = 0.0_wp |
| 1576 | |
| 1577 | |
| 1578 | IF ( limiter ) THEN |
| 1579 | ! |
| 1580 | !-- Compute monotone first-order fluxes which are required for mononteflux limitation. |
| 1581 | flux_t_1st(nzb) = 0.0_wp |
| 1582 | DO k = nzb+1, nzb_max_l |
| 1583 | flux_t_1st(k) = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & |
| 1584 | - ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) ) & |
| 1585 | * rho_air_zw(k) * adv_sca_1 |
| 1586 | ! |
| 1587 | !-- In flux limitation the total flux will be corrected. For the sake of cleariness the |
| 1588 | !-- higher-order advective and disspative fluxes will be merged onto flux_t. |
| 1589 | flux_t(k) = flux_t(k) + diss_t(k) |
| 1590 | diss_t(k) = 0.0_wp |
| 1591 | ENDDO |
| 1592 | ! |
| 1593 | !-- Flux limitation of vertical fluxes according to Skamarock (2006). |
| 1594 | !-- Please note, as flux limitation implies linear dependencies of fluxes, flux limitation is |
| 1595 | !-- only made for the vertical advection term. Limitation of the horizontal terms cannot be |
| 1596 | !-- parallelized. |
| 1597 | !-- Due to the linear dependency, the following loop will not be vectorized. |
| 1598 | !-- Further, note that the flux limiter is only applied within the urban layer, i.e up to the |
| 1599 | !-- topography top. |
| 1600 | DO k = nzb+1, nzb_max_l |
| 1601 | ! |
| 1602 | !-- Compute one-dimensional divergence along the vertical direction, which is used to correct |
| 1603 | !-- the advection discretization. This is necessary as in one-dimensional space the advection |
| 1604 | !-- velocity should actually be constant. |
| 1605 | div = ( w(k,j,i) * rho_air_zw(k) & |
| 1606 | - w(k-1,j,i) * rho_air_zw(k-1) & |
| 1607 | ) * drho_air(k) * ddzw(k) |
| 1608 | ! |
| 1609 | !-- Compute monotone solution of the advection equation from 1st-order fluxes. Please note, |
| 1610 | !-- the advection equation is corrected by the divergence term (in 1D the advective flow |
| 1611 | !-- should be divergence free). Moreover, please note, as time-increment the full timestep is |
| 1612 | !-- used, even though a Runge-Kutta scheme will be used. However, the length of the actual |
| 1613 | !-- time increment is not important at all since it cancels out later when the fluxes are |
| 1614 | !-- limited. |
| 1615 | mon = sk(k,j,i) + ( - ( flux_t_1st(k) - flux_t_1st(k-1) ) & |
| 1616 | * drho_air(k) * ddzw(k) & |
| 1617 | + div * sk(k,j,i) & |
| 1618 | ) * dt_3d |
| 1619 | ! |
| 1620 | !-- Determine minimum and maximum values along the numerical stencil. |
| 1621 | k_mmm = MAX( k - 3, nzb + 1 ) |
| 1622 | k_ppp = MIN( k + 3, nzt + 1 ) |
| 1623 | |
| 1624 | min_val = MINVAL( sk(k_mmm:k_ppp,j,i) ) |
| 1625 | max_val = MAXVAL( sk(k_mmm:k_ppp,j,i) ) |
| 1626 | ! |
| 1627 | !-- Compute difference between high- and low-order fluxes, which may act as correction fluxes |
| 1628 | f_corr_t = flux_t(k) - flux_t_1st(k) |
| 1629 | f_corr_d = flux_t(k-1) - flux_t_1st(k-1) |
| 1630 | ! |
| 1631 | !-- Determine outgoing fluxes, i.e. the part of the fluxes which can decrease the value within |
| 1632 | !-- the grid box |
| 1633 | f_corr_t_out = MAX( 0.0_wp, f_corr_t ) |
| 1634 | f_corr_d_out = MIN( 0.0_wp, f_corr_d ) |
| 1635 | ! |
| 1636 | !-- Determine ingoing fluxes, i.e. the part of the fluxes which can increase the value within |
| 1637 | !-- the grid box |
| 1638 | f_corr_t_in = MIN( 0.0_wp, f_corr_t) |
| 1639 | f_corr_d_in = MAX( 0.0_wp, f_corr_d) |
| 1640 | ! |
| 1641 | !-- Compute divergence of outgoing correction fluxes |
| 1642 | div_out = - ( f_corr_t_out - f_corr_d_out ) * drho_air(k) * ddzw(k) * dt_3d |
| 1643 | ! |
| 1644 | !-- Compute divergence of ingoing correction fluxes |
| 1645 | div_in = - ( f_corr_t_in - f_corr_d_in ) * drho_air(k) * ddzw(k) * dt_3d |
| 1646 | ! |
| 1647 | !-- Check if outgoing fluxes can lead to undershoots, i.e. values smaller than the minimum |
| 1648 | !-- value within the numerical stencil. If so, limit them. |
| 1649 | IF ( mon - min_val < - div_out .AND. ABS( div_out ) > 0.0_wp ) THEN |
| 1650 | fac_correction = ( mon - min_val ) / ( - div_out ) |
| 1651 | f_corr_t_out = f_corr_t_out * fac_correction |
| 1652 | f_corr_d_out = f_corr_d_out * fac_correction |
| 1653 | ENDIF |
| 1654 | ! |
| 1655 | !-- Check if ingoing fluxes can lead to overshoots, i.e. values larger than the maximum value |
| 1656 | !-- within the numerical stencil. If so, limit them. |
| 1657 | IF ( mon - max_val > - div_in .AND. ABS( div_in ) > 0.0_wp ) THEN |
| 1658 | fac_correction = ( mon - max_val ) / ( - div_in ) |
| 1659 | f_corr_t_in = f_corr_t_in * fac_correction |
| 1660 | f_corr_d_in = f_corr_d_in * fac_correction |
| 1661 | ENDIF |
| 1662 | ! |
| 1663 | !-- Finally add the limited fluxes to the original ones. If no flux limitation was done, the |
| 1664 | !-- fluxes equal the original ones. |
| 1665 | flux_t(k) = flux_t_1st(k) + f_corr_t_out + f_corr_t_in |
| 1666 | flux_t(k-1) = flux_t_1st(k-1) + f_corr_d_out + f_corr_d_in |
| 1667 | ENDDO |
| 1668 | ENDIF |
| 1669 | ! |
| 1670 | !-- Now compute the tendency term including divergence correction. |
| 1671 | DO k = nzb+1, nzb_max_l |
| 1672 | |
| 1673 | flux_d = flux_t(k-1) |
| 1674 | diss_d = diss_t(k-1) |
| 1675 | |
| 1676 | ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp ) |
| 1677 | ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp ) |
| 1678 | ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp ) |
| 1679 | |
| 1680 | ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp ) |
| 1681 | ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp ) |
| 1682 | ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp ) |
| 1683 | |
| 1684 | ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) |
| 1685 | ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) |
| 1686 | ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) |
| 1687 | ! |
| 1688 | !-- Calculate the divergence of the velocity field. A respective correction is needed to overcome |
| 1689 | !-- numerical instabilities introduced by an insufficient reduction of divergences near |
| 1690 | !-- topography. |
| 1691 | div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & |
| 1692 | - u(k,j,i) * ( & |
| 1693 | REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) & |
| 1694 | + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) & |
| 1695 | + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) & |
| 1696 | ) & |
| 1697 | ) * ddx & |
| 1698 | + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & |
| 1699 | - v(k,j,i) * ( & |
| 1700 | REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) & |
| 1701 | + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) & |
| 1702 | + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) & |
| 1703 | ) & |
| 1704 | ) * ddy & |
| 1705 | + ( w(k,j,i) * rho_air_zw(k) * ( ibit6 + ibit7 + ibit8 ) & |
| 1706 | - w(k-1,j,i) * rho_air_zw(k-1) * & |
| 1707 | ( & |
| 1708 | REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp ) & |
| 1709 | + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp ) & |
| 1710 | + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp ) & |
| 1711 | ) & |
| 1712 | ) * drho_air(k) * ddzw(k) |
| 1713 | |
| 1714 | tend(k,j,i) = tend(k,j,i) - ( & |
| 1715 | ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) & |
| 1716 | - swap_diss_x_local(k,j,tn) ) * ddx & |
| 1717 | + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) & |
| 1718 | - swap_diss_y_local(k,tn) ) * ddy & |
| 1719 | + ( ( flux_t(k) + diss_t(k) ) - ( flux_d + diss_d ) & |
| 1720 | ) * drho_air(k) * ddzw(k) & |
| 1721 | ) + sk(k,j,i) * div |
| 1722 | |
| 1723 | |
| 1724 | swap_flux_y_local(k,tn) = flux_n(k) |
| 1725 | swap_diss_y_local(k,tn) = diss_n(k) |
| 1726 | swap_flux_x_local(k,j,tn) = flux_r(k) |
| 1727 | swap_diss_x_local(k,j,tn) = diss_r(k) |
| 1728 | |
| 1729 | ENDDO |
| 1730 | |
| 1731 | DO k = nzb_max_l+1, nzt |
| 1732 | |
| 1733 | flux_d = flux_t(k-1) |
| 1734 | diss_d = diss_t(k-1) |
| 1735 | ! |
| 1736 | !-- Calculate the divergence of the velocity field. A respective correction is needed to overcome |
| 1737 | !-- numerical instabilities introduced by an insufficient reduction of divergences near |
| 1738 | !-- topography. |
| 1739 | div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & |
| 1740 | + ( v(k,j+1,i) - v(k,j,i) ) * ddy & |
| 1741 | + ( w(k,j,i) * rho_air_zw(k) & |
| 1742 | - w(k-1,j,i) * rho_air_zw(k-1) & |
| 1743 | ) * drho_air(k) * ddzw(k) |
| 1744 | |
| 1745 | tend(k,j,i) = tend(k,j,i) - ( & |
| 1746 | ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) & |
| 1747 | - swap_diss_x_local(k,j,tn) ) * ddx & |
| 1748 | + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) & |
| 1749 | - swap_diss_y_local(k,tn) ) * ddy & |
| 1750 | + ( ( flux_t(k) + diss_t(k) ) - ( flux_d + diss_d ) & |
| 1751 | ) * drho_air(k) * ddzw(k) & |
| 1752 | ) + sk(k,j,i) * div |
| 1753 | |
| 1754 | |
| 1755 | swap_flux_y_local(k,tn) = flux_n(k) |
| 1756 | swap_diss_y_local(k,tn) = diss_n(k) |
| 1757 | swap_flux_x_local(k,j,tn) = flux_r(k) |
| 1758 | swap_diss_x_local(k,j,tn) = diss_r(k) |
| 1759 | |
| 1760 | ENDDO |
| 1761 | |
| 1762 | ! |
| 1763 | !-- Evaluation of statistics. |
| 1764 | SELECT CASE ( sk_char ) |
| 1765 | |
| 1766 | CASE ( 'pt' ) |
| 1767 | |
| 1768 | DO k = nzb, nzt |
| 1769 | sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + & |
| 1770 | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
| 1771 | * ( w(k,j,i) - hom(k,1,3,0) ) & |
| 1772 | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
| 1773 | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
| 1774 | ) * weight_substep(intermediate_timestep_count) |
1453 | | ENDIF |
1454 | | ! |
1455 | | !-- Now compute the fluxes for the horizontal termns up to the highest |
1456 | | !-- topography. |
1457 | | DO k = nzb+1, nzb_max_l |
1458 | | |
1459 | | ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp ) |
1460 | | ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp ) |
1461 | | ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp ) |
1462 | | |
1463 | | u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) |
1464 | | flux_r(k) = u_comp * ( & |
1465 | | ( 37.0_wp * ibit2 * adv_sca_5 & |
1466 | | + 7.0_wp * ibit1 * adv_sca_3 & |
1467 | | + ibit0 * adv_sca_1 & |
1468 | | ) * & |
1469 | | ( sk(k,j,i+1) + sk(k,j,i) ) & |
1470 | | - ( 8.0_wp * ibit2 * adv_sca_5 & |
1471 | | + ibit1 * adv_sca_3 & |
1472 | | ) * & |
1473 | | ( sk(k,j,i+2) + sk(k,j,i-1) ) & |
1474 | | + ( ibit2 * adv_sca_5 & |
1475 | | ) * & |
1476 | | ( sk(k,j,i+3) + sk(k,j,i-2) ) & |
1477 | | ) |
1478 | | |
1479 | | diss_r(k) = -ABS( u_comp ) * ( & |
1480 | | ( 10.0_wp * ibit2 * adv_sca_5 & |
1481 | | + 3.0_wp * ibit1 * adv_sca_3 & |
1482 | | + ibit0 * adv_sca_1 & |
1483 | | ) * & |
1484 | | ( sk(k,j,i+1) - sk(k,j,i) ) & |
1485 | | - ( 5.0_wp * ibit2 * adv_sca_5 & |
1486 | | + ibit1 * adv_sca_3 & |
1487 | | ) * & |
1488 | | ( sk(k,j,i+2) - sk(k,j,i-1) ) & |
1489 | | + ( ibit2 * adv_sca_5 & |
1490 | | ) * & |
1491 | | ( sk(k,j,i+3) - sk(k,j,i-2) ) & |
1492 | | ) |
1493 | | |
1494 | | ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp ) |
1495 | | ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp ) |
1496 | | ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp ) |
1497 | | |
1498 | | v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) |
1499 | | flux_n(k) = v_comp * ( & |
1500 | | ( 37.0_wp * ibit5 * adv_sca_5 & |
1501 | | + 7.0_wp * ibit4 * adv_sca_3 & |
1502 | | + ibit3 * adv_sca_1 & |
1503 | | ) * & |
1504 | | ( sk(k,j+1,i) + sk(k,j,i) ) & |
1505 | | - ( 8.0_wp * ibit5 * adv_sca_5 & |
1506 | | + ibit4 * adv_sca_3 & |
1507 | | ) * & |
1508 | | ( sk(k,j+2,i) + sk(k,j-1,i) ) & |
1509 | | + ( ibit5 * adv_sca_5 & |
1510 | | ) * & |
1511 | | ( sk(k,j+3,i) + sk(k,j-2,i) ) & |
1512 | | ) |
1513 | | |
1514 | | diss_n(k) = -ABS( v_comp ) * ( & |
1515 | | ( 10.0_wp * ibit5 * adv_sca_5 & |
1516 | | + 3.0_wp * ibit4 * adv_sca_3 & |
1517 | | + ibit3 * adv_sca_1 & |
1518 | | ) * & |
1519 | | ( sk(k,j+1,i) - sk(k,j,i) ) & |
1520 | | - ( 5.0_wp * ibit5 * adv_sca_5 & |
1521 | | + ibit4 * adv_sca_3 & |
1522 | | ) * & |
1523 | | ( sk(k,j+2,i) - sk(k,j-1,i) ) & |
1524 | | + ( ibit5 * adv_sca_5 & |
1525 | | ) * & |
1526 | | ( sk(k,j+3,i) - sk(k,j-2,i) ) & |
1527 | | ) |
1528 | | ENDDO |
1529 | | ! |
1530 | | !-- Now compute the fluxes for the horizontal terms above the topography |
1531 | | !-- where no degradation along the horizontal parts is necessary (except |
1532 | | !-- for the non-cyclic lateral boundaries treated by nzb_max_l). |
1533 | | DO k = nzb_max_l+1, nzt |
1534 | | |
1535 | | u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) |
1536 | | flux_r(k) = u_comp * ( & |
1537 | | 37.0_wp * ( sk(k,j,i+1) + sk(k,j,i) ) & |
1538 | | - 8.0_wp * ( sk(k,j,i+2) + sk(k,j,i-1) ) & |
1539 | | + ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5 |
1540 | | diss_r(k) = -ABS( u_comp ) * ( & |
1541 | | 10.0_wp * ( sk(k,j,i+1) - sk(k,j,i) ) & |
1542 | | - 5.0_wp * ( sk(k,j,i+2) - sk(k,j,i-1) ) & |
1543 | | + ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5 |
1544 | | |
1545 | | v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) |
1546 | | flux_n(k) = v_comp * ( & |
1547 | | 37.0_wp * ( sk(k,j+1,i) + sk(k,j,i) ) & |
1548 | | - 8.0_wp * ( sk(k,j+2,i) + sk(k,j-1,i) ) & |
1549 | | + ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5 |
1550 | | diss_n(k) = -ABS( v_comp ) * ( & |
1551 | | 10.0_wp * ( sk(k,j+1,i) - sk(k,j,i) ) & |
1552 | | - 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) & |
1553 | | + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 |
1554 | | |
1555 | | ENDDO |
1556 | | ! |
1557 | | !-- Now, compute vertical fluxes. Split loop into a part treating the |
1558 | | !-- lowest grid points with indirect indexing, a main loop without |
1559 | | !-- indirect indexing, and a loop for the uppermost grip points with |
1560 | | !-- indirect indexing. This allows better vectorization for the main loop. |
1561 | | !-- First, compute the flux at model surface, which need has to be |
1562 | | !-- calculated explicetely for the tendency at |
1563 | | !-- the first w-level. For topography wall this is done implicitely by |
1564 | | !-- advc_flag. |
1565 | | flux_t(nzb) = 0.0_wp |
1566 | | diss_t(nzb) = 0.0_wp |
1567 | | |
1568 | | DO k = nzb+1, nzb+1 |
1569 | | ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) |
1570 | | ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) |
1571 | | ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) |
1572 | | ! |
1573 | | !-- k index has to be modified near bottom and top, else array |
1574 | | !-- subscripts will be exceeded. |
1575 | | k_ppp = k + 3 * ibit8 |
1576 | | k_pp = k + 2 * ( 1 - ibit6 ) |
1577 | | k_mm = k - 2 * ibit8 |
1578 | | |
1579 | | flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & |
1580 | | ( 37.0_wp * ibit8 * adv_sca_5 & |
1581 | | + 7.0_wp * ibit7 * adv_sca_3 & |
1582 | | + ibit6 * adv_sca_1 & |
1583 | | ) * & |
1584 | | ( sk(k+1,j,i) + sk(k,j,i) ) & |
1585 | | - ( 8.0_wp * ibit8 * adv_sca_5 & |
1586 | | + ibit7 * adv_sca_3 & |
1587 | | ) * & |
1588 | | ( sk(k_pp,j,i) + sk(k-1,j,i) ) & |
1589 | | + ( ibit8 * adv_sca_5 & |
1590 | | ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & |
1591 | | ) |
1592 | | |
1593 | | diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & |
1594 | | ( 10.0_wp * ibit8 * adv_sca_5 & |
1595 | | + 3.0_wp * ibit7 * adv_sca_3 & |
1596 | | + ibit6 * adv_sca_1 & |
1597 | | ) * & |
1598 | | ( sk(k+1,j,i) - sk(k,j,i) ) & |
1599 | | - ( 5.0_wp * ibit8 * adv_sca_5 & |
1600 | | + ibit7 * adv_sca_3 & |
1601 | | ) * & |
1602 | | ( sk(k_pp,j,i) - sk(k-1,j,i) ) & |
1603 | | + ( ibit8 * adv_sca_5 & |
1604 | | ) * & |
1605 | | ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & |
1606 | | ) |
1607 | | ENDDO |
1608 | | |
1609 | | DO k = nzb+2, nzt-2 |
1610 | | ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) |
1611 | | ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) |
1612 | | ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) |
1613 | | |
1614 | | flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & |
1615 | | ( 37.0_wp * ibit8 * adv_sca_5 & |
1616 | | + 7.0_wp * ibit7 * adv_sca_3 & |
1617 | | + ibit6 * adv_sca_1 & |
1618 | | ) * & |
1619 | | ( sk(k+1,j,i) + sk(k,j,i) ) & |
1620 | | - ( 8.0_wp * ibit8 * adv_sca_5 & |
1621 | | + ibit7 * adv_sca_3 & |
1622 | | ) * & |
1623 | | ( sk(k+2,j,i) + sk(k-1,j,i) ) & |
1624 | | + ( ibit8 * adv_sca_5 & |
1625 | | ) * ( sk(k+3,j,i)+ sk(k-2,j,i) ) & |
1626 | | ) |
1627 | | |
1628 | | diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & |
1629 | | ( 10.0_wp * ibit8 * adv_sca_5 & |
1630 | | + 3.0_wp * ibit7 * adv_sca_3 & |
1631 | | + ibit6 * adv_sca_1 & |
1632 | | ) * & |
1633 | | ( sk(k+1,j,i) - sk(k,j,i) ) & |
1634 | | - ( 5.0_wp * ibit8 * adv_sca_5 & |
1635 | | + ibit7 * adv_sca_3 & |
1636 | | ) * & |
1637 | | ( sk(k+2,j,i) - sk(k-1,j,i) ) & |
1638 | | + ( ibit8 * adv_sca_5 & |
1639 | | ) * & |
1640 | | ( sk(k+3,j,i) - sk(k-2,j,i) ) & |
1641 | | ) |
1642 | | ENDDO |
1643 | | |
1644 | | DO k = nzt-1, nzt-symmetry_flag |
1645 | | ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) |
1646 | | ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) |
1647 | | ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) |
1648 | | ! |
1649 | | !-- k index has to be modified near bottom and top, else array |
1650 | | !-- subscripts will be exceeded. |
1651 | | k_ppp = k + 3 * ibit8 |
1652 | | k_pp = k + 2 * ( 1 - ibit6 ) |
1653 | | k_mm = k - 2 * ibit8 |
1654 | | |
1655 | | |
1656 | | flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & |
1657 | | ( 37.0_wp * ibit8 * adv_sca_5 & |
1658 | | + 7.0_wp * ibit7 * adv_sca_3 & |
1659 | | + ibit6 * adv_sca_1 & |
1660 | | ) * & |
1661 | | ( sk(k+1,j,i) + sk(k,j,i) ) & |
1662 | | - ( 8.0_wp * ibit8 * adv_sca_5 & |
1663 | | + ibit7 * adv_sca_3 & |
1664 | | ) * & |
1665 | | ( sk(k_pp,j,i) + sk(k-1,j,i) ) & |
1666 | | + ( ibit8 * adv_sca_5 & |
1667 | | ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & |
1668 | | ) |
1669 | | |
1670 | | diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & |
1671 | | ( 10.0_wp * ibit8 * adv_sca_5 & |
1672 | | + 3.0_wp * ibit7 * adv_sca_3 & |
1673 | | + ibit6 * adv_sca_1 & |
1674 | | ) * & |
1675 | | ( sk(k+1,j,i) - sk(k,j,i) ) & |
1676 | | - ( 5.0_wp * ibit8 * adv_sca_5 & |
1677 | | + ibit7 * adv_sca_3 & |
1678 | | ) * & |
1679 | | ( sk(k_pp,j,i) - sk(k-1,j,i) ) & |
1680 | | + ( ibit8 * adv_sca_5 & |
1681 | | ) * & |
1682 | | ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & |
1683 | | ) |
1684 | | ENDDO |
1685 | | |
1686 | | ! |
1687 | | !-- Set resolved/turbulent flux at model top to zero (w-level). In case that |
1688 | | !-- a symmetric behavior between bottom and top shall be guaranteed (closed |
1689 | | !-- channel flow), the flux at nzt is also set to zero. |
1690 | | IF ( symmetry_flag == 1 ) THEN |
1691 | | flux_t(nzt) = 0.0_wp |
1692 | | diss_t(nzt) = 0.0_wp |
1693 | | ENDIF |
1694 | | flux_t(nzt+1) = 0.0_wp |
1695 | | diss_t(nzt+1) = 0.0_wp |
1696 | | |
1697 | | |
1698 | | IF ( limiter ) THEN |
1699 | | ! |
1700 | | !-- Compute monotone first-order fluxes which are required for mononte |
1701 | | !-- flux limitation. |
1702 | | flux_t_1st(nzb) = 0.0_wp |
1703 | | DO k = nzb+1, nzb_max_l |
1704 | | flux_t_1st(k) = ( w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) & |
1705 | | - ABS( w(k,j,i) ) * ( sk(k+1,j,i) - sk(k,j,i) ) ) & |
1706 | | * rho_air_zw(k) * adv_sca_1 |
1707 | | ! |
1708 | | !-- In flux limitation the total flux will be corrected. For the sake |
1709 | | !-- of cleariness the higher-order advective and disspative fluxes |
1710 | | !-- will be merged onto flux_t. |
1711 | | flux_t(k) = flux_t(k) + diss_t(k) |
1712 | | diss_t(k) = 0.0_wp |
| 1810 | |
| 1811 | CASE ( 'qi' ) |
| 1812 | |
| 1813 | DO k = nzb, nzt |
| 1814 | sums_wsqis_ws_l(k,tn) = sums_wsqis_ws_l(k,tn) + & |
| 1815 | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
| 1816 | * ( w(k,j,i) - hom(k,1,3,0) ) & |
| 1817 | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
| 1818 | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
| 1819 | ) * weight_substep(intermediate_timestep_count) |
1799 | | ENDIF |
1800 | | ! |
1801 | | !-- Now compute the tendency term including divergence correction. |
1802 | | DO k = nzb+1, nzb_max_l |
1803 | | |
1804 | | flux_d = flux_t(k-1) |
1805 | | diss_d = diss_t(k-1) |
1806 | | |
1807 | | ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp ) |
1808 | | ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp ) |
1809 | | ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp ) |
1810 | | |
1811 | | ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp ) |
1812 | | ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp ) |
1813 | | ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp ) |
1814 | | |
1815 | | ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) |
1816 | | ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) |
1817 | | ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) |
1818 | | ! |
1819 | | !-- Calculate the divergence of the velocity field. A respective |
1820 | | !-- correction is needed to overcome numerical instabilities introduced |
1821 | | !-- by a not sufficient reduction of divergences near topography. |
1822 | | div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & |
1823 | | - u(k,j,i) * ( & |
1824 | | REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) & |
1825 | | + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) & |
1826 | | + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) & |
1827 | | ) & |
1828 | | ) * ddx & |
1829 | | + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & |
1830 | | - v(k,j,i) * ( & |
1831 | | REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) & |
1832 | | + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) & |
1833 | | + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) & |
1834 | | ) & |
1835 | | ) * ddy & |
1836 | | + ( w(k,j,i) * rho_air_zw(k) * & |
1837 | | ( ibit6 + ibit7 + ibit8 ) & |
1838 | | - w(k-1,j,i) * rho_air_zw(k-1) * & |
1839 | | ( & |
1840 | | REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp ) & |
1841 | | + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp ) & |
1842 | | + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp ) & |
1843 | | ) & |
1844 | | ) * drho_air(k) * ddzw(k) |
1845 | | |
1846 | | tend(k,j,i) = tend(k,j,i) - ( & |
1847 | | ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & |
1848 | | swap_diss_x_local(k,j,tn) ) * ddx & |
1849 | | + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & |
1850 | | swap_diss_y_local(k,tn) ) * ddy & |
1851 | | + ( ( flux_t(k) + diss_t(k) ) - & |
1852 | | ( flux_d + diss_d ) & |
1853 | | ) * drho_air(k) * ddzw(k) & |
1854 | | ) + sk(k,j,i) * div |
1855 | | |
1856 | | |
1857 | | swap_flux_y_local(k,tn) = flux_n(k) |
1858 | | swap_diss_y_local(k,tn) = diss_n(k) |
1859 | | swap_flux_x_local(k,j,tn) = flux_r(k) |
1860 | | swap_diss_x_local(k,j,tn) = diss_r(k) |
1861 | | |
1862 | | ENDDO |
1863 | | |
1864 | | DO k = nzb_max_l+1, nzt |
1865 | | |
1866 | | flux_d = flux_t(k-1) |
1867 | | diss_d = diss_t(k-1) |
1868 | | ! |
1869 | | !-- Calculate the divergence of the velocity field. A respective |
1870 | | !-- correction is needed to overcome numerical instabilities introduced |
1871 | | !-- by a not sufficient reduction of divergences near topography. |
1872 | | div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & |
1873 | | + ( v(k,j+1,i) - v(k,j,i) ) * ddy & |
1874 | | + ( w(k,j,i) * rho_air_zw(k) & |
1875 | | - w(k-1,j,i) * rho_air_zw(k-1) & |
1876 | | ) * drho_air(k) * ddzw(k) |
1877 | | |
1878 | | tend(k,j,i) = tend(k,j,i) - ( & |
1879 | | ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - & |
1880 | | swap_diss_x_local(k,j,tn) ) * ddx & |
1881 | | + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - & |
1882 | | swap_diss_y_local(k,tn) ) * ddy & |
1883 | | + ( ( flux_t(k) + diss_t(k) ) - & |
1884 | | ( flux_d + diss_d ) & |
1885 | | ) * drho_air(k) * ddzw(k) & |
1886 | | ) + sk(k,j,i) * div |
1887 | | |
1888 | | |
1889 | | swap_flux_y_local(k,tn) = flux_n(k) |
1890 | | swap_diss_y_local(k,tn) = diss_n(k) |
1891 | | swap_flux_x_local(k,j,tn) = flux_r(k) |
1892 | | swap_diss_x_local(k,j,tn) = diss_r(k) |
1893 | | |
1894 | | ENDDO |
1895 | | |
1896 | | ! |
1897 | | !-- Evaluation of statistics. |
1898 | | SELECT CASE ( sk_char ) |
1899 | | |
1900 | | CASE ( 'pt' ) |
1901 | | |
1902 | | DO k = nzb, nzt |
1903 | | sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + & |
1904 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
1905 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
1906 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
1907 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
1908 | | ) * weight_substep(intermediate_timestep_count) |
1909 | | ENDDO |
1910 | | |
1911 | | CASE ( 'sa' ) |
1912 | | |
1913 | | DO k = nzb, nzt |
1914 | | sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + & |
1915 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
1916 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
1917 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
1918 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
1919 | | ) * weight_substep(intermediate_timestep_count) |
1920 | | ENDDO |
1921 | | |
1922 | | CASE ( 'q' ) |
1923 | | |
1924 | | DO k = nzb, nzt |
1925 | | sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + & |
1926 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
1927 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
1928 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
1929 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
1930 | | ) * weight_substep(intermediate_timestep_count) |
1931 | | ENDDO |
1932 | | |
1933 | | CASE ( 'qc' ) |
1934 | | |
1935 | | DO k = nzb, nzt |
1936 | | sums_wsqcs_ws_l(k,tn) = sums_wsqcs_ws_l(k,tn) + & |
1937 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
1938 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
1939 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
1940 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
1941 | | ) * weight_substep(intermediate_timestep_count) |
1942 | | ENDDO |
1943 | | |
1944 | | CASE ( 'qi' ) |
1945 | | |
1946 | | DO k = nzb, nzt |
1947 | | sums_wsqis_ws_l(k,tn) = sums_wsqis_ws_l(k,tn) + & |
1948 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
1949 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
1950 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
1951 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
1952 | | ) * weight_substep(intermediate_timestep_count) |
1953 | | ENDDO |
1954 | | |
1955 | | CASE ( 'qr' ) |
1956 | | |
1957 | | DO k = nzb, nzt |
1958 | | sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + & |
1959 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
1960 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
1961 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
1962 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
1963 | | ) * weight_substep(intermediate_timestep_count) |
1964 | | ENDDO |
1965 | | |
1966 | | CASE ( 'nc' ) |
1967 | | |
1968 | | DO k = nzb, nzt |
1969 | | sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) + & |
1970 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
1971 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
1972 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
1973 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
1974 | | ) * weight_substep(intermediate_timestep_count) |
1975 | | ENDDO |
1976 | | |
1977 | | CASE ( 'ni' ) |
1978 | | |
1979 | | DO k = nzb, nzt |
1980 | | sums_wsnis_ws_l(k,tn) = sums_wsnis_ws_l(k,tn) + & |
1981 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
1982 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
1983 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
1984 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
1985 | | ) * weight_substep(intermediate_timestep_count) |
1986 | | ENDDO |
1987 | | |
1988 | | CASE ( 'nr' ) |
1989 | | |
1990 | | DO k = nzb, nzt |
1991 | | sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + & |
1992 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
1993 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
1994 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
1995 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
1996 | | ) * weight_substep(intermediate_timestep_count) |
1997 | | ENDDO |
1998 | | |
1999 | | CASE ( 's' ) |
2000 | | |
2001 | | DO k = nzb, nzt |
2002 | | sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) + & |
2003 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
2004 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
2005 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
2006 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
2007 | | ) * weight_substep(intermediate_timestep_count) |
2008 | | ENDDO |
2009 | | |
2010 | | CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' ) |
2011 | | |
2012 | | DO k = nzb, nzt |
2013 | | sums_salsa_ws_l(k,tn) = sums_salsa_ws_l(k,tn) + & |
2014 | | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
2015 | | * ( w(k,j,i) - hom(k,1,3,0) ) & |
2016 | | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
2017 | | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
2018 | | ) * weight_substep(intermediate_timestep_count) |
2019 | | ENDDO |
2020 | | |
2021 | | ! CASE ( 'kc' ) |
2022 | | !kk Has to be implemented for kpp chemistry |
2023 | | |
2024 | | END SELECT |
2025 | | |
2026 | | END SUBROUTINE advec_s_ws_ij |
2027 | | |
2028 | | |
2029 | | |
2030 | | |
2031 | | !------------------------------------------------------------------------------! |
| 1832 | |
| 1833 | CASE ( 'nc' ) |
| 1834 | |
| 1835 | DO k = nzb, nzt |
| 1836 | sums_wsncs_ws_l(k,tn) = sums_wsncs_ws_l(k,tn) + & |
| 1837 | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
| 1838 | * ( w(k,j,i) - hom(k,1,3,0) ) & |
| 1839 | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
| 1840 | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
| 1841 | ) * weight_substep(intermediate_timestep_count) |
| 1842 | ENDDO |
| 1843 | |
| 1844 | CASE ( 'ni' ) |
| 1845 | |
| 1846 | DO k = nzb, nzt |
| 1847 | sums_wsnis_ws_l(k,tn) = sums_wsnis_ws_l(k,tn) + & |
| 1848 | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
| 1849 | * ( w(k,j,i) - hom(k,1,3,0) ) & |
| 1850 | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
| 1851 | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
| 1852 | ) * weight_substep(intermediate_timestep_count) |
| 1853 | ENDDO |
| 1854 | |
| 1855 | CASE ( 'nr' ) |
| 1856 | |
| 1857 | DO k = nzb, nzt |
| 1858 | sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + & |
| 1859 | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
| 1860 | * ( w(k,j,i) - hom(k,1,3,0) ) & |
| 1861 | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
| 1862 | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
| 1863 | ) * weight_substep(intermediate_timestep_count) |
| 1864 | ENDDO |
| 1865 | |
| 1866 | CASE ( 's' ) |
| 1867 | |
| 1868 | DO k = nzb, nzt |
| 1869 | sums_wsss_ws_l(k,tn) = sums_wsss_ws_l(k,tn) + & |
| 1870 | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
| 1871 | * ( w(k,j,i) - hom(k,1,3,0) ) & |
| 1872 | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
| 1873 | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
| 1874 | ) * weight_substep(intermediate_timestep_count) |
| 1875 | ENDDO |
| 1876 | |
| 1877 | CASE ( 'aerosol_mass', 'aerosol_number', 'salsa_gas' ) |
| 1878 | |
| 1879 | DO k = nzb, nzt |
| 1880 | sums_salsa_ws_l(k,tn) = sums_salsa_ws_l(k,tn) + & |
| 1881 | ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) ) & |
| 1882 | * ( w(k,j,i) - hom(k,1,3,0) ) & |
| 1883 | + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp ) & |
| 1884 | * ABS( w(k,j,i) - hom(k,1,3,0) ) & |
| 1885 | ) * weight_substep(intermediate_timestep_count) |
| 1886 | ENDDO |
| 1887 | |
| 1888 | ! CASE ( 'kc' ) |
| 1889 | !kk Has to be implemented for kpp chemistry |
| 1890 | |
| 1891 | END SELECT |
| 1892 | |
| 1893 | END SUBROUTINE advec_s_ws_ij |
| 1894 | |
| 1895 | |
| 1896 | |
| 1897 | |
| 1898 | !--------------------------------------------------------------------------------------------------! |
2035 | | !------------------------------------------------------------------------------! |
2036 | | SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn ) |
2037 | | |
2038 | | |
2039 | | INTEGER(iwp) :: i !< grid index along x-direction |
2040 | | INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread |
2041 | | INTEGER(iwp) :: j !< grid index along y-direction |
2042 | | INTEGER(iwp) :: k !< grid index along z-direction |
2043 | | INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults |
2044 | | INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults |
2045 | | INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults |
2046 | | INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms |
2047 | | INTEGER(iwp) :: tn !< number of OpenMP thread |
2048 | | |
2049 | | REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction |
2050 | | REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction |
2051 | | REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction |
2052 | | REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction |
2053 | | REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction |
2054 | | REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction |
2055 | | REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction |
2056 | | REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction |
2057 | | REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction |
2058 | | REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom |
2059 | | REAL(wp) :: div !< diverence on u-grid |
2060 | | REAL(wp) :: flux_d !< 6th-order flux at grid box bottom |
2061 | | REAL(wp) :: gu !< Galilei-transformation velocity along x |
2062 | | REAL(wp) :: gv !< Galilei-transformation velocity along y |
2063 | | REAL(wp) :: u_comp_l !< advection velocity along x at leftmost grid point on subdomain |
2064 | | |
2065 | | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box |
2066 | | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box |
2067 | | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top of the grid box |
2068 | | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side of the grid box |
2069 | | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box |
2070 | | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top of the grid box |
2071 | | REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !< advection velocity along x |
2072 | | REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< advection velocity along y |
2073 | | REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z |
2074 | | ! |
2075 | | !-- Used local modified copy of nzb_max (used to degrade order of |
2076 | | !-- discretization) at non-cyclic boundaries. Modify only at relevant points |
2077 | | !-- instead of the entire subdomain. This should lead to better |
2078 | | !-- load balance between boundary and non-boundary PEs. |
2079 | | IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & |
2080 | | ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR. & |
2081 | | ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR. & |
2082 | | ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j >= nyn - 2 ) THEN |
2083 | | nzb_max_l = nzt |
2084 | | ELSE |
2085 | | nzb_max_l = nzb_max |
2086 | | END IF |
2087 | | |
2088 | | gu = 2.0_wp * u_gtrans |
2089 | | gv = 2.0_wp * v_gtrans |
2090 | | ! |
2091 | | !-- Compute southside fluxes for the respective boundary of PE |
2092 | | IF ( j == nys ) THEN |
2093 | | |
| 1902 | !--------------------------------------------------------------------------------------------------! |
| 1903 | SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn ) |
| 1904 | |
| 1905 | |
| 1906 | INTEGER(iwp) :: i !< grid index along x-direction |
| 1907 | INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread |
| 1908 | INTEGER(iwp) :: j !< grid index along y-direction |
| 1909 | INTEGER(iwp) :: k !< grid index along z-direction |
| 1910 | INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults |
| 1911 | INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults |
| 1912 | INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults |
| 1913 | INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms |
| 1914 | INTEGER(iwp) :: tn !< number of OpenMP thread |
| 1915 | |
| 1916 | REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom |
| 1917 | REAL(wp) :: div !< diverence on u-grid |
| 1918 | REAL(wp) :: flux_d !< 6th-order flux at grid box bottom |
| 1919 | REAL(wp) :: gu !< Galilei-transformation velocity along x |
| 1920 | REAL(wp) :: gv !< Galilei-transformation velocity along y |
| 1921 | REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction |
| 1922 | REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction |
| 1923 | REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction |
| 1924 | REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction |
| 1925 | REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction |
| 1926 | REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction |
| 1927 | REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction |
| 1928 | REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction |
| 1929 | REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction |
| 1930 | REAL(wp) :: u_comp_l !< advection velocity along x at leftmost grid point on subdomain |
| 1931 | |
| 1932 | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box |
| 1933 | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box |
| 1934 | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top of the grid box |
| 1935 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side of the grid box |
| 1936 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box |
| 1937 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top of the grid box |
| 1938 | REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !< advection velocity along x |
| 1939 | REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< advection velocity along y |
| 1940 | REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z |
| 1941 | ! |
| 1942 | !-- Used local modified copy of nzb_max (used to degrade order of discretization) at non-cyclic |
| 1943 | !-- boundaries. Modify only at relevant points instead of the entire subdomain. This should lead to |
| 1944 | !-- better load balance between boundary and non-boundary PEs. |
| 1945 | IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & |
| 1946 | ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR. & |
| 1947 | ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR. & |
| 1948 | ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j >= nyn - 2 ) THEN |
| 1949 | nzb_max_l = nzt |
| 1950 | ELSE |
| 1951 | nzb_max_l = nzb_max |
| 1952 | END IF |
| 1953 | |
| 1954 | gu = 2.0_wp * u_gtrans |
| 1955 | gv = 2.0_wp * v_gtrans |
| 1956 | ! |
| 1957 | !-- Compute southside fluxes for the respective boundary of PE |
| 1958 | IF ( j == nys ) THEN |
| 1959 | |
| 1960 | DO k = nzb+1, nzb_max_l |
| 1961 | |
| 1962 | ibit5 = REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp ) |
| 1963 | ibit4 = REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp ) |
| 1964 | ibit3 = REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp ) |
| 1965 | |
| 1966 | v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv |
| 1967 | flux_s_u(k,tn) = v_comp(k) * ( & |
| 1968 | ( 37.0_wp * ibit5 * adv_mom_5 & |
| 1969 | + 7.0_wp * ibit4 * adv_mom_3 & |
| 1970 | + ibit3 * adv_mom_1 ) * ( u(k,j,i) + u(k,j-1,i) ) & |
| 1971 | - ( 8.0_wp * ibit5 * adv_mom_5 & |
| 1972 | + ibit4 * adv_mom_3 ) * ( u(k,j+1,i) + u(k,j-2,i) ) & |
| 1973 | + ( ibit5 * adv_mom_5 ) * ( u(k,j+2,i) + u(k,j-3,i) ) & |
| 1974 | ) |
| 1975 | |
| 1976 | diss_s_u(k,tn) = - ABS ( v_comp(k) ) * ( & |
| 1977 | ( 10.0_wp * ibit5 * adv_mom_5 & |
| 1978 | + 3.0_wp * ibit4 * adv_mom_3 & |
| 1979 | + ibit3 * adv_mom_1 ) * ( u(k,j,i) - u(k,j-1,i) ) & |
| 1980 | - ( 5.0_wp * ibit5 * adv_mom_5 & |
| 1981 | + ibit4 * adv_mom_3 ) * ( u(k,j+1,i) - u(k,j-2,i) ) & |
| 1982 | + ( ibit5 * adv_mom_5 ) * ( u(k,j+2,i) - u(k,j-3,i) ) & |
| 1983 | ) |
| 1984 | |
| 1985 | ENDDO |
| 1986 | |
| 1987 | DO k = nzb_max_l+1, nzt |
| 1988 | |
| 1989 | v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv |
| 1990 | flux_s_u(k,tn) = v_comp(k) * ( & |
| 1991 | 37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) & |
| 1992 | - 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) & |
| 1993 | + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 |
| 1994 | diss_s_u(k,tn) = - ABS(v_comp(k)) * ( & |
| 1995 | 10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) & |
| 1996 | - 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) & |
| 1997 | + ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5 |
| 1998 | |
| 1999 | ENDDO |
| 2000 | |
| 2001 | ENDIF |
| 2002 | ! |
| 2003 | !-- Compute leftside fluxes for the respective boundary of PE |
| 2004 | IF ( i == i_omp .OR. i == nxlu ) THEN |
| 2005 | |
| 2006 | DO k = nzb+1, nzb_max_l |
| 2007 | |
| 2008 | ibit2 = REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp ) |
| 2009 | ibit1 = REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp ) |
| 2010 | ibit0 = REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp ) |
| 2011 | |
| 2012 | u_comp_l = u(k,j,i) + u(k,j,i-1) - gu |
| 2013 | flux_l_u(k,j,tn) = u_comp_l * ( & |
| 2014 | ( 37.0_wp * ibit2 * adv_mom_5 & |
| 2015 | + 7.0_wp * ibit1 * adv_mom_3 & |
| 2016 | + ibit0 * adv_mom_1 ) * ( u(k,j,i) + u(k,j,i-1) ) & |
| 2017 | - ( 8.0_wp * ibit2 * adv_mom_5 & |
| 2018 | + ibit1 * adv_mom_3 ) * ( u(k,j,i+1) + u(k,j,i-2) ) & |
| 2019 | + ( ibit2 * adv_mom_5 ) * ( u(k,j,i+2) + u(k,j,i-3) ) & |
| 2020 | ) |
| 2021 | |
| 2022 | diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( & |
| 2023 | ( 10.0_wp * ibit2 * adv_mom_5 & |
| 2024 | + 3.0_wp * ibit1 * adv_mom_3 & |
| 2025 | + ibit0 * adv_mom_1 ) * ( u(k,j,i) - u(k,j,i-1) ) & |
| 2026 | - ( 5.0_wp * ibit2 * adv_mom_5 & |
| 2027 | + ibit1 * adv_mom_3 ) * ( u(k,j,i+1) - u(k,j,i-2) ) & |
| 2028 | + ( ibit2 * adv_mom_5 ) * ( u(k,j,i+2) - u(k,j,i-3) ) & |
| 2029 | ) |
| 2030 | |
| 2031 | ENDDO |
| 2032 | |
| 2033 | DO k = nzb_max_l+1, nzt |
| 2034 | |
| 2035 | u_comp_l = u(k,j,i) + u(k,j,i-1) - gu |
| 2036 | flux_l_u(k,j,tn) = u_comp_l * ( & |
| 2037 | 37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) & |
| 2038 | - 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) & |
| 2039 | + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 |
| 2040 | diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( & |
| 2041 | 10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) & |
| 2042 | - 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) & |
| 2043 | + ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5 |
| 2044 | |
| 2045 | ENDDO |
| 2046 | |
| 2047 | ENDIF |
| 2048 | ! |
| 2049 | !-- Now compute the fluxes tendency terms for the horizontal and vertical parts. |
| 2050 | DO k = nzb+1, nzb_max_l |
| 2051 | |
| 2052 | ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp ) |
| 2053 | ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp ) |
| 2054 | ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp ) |
| 2055 | |
| 2056 | u_comp(k) = u(k,j,i+1) + u(k,j,i) |
| 2057 | flux_r(k) = ( u_comp(k) - gu ) * ( & |
| 2058 | ( 37.0_wp * ibit2 * adv_mom_5 & |
| 2059 | + 7.0_wp * ibit1 * adv_mom_3 & |
| 2060 | + ibit0 * adv_mom_1 ) * ( u(k,j,i+1) + u(k,j,i) ) & |
| 2061 | - ( 8.0_wp * ibit2 * adv_mom_5 & |
| 2062 | + ibit1 * adv_mom_3 ) * ( u(k,j,i+2) + u(k,j,i-1) ) & |
| 2063 | + ( ibit2 * adv_mom_5 ) * ( u(k,j,i+3) + u(k,j,i-2) ) & |
| 2064 | ) |
| 2065 | |
| 2066 | diss_r(k) = - ABS( u_comp(k) - gu ) * ( & |
| 2067 | ( 10.0_wp * ibit2 * adv_mom_5 & |
| 2068 | + 3.0_wp * ibit1 * adv_mom_3 & |
| 2069 | + ibit0 * adv_mom_1 ) * ( u(k,j,i+1) - u(k,j,i) ) & |
| 2070 | - ( 5.0_wp * ibit2 * adv_mom_5 & |
| 2071 | + ibit1 * adv_mom_3 ) * ( u(k,j,i+2) - u(k,j,i-1) ) & |
| 2072 | + ( ibit2 * adv_mom_5 ) * ( u(k,j,i+3) - u(k,j,i-2) ) & |
| 2073 | ) |
| 2074 | |
| 2075 | ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp ) |
| 2076 | ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp ) |
| 2077 | ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp ) |
| 2078 | |
| 2079 | v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv |
| 2080 | flux_n(k) = v_comp(k) * ( & |
| 2081 | ( 37.0_wp * ibit5 * adv_mom_5 & |
| 2082 | + 7.0_wp * ibit4 * adv_mom_3 & |
| 2083 | + ibit3 * adv_mom_1 ) * ( u(k,j+1,i) + u(k,j,i) ) & |
| 2084 | - ( 8.0_wp * ibit5 * adv_mom_5 & |
| 2085 | + ibit4 * adv_mom_3 ) * ( u(k,j+2,i) + u(k,j-1,i) ) & |
| 2086 | + ( ibit5 * adv_mom_5 ) * ( u(k,j+3,i) + u(k,j-2,i) ) & |
| 2087 | ) |
| 2088 | |
| 2089 | diss_n(k) = - ABS ( v_comp(k) ) * ( & |
| 2090 | ( 10.0_wp * ibit5 * adv_mom_5 & |
| 2091 | + 3.0_wp * ibit4 * adv_mom_3 & |
| 2092 | + ibit3 * adv_mom_1 ) * ( u(k,j+1,i) - u(k,j,i) ) & |
| 2093 | - ( 5.0_wp * ibit5 * adv_mom_5 & |
| 2094 | + ibit4 * adv_mom_3 ) * ( u(k,j+2,i) - u(k,j-1,i) ) & |
| 2095 | + ( ibit5 * adv_mom_5 ) * ( u(k,j+3,i) - u(k,j-2,i) ) & |
| 2096 | ) |
| 2097 | ENDDO |
| 2098 | |
| 2099 | DO k = nzb_max_l+1, nzt |
| 2100 | |
| 2101 | u_comp(k) = u(k,j,i+1) + u(k,j,i) |
| 2102 | flux_r(k) = ( u_comp(k) - gu ) * ( & |
| 2103 | 37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) & |
| 2104 | - 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) & |
| 2105 | + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 |
| 2106 | diss_r(k) = - ABS( u_comp(k) - gu ) * ( & |
| 2107 | 10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) & |
| 2108 | - 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) & |
| 2109 | + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 |
| 2110 | |
| 2111 | v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv |
| 2112 | flux_n(k) = v_comp(k) * ( & |
| 2113 | 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) & |
| 2114 | - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) & |
| 2115 | + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 |
| 2116 | diss_n(k) = - ABS( v_comp(k) ) * ( & |
| 2117 | 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) & |
| 2118 | - 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) & |
| 2119 | + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 |
| 2120 | |
| 2121 | ENDDO |
| 2122 | ! |
| 2123 | !-- Now, compute vertical fluxes. Split loop into a part treating the lowest grid points with |
| 2124 | !-- indirect indexing, a main loop without indirect indexing, and a loop for the uppermost grid |
| 2125 | !-- points with indirect indexing. This allows better vectorization for the main loop. |
| 2126 | !-- First, compute the flux at model surface, which need has to be calculated explicitly for the |
| 2127 | !-- tendency at the first w-level. For topography wall this is done implicitely by advc_flags_m. |
| 2128 | flux_t(nzb) = 0.0_wp |
| 2129 | diss_t(nzb) = 0.0_wp |
| 2130 | w_comp(nzb) = 0.0_wp |
| 2131 | |
| 2132 | DO k = nzb+1, nzb+1 |
| 2133 | ! |
| 2134 | !-- k index has to be modified near bottom and top, else array subscripts will be exceeded. |
| 2135 | ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp ) |
| 2136 | ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp ) |
| 2137 | ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp ) |
| 2138 | |
| 2139 | k_ppp = k + 3 * ibit8 |
| 2140 | k_pp = k + 2 * ( 1 - ibit6 ) |
| 2141 | k_mm = k - 2 * ibit8 |
| 2142 | |
| 2143 | w_comp(k) = w(k,j,i) + w(k,j,i-1) |
| 2144 | flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & |
| 2145 | ( 37.0_wp * ibit8 * adv_mom_5 & |
| 2146 | + 7.0_wp * ibit7 * adv_mom_3 & |
| 2147 | + ibit6 * adv_mom_1 ) * ( u(k+1,j,i) + u(k,j,i) ) & |
| 2148 | - ( 8.0_wp * ibit8 * adv_mom_5 & |
| 2149 | + ibit7 * adv_mom_3 ) * ( u(k_pp,j,i) + u(k-1,j,i) ) & |
| 2150 | + ( ibit8 * adv_mom_5 ) * ( u(k_ppp,j,i) + u(k_mm,j,i) ) & |
| 2151 | ) |
| 2152 | |
| 2153 | diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & |
| 2154 | ( 10.0_wp * ibit8 * adv_mom_5 & |
| 2155 | + 3.0_wp * ibit7 * adv_mom_3 & |
| 2156 | + ibit6 * adv_mom_1 ) * ( u(k+1,j,i) - u(k,j,i) ) & |
| 2157 | - ( 5.0_wp * ibit8 * adv_mom_5 & |
| 2158 | + ibit7 * adv_mom_3 ) * ( u(k_pp,j,i) - u(k-1,j,i) ) & |
| 2159 | + ( ibit8 * adv_mom_5 ) * ( u(k_ppp,j,i) - u(k_mm,j,i) ) & |
| 2160 | ) |
| 2161 | ENDDO |
| 2162 | |
| 2163 | DO k = nzb+2, nzt-2 |
| 2164 | |
| 2165 | ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp ) |
| 2166 | ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp ) |
| 2167 | ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp ) |
| 2168 | |
| 2169 | w_comp(k) = w(k,j,i) + w(k,j,i-1) |
| 2170 | flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & |
| 2171 | ( 37.0_wp * ibit8 * adv_mom_5 & |
| 2172 | + 7.0_wp * ibit7 * adv_mom_3 & |
| 2173 | + ibit6 * adv_mom_1 ) * ( u(k+1,j,i) + u(k,j,i) ) & |
| 2174 | - ( 8.0_wp * ibit8 * adv_mom_5 & |
| 2175 | + ibit7 * adv_mom_3 ) * ( u(k+2,j,i) + u(k-1,j,i) ) & |
| 2176 | + ( ibit8 * adv_mom_5 ) * ( u(k+3,j,i) + u(k-2,j,i) ) & |
| 2177 | ) |
| 2178 | |
| 2179 | diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & |
| 2180 | ( 10.0_wp * ibit8 * adv_mom_5 & |
| 2181 | + 3.0_wp * ibit7 * adv_mom_3 & |
| 2182 | + ibit6 * adv_mom_1 ) * ( u(k+1,j,i) - u(k,j,i) ) & |
| 2183 | - ( 5.0_wp * ibit8 * adv_mom_5 & |
| 2184 | + ibit7 * adv_mom_3 ) * ( u(k+2,j,i) - u(k-1,j,i) ) & |
| 2185 | + ( ibit8 * adv_mom_5 ) * ( u(k+3,j,i) - u(k-2,j,i) ) & |
| 2186 | ) |
| 2187 | ENDDO |
| 2188 | |
| 2189 | DO k = nzt-1, nzt-symmetry_flag |
| 2190 | ! |
| 2191 | !-- k index has to be modified near bottom and top, else array subscripts will be exceeded. |
| 2192 | ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp ) |
| 2193 | ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp ) |
| 2194 | ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp ) |
| 2195 | |
| 2196 | k_ppp = k + 3 * ibit8 |
| 2197 | k_pp = k + 2 * ( 1 - ibit6 ) |
| 2198 | k_mm = k - 2 * ibit8 |
| 2199 | |
| 2200 | w_comp(k) = w(k,j,i) + w(k,j,i-1) |
| 2201 | flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & |
| 2202 | ( 37.0_wp * ibit8 * adv_mom_5 & |
| 2203 | + 7.0_wp * ibit7 * adv_mom_3 & |
| 2204 | + ibit6 * adv_mom_1 ) * ( u(k+1,j,i) + u(k,j,i) ) & |
| 2205 | - ( 8.0_wp * ibit8 * adv_mom_5 & |
| 2206 | + ibit7 * adv_mom_3 ) * ( u(k_pp,j,i) + u(k-1,j,i) ) & |
| 2207 | + ( ibit8 * adv_mom_5 ) * ( u(k_ppp,j,i) + u(k_mm,j,i) ) & |
| 2208 | ) |
| 2209 | |
| 2210 | diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & |
| 2211 | ( 10.0_wp * ibit8 * adv_mom_5 & |
| 2212 | + 3.0_wp * ibit7 * adv_mom_3 & |
| 2213 | + ibit6 * adv_mom_1 ) * ( u(k+1,j,i) - u(k,j,i) ) & |
| 2214 | - ( 5.0_wp * ibit8 * adv_mom_5 & |
| 2215 | + ibit7 * adv_mom_3 ) * ( u(k_pp,j,i) - u(k-1,j,i) ) & |
| 2216 | + ( ibit8 * adv_mom_5 ) * ( u(k_ppp,j,i) - u(k_mm,j,i) ) & |
| 2217 | ) |
| 2218 | ENDDO |
| 2219 | |
| 2220 | ! |
| 2221 | !-- Set resolved/turbulent flux at model top to zero (w-level). In case that a symmetric behavior |
| 2222 | !-- between bottom and top shall be guaranteed (closed channel flow), the flux at nzt is also set to |
| 2223 | !-- zero. |
| 2224 | IF ( symmetry_flag == 1 ) THEN |
| 2225 | flux_t(nzt) = 0.0_wp |
| 2226 | diss_t(nzt) = 0.0_wp |
| 2227 | w_comp(nzt) = 0.0_wp |
| 2228 | ENDIF |
| 2229 | flux_t(nzt+1) = 0.0_wp |
| 2230 | diss_t(nzt+1) = 0.0_wp |
| 2231 | w_comp(nzt+1) = 0.0_wp |
| 2232 | |
| 2233 | DO k = nzb+1, nzb_max_l |
| 2234 | |
| 2235 | flux_d = flux_t(k-1) |
| 2236 | diss_d = diss_t(k-1) |
| 2237 | |
| 2238 | ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp ) |
| 2239 | ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp ) |
| 2240 | ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp ) |
| 2241 | |
| 2242 | ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp ) |
| 2243 | ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp ) |
| 2244 | ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp ) |
| 2245 | |
| 2246 | ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp ) |
| 2247 | ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp ) |
| 2248 | ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp ) |
| 2249 | ! |
| 2250 | !-- Calculate the divergence of the velocity field. A respective correction is needed to overcome |
| 2251 | !-- numerical instabilities introduced by an insufficient reduction of divergences near |
| 2252 | !-- topography. |
| 2253 | div = ( ( u_comp(k) * ( ibit0 + ibit1 + ibit2 ) & |
| 2254 | - ( u(k,j,i) + u(k,j,i-1) ) & |
| 2255 | * ( & |
| 2256 | REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp ) & |
| 2257 | + REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp ) & |
| 2258 | + REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp ) & |
| 2259 | ) & |
| 2260 | ) * ddx & |
| 2261 | + ( ( v_comp(k) + gv ) * ( ibit3 + ibit4 + ibit5 ) & |
| 2262 | - ( v(k,j,i) + v(k,j,i-1 ) ) & |
| 2263 | * ( & |
| 2264 | REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp ) & |
| 2265 | + REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp ) & |
| 2266 | + REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp ) & |
| 2267 | ) & |
| 2268 | ) * ddy & |
| 2269 | + ( w_comp(k) * rho_air_zw(k) * ( ibit6 + ibit7 + ibit8 ) & |
| 2270 | - w_comp(k-1) * rho_air_zw(k-1) & |
| 2271 | * ( & |
| 2272 | REAL( IBITS(advc_flags_m(k-1,j,i),6,1), KIND = wp ) & |
| 2273 | + REAL( IBITS(advc_flags_m(k-1,j,i),7,1), KIND = wp ) & |
| 2274 | + REAL( IBITS(advc_flags_m(k-1,j,i),8,1), KIND = wp ) & |
| 2275 | ) & |
| 2276 | ) * drho_air(k) * ddzw(k) & |
| 2277 | ) * 0.5_wp |
| 2278 | |
| 2279 | tend(k,j,i) = tend(k,j,i) - ( & |
| 2280 | ( flux_r(k) + diss_r(k) & |
| 2281 | - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & |
| 2282 | + ( flux_n(k) + diss_n(k) & |
| 2283 | - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & |
| 2284 | + ( ( flux_t(k) + diss_t(k) ) & |
| 2285 | - ( flux_d + diss_d ) & |
| 2286 | ) * drho_air(k) * ddzw(k) & |
| 2287 | ) + div * u(k,j,i) |
| 2288 | |
| 2289 | flux_l_u(k,j,tn) = flux_r(k) |
| 2290 | diss_l_u(k,j,tn) = diss_r(k) |
| 2291 | flux_s_u(k,tn) = flux_n(k) |
| 2292 | diss_s_u(k,tn) = diss_n(k) |
| 2293 | ! |
| 2294 | !-- Statistical Evaluation of u'u'. The factor has to be applied for right evaluation when |
| 2295 | !-- gallilei_trans = .T. . |
| 2296 | sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) + & |
| 2297 | ( flux_r(k) & |
| 2298 | * ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & |
| 2299 | / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) & |
| 2300 | + diss_r(k) & |
| 2301 | * ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & |
| 2302 | / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) & |
| 2303 | ) * weight_substep(intermediate_timestep_count) |
| 2304 | ! |
| 2305 | !-- Statistical Evaluation of w'u'. |
| 2306 | sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) + & |
| 2307 | ( flux_t(k) & |
| 2308 | * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & |
| 2309 | / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & |
| 2310 | + diss_t(k) & |
| 2311 | * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & |
| 2312 | / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & |
| 2313 | ) * weight_substep(intermediate_timestep_count) |
| 2314 | ENDDO |
| 2315 | |
| 2316 | DO k = nzb_max_l+1, nzt |
| 2317 | |
| 2318 | flux_d = flux_t(k-1) |
| 2319 | diss_d = diss_t(k-1) |
| 2320 | ! |
| 2321 | !-- Calculate the divergence of the velocity field. A respective correction is needed to overcome |
| 2322 | !-- numerical instabilities introduced by an insufficient reduction of divergences near |
| 2323 | !-- topography. |
| 2324 | div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx & |
| 2325 | + ( v_comp(k) + gv - ( v(k,j,i) + v(k,j,i-1) ) ) * ddy & |
| 2326 | + ( w_comp(k) * rho_air_zw(k) & |
| 2327 | - w_comp(k-1) * rho_air_zw(k-1) & |
| 2328 | ) * drho_air(k) * ddzw(k) & |
| 2329 | ) * 0.5_wp |
| 2330 | |
| 2331 | tend(k,j,i) = tend(k,j,i) - ( & |
| 2332 | ( flux_r(k) + diss_r(k) & |
| 2333 | - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & |
| 2334 | + ( flux_n(k) + diss_n(k) & |
| 2335 | - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & |
| 2336 | + ( ( flux_t(k) + diss_t(k) ) & |
| 2337 | - ( flux_d + diss_d ) & |
| 2338 | ) * drho_air(k) * ddzw(k) & |
| 2339 | ) + div * u(k,j,i) |
| 2340 | |
| 2341 | flux_l_u(k,j,tn) = flux_r(k) |
| 2342 | diss_l_u(k,j,tn) = diss_r(k) |
| 2343 | flux_s_u(k,tn) = flux_n(k) |
| 2344 | diss_s_u(k,tn) = diss_n(k) |
| 2345 | ! |
| 2346 | !-- Statistical Evaluation of u'u'. The factor has to be applied for right evaluation when |
| 2347 | !-- gallilei_trans = .T. . |
| 2348 | sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) + & |
| 2349 | ( flux_r(k) & |
| 2350 | * ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & |
| 2351 | / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) & |
| 2352 | + diss_r(k) & |
| 2353 | * ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) & |
| 2354 | / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) & |
| 2355 | ) * weight_substep(intermediate_timestep_count) |
| 2356 | ! |
| 2357 | !-- Statistical Evaluation of w'u'. |
| 2358 | sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) + & |
| 2359 | ( flux_t(k) & |
| 2360 | * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & |
| 2361 | / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & |
| 2362 | + diss_t(k) & |
| 2363 | * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & |
| 2364 | / ( ABS( w_comp(k) ) + 1.0E-20_wp ) & |
| 2365 | ) * weight_substep(intermediate_timestep_count) |
| 2366 | ENDDO |
| 2367 | |
| 2368 | |
| 2369 | |
| 2370 | END SUBROUTINE advec_u_ws_ij |
| 2371 | |
| 2372 | |
| 2373 | |
| 2374 | !--------------------------------------------------------------------------------------------------! |
| 2375 | ! Description: |
| 2376 | ! ------------ |
| 2377 | !> Advection of v-component - Call for grid point i,j |
| 2378 | !--------------------------------------------------------------------------------------------------! |
| 2379 | SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn ) |
| 2380 | |
| 2381 | |
| 2382 | INTEGER(iwp) :: i !< grid index along x-direction |
| 2383 | INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread |
| 2384 | INTEGER(iwp) :: j !< grid index along y-direction |
| 2385 | INTEGER(iwp) :: k !< grid index along z-direction |
| 2386 | INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults |
| 2387 | INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults |
| 2388 | INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults |
| 2389 | INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms |
| 2390 | INTEGER(iwp) :: tn !< number of OpenMP thread |
| 2391 | |
| 2392 | REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom |
| 2393 | REAL(wp) :: div !< divergence on v-grid |
| 2394 | REAL(wp) :: flux_d !< 6th-order flux at grid box bottom |
| 2395 | REAL(wp) :: gu !< Galilei-transformation velocity along x |
| 2396 | REAL(wp) :: gv !< Galilei-transformation velocity along y |
| 2397 | REAL(wp) :: ibit9 !< flag indicating 1st-order scheme along x-direction |
| 2398 | REAL(wp) :: ibit10 !< flag indicating 3rd-order scheme along x-direction |
| 2399 | REAL(wp) :: ibit11 !< flag indicating 5th-order scheme along x-direction |
| 2400 | REAL(wp) :: ibit12 !< flag indicating 1st-order scheme along y-direction |
| 2401 | REAL(wp) :: ibit13 !< flag indicating 3rd-order scheme along y-direction |
| 2402 | REAL(wp) :: ibit14 !< flag indicating 3rd-order scheme along y-direction |
| 2403 | REAL(wp) :: ibit15 !< flag indicating 1st-order scheme along z-direction |
| 2404 | REAL(wp) :: ibit16 !< flag indicating 3rd-order scheme along z-direction |
| 2405 | REAL(wp) :: ibit17 !< flag indicating 3rd-order scheme along z-direction |
| 2406 | REAL(wp) :: v_comp_l !< advection velocity along y on leftmost grid point on subdomain |
| 2407 | |
| 2408 | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box |
| 2409 | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box |
| 2410 | REAL(wp), DIMENSION(nzb:nzt+1) :: diss_t !< discretized artificial dissipation at top of the grid box |
| 2411 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_n !< discretized 6th-order flux at northward-side of the grid box |
| 2412 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_r !< discretized 6th-order flux at rightward-side of the grid box |
| 2413 | REAL(wp), DIMENSION(nzb:nzt+1) :: flux_t !< discretized 6th-order flux at top of the grid box |
| 2414 | REAL(wp), DIMENSION(nzb:nzt+1) :: u_comp !< advection velocity along x |
| 2415 | REAL(wp), DIMENSION(nzb:nzt+1) :: v_comp !< advection velocity along y |
| 2416 | REAL(wp), DIMENSION(nzb:nzt+1) :: w_comp !< advection velocity along z |
| 2417 | ! |
| 2418 | !-- Used local modified copy of nzb_max (used to degrade order of discretization) at non-cyclic |
| 2419 | !-- boundaries. Modify only at relevant points instead of the entire subdomain. This should lead to |
| 2420 | !-- better load balance between boundary and non-boundary PEs. |
| 2421 | IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR. & |
| 2422 | ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR. & |
| 2423 | ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR. & |
| 2424 | ( bc_dirichlet_n .OR. bc_radiation_n ) .AND. j >= nyn - 2 ) THEN |
| 2425 | nzb_max_l = nzt |
| 2426 | ELSE |
| 2427 | nzb_max_l = nzb_max |
| 2428 | END IF |
| 2429 | |
| 2430 | gu = 2.0_wp * u_gtrans |
| 2431 | gv = 2.0_wp * v_gtrans |
| 2432 | |
| 2433 | ! |
| 2434 | !-- Compute leftside fluxes for the respective boundary. |
| 2435 | IF ( i == i_omp ) THEN |
| 2436 | |
| 2437 | DO k = nzb+1, nzb_max_l |
| 2438 | |
| 2439 | ibit11 = REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp ) |
| 2440 | ibit10 = REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp ) |
| 2441 | ibit9 = REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp ) |
| 2442 | |
| 2443 | u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu |
| 2444 | flux_l_v(k,j,tn) = u_comp(k) * ( & |
| 2445 | ( 37.0_wp * ibit11 * adv_mom_5 & |
| 2446 | + 7.0_wp * ibit10 * adv_mom_3 & |
| 2447 | + ibit9 * adv_mom_1 ) * ( v(k,j,i) + v(k,j,i-1) ) & |
| 2448 | - ( 8.0_wp * ibit11 * adv_mom_5 & |
| 2449 | + ibit10 * adv_mom_3 ) * ( v(k,j,i+1) + v(k,j,i-2) ) & |
| 2450 | + ( ibit11 * adv_mom_5 ) * ( v(k,j,i+2) + v(k,j,i-3) ) & |
| 2451 | ) |
| 2452 | |
| 2453 | diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( & |
| 2454 | ( 10.0_wp * ibit11 * adv_mom_5 & |
| 2455 | + 3.0_wp * ibit10 * adv_mom_3 & |
| 2456 | + ibit9 * adv_mom_1 ) * ( v(k,j,i) - v(k,j,i-1) ) & |
| 2457 | - ( 5.0_wp * ibit11 * adv_mom_5 & |
| 2458 | + ibit10 * adv_mom_3 ) * ( v(k,j,i+1) - v(k,j,i-2) ) & |
| 2459 | + ( ibit11 * adv_mom_5 ) * ( v(k,j,i+2) - v(k,j,i-3) ) & |
| 2460 | ) |
| 2461 | |
| 2462 | ENDDO |
| 2463 | |
| 2464 | DO k = nzb_max_l+1, nzt |
| 2465 | |
| 2466 | u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu |
| 2467 | flux_l_v(k,j,tn) = u_comp(k) * ( & |
| 2468 | 37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) & |
| 2469 | - 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) & |
| 2470 | + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 |
| 2471 | diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( & |
| 2472 | 10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) & |
| 2473 | - 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) & |
| 2474 | + ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5 |
| 2475 | |
| 2476 | ENDDO |
| 2477 | |
| 2478 | ENDIF |
| 2479 | ! |
| 2480 | !-- Compute southside fluxes for the respective boundary. |
| 2481 | IF ( j == nysv ) THEN |
| 2482 | |
| 2483 | DO k = nzb+1, nzb_max_l |
| 2484 | |
| 2485 | ibit14 = REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp ) |
| 2486 | ibit13 = REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp ) |
| 2487 | ibit12 = REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp ) |
| 2488 | |
| 2489 | v_comp_l = v(k,j,i) + v(k,j-1,i) - gv |
| 2490 | flux_s_v(k,tn) = v_comp_l * ( & |
| 2491 | ( 37.0_wp * ibit14 * adv_mom_5 & |
| 2492 | + 7.0_wp * ibit13 * adv_mom_3 & |
| 2493 | + ibit12 * adv_mom_1 ) * ( v(k,j,i) + v(k,j-1,i) ) & |
| 2494 | - ( 8.0_wp * ibit14 * adv_mom_5 & |
| 2495 | + ibit13 * adv_mom_3 ) * ( v(k,j+1,i) + v(k,j-2,i) ) & |
| 2496 | + ( ibit14 * adv_mom_5 ) * ( v(k,j+2,i) + v(k,j-3,i) ) & |
| 2497 | ) |
| 2498 | |
| 2499 | diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & |
| 2500 | ( 10.0_wp * ibit14 * adv_mom_5 & |
| 2501 | + 3.0_wp * ibit13 * adv_mom_3 & |
| 2502 | + ibit12 * adv_mom_1 ) * ( v(k,j,i) - v(k,j-1,i) ) & |
| 2503 | - ( 5.0_wp * ibit14 * adv_mom_5 & |
| 2504 | + ibit13 * adv_mom_3 ) * ( v(k,j+1,i) - v(k,j-2,i) ) & |
| 2505 | + ( ibit14 * adv_mom_5 ) * ( v(k,j+2,i) - v(k,j-3,i) ) & |
| 2506 | ) |
| 2507 | |
| 2508 | ENDDO |
| 2509 | |
| 2510 | DO k = nzb_max_l+1, nzt |
| 2511 | |
| 2512 | v_comp_l = v(k,j,i) + v(k,j-1,i) - gv |
| 2513 | flux_s_v(k,tn) = v_comp_l * ( & |
| 2514 | 37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) & |
| 2515 | - 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) & |
| 2516 | + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 |
| 2517 | diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & |
| 2518 | 10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) & |
| 2519 | - 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) & |
| 2520 | + ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5 |
| 2521 | |
| 2522 | ENDDO |
| 2523 | |
| 2524 | ENDIF |
| 2525 | ! |
| 2526 | !-- Now compute the fluxes and tendency terms for the horizontal and |
| 2527 | !-- verical parts. |
| 2528 | DO k = nzb+1, nzb_max_l |
| 2529 | |
| 2530 | ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp ) |
| 2531 | ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp ) |
| 2532 | ibit9 = REAL( IBITS(advc_flags_m(k,j,i),9,1), KIND = wp ) |
| 2533 | |
| 2534 | u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu |
| 2535 | flux_r(k) = u_comp(k) * ( & |
| 2536 | ( 37.0_wp * ibit11 * adv_mom_5 & |
| 2537 | + 7.0_wp * ibit10 * adv_mom_3 & |
| 2538 | + ibit9 * adv_mom_1 ) * ( v(k,j,i+1) + v(k,j,i) ) & |
| 2539 | - ( 8.0_wp * ibit11 * adv_mom_5 & |
| 2540 | + ibit10 * adv_mom_3 ) * ( v(k,j,i+2) + v(k,j,i-1) ) & |
| 2541 | + ( ibit11 * adv_mom_5 ) * ( v(k,j,i+3) + v(k,j,i-2) ) & |
| 2542 | ) |
|