Changeset 4109
- Timestamp:
- Jul 22, 2019 5:00:34 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r4079 r4109 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! - Separate initialization of advection flags for momentum and scalars. In this 23 ! context, resort the bits and do some minor formatting. 24 ! - Make flag initialization for scalars more flexible, introduce an 25 ! arguemnt list to indicate non-cyclic boundaries (required for decycled 26 ! scalars such as chemical species or aerosols) 27 ! - Introduce extended 'degradation zones', where horizontal advection of 28 ! passive scalars is discretized by first-order scheme at all grid points 29 ! that in the vicinity of buildings (<= 3 grid points). Even though no 30 ! building is within the numerical stencil, first-order scheme is used. 31 ! At fourth and fifth grid point the order of the horizontal advection scheme 32 ! is successively upgraded. 33 ! These extended degradation zones are used to avoid stationary numerical 34 ! oscillations, which are responsible for high concentration maxima that may 35 ! appear under shear-free stable conditions. 36 ! - Change interface for scalar advection routine. 37 ! - Bugfix, avoid uninitialized value sk_num in vector version of scalar 38 ! advection 23 39 ! 24 40 ! Former revisions: … … 106 122 ! Change in file header (GPL part) 107 123 ! Implement advection for TKE-dissipation in case of RANS-TKE-e closure (TG) 108 ! Allocate advc_flags_ 1/2 within ws_init_flags instead of init_grid124 ! Allocate advc_flags_m/2 within ws_init_flags instead of init_grid 109 125 ! Change argument list for exchange_horiz_2d_int (MS) 110 126 ! … … 120 136 ! 121 137 ! 2232 2017-05-30 17:47:52Z suehring 122 ! Rename wall_flags_0 and wall_flags_00 into advc_flags_ 1 and advc_flags_2,138 ! Rename wall_flags_0 and wall_flags_00 into advc_flags_m and advc_flags_m, 123 139 ! respectively. 124 ! Set advc_flags_ 1/2 on basis of wall_flags_0/00 instead of nzb_s/u/v/w_inner.125 ! Setting advc_flags_ 1/2 also for downward-facing walls140 ! Set advc_flags_m/2 on basis of wall_flags_0/00 instead of nzb_s/u/v/w_inner. 141 ! Setting advc_flags_m/2 also for downward-facing walls 126 142 ! 127 143 ! 2200 2017-04-11 11:37:51Z suehring … … 246 262 ! vector version. 247 263 ! Degradation of the applied order of scheme is now steered by multiplying with 248 ! Integer advc_flags_ 1. 2nd order scheme, WS3 and WS5 are calculated on each264 ! Integer advc_flags_m. 2nd order scheme, WS3 and WS5 are calculated on each 249 265 ! grid point and mulitplied with the appropriate flag. 250 266 ! 2nd order numerical dissipation term changed. Now the appropriate 2nd order … … 267 283 ! 268 284 ! 411 2009-12-11 12:31:43 Z suehring 285 ! 286 ! 287 ! 288 ! @author Matthias Suehring 289 ! 269 290 ! 270 291 ! Description: … … 282 303 !> 283 304 !> @todo Implement monotonic flux limiter also for vector version. 305 !> @todo Move 3d arrays advc_flag, advc_flags_m from modules to advec_ws 306 !> @todo Move arrays flux_l_x from modules to advec_ws 284 307 !------------------------------------------------------------------------------! 285 308 MODULE advec_ws … … 299 322 300 323 USE control_parameters, & 301 ONLY: humidity, loop_optimization, passive_scalar, & 302 rans_tke_e, ws_scheme_mom, ws_scheme_sca, & 303 momentum_advec, scalar_advec, & 304 bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, & 305 bc_dirichlet_s, bc_radiation_l, bc_radiation_n, & 306 bc_radiation_r, bc_radiation_s, intermediate_timestep_count, & 307 u_gtrans, v_gtrans, dt_3d 324 ONLY: air_chemistry, & 325 bc_dirichlet_l, & 326 bc_dirichlet_n, & 327 bc_dirichlet_r, & 328 bc_dirichlet_s, & 329 bc_radiation_l, & 330 bc_radiation_n, & 331 bc_radiation_r, & 332 bc_radiation_s, & 333 humidity, & 334 loop_optimization, & 335 passive_scalar, & 336 rans_tke_e, & 337 momentum_advec, & 338 salsa, & 339 scalar_advec, & 340 intermediate_timestep_count, & 341 u_gtrans, & 342 v_gtrans, & 343 ws_scheme_mom, & 344 ws_scheme_sca, & 345 dt_3d 308 346 309 347 USE indices, & 310 ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, & 311 nzb, nzb_max, nzt, advc_flags_1, advc_flags_2, wall_flags_0 348 ONLY: advc_flags_m, & 349 advc_flags_s, & 350 nbgp, & 351 nx, & 352 nxl, & 353 nxlg, & 354 nxlu, & 355 nxr, & 356 nxrg, & 357 ny, & 358 nyn, & 359 nyng, & 360 nys, & 361 nysg, & 362 nysv, & 363 nzb, & 364 nzb_max, & 365 nzt, & 366 wall_flags_0 312 367 313 368 USE grid_variables, & … … 338 393 PRIVATE 339 394 PUBLIC advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, ws_init, & 340 ws_init_flags , ws_statistics395 ws_init_flags_momentum, ws_init_flags_scalar, ws_statistics 341 396 342 397 INTERFACE ws_init 343 398 MODULE PROCEDURE ws_init 344 END INTERFACE ws_init 345 346 INTERFACE ws_init_flags 347 MODULE PROCEDURE ws_init_flags 348 END INTERFACE ws_init_flags 399 END INTERFACE ws_init 400 401 INTERFACE ws_init_flags_momentum 402 MODULE PROCEDURE ws_init_flags_momentum 403 END INTERFACE ws_init_flags_momentum 404 405 INTERFACE ws_init_flags_scalar 406 MODULE PROCEDURE ws_init_flags_scalar 407 END INTERFACE ws_init_flags_scalar 349 408 350 409 INTERFACE ws_statistics … … 490 549 ! Description: 491 550 ! ------------ 492 !> Initialization of flags for WS-scheme used to degrade the order of the scheme 493 !> near walls. 551 !> Initialization of flags to control the order of the advection scheme near 552 !> solid walls and non-cyclic inflow boundaries, where the order is sucessively 553 !> degraded. 494 554 !------------------------------------------------------------------------------! 495 SUBROUTINE ws_init_flags 555 SUBROUTINE ws_init_flags_momentum 496 556 497 557 … … 505 565 LOGICAL :: flag_set !< steering variable for advection flags 506 566 507 ALLOCATE( advc_flags_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 508 advc_flags_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 509 advc_flags_1 = 0 510 advc_flags_2 = 0 511 IF ( scalar_advec == 'ws-scheme' ) THEN 512 ! 513 !-- Set flags to steer the degradation of the advection scheme in advec_ws 514 !-- near topography, inflow- and outflow boundaries as well as bottom and 515 !-- top of model domain. advc_flags_1 remains zero for all non-prognostic 516 !-- grid points. 517 DO i = nxl, nxr 518 DO j = nys, nyn 519 DO k = nzb+1, nzt 520 ! 521 !-- scalar - x-direction 522 !-- WS1 (0), WS3 (1), WS5 (2) 523 IF ( ( .NOT. BTEST(wall_flags_0(k,j,i+1),0) & 524 .OR. .NOT. BTEST(wall_flags_0(k,j,i+2),0) & 525 .OR. .NOT. BTEST(wall_flags_0(k,j,i-1),0) ) & 526 .OR. ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 527 .AND. i == nxl ) & 528 .OR. ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 529 .AND. i == nxr ) ) & 530 THEN 531 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 0 ) 532 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+3),0) & 533 .AND. BTEST(wall_flags_0(k,j,i+1),0) & 534 .AND. BTEST(wall_flags_0(k,j,i+2),0) & 535 .AND. BTEST(wall_flags_0(k,j,i-1),0) & 536 ) .OR. & 537 ( .NOT. BTEST(wall_flags_0(k,j,i-2),0) & 538 .AND. BTEST(wall_flags_0(k,j,i+1),0) & 539 .AND. BTEST(wall_flags_0(k,j,i+2),0) & 540 .AND. BTEST(wall_flags_0(k,j,i-1),0) & 541 ) & 542 .OR. & 543 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 544 .AND. i == nxr-1 ) .OR. & 545 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 546 .AND. i == nxlu ) ) & ! why not nxl+1 547 THEN 548 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 1 ) 549 ELSEIF ( BTEST(wall_flags_0(k,j,i+1),0) & 550 .AND. BTEST(wall_flags_0(k,j,i+2),0) & 551 .AND. BTEST(wall_flags_0(k,j,i+3),0) & 552 .AND. BTEST(wall_flags_0(k,j,i-1),0) & 553 .AND. BTEST(wall_flags_0(k,j,i-2),0) ) & 554 THEN 555 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 2 ) 556 ENDIF 557 ! 558 !-- scalar - y-direction 559 !-- WS1 (3), WS3 (4), WS5 (5) 560 IF ( ( .NOT. BTEST(wall_flags_0(k,j+1,i),0) & 561 .OR. .NOT. BTEST(wall_flags_0(k,j+2,i),0) & 562 .OR. .NOT. BTEST(wall_flags_0(k,j-1,i),0)) & 563 .OR. ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 564 .AND. j == nys ) & 565 .OR. ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 566 .AND. j == nyn ) ) & 567 THEN 568 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 3 ) 569 ! 570 !-- WS3 571 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+3,i),0) & 572 .AND. BTEST(wall_flags_0(k,j+1,i),0) & 573 .AND. BTEST(wall_flags_0(k,j+2,i),0) & 574 .AND. BTEST(wall_flags_0(k,j-1,i),0) & 575 ) .OR. & 576 ( .NOT. BTEST(wall_flags_0(k,j-2,i),0) & 577 .AND. BTEST(wall_flags_0(k,j+1,i),0) & 578 .AND. BTEST(wall_flags_0(k,j+2,i),0) & 579 .AND. BTEST(wall_flags_0(k,j-1,i),0) & 580 ) & 581 .OR. & 582 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 583 .AND. j == nysv ) .OR. & ! why not nys+1 584 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 585 .AND. j == nyn-1 ) ) & 586 THEN 587 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 4 ) 588 ! 589 !-- WS5 590 ELSEIF ( BTEST(wall_flags_0(k,j+1,i),0) & 591 .AND. BTEST(wall_flags_0(k,j+2,i),0) & 592 .AND. BTEST(wall_flags_0(k,j+3,i),0) & 593 .AND. BTEST(wall_flags_0(k,j-1,i),0) & 594 .AND. BTEST(wall_flags_0(k,j-2,i),0) ) & 595 THEN 596 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 5 ) 597 ENDIF 598 ! 599 !-- scalar - z-direction 600 !-- WS1 (6), WS3 (7), WS5 (8) 601 IF ( k == nzb+1 ) THEN 602 k_mm = nzb 603 ELSE 604 k_mm = k - 2 605 ENDIF 606 IF ( k > nzt-1 ) THEN 607 k_pp = nzt+1 608 ELSE 609 k_pp = k + 2 610 ENDIF 611 IF ( k > nzt-2 ) THEN 612 k_ppp = nzt+1 613 ELSE 614 k_ppp = k + 3 615 ENDIF 616 617 flag_set = .FALSE. 618 IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),0) .AND. & 619 BTEST(wall_flags_0(k,j,i),0) .OR. & 620 .NOT. BTEST(wall_flags_0(k_pp,j,i),0) .AND. & 621 BTEST(wall_flags_0(k,j,i),0) .OR. & 622 k == nzt ) & 623 THEN 624 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 6 ) 625 flag_set = .TRUE. 626 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),0) .OR. & 627 .NOT. BTEST(wall_flags_0(k_ppp,j,i),0) ) .AND. & 628 BTEST(wall_flags_0(k-1,j,i),0) .AND. & 629 BTEST(wall_flags_0(k,j,i),0) .AND. & 630 BTEST(wall_flags_0(k+1,j,i),0) .AND. & 631 BTEST(wall_flags_0(k_pp,j,i),0) .OR. & 632 k == nzt - 1 ) & 633 THEN 634 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 7 ) 635 flag_set = .TRUE. 636 ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0) & 637 .AND. BTEST(wall_flags_0(k-1,j,i),0) & 638 .AND. BTEST(wall_flags_0(k,j,i),0) & 639 .AND. BTEST(wall_flags_0(k+1,j,i),0) & 640 .AND. BTEST(wall_flags_0(k_pp,j,i),0) & 641 .AND. BTEST(wall_flags_0(k_ppp,j,i),0) & 642 .AND. .NOT. flag_set ) & 643 THEN 644 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 8 ) 645 ENDIF 646 647 ENDDO 567 advc_flags_m = 0 568 569 ! 570 !-- Set advc_flags_m to steer the degradation of the advection scheme in advec_ws 571 !-- near topography, inflow- and outflow boundaries as well as bottom and 572 !-- top of model domain. advc_flags_m remains zero for all non-prognostic 573 !-- grid points. 574 !-- u-component 575 DO i = nxl, nxr 576 DO j = nys, nyn 577 DO k = nzb+1, nzt 578 ! 579 !-- At first, set flags to WS1. 580 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 581 !-- in order to handle the left/south flux. 582 !-- near vertical walls. 583 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 ) 584 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 ) 585 ! 586 !-- u component - x-direction 587 !-- WS1 (0), WS3 (1), WS5 (2) 588 IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),1) .OR. & 589 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 590 .AND. i <= nxlu ) .OR. & 591 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 592 .AND. i == nxr ) ) & 593 THEN 594 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 0 ) 595 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),1) .AND. & 596 BTEST(wall_flags_0(k,j,i+1),1) .OR. & 597 .NOT. BTEST(wall_flags_0(k,j,i-1),1) ) & 598 .OR. & 599 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 600 .AND. i == nxr-1 ) .OR. & 601 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 602 .AND. i == nxlu+1) ) & 603 THEN 604 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 1 ) 605 ! 606 !-- Clear flag for WS1 607 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 ) 608 ELSEIF ( BTEST(wall_flags_0(k,j,i+1),1) .AND. & 609 BTEST(wall_flags_0(k,j,i+2),1) .AND. & 610 BTEST(wall_flags_0(k,j,i-1),1) ) & 611 THEN 612 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 2 ) 613 ! 614 !-- Clear flag for WS1 615 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 0 ) 616 ENDIF 617 ! 618 !-- u component - y-direction 619 !-- WS1 (3), WS3 (4), WS5 (5) 620 IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),1) .OR. & 621 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 622 .AND. j == nys ) .OR. & 623 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 624 .AND. j == nyn ) ) & 625 THEN 626 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 3 ) 627 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),1) .AND. & 628 BTEST(wall_flags_0(k,j+1,i),1) .OR. & 629 .NOT. BTEST(wall_flags_0(k,j-1,i),1) ) & 630 .OR. & 631 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 632 .AND. j == nysv ) .OR. & 633 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 634 .AND. j == nyn-1 ) ) & 635 THEN 636 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 4 ) 637 ! 638 !-- Clear flag for WS1 639 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) 640 ELSEIF ( BTEST(wall_flags_0(k,j+1,i),1) .AND. & 641 BTEST(wall_flags_0(k,j+2,i),1) .AND. & 642 BTEST(wall_flags_0(k,j-1,i),1) ) & 643 THEN 644 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 5 ) 645 ! 646 !-- Clear flag for WS1 647 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 3 ) 648 ENDIF 649 ! 650 !-- u component - z-direction 651 !-- WS1 (6), WS3 (7), WS5 (8) 652 IF ( k == nzb+1 ) THEN 653 k_mm = nzb 654 ELSE 655 k_mm = k - 2 656 ENDIF 657 IF ( k > nzt-1 ) THEN 658 k_pp = nzt+1 659 ELSE 660 k_pp = k + 2 661 ENDIF 662 IF ( k > nzt-2 ) THEN 663 k_ppp = nzt+1 664 ELSE 665 k_ppp = k + 3 666 ENDIF 667 668 flag_set = .FALSE. 669 IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),1) .AND. & 670 BTEST(wall_flags_0(k,j,i),1) .OR. & 671 .NOT. BTEST(wall_flags_0(k_pp,j,i),1) .AND. & 672 BTEST(wall_flags_0(k,j,i),1) .OR. & 673 k == nzt ) & 674 THEN 675 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 6 ) 676 flag_set = .TRUE. 677 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),1) .OR. & 678 .NOT. BTEST(wall_flags_0(k_ppp,j,i),1) ) .AND. & 679 BTEST(wall_flags_0(k-1,j,i),1) .AND. & 680 BTEST(wall_flags_0(k,j,i),1) .AND. & 681 BTEST(wall_flags_0(k+1,j,i),1) .AND. & 682 BTEST(wall_flags_0(k_pp,j,i),1) .OR. & 683 k == nzt - 1 ) & 684 THEN 685 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 ) 686 flag_set = .TRUE. 687 ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),1) .AND. & 688 BTEST(wall_flags_0(k-1,j,i),1) .AND. & 689 BTEST(wall_flags_0(k,j,i),1) .AND. & 690 BTEST(wall_flags_0(k+1,j,i),1) .AND. & 691 BTEST(wall_flags_0(k_pp,j,i),1) .AND. & 692 BTEST(wall_flags_0(k_ppp,j,i),1) .AND. & 693 .NOT. flag_set ) & 694 THEN 695 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 8 ) 696 ENDIF 697 648 698 ENDDO 649 699 ENDDO 650 ENDIF 651 652 IF ( momentum_advec == 'ws-scheme' ) THEN 653 ! 654 !-- Set advc_flags_1 to steer the degradation of the advection scheme in advec_ws 655 !-- near topography, inflow- and outflow boundaries as well as bottom and 656 !-- top of model domain. advc_flags_1 remains zero for all non-prognostic 657 !-- grid points. 658 DO i = nxl, nxr 659 DO j = nys, nyn 660 DO k = nzb+1, nzt 661 ! 662 !-- At first, set flags to WS1. 663 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 664 !-- in order to handle the left/south flux. 665 !-- near vertical walls. 666 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 ) 667 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 ) 668 ! 669 !-- u component - x-direction 670 !-- WS1 (9), WS3 (10), WS5 (11) 671 IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),1) .OR. & 672 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 673 .AND. i <= nxlu ) .OR. & 674 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 675 .AND. i == nxr ) ) & 676 THEN 677 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 9 ) 678 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),1) .AND. & 679 BTEST(wall_flags_0(k,j,i+1),1) .OR. & 680 .NOT. BTEST(wall_flags_0(k,j,i-1),1) ) & 681 .OR. & 682 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 683 .AND. i == nxr-1 ) .OR. & 684 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 685 .AND. i == nxlu+1) ) & 686 THEN 687 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 10 ) 688 ! 689 !-- Clear flag for WS1 690 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 ) 691 ELSEIF ( BTEST(wall_flags_0(k,j,i+1),1) .AND. & 692 BTEST(wall_flags_0(k,j,i+2),1) .AND. & 693 BTEST(wall_flags_0(k,j,i-1),1) ) & 694 THEN 695 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 11 ) 696 ! 697 !-- Clear flag for WS1 698 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 9 ) 699 ENDIF 700 ! 701 !-- u component - y-direction 702 !-- WS1 (12), WS3 (13), WS5 (14) 703 IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),1) .OR. & 704 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 705 .AND. j == nys ) .OR. & 706 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 707 .AND. j == nyn ) ) & 708 THEN 709 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 12 ) 710 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),1) .AND. & 711 BTEST(wall_flags_0(k,j+1,i),1) .OR. & 712 .NOT. BTEST(wall_flags_0(k,j-1,i),1) ) & 713 .OR. & 714 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 715 .AND. j == nysv ) .OR. & 716 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 717 .AND. j == nyn-1 ) ) & 718 THEN 719 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 13 ) 720 ! 721 !-- Clear flag for WS1 722 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 ) 723 ELSEIF ( BTEST(wall_flags_0(k,j+1,i),1) .AND. & 724 BTEST(wall_flags_0(k,j+2,i),1) .AND. & 725 BTEST(wall_flags_0(k,j-1,i),1) ) & 726 THEN 727 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 14 ) 728 ! 729 !-- Clear flag for WS1 730 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 12 ) 731 ENDIF 732 ! 733 !-- u component - z-direction 734 !-- WS1 (15), WS3 (16), WS5 (17) 735 IF ( k == nzb+1 ) THEN 736 k_mm = nzb 737 ELSE 738 k_mm = k - 2 739 ENDIF 740 IF ( k > nzt-1 ) THEN 741 k_pp = nzt+1 742 ELSE 743 k_pp = k + 2 744 ENDIF 745 IF ( k > nzt-2 ) THEN 746 k_ppp = nzt+1 747 ELSE 748 k_ppp = k + 3 749 ENDIF 750 751 flag_set = .FALSE. 752 IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),1) .AND. & 753 BTEST(wall_flags_0(k,j,i),1) .OR. & 754 .NOT. BTEST(wall_flags_0(k_pp,j,i),1) .AND. & 755 BTEST(wall_flags_0(k,j,i),1) .OR. & 756 k == nzt ) & 757 THEN 758 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 15 ) 759 flag_set = .TRUE. 760 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),1) .OR. & 761 .NOT. BTEST(wall_flags_0(k_ppp,j,i),1) ) .AND. & 762 BTEST(wall_flags_0(k-1,j,i),1) .AND. & 763 BTEST(wall_flags_0(k,j,i),1) .AND. & 764 BTEST(wall_flags_0(k+1,j,i),1) .AND. & 765 BTEST(wall_flags_0(k_pp,j,i),1) .OR. & 766 k == nzt - 1 ) & 767 THEN 768 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 16 ) 769 flag_set = .TRUE. 770 ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),1) .AND. & 771 BTEST(wall_flags_0(k-1,j,i),1) .AND. & 772 BTEST(wall_flags_0(k,j,i),1) .AND. & 773 BTEST(wall_flags_0(k+1,j,i),1) .AND. & 774 BTEST(wall_flags_0(k_pp,j,i),1) .AND. & 775 BTEST(wall_flags_0(k_ppp,j,i),1) .AND. & 776 .NOT. flag_set ) & 777 THEN 778 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 17 ) 779 ENDIF 780 781 ENDDO 700 ENDDO 701 ! 702 !-- v-component 703 DO i = nxl, nxr 704 DO j = nys, nyn 705 DO k = nzb+1, nzt 706 ! 707 !-- At first, set flags to WS1. 708 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 709 !-- in order to handle the left/south flux. 710 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 ) 711 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 ) 712 ! 713 !-- v component - x-direction 714 !-- WS1 (9), WS3 (10), WS5 (11) 715 IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),2) .OR. & 716 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 717 .AND. i == nxl ) .OR. & 718 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 719 .AND. i == nxr ) ) & 720 THEN 721 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 9 ) 722 ! 723 !-- WS3 724 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),2) .AND. & 725 BTEST(wall_flags_0(k,j,i+1),2) ) .OR. & 726 .NOT. BTEST(wall_flags_0(k,j,i-1),2) & 727 .OR. & 728 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 729 .AND. i == nxr-1 ) .OR. & 730 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 731 .AND. i == nxlu ) ) & 732 THEN 733 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 10 ) 734 ! 735 !-- Clear flag for WS1 736 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 ) 737 ELSEIF ( BTEST(wall_flags_0(k,j,i+1),2) .AND. & 738 BTEST(wall_flags_0(k,j,i+2),2) .AND. & 739 BTEST(wall_flags_0(k,j,i-1),2) ) & 740 THEN 741 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 11 ) 742 ! 743 !-- Clear flag for WS1 744 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 9 ) 745 ENDIF 746 ! 747 !-- v component - y-direction 748 !-- WS1 (12), WS3 (13), WS5 (14) 749 IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),2) .OR. & 750 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 751 .AND. j <= nysv ) .OR. & 752 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 753 .AND. j == nyn ) ) & 754 THEN 755 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 12 ) 756 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),2) .AND. & 757 BTEST(wall_flags_0(k,j+1,i),2) .OR. & 758 .NOT. BTEST(wall_flags_0(k,j-1,i),2) ) & 759 .OR. & 760 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 761 .AND. j == nysv+1) .OR. & 762 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 763 .AND. j == nyn-1 ) ) & 764 THEN 765 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 13 ) 766 ! 767 !-- Clear flag for WS1 768 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 ) 769 ELSEIF ( BTEST(wall_flags_0(k,j+1,i),2) .AND. & 770 BTEST(wall_flags_0(k,j+2,i),2) .AND. & 771 BTEST(wall_flags_0(k,j-1,i),2) ) & 772 THEN 773 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 14 ) 774 ! 775 !-- Clear flag for WS1 776 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 12 ) 777 ENDIF 778 ! 779 !-- v component - z-direction 780 !-- WS1 (15), WS3 (16), WS5 (17) 781 IF ( k == nzb+1 ) THEN 782 k_mm = nzb 783 ELSE 784 k_mm = k - 2 785 ENDIF 786 IF ( k > nzt-1 ) THEN 787 k_pp = nzt+1 788 ELSE 789 k_pp = k + 2 790 ENDIF 791 IF ( k > nzt-2 ) THEN 792 k_ppp = nzt+1 793 ELSE 794 k_ppp = k + 3 795 ENDIF 796 797 flag_set = .FALSE. 798 IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),2) .AND. & 799 BTEST(wall_flags_0(k,j,i),2) .OR. & 800 .NOT. BTEST(wall_flags_0(k_pp,j,i),2) .AND. & 801 BTEST(wall_flags_0(k,j,i),2) .OR. & 802 k == nzt ) & 803 THEN 804 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 15 ) 805 flag_set = .TRUE. 806 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),2) .OR. & 807 .NOT. BTEST(wall_flags_0(k_ppp,j,i),2) ) .AND. & 808 BTEST(wall_flags_0(k-1,j,i),2) .AND. & 809 BTEST(wall_flags_0(k,j,i),2) .AND. & 810 BTEST(wall_flags_0(k+1,j,i),2) .AND. & 811 BTEST(wall_flags_0(k_pp,j,i),2) .OR. & 812 k == nzt - 1 ) & 813 THEN 814 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 ) 815 flag_set = .TRUE. 816 ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),2) .AND. & 817 BTEST(wall_flags_0(k-1,j,i),2) .AND. & 818 BTEST(wall_flags_0(k,j,i),2) .AND. & 819 BTEST(wall_flags_0(k+1,j,i),2) .AND. & 820 BTEST(wall_flags_0(k_pp,j,i),2) .AND. & 821 BTEST(wall_flags_0(k_ppp,j,i),2) .AND. & 822 .NOT. flag_set ) & 823 THEN 824 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 17 ) 825 ENDIF 826 782 827 ENDDO 783 828 ENDDO 784 785 DO i = nxl, nxr 786 DO j = nys, nyn 787 DO k = nzb+1, nzt 788 ! 789 !-- At first, set flags to WS1. 790 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 791 !-- in order to handle the left/south flux. 792 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 ) 793 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 ) 794 ! 795 !-- v component - x-direction 796 !-- WS1 (18), WS3 (19), WS5 (20) 797 IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),2) .OR. & 798 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 799 .AND. i == nxl ) .OR. & 800 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 801 .AND. i == nxr ) ) & 802 THEN 803 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 18 ) 804 ! 805 !-- WS3 806 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),2) .AND. & 807 BTEST(wall_flags_0(k,j,i+1),2) ) .OR. & 808 .NOT. BTEST(wall_flags_0(k,j,i-1),2) & 809 .OR. & 810 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 811 .AND. i == nxr-1 ) .OR. & 812 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 813 .AND. i == nxlu ) ) & 814 THEN 815 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 19 ) 816 ! 817 !-- Clear flag for WS1 818 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 ) 819 ELSEIF ( BTEST(wall_flags_0(k,j,i+1),2) .AND. & 820 BTEST(wall_flags_0(k,j,i+2),2) .AND. & 821 BTEST(wall_flags_0(k,j,i-1),2) ) & 822 THEN 823 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 20 ) 824 ! 825 !-- Clear flag for WS1 826 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 18 ) 827 ENDIF 828 ! 829 !-- v component - y-direction 830 !-- WS1 (21), WS3 (22), WS5 (23) 831 IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),2) .OR. & 832 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 833 .AND. j <= nysv ) .OR. & 834 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 835 .AND. j == nyn ) ) & 836 THEN 837 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 21 ) 838 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),2) .AND. & 839 BTEST(wall_flags_0(k,j+1,i),2) .OR. & 840 .NOT. BTEST(wall_flags_0(k,j-1,i),2) ) & 841 .OR. & 842 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 843 .AND. j == nysv+1) .OR. & 844 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 845 .AND. j == nyn-1 ) ) & 846 THEN 847 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 22 ) 848 ! 849 !-- Clear flag for WS1 850 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 ) 851 ELSEIF ( BTEST(wall_flags_0(k,j+1,i),2) .AND. & 852 BTEST(wall_flags_0(k,j+2,i),2) .AND. & 853 BTEST(wall_flags_0(k,j-1,i),2) ) & 854 THEN 855 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 23 ) 856 ! 857 !-- Clear flag for WS1 858 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 21 ) 859 ENDIF 860 ! 861 !-- v component - z-direction 862 !-- WS1 (24), WS3 (25), WS5 (26) 863 IF ( k == nzb+1 ) THEN 864 k_mm = nzb 865 ELSE 866 k_mm = k - 2 867 ENDIF 868 IF ( k > nzt-1 ) THEN 869 k_pp = nzt+1 870 ELSE 871 k_pp = k + 2 872 ENDIF 873 IF ( k > nzt-2 ) THEN 874 k_ppp = nzt+1 875 ELSE 876 k_ppp = k + 3 877 ENDIF 878 879 flag_set = .FALSE. 880 IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),2) .AND. & 881 BTEST(wall_flags_0(k,j,i),2) .OR. & 882 .NOT. BTEST(wall_flags_0(k_pp,j,i),2) .AND. & 883 BTEST(wall_flags_0(k,j,i),2) .OR. & 884 k == nzt ) & 885 THEN 886 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 24 ) 887 flag_set = .TRUE. 888 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),2) .OR. & 889 .NOT. BTEST(wall_flags_0(k_ppp,j,i),2) ) .AND. & 890 BTEST(wall_flags_0(k-1,j,i),2) .AND. & 891 BTEST(wall_flags_0(k,j,i),2) .AND. & 892 BTEST(wall_flags_0(k+1,j,i),2) .AND. & 893 BTEST(wall_flags_0(k_pp,j,i),2) .OR. & 894 k == nzt - 1 ) & 895 THEN 896 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 25 ) 897 flag_set = .TRUE. 898 ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),2) .AND. & 899 BTEST(wall_flags_0(k-1,j,i),2) .AND. & 900 BTEST(wall_flags_0(k,j,i),2) .AND. & 901 BTEST(wall_flags_0(k+1,j,i),2) .AND. & 902 BTEST(wall_flags_0(k_pp,j,i),2) .AND. & 903 BTEST(wall_flags_0(k_ppp,j,i),2) .AND. & 904 .NOT. flag_set ) & 905 THEN 906 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 26 ) 907 ENDIF 908 909 ENDDO 829 ENDDO 830 ! 831 !-- w - component 832 DO i = nxl, nxr 833 DO j = nys, nyn 834 DO k = nzb+1, nzt 835 ! 836 !-- At first, set flags to WS1. 837 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 838 !-- in order to handle the left/south flux. 839 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 ) 840 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 ) 841 ! 842 !-- w component - x-direction 843 !-- WS1 (18), WS3 (19), WS5 (20) 844 IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),3) .OR. & 845 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 846 .AND. i == nxl ) .OR. & 847 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 848 .AND. i == nxr ) ) & 849 THEN 850 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 18 ) 851 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),3) .AND. & 852 BTEST(wall_flags_0(k,j,i+1),3) .OR. & 853 .NOT. BTEST(wall_flags_0(k,j,i-1),3) ) & 854 .OR. & 855 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 856 .AND. i == nxr-1 ) .OR. & 857 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 858 .AND. i == nxlu ) ) & 859 THEN 860 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 19 ) 861 ! 862 !-- Clear flag for WS1 863 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 ) 864 ELSEIF ( BTEST(wall_flags_0(k,j,i+1),3) .AND. & 865 BTEST(wall_flags_0(k,j,i+2),3) .AND. & 866 BTEST(wall_flags_0(k,j,i-1),3) ) & 867 THEN 868 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i),20 ) 869 ! 870 !-- Clear flag for WS1 871 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 18 ) 872 ENDIF 873 ! 874 !-- w component - y-direction 875 !-- WS1 (21), WS3 (22), WS5 (23) 876 IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),3) .OR. & 877 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 878 .AND. j == nys ) .OR. & 879 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 880 .AND. j == nyn ) ) & 881 THEN 882 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 21 ) 883 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),3) .AND. & 884 BTEST(wall_flags_0(k,j+1,i),3) .OR. & 885 .NOT. BTEST(wall_flags_0(k,j-1,i),3) ) & 886 .OR. & 887 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 888 .AND. j == nysv ) .OR. & 889 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 890 .AND. j == nyn-1 ) ) & 891 THEN 892 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 22 ) 893 ! 894 !-- Clear flag for WS1 895 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 ) 896 ELSEIF ( BTEST(wall_flags_0(k,j+1,i),3) .AND. & 897 BTEST(wall_flags_0(k,j+2,i),3) .AND. & 898 BTEST(wall_flags_0(k,j-1,i),3) ) & 899 THEN 900 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 23 ) 901 ! 902 !-- Clear flag for WS1 903 advc_flags_m(k,j,i) = IBCLR( advc_flags_m(k,j,i), 21 ) 904 ENDIF 905 ! 906 !-- w component - z-direction 907 !-- WS1 (24), WS3 (25), WS5 (26) 908 flag_set = .FALSE. 909 IF ( k == nzb+1 ) THEN 910 k_mm = nzb 911 ELSE 912 k_mm = k - 2 913 ENDIF 914 IF ( k > nzt-1 ) THEN 915 k_pp = nzt+1 916 ELSE 917 k_pp = k + 2 918 ENDIF 919 IF ( k > nzt-2 ) THEN 920 k_ppp = nzt+1 921 ELSE 922 k_ppp = k + 3 923 ENDIF 924 925 IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. & 926 .NOT. BTEST(wall_flags_0(k,j,i),3) .AND. & 927 BTEST(wall_flags_0(k+1,j,i),3) ) .OR. & 928 ( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. & 929 BTEST(wall_flags_0(k,j,i),3) ) .OR. & 930 ( .NOT. BTEST(wall_flags_0(k+1,j,i),3) .AND. & 931 BTEST(wall_flags_0(k,j,i),3) ) .OR. & 932 k == nzt ) & 933 THEN 934 ! 935 !-- Please note, at k == nzb_w_inner(j,i) a flag is explictely 936 !-- set, although this is not a prognostic level. However, 937 !-- contrary to the advection of u,v and s this is necessary 938 !-- because flux_t(nzb_w_inner(j,i)) is used for the tendency 939 !-- at k == nzb_w_inner(j,i)+1. 940 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 ) 941 flag_set = .TRUE. 942 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3) .OR. & 943 .NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND. & 944 BTEST(wall_flags_0(k-1,j,i),3) .AND. & 945 BTEST(wall_flags_0(k,j,i),3) .AND. & 946 BTEST(wall_flags_0(k+1,j,i),3) .OR. & 947 k == nzt - 1 ) & 948 THEN 949 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 25 ) 950 flag_set = .TRUE. 951 ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),3) .AND. & 952 BTEST(wall_flags_0(k-1,j,i),3) .AND. & 953 BTEST(wall_flags_0(k,j,i),3) .AND. & 954 BTEST(wall_flags_0(k+1,j,i),3) .AND. & 955 BTEST(wall_flags_0(k_pp,j,i),3) .AND. & 956 BTEST(wall_flags_0(k_ppp,j,i),3) .AND. & 957 .NOT. flag_set ) & 958 THEN 959 advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 26 ) 960 ENDIF 961 910 962 ENDDO 911 963 ENDDO 912 DO i = nxl, nxr 913 DO j = nys, nyn 914 DO k = nzb+1, nzt 915 ! 916 !-- At first, set flags to WS1. 917 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 918 !-- in order to handle the left/south flux. 919 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 ) 920 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 ) 921 ! 922 !-- w component - x-direction 923 !-- WS1 (27), WS3 (28), WS5 (29) 924 IF ( .NOT. BTEST(wall_flags_0(k,j,i+1),3) .OR. & 925 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 926 .AND. i == nxl ) .OR. & 927 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 928 .AND. i == nxr ) ) & 929 THEN 930 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 27 ) 931 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+2),3) .AND. & 932 BTEST(wall_flags_0(k,j,i+1),3) .OR. & 933 .NOT. BTEST(wall_flags_0(k,j,i-1),3) ) & 934 .OR. & 935 ( ( bc_dirichlet_r .OR. bc_radiation_r ) & 936 .AND. i == nxr-1 ) .OR. & 937 ( ( bc_dirichlet_l .OR. bc_radiation_l ) & 938 .AND. i == nxlu ) ) & 939 THEN 940 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 28 ) 941 ! 942 !-- Clear flag for WS1 943 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 ) 944 ELSEIF ( BTEST(wall_flags_0(k,j,i+1),3) .AND. & 945 BTEST(wall_flags_0(k,j,i+2),3) .AND. & 946 BTEST(wall_flags_0(k,j,i-1),3) ) & 947 THEN 948 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i),29 ) 949 ! 950 !-- Clear flag for WS1 951 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 27 ) 964 ENDDO 965 ! 966 !-- Exchange ghost points for advection flags 967 CALL exchange_horiz_int( advc_flags_m, nys, nyn, nxl, nxr, nzt, nbgp ) 968 ! 969 !-- Set boundary flags at inflow and outflow boundary in case of 970 !-- non-cyclic boundary conditions. 971 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 972 advc_flags_m(:,:,nxl-1) = advc_flags_m(:,:,nxl) 973 ENDIF 974 975 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 976 advc_flags_m(:,:,nxr+1) = advc_flags_m(:,:,nxr) 977 ENDIF 978 979 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 980 advc_flags_m(:,nyn+1,:) = advc_flags_m(:,nyn,:) 981 ENDIF 982 983 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 984 advc_flags_m(:,nys-1,:) = advc_flags_m(:,nys,:) 985 ENDIF 986 987 END SUBROUTINE ws_init_flags_momentum 988 989 990 !------------------------------------------------------------------------------! 991 ! Description: 992 ! ------------ 993 !> Initialization of flags to control the order of the advection scheme near 994 !> solid walls and non-cyclic inflow boundaries, where the order is sucessively 995 !> degraded. 996 !------------------------------------------------------------------------------! 997 SUBROUTINE ws_init_flags_scalar( non_cyclic_l, non_cyclic_n, non_cyclic_r, & 998 non_cyclic_s, advc_flag, extensive_degrad ) 999 1000 1001 INTEGER(iwp) :: i !< index variable along x 1002 INTEGER(iwp) :: j !< index variable along y 1003 INTEGER(iwp) :: k !< index variable along z 1004 INTEGER(iwp) :: k_mm !< dummy index along z 1005 INTEGER(iwp) :: k_pp !< dummy index along z 1006 INTEGER(iwp) :: k_ppp !< dummy index along z 1007 1008 INTEGER(iwp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::& 1009 advc_flag !< flag array to control order of scalar advection 1010 1011 LOGICAL :: flag_set !< steering variable for advection flags 1012 LOGICAL :: non_cyclic_l !< flag that indicates non-cyclic boundary on the left 1013 LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north 1014 LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right 1015 LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south 1016 1017 LOGICAL, OPTIONAL :: extensive_degrad !< flag indicating that extensive degradation is required, e.g. for 1018 !< passive scalars nearby topography along the horizontal directions, 1019 !< as no monotonic limiter can be applied there 1020 ! 1021 !-- Set flags to steer the degradation of the advection scheme in advec_ws 1022 !-- near topography, inflow- and outflow boundaries as well as bottom and 1023 !-- top of model domain. advc_flags_m remains zero for all non-prognostic 1024 !-- grid points. 1025 DO i = nxl, nxr 1026 DO j = nys, nyn 1027 DO k = nzb+1, nzt 1028 IF ( .NOT. BTEST(wall_flags_0(k,j,i),0) ) CYCLE 1029 ! 1030 !-- scalar - x-direction 1031 !-- WS1 (0), WS3 (1), WS5 (2) 1032 IF ( ( .NOT. BTEST(wall_flags_0(k,j,i+1),0) & 1033 .OR. .NOT. BTEST(wall_flags_0(k,j,i+2),0) & 1034 .OR. .NOT. BTEST(wall_flags_0(k,j,i-1),0) ) & 1035 .OR. ( non_cyclic_l .AND. i == 0 ) & 1036 .OR. ( non_cyclic_r .AND. i == nx ) ) THEN 1037 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) 1038 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j,i+3),0) & 1039 .AND. BTEST(wall_flags_0(k,j,i+1),0) & 1040 .AND. BTEST(wall_flags_0(k,j,i+2),0) & 1041 .AND. BTEST(wall_flags_0(k,j,i-1),0) & 1042 ) .OR. & 1043 ( .NOT. BTEST(wall_flags_0(k,j,i-2),0) & 1044 .AND. BTEST(wall_flags_0(k,j,i+1),0) & 1045 .AND. BTEST(wall_flags_0(k,j,i+2),0) & 1046 .AND. BTEST(wall_flags_0(k,j,i-1),0) & 1047 ) & 1048 .OR. & 1049 ( non_cyclic_r .AND. i == nx-1 ) .OR. & 1050 ( non_cyclic_l .AND. i == 1 ) ) THEN 1051 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 ) 1052 ELSEIF ( BTEST(wall_flags_0(k,j,i+1),0) & 1053 .AND. BTEST(wall_flags_0(k,j,i+2),0) & 1054 .AND. BTEST(wall_flags_0(k,j,i+3),0) & 1055 .AND. BTEST(wall_flags_0(k,j,i-1),0) & 1056 .AND. BTEST(wall_flags_0(k,j,i-2),0) ) & 1057 THEN 1058 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 ) 1059 ENDIF 1060 ! 1061 !-- scalar - y-direction 1062 !-- WS1 (3), WS3 (4), WS5 (5) 1063 IF ( ( .NOT. BTEST(wall_flags_0(k,j+1,i),0) & 1064 .OR. .NOT. BTEST(wall_flags_0(k,j+2,i),0) & 1065 .OR. .NOT. BTEST(wall_flags_0(k,j-1,i),0)) & 1066 .OR. ( non_cyclic_s .AND. j == 0 ) & 1067 .OR. ( non_cyclic_n .AND. j == ny ) ) THEN 1068 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 1069 ! 1070 !-- WS3 1071 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+3,i),0) & 1072 .AND. BTEST(wall_flags_0(k,j+1,i),0) & 1073 .AND. BTEST(wall_flags_0(k,j+2,i),0) & 1074 .AND. BTEST(wall_flags_0(k,j-1,i),0) & 1075 ) .OR. & 1076 ( .NOT. BTEST(wall_flags_0(k,j-2,i),0) & 1077 .AND. BTEST(wall_flags_0(k,j+1,i),0) & 1078 .AND. BTEST(wall_flags_0(k,j+2,i),0) & 1079 .AND. BTEST(wall_flags_0(k,j-1,i),0) & 1080 ) & 1081 .OR. & 1082 ( non_cyclic_s .AND. j == 1 ) .OR. & 1083 ( non_cyclic_n .AND. j == ny-1 ) ) THEN 1084 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 ) 1085 ! 1086 !-- WS5 1087 ELSEIF ( BTEST(wall_flags_0(k,j+1,i),0) & 1088 .AND. BTEST(wall_flags_0(k,j+2,i),0) & 1089 .AND. BTEST(wall_flags_0(k,j+3,i),0) & 1090 .AND. BTEST(wall_flags_0(k,j-1,i),0) & 1091 .AND. BTEST(wall_flags_0(k,j-2,i),0) ) & 1092 THEN 1093 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 ) 1094 ENDIF 1095 ! 1096 !-- Near topography, set horizontal advection scheme to 1st order 1097 !-- for passive scalars, even if only one direction may be 1098 !-- blocked by topography. These locations will be identified 1099 !-- by wall_flags_0 bit 31. Note, since several modules define 1100 !-- advection flags but may apply different scalar boundary 1101 !-- conditions, bit 31 is temporarily stored on advc_flags. 1102 !-- Moreover, note that this extended degradtion for passive 1103 !-- scalars is not required for the vertical direction as there 1104 !-- the monotonic limiter can be applied. 1105 IF ( PRESENT( extensive_degrad ) ) THEN 1106 IF ( extensive_degrad ) THEN 1107 ! 1108 !-- At all grid points that are within a three-grid point 1109 !-- range to topography, set 1st-order scheme. 1110 IF( BTEST( advc_flag(k,j,i), 31 ) ) THEN 1111 ! 1112 !-- Clear flags that might indicate higher-order 1113 !-- advection along x- and y-direction. 1114 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) 1115 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) 1116 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) 1117 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) 1118 ! 1119 !-- Set flags that indicate 1st-order advection along 1120 !-- x- and y-direction. 1121 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) 1122 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 1123 ENDIF 1124 ! 1125 !-- Adjacent to this extended degradation zone, successively 1126 !-- upgrade the order of the scheme if this grid point isn't 1127 !-- flagged with bit 31 (indicating extended degradation 1128 !-- zone). 1129 IF ( .NOT. BTEST( advc_flag(k,j,i), 31 ) ) THEN 1130 ! 1131 !-- x-direction. First, clear all previous settings, than 1132 !-- set flag for 3rd-order scheme. 1133 IF ( BTEST( advc_flag(k,j,i-1), 31 ) .AND. & 1134 BTEST( advc_flag(k,j,i+1), 31 ) ) THEN 1135 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 ) 1136 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) 1137 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) 1138 1139 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 ) 1140 ENDIF 1141 ! 1142 !-- x-direction. First, clear all previous settings, than 1143 !-- set flag for 5rd-order scheme. 1144 IF ( .NOT. BTEST( advc_flag(k,j,i-1), 31 ) .AND. & 1145 BTEST( advc_flag(k,j,i-2), 31 ) .AND. & 1146 .NOT. BTEST( advc_flag(k,j,i+1), 31 ) .AND. & 1147 BTEST( advc_flag(k,j,i+2), 31 ) ) THEN 1148 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 ) 1149 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) 1150 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) 1151 1152 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 2 ) 1153 ENDIF 1154 ! 1155 !-- y-direction. First, clear all previous settings, than 1156 !-- set flag for 3rd-order scheme. 1157 IF ( BTEST( advc_flag(k,j-1,i), 31 ) .AND. & 1158 BTEST( advc_flag(k,j+1,i), 31 ) ) THEN 1159 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 ) 1160 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) 1161 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) 1162 1163 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 ) 1164 ENDIF 1165 ! 1166 !-- y-direction. First, clear all previous settings, than 1167 !-- set flag for 5rd-order scheme. 1168 IF ( .NOT. BTEST( advc_flag(k,j-1,i), 31 ) .AND. & 1169 BTEST( advc_flag(k,j-2,i), 31 ) .AND. & 1170 .NOT. BTEST( advc_flag(k,j+1,i), 31 ) .AND. & 1171 BTEST( advc_flag(k,j+2,i), 31 ) ) THEN 1172 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 ) 1173 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) 1174 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) 1175 1176 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 5 ) 1177 ENDIF 1178 ENDIF 1179 952 1180 ENDIF 953 ! 954 !-- w component - y-direction 955 !-- WS1 (30), WS3 (31), WS5 (32) 956 IF ( .NOT. BTEST(wall_flags_0(k,j+1,i),3) .OR. & 957 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 958 .AND. j == nys ) .OR. & 959 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 960 .AND. j == nyn ) ) & 961 THEN 962 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 30 ) 963 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k,j+2,i),3) .AND. & 964 BTEST(wall_flags_0(k,j+1,i),3) .OR. & 965 .NOT. BTEST(wall_flags_0(k,j-1,i),3) ) & 966 .OR. & 967 ( ( bc_dirichlet_s .OR. bc_radiation_s ) & 968 .AND. j == nysv ) .OR. & 969 ( ( bc_dirichlet_n .OR. bc_radiation_n ) & 970 .AND. j == nyn-1 ) ) & 971 THEN 972 advc_flags_1(k,j,i) = IBSET( advc_flags_1(k,j,i), 31 ) 973 ! 974 !-- Clear flag for WS1 975 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 ) 976 ELSEIF ( BTEST(wall_flags_0(k,j+1,i),3) .AND. & 977 BTEST(wall_flags_0(k,j+2,i),3) .AND. & 978 BTEST(wall_flags_0(k,j-1,i),3) ) & 979 THEN 980 advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 0 ) 981 ! 982 !-- Clear flag for WS1 983 advc_flags_1(k,j,i) = IBCLR( advc_flags_1(k,j,i), 30 ) 1181 1182 ! 1183 !-- Near lateral boundary flags might be overwritten. Set 1184 !-- them again. 1185 !-- x-direction 1186 IF ( ( non_cyclic_l .AND. i == 0 ) .OR. & 1187 ( non_cyclic_r .AND. i == nx ) ) THEN 1188 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 ) 1189 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) 1190 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) 1191 1192 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 0 ) 984 1193 ENDIF 985 ! 986 !-- w component - z-direction 987 !-- WS1 (33), WS3 (34), WS5 (35) 988 flag_set = .FALSE.989 IF ( k == nzb+1 ) THEN990 k_mm = nzb991 ELSE992 k_mm = k - 21194 1195 IF ( ( non_cyclic_l .AND. i == 1 ) .OR. & 1196 ( non_cyclic_r .AND. i == nx-1 ) ) THEN 1197 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 0 ) 1198 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 1 ) 1199 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 2 ) 1200 1201 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 1 ) 993 1202 ENDIF 994 IF ( k > nzt-1 ) THEN 995 k_pp = nzt+1 996 ELSE 997 k_pp = k + 2 1203 ! 1204 !-- y-direction 1205 IF ( ( non_cyclic_n .AND. j == 0 ) .OR. & 1206 ( non_cyclic_s .AND. j == ny ) ) THEN 1207 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 ) 1208 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) 1209 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) 1210 1211 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 3 ) 998 1212 ENDIF 999 IF ( k > nzt-2 ) THEN1000 k_ppp = nzt+11001 ELSE1002 k_ppp = k + 31003 ENDIF1004 1213 1005 IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. & 1006 .NOT. BTEST(wall_flags_0(k,j,i),3) .AND. & 1007 BTEST(wall_flags_0(k+1,j,i),3) ) .OR. & 1008 ( .NOT. BTEST(wall_flags_0(k-1,j,i),3) .AND. & 1009 BTEST(wall_flags_0(k,j,i),3) ) .OR. & 1010 ( .NOT. BTEST(wall_flags_0(k+1,j,i),3) .AND. & 1011 BTEST(wall_flags_0(k,j,i),3) ) .OR. & 1012 k == nzt ) & 1013 THEN 1014 ! 1015 !-- Please note, at k == nzb_w_inner(j,i) a flag is explictely 1016 !-- set, although this is not a prognostic level. However, 1017 !-- contrary to the advection of u,v and s this is necessary 1018 !-- because flux_t(nzb_w_inner(j,i)) is used for the tendency 1019 !-- at k == nzb_w_inner(j,i)+1. 1020 advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 1 ) 1021 flag_set = .TRUE. 1022 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3) .OR. & 1023 .NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND. & 1024 BTEST(wall_flags_0(k-1,j,i),3) .AND. & 1025 BTEST(wall_flags_0(k,j,i),3) .AND. & 1026 BTEST(wall_flags_0(k+1,j,i),3) .OR. & 1027 k == nzt - 1 ) & 1028 THEN 1029 advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 2 ) 1030 flag_set = .TRUE. 1031 ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),3) .AND. & 1032 BTEST(wall_flags_0(k-1,j,i),3) .AND. & 1033 BTEST(wall_flags_0(k,j,i),3) .AND. & 1034 BTEST(wall_flags_0(k+1,j,i),3) .AND. & 1035 BTEST(wall_flags_0(k_pp,j,i),3) .AND. & 1036 BTEST(wall_flags_0(k_ppp,j,i),3) .AND. & 1037 .NOT. flag_set ) & 1038 THEN 1039 advc_flags_2(k,j,i) = IBSET( advc_flags_2(k,j,i), 3 ) 1214 IF ( ( non_cyclic_n .AND. j == 1 ) .OR. & 1215 ( non_cyclic_s .AND. j == ny-1 ) ) THEN 1216 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 3 ) 1217 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 4 ) 1218 advc_flag(k,j,i) = IBCLR( advc_flag(k,j,i), 5 ) 1219 1220 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 4 ) 1040 1221 ENDIF 1041 1042 ENDDO 1222 1223 ENDIF 1224 1225 1226 ! 1227 !-- scalar - z-direction 1228 !-- WS1 (6), WS3 (7), WS5 (8) 1229 IF ( k == nzb+1 ) THEN 1230 k_mm = nzb 1231 ELSE 1232 k_mm = k - 2 1233 ENDIF 1234 IF ( k > nzt-1 ) THEN 1235 k_pp = nzt+1 1236 ELSE 1237 k_pp = k + 2 1238 ENDIF 1239 IF ( k > nzt-2 ) THEN 1240 k_ppp = nzt+1 1241 ELSE 1242 k_ppp = k + 3 1243 ENDIF 1244 1245 flag_set = .FALSE. 1246 IF ( .NOT. BTEST(wall_flags_0(k-1,j,i),0) .AND. & 1247 BTEST(wall_flags_0(k,j,i),0) .OR. & 1248 .NOT. BTEST(wall_flags_0(k_pp,j,i),0) .AND. & 1249 BTEST(wall_flags_0(k,j,i),0) .OR. & 1250 k == nzt ) & 1251 THEN 1252 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 6 ) 1253 flag_set = .TRUE. 1254 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),0) .OR. & 1255 .NOT. BTEST(wall_flags_0(k_ppp,j,i),0) ) .AND. & 1256 BTEST(wall_flags_0(k-1,j,i),0) .AND. & 1257 BTEST(wall_flags_0(k,j,i),0) .AND. & 1258 BTEST(wall_flags_0(k+1,j,i),0) .AND. & 1259 BTEST(wall_flags_0(k_pp,j,i),0) .OR. & 1260 k == nzt - 1 ) & 1261 THEN 1262 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 ) 1263 flag_set = .TRUE. 1264 ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0) & 1265 .AND. BTEST(wall_flags_0(k-1,j,i),0) & 1266 .AND. BTEST(wall_flags_0(k,j,i),0) & 1267 .AND. BTEST(wall_flags_0(k+1,j,i),0) & 1268 .AND. BTEST(wall_flags_0(k_pp,j,i),0) & 1269 .AND. BTEST(wall_flags_0(k_ppp,j,i),0) & 1270 .AND. .NOT. flag_set ) & 1271 THEN 1272 advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 8 ) 1273 ENDIF 1274 1043 1275 ENDDO 1044 1276 ENDDO 1045 1277 ENDDO 1278 ! 1279 !-- Exchange 3D integer wall_flags. 1280 ! 1281 !-- Exchange ghost points for advection flags 1282 CALL exchange_horiz_int( advc_flag, nys, nyn, nxl, nxr, nzt, nbgp ) 1283 ! 1284 !-- Set boundary flags at inflow and outflow boundary in case of 1285 !-- non-cyclic boundary conditions. 1286 IF ( non_cyclic_l ) THEN 1287 advc_flag(:,:,nxl-1) = advc_flag(:,:,nxl) 1046 1288 ENDIF 1047 1289 1048 1049 ! 1050 !-- Exchange 3D integer wall_flags. 1051 IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' & 1052 ) THEN 1053 ! 1054 !-- Exchange ghost points for advection flags 1055 CALL exchange_horiz_int( advc_flags_1, nys, nyn, nxl, nxr, nzt, nbgp ) 1056 CALL exchange_horiz_int( advc_flags_2, nys, nyn, nxl, nxr, nzt, nbgp ) 1057 ! 1058 !-- Set boundary flags at inflow and outflow boundary in case of 1059 !-- non-cyclic boundary conditions. 1060 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 1061 advc_flags_1(:,:,nxl-1) = advc_flags_1(:,:,nxl) 1062 advc_flags_2(:,:,nxl-1) = advc_flags_2(:,:,nxl) 1063 ENDIF 1064 1065 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 1066 advc_flags_1(:,:,nxr+1) = advc_flags_1(:,:,nxr) 1067 advc_flags_2(:,:,nxr+1) = advc_flags_2(:,:,nxr) 1068 ENDIF 1069 1070 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 1071 advc_flags_1(:,nyn+1,:) = advc_flags_1(:,nyn,:) 1072 advc_flags_2(:,nyn+1,:) = advc_flags_2(:,nyn,:) 1073 ENDIF 1074 1075 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 1076 advc_flags_1(:,nys-1,:) = advc_flags_1(:,nys,:) 1077 advc_flags_2(:,nys-1,:) = advc_flags_2(:,nys,:) 1078 ENDIF 1290 IF ( non_cyclic_r ) THEN 1291 advc_flag(:,:,nxr+1) = advc_flag(:,:,nxr) 1292 ENDIF 1293 1294 IF ( non_cyclic_n ) THEN 1295 advc_flag(:,nyn+1,:) = advc_flag(:,nyn,:) 1296 ENDIF 1297 1298 IF ( non_cyclic_s ) THEN 1299 advc_flag(:,nys-1,:) = advc_flag(:,nys,:) 1300 ENDIF 1079 1301 1080 ENDIF 1081 1082 1083 END SUBROUTINE ws_init_flags 1084 1085 1302 1303 1304 END SUBROUTINE ws_init_flags_scalar 1305 1086 1306 !------------------------------------------------------------------------------! 1087 1307 ! Description: … … 1123 1343 !> Scalar advection - Call for grid point i,j 1124 1344 !------------------------------------------------------------------------------! 1125 SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local,&1345 SUBROUTINE advec_s_ws_ij( advc_flag, i, j, sk, sk_char, swap_flux_y_local, & 1126 1346 swap_diss_y_local, swap_flux_x_local, & 1127 swap_diss_x_local, i_omp, tn, flux_limitation ) 1347 swap_diss_x_local, i_omp, tn, & 1348 non_cyclic_l, non_cyclic_n, & 1349 non_cyclic_r, non_cyclic_s, & 1350 flux_limitation ) 1128 1351 1129 1352 … … 1140 1363 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 1141 1364 INTEGER(iwp) :: tn !< number of OpenMP thread 1142 1365 1366 INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: & 1367 advc_flag !< flag array to control order of scalar advection 1368 1369 LOGICAL :: non_cyclic_l !< flag that indicates non-cyclic boundary on the left 1370 LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north 1371 LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right 1372 LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south 1143 1373 LOGICAL, OPTIONAL :: flux_limitation !< flag indicating flux limitation of the vertical advection 1144 1374 LOGICAL :: limiter !< control flag indicating the application of flux limitation … … 1196 1426 !-- instead of the entire subdomain. This should lead to better 1197 1427 !-- load balance between boundary and non-boundary PEs. 1198 IF( ( bc_dirichlet_l .OR. bc_radiation_l ) .AND. i <= nxl + 2 .OR.&1199 ( bc_dirichlet_r .OR. bc_radiation_r ) .AND. i >= nxr - 2 .OR.&1200 ( bc_dirichlet_s .OR. bc_radiation_s ) .AND. j <= nys + 2 .OR.&1201 ( bc_dirichlet_n .OR. bc_radiation_n ).AND. j >= nyn - 2 ) THEN1428 IF( non_cyclic_l .AND. i <= nxl + 2 .OR. & 1429 non_cyclic_r .AND. i >= nxr - 2 .OR. & 1430 non_cyclic_s .AND. j <= nys + 2 .OR. & 1431 non_cyclic_n .AND. j >= nyn - 2 ) THEN 1202 1432 nzb_max_l = nzt 1203 1433 ELSE … … 1215 1445 DO k = nzb+1, nzb_max_l 1216 1446 1217 ibit5 = REAL( IBITS(advc_flag s_1(k,j-1,i),5,1), KIND = wp )1218 ibit4 = REAL( IBITS(advc_flag s_1(k,j-1,i),4,1), KIND = wp )1219 ibit3 = REAL( IBITS(advc_flag s_1(k,j-1,i),3,1), KIND = wp )1447 ibit5 = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) 1448 ibit4 = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) 1449 ibit3 = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) 1220 1450 1221 1451 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) … … 1276 1506 DO k = nzb+1, nzb_max_l 1277 1507 1278 ibit2 = REAL( IBITS(advc_flag s_1(k,j,i-1),2,1), KIND = wp )1279 ibit1 = REAL( IBITS(advc_flag s_1(k,j,i-1),1,1), KIND = wp )1280 ibit0 = REAL( IBITS(advc_flag s_1(k,j,i-1),0,1), KIND = wp )1508 ibit2 = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) 1509 ibit1 = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) 1510 ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) 1281 1511 1282 1512 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) … … 1340 1570 !-- flux at the end. 1341 1571 1342 ibit2 = REAL( IBITS(advc_flag s_1(k,j,i),2,1), KIND = wp )1343 ibit1 = REAL( IBITS(advc_flag s_1(k,j,i),1,1), KIND = wp )1344 ibit0 = REAL( IBITS(advc_flag s_1(k,j,i),0,1), KIND = wp )1572 ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp ) 1573 ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp ) 1574 ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp ) 1345 1575 1346 1576 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) … … 1375 1605 ) 1376 1606 1377 ibit5 = REAL( IBITS(advc_flag s_1(k,j,i),5,1), KIND = wp )1378 ibit4 = REAL( IBITS(advc_flag s_1(k,j,i),4,1), KIND = wp )1379 ibit3 = REAL( IBITS(advc_flag s_1(k,j,i),3,1), KIND = wp )1607 ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp ) 1608 ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp ) 1609 ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp ) 1380 1610 1381 1611 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) … … 1445 1675 !-- calculated explicetely for the tendency at 1446 1676 !-- the first w-level. For topography wall this is done implicitely by 1447 !-- advc_flag s_1.1677 !-- advc_flag. 1448 1678 flux_t(nzb) = 0.0_wp 1449 1679 diss_t(nzb) = 0.0_wp 1450 1680 DO k = nzb+1, nzb+2 1451 ibit8 = REAL( IBITS(advc_flag s_1(k,j,i),8,1), KIND = wp )1452 ibit7 = REAL( IBITS(advc_flag s_1(k,j,i),7,1), KIND = wp )1453 ibit6 = REAL( IBITS(advc_flag s_1(k,j,i),6,1), KIND = wp )1681 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) 1682 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) 1683 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) 1454 1684 ! 1455 1685 !-- k index has to be modified near bottom and top, else array … … 1491 1721 1492 1722 DO k = nzb+3, nzt-2 1493 ibit8 = REAL( IBITS(advc_flag s_1(k,j,i),8,1), KIND = wp )1494 ibit7 = REAL( IBITS(advc_flag s_1(k,j,i),7,1), KIND = wp )1495 ibit6 = REAL( IBITS(advc_flag s_1(k,j,i),6,1), KIND = wp )1723 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) 1724 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) 1725 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) 1496 1726 1497 1727 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & … … 1526 1756 1527 1757 DO k = nzt-1, nzt 1528 ibit8 = REAL( IBITS(advc_flag s_1(k,j,i),8,1), KIND = wp )1529 ibit7 = REAL( IBITS(advc_flag s_1(k,j,i),7,1), KIND = wp )1530 ibit6 = REAL( IBITS(advc_flag s_1(k,j,i),6,1), KIND = wp )1758 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) 1759 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) 1760 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) 1531 1761 ! 1532 1762 !-- k index has to be modified near bottom and top, else array … … 1675 1905 diss_d = diss_t(k-1) 1676 1906 1677 ibit2 = REAL( IBITS(advc_flag s_1(k,j,i),2,1), KIND = wp )1678 ibit1 = REAL( IBITS(advc_flag s_1(k,j,i),1,1), KIND = wp )1679 ibit0 = REAL( IBITS(advc_flag s_1(k,j,i),0,1), KIND = wp )1907 ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp ) 1908 ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp ) 1909 ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp ) 1680 1910 1681 ibit5 = REAL( IBITS(advc_flag s_1(k,j,i),5,1), KIND = wp )1682 ibit4 = REAL( IBITS(advc_flag s_1(k,j,i),4,1), KIND = wp )1683 ibit3 = REAL( IBITS(advc_flag s_1(k,j,i),3,1), KIND = wp )1911 ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp ) 1912 ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp ) 1913 ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp ) 1684 1914 1685 ibit8 = REAL( IBITS(advc_flag s_1(k,j,i),8,1), KIND = wp )1686 ibit7 = REAL( IBITS(advc_flag s_1(k,j,i),7,1), KIND = wp )1687 ibit6 = REAL( IBITS(advc_flag s_1(k,j,i),6,1), KIND = wp )1915 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) 1916 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) 1917 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) 1688 1918 ! 1689 1919 !-- Calculate the divergence of the velocity field. A respective … … 1692 1922 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & 1693 1923 - u(k,j,i) * ( & 1694 REAL( IBITS(advc_flag s_1(k,j,i-1),0,1), KIND = wp )&1695 + REAL( IBITS(advc_flag s_1(k,j,i-1),1,1), KIND = wp )&1696 + REAL( IBITS(advc_flag s_1(k,j,i-1),2,1), KIND = wp )&1924 REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) & 1925 + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) & 1926 + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) & 1697 1927 ) & 1698 1928 ) * ddx & 1699 1929 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 1700 1930 - v(k,j,i) * ( & 1701 REAL( IBITS(advc_flag s_1(k,j-1,i),3,1), KIND = wp )&1702 + REAL( IBITS(advc_flag s_1(k,j-1,i),4,1), KIND = wp )&1703 + REAL( IBITS(advc_flag s_1(k,j-1,i),5,1), KIND = wp )&1931 REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) & 1932 + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) & 1933 + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) & 1704 1934 ) & 1705 1935 ) * ddy & … … 1708 1938 - w(k-1,j,i) * rho_air_zw(k-1) * & 1709 1939 ( & 1710 REAL( IBITS(advc_flag s_1(k-1,j,i),6,1), KIND = wp )&1711 + REAL( IBITS(advc_flag s_1(k-1,j,i),7,1), KIND = wp )&1712 + REAL( IBITS(advc_flag s_1(k-1,j,i),8,1), KIND = wp )&1940 REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp ) & 1941 + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp ) & 1942 + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp ) & 1713 1943 ) & 1714 1944 ) * drho_air(k) * ddzw(k) … … 1897 2127 INTEGER(iwp) :: tn !< number of OpenMP thread 1898 2128 1899 REAL(wp) :: ibit 9!< flag indicating 1st-order scheme along x-direction1900 REAL(wp) :: ibit1 0!< flag indicating 3rd-order scheme along x-direction1901 REAL(wp) :: ibit 11!< flag indicating 5th-order scheme along x-direction1902 REAL(wp) :: ibit 12!< flag indicating 1st-order scheme along y-direction1903 REAL(wp) :: ibit 13!< flag indicating 3rd-order scheme along y-direction1904 REAL(wp) :: ibit 14!< flag indicating 5th-order scheme along y-direction1905 REAL(wp) :: ibit 15!< flag indicating 1st-order scheme along z-direction1906 REAL(wp) :: ibit 16!< flag indicating 3rd-order scheme along z-direction1907 REAL(wp) :: ibit 17!< flag indicating 5th-order scheme along z-direction2129 REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction 2130 REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction 2131 REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction 2132 REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction 2133 REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction 2134 REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction 2135 REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction 2136 REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction 2137 REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction 1908 2138 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 1909 2139 REAL(wp) :: div !< diverence on u-grid … … 1944 2174 DO k = nzb+1, nzb_max_l 1945 2175 1946 ibit 14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp )1947 ibit 13 = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp )1948 ibit 12 = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp )2176 ibit5 = REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp ) 2177 ibit4 = REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp ) 2178 ibit3 = REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp ) 1949 2179 1950 2180 v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv 1951 flux_s_u(k,tn) = v_comp(k) * ( &1952 ( 37.0_wp * ibit 14 * adv_mom_5&1953 + 7.0_wp * ibit 13 * adv_mom_3&1954 + ibit 12 * adv_mom_1&1955 ) * &1956 ( u(k,j,i) + u(k,j-1,i) ) &1957 - ( 8.0_wp * ibit 14 * adv_mom_5&1958 + ibit 13 * adv_mom_3&1959 ) * &1960 ( u(k,j+1,i) + u(k,j-2,i) ) &1961 + ( ibit 14 * adv_mom_5&1962 ) * &1963 ( u(k,j+2,i) + u(k,j-3,i) ) &2181 flux_s_u(k,tn) = v_comp(k) * ( & 2182 ( 37.0_wp * ibit5 * adv_mom_5 & 2183 + 7.0_wp * ibit4 * adv_mom_3 & 2184 + ibit3 * adv_mom_1 & 2185 ) * & 2186 ( u(k,j,i) + u(k,j-1,i) ) & 2187 - ( 8.0_wp * ibit5 * adv_mom_5 & 2188 + ibit4 * adv_mom_3 & 2189 ) * & 2190 ( u(k,j+1,i) + u(k,j-2,i) ) & 2191 + ( ibit5 * adv_mom_5 & 2192 ) * & 2193 ( u(k,j+2,i) + u(k,j-3,i) ) & 1964 2194 ) 1965 2195 1966 diss_s_u(k,tn) = - ABS ( v_comp(k) ) * ( &1967 ( 10.0_wp * ibit 14 * adv_mom_5&1968 + 3.0_wp * ibit 13 * adv_mom_3&1969 + ibit 12 * adv_mom_1&1970 ) * &1971 ( u(k,j,i) - u(k,j-1,i) ) &1972 - ( 5.0_wp * ibit 14 * adv_mom_5&1973 + ibit 13 * adv_mom_3&1974 ) * &1975 ( u(k,j+1,i) - u(k,j-2,i) ) &1976 + ( ibit 14 * adv_mom_5&1977 ) * &1978 ( u(k,j+2,i) - u(k,j-3,i) ) &2196 diss_s_u(k,tn) = - ABS ( v_comp(k) ) * ( & 2197 ( 10.0_wp * ibit5 * adv_mom_5 & 2198 + 3.0_wp * ibit4 * adv_mom_3 & 2199 + ibit3 * adv_mom_1 & 2200 ) * & 2201 ( u(k,j,i) - u(k,j-1,i) ) & 2202 - ( 5.0_wp * ibit5 * adv_mom_5 & 2203 + ibit4 * adv_mom_3 & 2204 ) * & 2205 ( u(k,j+1,i) - u(k,j-2,i) ) & 2206 + ( ibit5 * adv_mom_5 & 2207 ) * & 2208 ( u(k,j+2,i) - u(k,j-3,i) ) & 1979 2209 ) 1980 2210 … … 1984 2214 1985 2215 v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv 1986 flux_s_u(k,tn) = v_comp(k) * ( &1987 37.0_wp * ( u(k,j,i) + u(k,j-1,i) )&1988 - 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) &1989 + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 1990 diss_s_u(k,tn) = - ABS(v_comp(k)) * ( &1991 10.0_wp * ( u(k,j,i) - u(k,j-1,i) )&1992 - 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) &2216 flux_s_u(k,tn) = v_comp(k) * ( & 2217 37.0_wp * ( u(k,j,i) + u(k,j-1,i) ) & 2218 - 8.0_wp * ( u(k,j+1,i) + u(k,j-2,i) ) & 2219 + ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5 2220 diss_s_u(k,tn) = - ABS(v_comp(k)) * ( & 2221 10.0_wp * ( u(k,j,i) - u(k,j-1,i) ) & 2222 - 5.0_wp * ( u(k,j+1,i) - u(k,j-2,i) ) & 1993 2223 + ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5 1994 2224 … … 2002 2232 DO k = nzb+1, nzb_max_l 2003 2233 2004 ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp ) 2005 ibit10 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp ) 2006 ibit9 = REAL( IBITS(advc_flags_1(k,j,i-1),9,1), KIND = wp ) 2007 2008 u_comp_l = u(k,j,i) + u(k,j,i-1) - gu 2009 flux_l_u(k,j,tn) = u_comp_l * ( & 2010 ( 37.0_wp * ibit11 * adv_mom_5 & 2011 + 7.0_wp * ibit10 * adv_mom_3 & 2012 + ibit9 * adv_mom_1 & 2013 ) * & 2014 ( u(k,j,i) + u(k,j,i-1) ) & 2015 - ( 8.0_wp * ibit11 * adv_mom_5 & 2016 + ibit10 * adv_mom_3 & 2017 ) * & 2018 ( u(k,j,i+1) + u(k,j,i-2) ) & 2019 + ( ibit11 * adv_mom_5 & 2020 ) * & 2021 ( u(k,j,i+2) + u(k,j,i-3) ) & 2022 ) 2023 2024 diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( & 2025 ( 10.0_wp * ibit11 * adv_mom_5 & 2026 + 3.0_wp * ibit10 * adv_mom_3 & 2027 + ibit9 * adv_mom_1 & 2028 ) * & 2029 ( u(k,j,i) - u(k,j,i-1) ) & 2030 - ( 5.0_wp * ibit11 * adv_mom_5 & 2031 + ibit10 * adv_mom_3 & 2032 ) * & 2033 ( u(k,j,i+1) - u(k,j,i-2) ) & 2034 + ( ibit11 * adv_mom_5 & 2035 ) * & 2036 ( u(k,j,i+2) - u(k,j,i-3) ) & 2037 ) 2038 2039 ENDDO 2040 2041 DO k = nzb_max_l+1, nzt 2234 ibit2 = REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp ) 2235 ibit1 = REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp ) 2236 ibit0 = REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp ) 2042 2237 2043 2238 u_comp_l = u(k,j,i) + u(k,j,i-1) - gu 2044 2239 flux_l_u(k,j,tn) = u_comp_l * ( & 2045 37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) & 2240 ( 37.0_wp * ibit2 * adv_mom_5 & 2241 + 7.0_wp * ibit1 * adv_mom_3 & 2242 + ibit0 * adv_mom_1 & 2243 ) * & 2244 ( u(k,j,i) + u(k,j,i-1) ) & 2245 - ( 8.0_wp * ibit2 * adv_mom_5 & 2246 + ibit1 * adv_mom_3 & 2247 ) * & 2248 ( u(k,j,i+1) + u(k,j,i-2) ) & 2249 + ( ibit2 * adv_mom_5 & 2250 ) * & 2251 ( u(k,j,i+2) + u(k,j,i-3) ) & 2252 ) 2253 2254 diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( & 2255 ( 10.0_wp * ibit2 * adv_mom_5 & 2256 + 3.0_wp * ibit1 * adv_mom_3 & 2257 + ibit0 * adv_mom_1 & 2258 ) * & 2259 ( u(k,j,i) - u(k,j,i-1) ) & 2260 - ( 5.0_wp * ibit2 * adv_mom_5 & 2261 + ibit1 * adv_mom_3 & 2262 ) * & 2263 ( u(k,j,i+1) - u(k,j,i-2) ) & 2264 + ( ibit2 * adv_mom_5 & 2265 ) * & 2266 ( u(k,j,i+2) - u(k,j,i-3) ) & 2267 ) 2268 2269 ENDDO 2270 2271 DO k = nzb_max_l+1, nzt 2272 2273 u_comp_l = u(k,j,i) + u(k,j,i-1) - gu 2274 flux_l_u(k,j,tn) = u_comp_l * ( & 2275 37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) & 2046 2276 - 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) & 2047 2277 + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 2048 2278 diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( & 2049 10.0_wp * ( u(k,j,i) - u(k,j,i-1) )&2279 10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) & 2050 2280 - 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) & 2051 2281 + ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5 … … 2059 2289 DO k = nzb+1, nzb_max_l 2060 2290 2061 ibit 11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )2062 ibit1 0 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp )2063 ibit 9 = REAL( IBITS(advc_flags_1(k,j,i),9,1),KIND = wp )2291 ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp ) 2292 ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp ) 2293 ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp ) 2064 2294 2065 2295 u_comp(k) = u(k,j,i+1) + u(k,j,i) 2066 flux_r(k) = ( u_comp(k) - gu ) * ( &2067 ( 37.0_wp * ibit 11 * adv_mom_5&2068 + 7.0_wp * ibit1 0 * adv_mom_3&2069 + ibit 9 * adv_mom_1&2070 ) * &2071 ( u(k,j,i+1) + u(k,j,i) ) &2072 - ( 8.0_wp * ibit 11 * adv_mom_5&2073 + ibit1 0 * adv_mom_3&2074 ) * &2075 ( u(k,j,i+2) + u(k,j,i-1) ) &2076 + ( ibit 11 * adv_mom_5&2077 ) * &2078 ( u(k,j,i+3) + u(k,j,i-2) ) &2079 ) 2080 2081 diss_r(k) = - ABS( u_comp(k) - gu ) * ( &2082 ( 10.0_wp * ibit 11 * adv_mom_5&2083 + 3.0_wp * ibit1 0 * adv_mom_3&2084 + ibit 9 * adv_mom_1&2085 ) * &2086 ( u(k,j,i+1) - u(k,j,i) ) &2087 - ( 5.0_wp * ibit 11 * adv_mom_5&2088 + ibit1 0 * adv_mom_3&2089 ) * &2090 ( u(k,j,i+2) - u(k,j,i-1) ) &2091 + ( ibit 11 * adv_mom_5&2092 ) * &2093 ( u(k,j,i+3) - u(k,j,i-2) ) &2296 flux_r(k) = ( u_comp(k) - gu ) * ( & 2297 ( 37.0_wp * ibit2 * adv_mom_5 & 2298 + 7.0_wp * ibit1 * adv_mom_3 & 2299 + ibit0 * adv_mom_1 & 2300 ) * & 2301 ( u(k,j,i+1) + u(k,j,i) ) & 2302 - ( 8.0_wp * ibit2 * adv_mom_5 & 2303 + ibit1 * adv_mom_3 & 2304 ) * & 2305 ( u(k,j,i+2) + u(k,j,i-1) ) & 2306 + ( ibit2 * adv_mom_5 & 2307 ) * & 2308 ( u(k,j,i+3) + u(k,j,i-2) ) & 2309 ) 2310 2311 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 2312 ( 10.0_wp * ibit2 * adv_mom_5 & 2313 + 3.0_wp * ibit1 * adv_mom_3 & 2314 + ibit0 * adv_mom_1 & 2315 ) * & 2316 ( u(k,j,i+1) - u(k,j,i) ) & 2317 - ( 5.0_wp * ibit2 * adv_mom_5 & 2318 + ibit1 * adv_mom_3 & 2319 ) * & 2320 ( u(k,j,i+2) - u(k,j,i-1) ) & 2321 + ( ibit2 * adv_mom_5 & 2322 ) * & 2323 ( u(k,j,i+3) - u(k,j,i-2) ) & 2094 2324 ) 2095 2325 2096 ibit 14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )2097 ibit 13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )2098 ibit 12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )2326 ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp ) 2327 ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp ) 2328 ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp ) 2099 2329 2100 2330 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv 2101 flux_n(k) = v_comp(k) * ( &2102 ( 37.0_wp * ibit 14 * adv_mom_5&2103 + 7.0_wp * ibit 13 * adv_mom_3&2104 + ibit 12 * adv_mom_1&2105 ) * &2106 ( u(k,j+1,i) + u(k,j,i) ) &2107 - ( 8.0_wp * ibit 14 * adv_mom_5&2108 + ibit 13 * adv_mom_3&2109 ) * &2110 ( u(k,j+2,i) + u(k,j-1,i) ) &2111 + ( ibit 14 * adv_mom_5&2112 ) * &2113 ( u(k,j+3,i) + u(k,j-2,i) ) &2114 ) 2115 2116 diss_n(k) = - ABS ( v_comp(k) ) * ( &2117 ( 10.0_wp * ibit 14 * adv_mom_5&2118 + 3.0_wp * ibit 13 * adv_mom_3&2119 + ibit 12 * adv_mom_1&2120 ) * &2121 ( u(k,j+1,i) - u(k,j,i) ) &2122 - ( 5.0_wp * ibit 14 * adv_mom_5&2123 + ibit 13 * adv_mom_3&2124 ) * &2125 ( u(k,j+2,i) - u(k,j-1,i) ) &2126 + ( ibit 14 * adv_mom_5&2127 ) * &2128 ( u(k,j+3,i) - u(k,j-2,i) ) &2331 flux_n(k) = v_comp(k) * ( & 2332 ( 37.0_wp * ibit5 * adv_mom_5 & 2333 + 7.0_wp * ibit4 * adv_mom_3 & 2334 + ibit3 * adv_mom_1 & 2335 ) * & 2336 ( u(k,j+1,i) + u(k,j,i) ) & 2337 - ( 8.0_wp * ibit5 * adv_mom_5 & 2338 + ibit4 * adv_mom_3 & 2339 ) * & 2340 ( u(k,j+2,i) + u(k,j-1,i) ) & 2341 + ( ibit5 * adv_mom_5 & 2342 ) * & 2343 ( u(k,j+3,i) + u(k,j-2,i) ) & 2344 ) 2345 2346 diss_n(k) = - ABS ( v_comp(k) ) * ( & 2347 ( 10.0_wp * ibit5 * adv_mom_5 & 2348 + 3.0_wp * ibit4 * adv_mom_3 & 2349 + ibit3 * adv_mom_1 & 2350 ) * & 2351 ( u(k,j+1,i) - u(k,j,i) ) & 2352 - ( 5.0_wp * ibit5 * adv_mom_5 & 2353 + ibit4 * adv_mom_3 & 2354 ) * & 2355 ( u(k,j+2,i) - u(k,j-1,i) ) & 2356 + ( ibit5 * adv_mom_5 & 2357 ) * & 2358 ( u(k,j+3,i) - u(k,j-2,i) ) & 2129 2359 ) 2130 2360 ENDDO … … 2133 2363 2134 2364 u_comp(k) = u(k,j,i+1) + u(k,j,i) 2135 flux_r(k) = ( u_comp(k) - gu ) * ( &2136 37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) &2137 - 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) &2138 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 2139 diss_r(k) = - ABS( u_comp(k) - gu ) * ( &2140 10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) &2141 - 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) &2142 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 2143 2144 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv 2145 flux_n(k) = v_comp(k) * ( &2146 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) &2147 - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) &2148 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 2149 diss_n(k) = - ABS( v_comp(k) ) * ( &2150 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) &2151 - 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) &2365 flux_r(k) = ( u_comp(k) - gu ) * ( & 2366 37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) & 2367 - 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) & 2368 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 2369 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 2370 10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) & 2371 - 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) & 2372 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 2373 2374 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv 2375 flux_n(k) = v_comp(k) * ( & 2376 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) & 2377 - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) & 2378 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 2379 diss_n(k) = - ABS( v_comp(k) ) * ( & 2380 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) & 2381 - 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) & 2152 2382 + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 2153 2383 … … 2161 2391 !-- calculated explicetely for the tendency at 2162 2392 !-- the first w-level. For topography wall this is done implicitely by 2163 !-- advc_flags_ 1.2393 !-- advc_flags_m. 2164 2394 flux_t(nzb) = 0.0_wp 2165 2395 diss_t(nzb) = 0.0_wp … … 2169 2399 !-- k index has to be modified near bottom and top, else array 2170 2400 !-- subscripts will be exceeded. 2171 ibit 17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )2172 ibit 16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )2173 ibit 15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )2174 2175 k_ppp = k + 3 * ibit 172176 k_pp = k + 2 * ( 1 - ibit 15)2177 k_mm = k - 2 * ibit 172401 ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp ) 2402 ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp ) 2403 ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp ) 2404 2405 k_ppp = k + 3 * ibit8 2406 k_pp = k + 2 * ( 1 - ibit6 ) 2407 k_mm = k - 2 * ibit8 2178 2408 2179 2409 w_comp(k) = w(k,j,i) + w(k,j,i-1) 2180 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &2181 ( 37.0_wp * ibit 17 * adv_mom_5&2182 + 7.0_wp * ibit 16 * adv_mom_3&2183 + ibit 15 * adv_mom_1&2184 ) * &2185 ( u(k+1,j,i) + u(k,j,i) )&2186 - ( 8.0_wp * ibit 17 * adv_mom_5&2187 + ibit 16 * adv_mom_3&2188 ) * &2189 ( u(k_pp,j,i) + u(k-1,j,i) )&2190 + ( ibit 17 * adv_mom_5&2191 ) * &2192 ( u(k_ppp,j,i) + u(k_mm,j,i) ) &2193 ) 2194 2195 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &2196 ( 10.0_wp * ibit 17 * adv_mom_5&2197 + 3.0_wp * ibit 16 * adv_mom_3&2198 + ibit 15 * adv_mom_1&2199 ) * &2200 ( u(k+1,j,i) - u(k,j,i) ) &2201 - ( 5.0_wp * ibit 17 * adv_mom_5&2202 + ibit 16 * adv_mom_3&2203 ) * &2204 ( u(k_pp,j,i) - u(k-1,j,i) ) &2205 + ( ibit 17 * adv_mom_5&2206 ) * &2207 ( u(k_ppp,j,i) - u(k_mm,j,i) ) &2410 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2411 ( 37.0_wp * ibit8 * adv_mom_5 & 2412 + 7.0_wp * ibit7 * adv_mom_3 & 2413 + ibit6 * adv_mom_1 & 2414 ) * & 2415 ( u(k+1,j,i) + u(k,j,i) ) & 2416 - ( 8.0_wp * ibit8 * adv_mom_5 & 2417 + ibit7 * adv_mom_3 & 2418 ) * & 2419 ( u(k_pp,j,i) + u(k-1,j,i) ) & 2420 + ( ibit8 * adv_mom_5 & 2421 ) * & 2422 ( u(k_ppp,j,i) + u(k_mm,j,i) ) & 2423 ) 2424 2425 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2426 ( 10.0_wp * ibit8 * adv_mom_5 & 2427 + 3.0_wp * ibit7 * adv_mom_3 & 2428 + ibit6 * adv_mom_1 & 2429 ) * & 2430 ( u(k+1,j,i) - u(k,j,i) ) & 2431 - ( 5.0_wp * ibit8 * adv_mom_5 & 2432 + ibit7 * adv_mom_3 & 2433 ) * & 2434 ( u(k_pp,j,i) - u(k-1,j,i) ) & 2435 + ( ibit8 * adv_mom_5 & 2436 ) * & 2437 ( u(k_ppp,j,i) - u(k_mm,j,i) ) & 2208 2438 ) 2209 2439 ENDDO … … 2211 2441 DO k = nzb+3, nzt-2 2212 2442 2213 ibit 17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )2214 ibit 16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )2215 ibit 15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )2443 ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp ) 2444 ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp ) 2445 ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp ) 2216 2446 2217 2447 w_comp(k) = w(k,j,i) + w(k,j,i-1) 2218 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &2219 ( 37.0_wp * ibit 17 * adv_mom_5&2220 + 7.0_wp * ibit 16 * adv_mom_3&2221 + ibit 15 * adv_mom_1&2222 ) * &2223 ( u(k+1,j,i) + u(k,j,i) ) &2224 - ( 8.0_wp * ibit 17 * adv_mom_5&2225 + ibit 16 * adv_mom_3&2226 ) * &2227 ( u(k+2,j,i) + u(k-1,j,i) ) &2228 + ( ibit 17 * adv_mom_5 &2229 ) * &2230 ( u(k+3,j,i) + u(k-2,j,i) ) &2448 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2449 ( 37.0_wp * ibit8 * adv_mom_5 & 2450 + 7.0_wp * ibit7 * adv_mom_3 & 2451 + ibit6 * adv_mom_1 & 2452 ) * & 2453 ( u(k+1,j,i) + u(k,j,i) ) & 2454 - ( 8.0_wp * ibit8 * adv_mom_5 & 2455 + ibit7 * adv_mom_3 & 2456 ) * & 2457 ( u(k+2,j,i) + u(k-1,j,i) ) & 2458 + ( ibit8 * adv_mom_5 & 2459 ) * & 2460 ( u(k+3,j,i) + u(k-2,j,i) ) & 2231 2461 ) 2232 2462 2233 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &2234 ( 10.0_wp * ibit 17 * adv_mom_5&2235 + 3.0_wp * ibit 16 * adv_mom_3&2236 + ibit 15 * adv_mom_1&2237 ) * &2238 ( u(k+1,j,i) - u(k,j,i) )&2239 - ( 5.0_wp * ibit 17 * adv_mom_5&2240 + ibit 16 * adv_mom_3&2241 ) * &2242 ( u(k+2,j,i) - u(k-1,j,i) ) &2243 + ( ibit 17 * adv_mom_5&2244 ) * &2245 ( u(k+3,j,i) - u(k-2,j,i) )&2463 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2464 ( 10.0_wp * ibit8 * adv_mom_5 & 2465 + 3.0_wp * ibit7 * adv_mom_3 & 2466 + ibit6 * adv_mom_1 & 2467 ) * & 2468 ( u(k+1,j,i) - u(k,j,i) ) & 2469 - ( 5.0_wp * ibit8 * adv_mom_5 & 2470 + ibit7 * adv_mom_3 & 2471 ) * & 2472 ( u(k+2,j,i) - u(k-1,j,i) ) & 2473 + ( ibit8 * adv_mom_5 & 2474 ) * & 2475 ( u(k+3,j,i) - u(k-2,j,i) ) & 2246 2476 ) 2247 2477 ENDDO … … 2251 2481 !-- k index has to be modified near bottom and top, else array 2252 2482 !-- subscripts will be exceeded. 2253 ibit 17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )2254 ibit 16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )2255 ibit 15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )2256 2257 k_ppp = k + 3 * ibit 172258 k_pp = k + 2 * ( 1 - ibit 15)2259 k_mm = k - 2 * ibit 172483 ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp ) 2484 ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp ) 2485 ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp ) 2486 2487 k_ppp = k + 3 * ibit8 2488 k_pp = k + 2 * ( 1 - ibit6 ) 2489 k_mm = k - 2 * ibit8 2260 2490 2261 2491 w_comp(k) = w(k,j,i) + w(k,j,i-1) 2262 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &2263 ( 37.0_wp * ibit 17 * adv_mom_5&2264 + 7.0_wp * ibit 16 * adv_mom_3&2265 + ibit 15 * adv_mom_1&2266 ) * &2267 ( u(k+1,j,i) + u(k,j,i) )&2268 - ( 8.0_wp * ibit 17 * adv_mom_5&2269 + ibit 16 * adv_mom_3&2270 ) * &2271 ( u(k_pp,j,i) + u(k-1,j,i) )&2272 + ( ibit 17 * adv_mom_5&2273 ) * &2274 ( u(k_ppp,j,i) + u(k_mm,j,i) ) &2492 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2493 ( 37.0_wp * ibit8 * adv_mom_5 & 2494 + 7.0_wp * ibit7 * adv_mom_3 & 2495 + ibit6 * adv_mom_1 & 2496 ) * & 2497 ( u(k+1,j,i) + u(k,j,i) ) & 2498 - ( 8.0_wp * ibit8 * adv_mom_5 & 2499 + ibit7 * adv_mom_3 & 2500 ) * & 2501 ( u(k_pp,j,i) + u(k-1,j,i) ) & 2502 + ( ibit8 * adv_mom_5 & 2503 ) * & 2504 ( u(k_ppp,j,i) + u(k_mm,j,i) ) & 2275 2505 ) 2276 2506 2277 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &2278 ( 10.0_wp * ibit 17 * adv_mom_5&2279 + 3.0_wp * ibit 16 * adv_mom_3&2280 + ibit 15 * adv_mom_1&2281 ) * &2282 ( u(k+1,j,i) - u(k,j,i) ) &2283 - ( 5.0_wp * ibit 17 * adv_mom_5&2284 + ibit 16 * adv_mom_3&2285 ) * &2286 ( u(k_pp,j,i) - u(k-1,j,i) ) &2287 + ( ibit 17 * adv_mom_5&2288 ) * &2289 ( u(k_ppp,j,i) - u(k_mm,j,i) ) &2507 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2508 ( 10.0_wp * ibit8 * adv_mom_5 & 2509 + 3.0_wp * ibit7 * adv_mom_3 & 2510 + ibit6 * adv_mom_1 & 2511 ) * & 2512 ( u(k+1,j,i) - u(k,j,i) ) & 2513 - ( 5.0_wp * ibit8 * adv_mom_5 & 2514 + ibit7 * adv_mom_3 & 2515 ) * & 2516 ( u(k_pp,j,i) - u(k-1,j,i) ) & 2517 + ( ibit8 * adv_mom_5 & 2518 ) * & 2519 ( u(k_ppp,j,i) - u(k_mm,j,i) ) & 2290 2520 ) 2291 2521 ENDDO … … 2296 2526 diss_d = diss_t(k-1) 2297 2527 2298 ibit 11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp )2299 ibit1 0 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp )2300 ibit 9 = REAL( IBITS(advc_flags_1(k,j,i),9,1),KIND = wp )2528 ibit2 = REAL( IBITS(advc_flags_m(k,j,i),2,1), KIND = wp ) 2529 ibit1 = REAL( IBITS(advc_flags_m(k,j,i),1,1), KIND = wp ) 2530 ibit0 = REAL( IBITS(advc_flags_m(k,j,i),0,1), KIND = wp ) 2301 2531 2302 ibit 14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp )2303 ibit 13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp )2304 ibit 12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp )2532 ibit5 = REAL( IBITS(advc_flags_m(k,j,i),5,1), KIND = wp ) 2533 ibit4 = REAL( IBITS(advc_flags_m(k,j,i),4,1), KIND = wp ) 2534 ibit3 = REAL( IBITS(advc_flags_m(k,j,i),3,1), KIND = wp ) 2305 2535 2306 ibit 17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp )2307 ibit 16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp )2308 ibit 15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp )2536 ibit8 = REAL( IBITS(advc_flags_m(k,j,i),8,1), KIND = wp ) 2537 ibit7 = REAL( IBITS(advc_flags_m(k,j,i),7,1), KIND = wp ) 2538 ibit6 = REAL( IBITS(advc_flags_m(k,j,i),6,1), KIND = wp ) 2309 2539 ! 2310 2540 !-- Calculate the divergence of the velocity field. A respective 2311 2541 !-- correction is needed to overcome numerical instabilities introduced 2312 2542 !-- by a not sufficient reduction of divergences near topography. 2313 div = ( ( u_comp(k) * ( ibit 9 + ibit10 + ibit11 )&2314 - ( u(k,j,i) + u(k,j,i-1) ) &2315 * ( &2316 REAL( IBITS(advc_flags_ 1(k,j,i-1),9,1), KIND = wp )&2317 + REAL( IBITS(advc_flags_ 1(k,j,i-1),10,1), KIND = wp )&2318 + REAL( IBITS(advc_flags_ 1(k,j,i-1),11,1), KIND = wp )&2319 ) &2320 ) * ddx &2321 + ( ( v_comp(k) + gv ) * ( ibit 12 + ibit13 + ibit14 )&2322 - ( v(k,j,i) + v(k,j,i-1 ) ) &2323 * ( &2324 REAL( IBITS(advc_flags_ 1(k,j-1,i),12,1), KIND = wp )&2325 + REAL( IBITS(advc_flags_ 1(k,j-1,i),13,1), KIND = wp )&2326 + REAL( IBITS(advc_flags_ 1(k,j-1,i),14,1), KIND = wp )&2327 ) &2328 ) * ddy &2329 + ( w_comp(k) * rho_air_zw(k) * ( ibit 15 + ibit16 + ibit17 )&2330 - w_comp(k-1) * rho_air_zw(k-1) &2331 * ( &2332 REAL( IBITS(advc_flags_ 1(k-1,j,i),15,1), KIND = wp )&2333 + REAL( IBITS(advc_flags_ 1(k-1,j,i),16,1), KIND = wp )&2334 + REAL( IBITS(advc_flags_ 1(k-1,j,i),17,1), KIND = wp )&2335 ) &2336 ) * drho_air(k) * ddzw(k) &2337 ) * 0.5_wp 2338 2339 tend(k,j,i) = tend(k,j,i) - ( &2340 ( flux_r(k) + diss_r(k) &2341 - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &2342 + ( flux_n(k) + diss_n(k) &2343 - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &2344 + ( ( flux_t(k) + diss_t(k) ) &2345 - ( flux_d + diss_d ) &2346 ) * drho_air(k) * ddzw(k) &2543 div = ( ( u_comp(k) * ( ibit0 + ibit1 + ibit2 ) & 2544 - ( u(k,j,i) + u(k,j,i-1) ) & 2545 * ( & 2546 REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp ) & 2547 + REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp ) & 2548 + REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp ) & 2549 ) & 2550 ) * ddx & 2551 + ( ( v_comp(k) + gv ) * ( ibit3 + ibit4 + ibit5 ) & 2552 - ( v(k,j,i) + v(k,j,i-1 ) ) & 2553 * ( & 2554 REAL( IBITS(advc_flags_m(k,j-1,i),3,1), KIND = wp ) & 2555 + REAL( IBITS(advc_flags_m(k,j-1,i),4,1), KIND = wp ) & 2556 + REAL( IBITS(advc_flags_m(k,j-1,i),5,1), KIND = wp ) & 2557 ) & 2558 ) * ddy & 2559 + ( w_comp(k) * rho_air_zw(k) * ( ibit6 + ibit7 + ibit8 ) & 2560 - w_comp(k-1) * rho_air_zw(k-1) & 2561 * ( & 2562 REAL( IBITS(advc_flags_m(k-1,j,i),6,1), KIND = wp ) & 2563 + REAL( IBITS(advc_flags_m(k-1,j,i),7,1), KIND = wp ) & 2564 + REAL( IBITS(advc_flags_m(k-1,j,i),8,1), KIND = wp ) & 2565 ) & 2566 ) * drho_air(k) * ddzw(k) & 2567 ) * 0.5_wp 2568 2569 tend(k,j,i) = tend(k,j,i) - ( & 2570 ( flux_r(k) + diss_r(k) & 2571 - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & 2572 + ( flux_n(k) + diss_n(k) & 2573 - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & 2574 + ( ( flux_t(k) + diss_t(k) ) & 2575 - ( flux_d + diss_d ) & 2576 ) * drho_air(k) * ddzw(k) & 2347 2577 ) + div * u(k,j,i) 2348 2578 … … 2450 2680 INTEGER(iwp) :: tn !< number of OpenMP thread 2451 2681 2452 REAL(wp) :: ibit 18!< flag indicating 1st-order scheme along x-direction2453 REAL(wp) :: ibit1 9!< flag indicating 3rd-order scheme along x-direction2454 REAL(wp) :: ibit 20!< flag indicating 5th-order scheme along x-direction2455 REAL(wp) :: ibit 21!< flag indicating 1st-order scheme along y-direction2456 REAL(wp) :: ibit 22!< flag indicating 3rd-order scheme along y-direction2457 REAL(wp) :: ibit 23!< flag indicating 3rd-order scheme along y-direction2458 REAL(wp) :: ibit 24!< flag indicating 1st-order scheme along z-direction2459 REAL(wp) :: ibit 25!< flag indicating 3rd-order scheme along z-direction2460 REAL(wp) :: ibit 26!< flag indicating 3rd-order scheme along z-direction2682 REAL(wp) :: ibit9 !< flag indicating 1st-order scheme along x-direction 2683 REAL(wp) :: ibit10 !< flag indicating 3rd-order scheme along x-direction 2684 REAL(wp) :: ibit11 !< flag indicating 5th-order scheme along x-direction 2685 REAL(wp) :: ibit12 !< flag indicating 1st-order scheme along y-direction 2686 REAL(wp) :: ibit13 !< flag indicating 3rd-order scheme along y-direction 2687 REAL(wp) :: ibit14 !< flag indicating 3rd-order scheme along y-direction 2688 REAL(wp) :: ibit15 !< flag indicating 1st-order scheme along z-direction 2689 REAL(wp) :: ibit16 !< flag indicating 3rd-order scheme along z-direction 2690 REAL(wp) :: ibit17 !< flag indicating 3rd-order scheme along z-direction 2461 2691 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 2462 2692 REAL(wp) :: div !< divergence on v-grid … … 2498 2728 DO k = nzb+1, nzb_max_l 2499 2729 2500 ibit 20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp )2501 ibit1 9 = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp )2502 ibit 18 = REAL( IBITS(advc_flags_1(k,j,i-1),18,1),KIND = wp )2730 ibit11 = REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp ) 2731 ibit10 = REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp ) 2732 ibit9 = REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp ) 2503 2733 2504 2734 u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu 2505 2735 flux_l_v(k,j,tn) = u_comp(k) * ( & 2506 ( 37.0_wp * ibit 20* adv_mom_5 &2507 + 7.0_wp * ibit1 9* adv_mom_3 &2508 + ibit 18* adv_mom_1 &2736 ( 37.0_wp * ibit11 * adv_mom_5 & 2737 + 7.0_wp * ibit10 * adv_mom_3 & 2738 + ibit9 * adv_mom_1 & 2509 2739 ) * & 2510 2740 ( v(k,j,i) + v(k,j,i-1) ) & 2511 - ( 8.0_wp * ibit 20* adv_mom_5 &2512 + ibit1 9* adv_mom_3 &2741 - ( 8.0_wp * ibit11 * adv_mom_5 & 2742 + ibit10 * adv_mom_3 & 2513 2743 ) * & 2514 2744 ( v(k,j,i+1) + v(k,j,i-2) ) & 2515 + ( ibit 20* adv_mom_5 &2745 + ( ibit11 * adv_mom_5 & 2516 2746 ) * & 2517 2747 ( v(k,j,i+2) + v(k,j,i-3) ) & … … 2519 2749 2520 2750 diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( & 2521 ( 10.0_wp * ibit 20* adv_mom_5 &2522 + 3.0_wp * ibit1 9* adv_mom_3 &2523 + ibit 18* adv_mom_1 &2751 ( 10.0_wp * ibit11 * adv_mom_5 & 2752 + 3.0_wp * ibit10 * adv_mom_3 & 2753 + ibit9 * adv_mom_1 & 2524 2754 ) * & 2525 2755 ( v(k,j,i) - v(k,j,i-1) ) & 2526 - ( 5.0_wp * ibit 20* adv_mom_5 &2527 + ibit1 9* adv_mom_3 &2756 - ( 5.0_wp * ibit11 * adv_mom_5 & 2757 + ibit10 * adv_mom_3 & 2528 2758 ) * & 2529 2759 ( v(k,j,i+1) - v(k,j,i-2) ) & 2530 + ( ibit 20* adv_mom_5 &2760 + ( ibit11 * adv_mom_5 & 2531 2761 ) * & 2532 2762 ( v(k,j,i+2) - v(k,j,i-3) ) & … … 2537 2767 DO k = nzb_max_l+1, nzt 2538 2768 2539 u_comp(k) 2769 u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu 2540 2770 flux_l_v(k,j,tn) = u_comp(k) * ( & 2541 37.0_wp * ( v(k,j,i) + v(k,j,i-1) )&2771 37.0_wp * ( v(k,j,i) + v(k,j,i-1) ) & 2542 2772 - 8.0_wp * ( v(k,j,i+1) + v(k,j,i-2) ) & 2543 2773 + ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5 2544 2774 diss_l_v(k,j,tn) = - ABS( u_comp(k) ) * ( & 2545 10.0_wp * ( v(k,j,i) - v(k,j,i-1) )&2775 10.0_wp * ( v(k,j,i) - v(k,j,i-1) ) & 2546 2776 - 5.0_wp * ( v(k,j,i+1) - v(k,j,i-2) ) & 2547 2777 + ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5 … … 2556 2786 DO k = nzb+1, nzb_max_l 2557 2787 2558 ibit 23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp )2559 ibit 22 = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp )2560 ibit 21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp )2788 ibit14 = REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp ) 2789 ibit13 = REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp ) 2790 ibit12 = REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp ) 2561 2791 2562 2792 v_comp_l = v(k,j,i) + v(k,j-1,i) - gv 2563 2793 flux_s_v(k,tn) = v_comp_l * ( & 2564 ( 37.0_wp * ibit 23* adv_mom_5 &2565 + 7.0_wp * ibit 22* adv_mom_3 &2566 + ibit 21* adv_mom_1 &2794 ( 37.0_wp * ibit14 * adv_mom_5 & 2795 + 7.0_wp * ibit13 * adv_mom_3 & 2796 + ibit12 * adv_mom_1 & 2567 2797 ) * & 2568 2798 ( v(k,j,i) + v(k,j-1,i) ) & 2569 - ( 8.0_wp * ibit 23* adv_mom_5 &2570 + ibit 22* adv_mom_3 &2799 - ( 8.0_wp * ibit14 * adv_mom_5 & 2800 + ibit13 * adv_mom_3 & 2571 2801 ) * & 2572 2802 ( v(k,j+1,i) + v(k,j-2,i) ) & 2573 + ( ibit 23* adv_mom_5 &2803 + ( ibit14 * adv_mom_5 & 2574 2804 ) * & 2575 2805 ( v(k,j+2,i) + v(k,j-3,i) ) & … … 2577 2807 2578 2808 diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & 2579 ( 10.0_wp * ibit 23* adv_mom_5 &2580 + 3.0_wp * ibit 22* adv_mom_3 &2581 + ibit 21* adv_mom_1 &2809 ( 10.0_wp * ibit14 * adv_mom_5 & 2810 + 3.0_wp * ibit13 * adv_mom_3 & 2811 + ibit12 * adv_mom_1 & 2582 2812 ) * & 2583 2813 ( v(k,j,i) - v(k,j-1,i) ) & 2584 - ( 5.0_wp * ibit 23* adv_mom_5 &2585 + ibit 22* adv_mom_3 &2814 - ( 5.0_wp * ibit14 * adv_mom_5 & 2815 + ibit13 * adv_mom_3 & 2586 2816 ) * & 2587 2817 ( v(k,j+1,i) - v(k,j-2,i) ) & 2588 + ( ibit 23* adv_mom_5 &2818 + ( ibit14 * adv_mom_5 & 2589 2819 ) * & 2590 2820 ( v(k,j+2,i) - v(k,j-3,i) ) & … … 2597 2827 v_comp_l = v(k,j,i) + v(k,j-1,i) - gv 2598 2828 flux_s_v(k,tn) = v_comp_l * ( & 2599 37.0_wp * ( v(k,j,i) + v(k,j-1,i) )&2829 37.0_wp * ( v(k,j,i) + v(k,j-1,i) ) & 2600 2830 - 8.0_wp * ( v(k,j+1,i) + v(k,j-2,i) ) & 2601 2831 + ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5 2602 2832 diss_s_v(k,tn) = - ABS( v_comp_l ) * ( & 2603 10.0_wp * ( v(k,j,i) - v(k,j-1,i) )&2833 10.0_wp * ( v(k,j,i) - v(k,j-1,i) ) & 2604 2834 - 5.0_wp * ( v(k,j+1,i) - v(k,j-2,i) ) & 2605 2835 + ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5 … … 2613 2843 DO k = nzb+1, nzb_max_l 2614 2844 2615 ibit 20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )2616 ibit1 9 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp )2617 ibit 18 = REAL( IBITS(advc_flags_1(k,j,i),18,1),KIND = wp )2845 ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp ) 2846 ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp ) 2847 ibit9 = REAL( IBITS(advc_flags_m(k,j,i),9,1), KIND = wp ) 2618 2848 2619 2849 u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu 2620 2850 flux_r(k) = u_comp(k) * ( & 2621 ( 37.0_wp * ibit 20* adv_mom_5 &2622 + 7.0_wp * ibit1 9* adv_mom_3 &2623 + ibit 18* adv_mom_1 &2851 ( 37.0_wp * ibit11 * adv_mom_5 & 2852 + 7.0_wp * ibit10 * adv_mom_3 & 2853 + ibit9 * adv_mom_1 & 2624 2854 ) * & 2625 2855 ( v(k,j,i+1) + v(k,j,i) ) & 2626 - ( 8.0_wp * ibit 20* adv_mom_5 &2627 + ibit1 9* adv_mom_3 &2856 - ( 8.0_wp * ibit11 * adv_mom_5 & 2857 + ibit10 * adv_mom_3 & 2628 2858 ) * & 2629 2859 ( v(k,j,i+2) + v(k,j,i-1) ) & 2630 + ( ibit 20* adv_mom_5 &2860 + ( ibit11 * adv_mom_5 & 2631 2861 ) * & 2632 2862 ( v(k,j,i+3) + v(k,j,i-2) ) & … … 2634 2864 2635 2865 diss_r(k) = - ABS( u_comp(k) ) * ( & 2636 ( 10.0_wp * ibit 20* adv_mom_5 &2637 + 3.0_wp * ibit1 9* adv_mom_3 &2638 + ibit 18* adv_mom_1 &2866 ( 10.0_wp * ibit11 * adv_mom_5 & 2867 + 3.0_wp * ibit10 * adv_mom_3 & 2868 + ibit9 * adv_mom_1 & 2639 2869 ) * & 2640 2870 ( v(k,j,i+1) - v(k,j,i) ) & 2641 - ( 5.0_wp * ibit 20* adv_mom_5 &2642 + ibit1 9* adv_mom_3 &2871 - ( 5.0_wp * ibit11 * adv_mom_5 & 2872 + ibit10 * adv_mom_3 & 2643 2873 ) * & 2644 2874 ( v(k,j,i+2) - v(k,j,i-1) ) & 2645 + ( ibit 20* adv_mom_5 &2875 + ( ibit11 * adv_mom_5 & 2646 2876 ) * & 2647 2877 ( v(k,j,i+3) - v(k,j,i-2) ) & 2648 2878 ) 2649 2879 2650 ibit 23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )2651 ibit 22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )2652 ibit 21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )2880 ibit14 = REAL( IBITS(advc_flags_m(k,j,i),14,1), KIND = wp ) 2881 ibit13 = REAL( IBITS(advc_flags_m(k,j,i),13,1), KIND = wp ) 2882 ibit12 = REAL( IBITS(advc_flags_m(k,j,i),12,1), KIND = wp ) 2653 2883 2654 2884 2655 2885 v_comp(k) = v(k,j+1,i) + v(k,j,i) 2656 2886 flux_n(k) = ( v_comp(k) - gv ) * ( & 2657 ( 37.0_wp * ibit 23* adv_mom_5 &2658 + 7.0_wp * ibit 22* adv_mom_3 &2659 + ibit 21* adv_mom_1 &2887 ( 37.0_wp * ibit14 * adv_mom_5 & 2888 + 7.0_wp * ibit13 * adv_mom_3 & 2889 + ibit12 * adv_mom_1 & 2660 2890 ) * & 2661 2891 ( v(k,j+1,i) + v(k,j,i) ) & 2662 - ( 8.0_wp * ibit 23* adv_mom_5 &2663 + ibit 22* adv_mom_3 &2892 - ( 8.0_wp * ibit14 * adv_mom_5 & 2893 + ibit13 * adv_mom_3 & 2664 2894 ) * & 2665 2895 ( v(k,j+2,i) + v(k,j-1,i) ) & 2666 + ( ibit 23* adv_mom_5 &2896 + ( ibit14 * adv_mom_5 & 2667 2897 ) * & 2668 2898 ( v(k,j+3,i) + v(k,j-2,i) ) & … … 2670 2900 2671 2901 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 2672 ( 10.0_wp * ibit 23* adv_mom_5 &2673 + 3.0_wp * ibit 22* adv_mom_3 &2674 + ibit 21* adv_mom_1 &2902 ( 10.0_wp * ibit14 * adv_mom_5 & 2903 + 3.0_wp * ibit13 * adv_mom_3 & 2904 + ibit12 * adv_mom_1 & 2675 2905 ) * & 2676 2906 ( v(k,j+1,i) - v(k,j,i) ) & 2677 - ( 5.0_wp * ibit 23* adv_mom_5 &2678 + ibit 22* adv_mom_3 &2907 - ( 5.0_wp * ibit14 * adv_mom_5 & 2908 + ibit13 * adv_mom_3 & 2679 2909 ) * & 2680 2910 ( v(k,j+2,i) - v(k,j-1,i) ) & 2681 + ( ibit 23* adv_mom_5 &2911 + ( ibit14 * adv_mom_5 & 2682 2912 ) * & 2683 2913 ( v(k,j+3,i) - v(k,j-2,i) ) & … … 2718 2948 !-- calculated explicetely for the tendency at 2719 2949 !-- the first w-level. For topography wall this is done implicitely by 2720 !-- advc_flags_ 1.2950 !-- advc_flags_m. 2721 2951 flux_t(nzb) = 0.0_wp 2722 2952 diss_t(nzb) = 0.0_wp … … 2727 2957 !-- k index has to be modified near bottom and top, else array 2728 2958 !-- subscripts will be exceeded. 2729 ibit 26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )2730 ibit 25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )2731 ibit 24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )2732 2733 k_ppp = k + 3 * ibit 262734 k_pp = k + 2 * ( 1 - ibit 24)2735 k_mm = k - 2 * ibit 262959 ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp ) 2960 ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp ) 2961 ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp ) 2962 2963 k_ppp = k + 3 * ibit17 2964 k_pp = k + 2 * ( 1 - ibit15 ) 2965 k_mm = k - 2 * ibit17 2736 2966 2737 2967 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2738 2968 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2739 ( 37.0_wp * ibit 26* adv_mom_5 &2740 + 7.0_wp * ibit 25* adv_mom_3 &2741 + ibit 24* adv_mom_1 &2969 ( 37.0_wp * ibit17 * adv_mom_5 & 2970 + 7.0_wp * ibit16 * adv_mom_3 & 2971 + ibit15 * adv_mom_1 & 2742 2972 ) * & 2743 2973 ( v(k+1,j,i) + v(k,j,i) ) & 2744 - ( 8.0_wp * ibit 26* adv_mom_5 &2745 + ibit 25* adv_mom_3 &2974 - ( 8.0_wp * ibit17 * adv_mom_5 & 2975 + ibit16 * adv_mom_3 & 2746 2976 ) * & 2747 2977 ( v(k_pp,j,i) + v(k-1,j,i) ) & 2748 + ( ibit 26* adv_mom_5 &2978 + ( ibit17 * adv_mom_5 & 2749 2979 ) * & 2750 2980 ( v(k_ppp,j,i) + v(k_mm,j,i) ) & … … 2752 2982 2753 2983 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2754 ( 10.0_wp * ibit 26* adv_mom_5 &2755 + 3.0_wp * ibit 25* adv_mom_3 &2756 + ibit 24* adv_mom_1 &2984 ( 10.0_wp * ibit17 * adv_mom_5 & 2985 + 3.0_wp * ibit16 * adv_mom_3 & 2986 + ibit15 * adv_mom_1 & 2757 2987 ) * & 2758 2988 ( v(k+1,j,i) - v(k,j,i) ) & 2759 - ( 5.0_wp * ibit 26* adv_mom_5 &2760 + ibit 25* adv_mom_3 &2989 - ( 5.0_wp * ibit17 * adv_mom_5 & 2990 + ibit16 * adv_mom_3 & 2761 2991 ) * & 2762 2992 ( v(k_pp,j,i) - v(k-1,j,i) ) & 2763 + ( ibit 26* adv_mom_5 &2993 + ( ibit17 * adv_mom_5 & 2764 2994 ) * & 2765 2995 ( v(k_ppp,j,i) - v(k_mm,j,i) ) & … … 2769 2999 DO k = nzb+3, nzt-2 2770 3000 2771 ibit 26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )2772 ibit 25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )2773 ibit 24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )3001 ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp ) 3002 ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp ) 3003 ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp ) 2774 3004 2775 3005 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2776 3006 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2777 ( 37.0_wp * ibit 26* adv_mom_5 &2778 + 7.0_wp * ibit 25* adv_mom_3 &2779 + ibit 24* adv_mom_1 &2780 ) * & 2781 ( v(k+1,j,i) + v(k,j,i) )&2782 - ( 8.0_wp * ibit 26* adv_mom_5 &2783 + ibit 25* adv_mom_3 &2784 ) * & 2785 ( v(k+2,j,i) + v(k-1,j,i) )&2786 + ( ibit 26* adv_mom_5 &3007 ( 37.0_wp * ibit17 * adv_mom_5 & 3008 + 7.0_wp * ibit16 * adv_mom_3 & 3009 + ibit15 * adv_mom_1 & 3010 ) * & 3011 ( v(k+1,j,i) + v(k,j,i) ) & 3012 - ( 8.0_wp * ibit17 * adv_mom_5 & 3013 + ibit16 * adv_mom_3 & 3014 ) * & 3015 ( v(k+2,j,i) + v(k-1,j,i) ) & 3016 + ( ibit17 * adv_mom_5 & 2787 3017 ) * & 2788 3018 ( v(k+3,j,i) + v(k-2,j,i) ) & … … 2790 3020 2791 3021 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2792 ( 10.0_wp * ibit 26* adv_mom_5 &2793 + 3.0_wp * ibit 25* adv_mom_3 &2794 + ibit 24* adv_mom_1 &2795 ) * & 2796 ( v(k+1,j,i) - v(k,j,i) )&2797 - ( 5.0_wp * ibit 26* adv_mom_5 &2798 + ibit 25* adv_mom_3 &2799 ) * & 2800 ( v(k+2,j,i) - v(k-1,j,i) )&2801 + ( ibit 26* adv_mom_5 &3022 ( 10.0_wp * ibit17 * adv_mom_5 & 3023 + 3.0_wp * ibit16 * adv_mom_3 & 3024 + ibit15 * adv_mom_1 & 3025 ) * & 3026 ( v(k+1,j,i) - v(k,j,i) ) & 3027 - ( 5.0_wp * ibit17 * adv_mom_5 & 3028 + ibit16 * adv_mom_3 & 3029 ) * & 3030 ( v(k+2,j,i) - v(k-1,j,i) ) & 3031 + ( ibit17 * adv_mom_5 & 2802 3032 ) * & 2803 3033 ( v(k+3,j,i) - v(k-2,j,i) ) & … … 2809 3039 !-- k index has to be modified near bottom and top, else array 2810 3040 !-- subscripts will be exceeded. 2811 ibit 26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )2812 ibit 25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )2813 ibit 24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )2814 2815 k_ppp = k + 3 * ibit 262816 k_pp = k + 2 * ( 1 - ibit 24)2817 k_mm = k - 2 * ibit 263041 ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp ) 3042 ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp ) 3043 ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp ) 3044 3045 k_ppp = k + 3 * ibit17 3046 k_pp = k + 2 * ( 1 - ibit15 ) 3047 k_mm = k - 2 * ibit17 2818 3048 2819 3049 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2820 3050 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2821 ( 37.0_wp * ibit 26* adv_mom_5 &2822 + 7.0_wp * ibit 25* adv_mom_3 &2823 + ibit 24* adv_mom_1 &3051 ( 37.0_wp * ibit17 * adv_mom_5 & 3052 + 7.0_wp * ibit16 * adv_mom_3 & 3053 + ibit15 * adv_mom_1 & 2824 3054 ) * & 2825 3055 ( v(k+1,j,i) + v(k,j,i) ) & 2826 - ( 8.0_wp * ibit 26* adv_mom_5 &2827 + ibit 25* adv_mom_3 &3056 - ( 8.0_wp * ibit17 * adv_mom_5 & 3057 + ibit16 * adv_mom_3 & 2828 3058 ) * & 2829 3059 ( v(k_pp,j,i) + v(k-1,j,i) ) & 2830 + ( ibit 26* adv_mom_5 &3060 + ( ibit17 * adv_mom_5 & 2831 3061 ) * & 2832 3062 ( v(k_ppp,j,i) + v(k_mm,j,i) ) & … … 2834 3064 2835 3065 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2836 ( 10.0_wp * ibit 26* adv_mom_5 &2837 + 3.0_wp * ibit 25* adv_mom_3 &2838 + ibit 24* adv_mom_1 &3066 ( 10.0_wp * ibit17 * adv_mom_5 & 3067 + 3.0_wp * ibit16 * adv_mom_3 & 3068 + ibit15 * adv_mom_1 & 2839 3069 ) * & 2840 3070 ( v(k+1,j,i) - v(k,j,i) ) & 2841 - ( 5.0_wp * ibit 26* adv_mom_5 &2842 + ibit 25* adv_mom_3 &3071 - ( 5.0_wp * ibit17 * adv_mom_5 & 3072 + ibit16 * adv_mom_3 & 2843 3073 ) * & 2844 3074 ( v(k_pp,j,i) - v(k-1,j,i) ) & 2845 + ( ibit 26* adv_mom_5 &3075 + ( ibit17 * adv_mom_5 & 2846 3076 ) * & 2847 3077 ( v(k_ppp,j,i) - v(k_mm,j,i) ) & … … 2854 3084 diss_d = diss_t(k-1) 2855 3085 2856 ibit 20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp )2857 ibit1 9 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp )2858 ibit 18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp )3086 ibit11 = REAL( IBITS(advc_flags_m(k,j,i),11,1), KIND = wp ) 3087 ibit10 = REAL( IBITS(advc_flags_m(k,j,i),10,1), KIND = wp ) 3088 ibit9 = REAL( IBITS(advc_flags_m(k,j,i),9,1), KIND = wp ) 2859 3089 2860 ibit 23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp )2861 ibit 22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp )2862 ibit 21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp )3090 ibit14 = REAL( IBITS(advc_flags_m(k,j,i),14,1), KIND = wp ) 3091 ibit13 = REAL( IBITS(advc_flags_m(k,j,i),13,1), KIND = wp ) 3092 ibit12 = REAL( IBITS(advc_flags_m(k,j,i),12,1), KIND = wp ) 2863 3093 2864 ibit 26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp )2865 ibit 25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp )2866 ibit 24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp )3094 ibit17 = REAL( IBITS(advc_flags_m(k,j,i),17,1), KIND = wp ) 3095 ibit16 = REAL( IBITS(advc_flags_m(k,j,i),16,1), KIND = wp ) 3096 ibit15 = REAL( IBITS(advc_flags_m(k,j,i),15,1), KIND = wp ) 2867 3097 ! 2868 3098 !-- Calculate the divergence of the velocity field. A respective … … 2870 3100 !-- by a not sufficient reduction of divergences near topography. 2871 3101 div = ( ( ( u_comp(k) + gu ) & 2872 * ( ibit 18 + ibit19 + ibit20 )&3102 * ( ibit9 + ibit10 + ibit11 ) & 2873 3103 - ( u(k,j-1,i) + u(k,j,i) ) & 2874 3104 * ( & 2875 REAL( IBITS(advc_flags_ 1(k,j,i-1),18,1),KIND = wp ) &2876 + REAL( IBITS(advc_flags_ 1(k,j,i-1),19,1), KIND = wp ) &2877 + REAL( IBITS(advc_flags_ 1(k,j,i-1),20,1), KIND = wp ) &3105 REAL( IBITS(advc_flags_m(k,j,i-1),9,1), KIND = wp ) & 3106 + REAL( IBITS(advc_flags_m(k,j,i-1),10,1), KIND = wp ) & 3107 + REAL( IBITS(advc_flags_m(k,j,i-1),11,1), KIND = wp ) & 2878 3108 ) & 2879 3109 ) * ddx & 2880 3110 + ( v_comp(k) & 2881 * ( ibit 21 + ibit22 + ibit23) &3111 * ( ibit12 + ibit13 + ibit14 ) & 2882 3112 - ( v(k,j,i) + v(k,j-1,i) ) & 2883 3113 * ( & 2884 REAL( IBITS(advc_flags_ 1(k,j-1,i),21,1), KIND = wp ) &2885 + REAL( IBITS(advc_flags_ 1(k,j-1,i),22,1), KIND = wp ) &2886 + REAL( IBITS(advc_flags_ 1(k,j-1,i),23,1), KIND = wp ) &3114 REAL( IBITS(advc_flags_m(k,j-1,i),12,1), KIND = wp ) & 3115 + REAL( IBITS(advc_flags_m(k,j-1,i),13,1), KIND = wp ) & 3116 + REAL( IBITS(advc_flags_m(k,j-1,i),14,1), KIND = wp ) & 2887 3117 ) & 2888 3118 ) * ddy & 2889 + ( w_comp(k) * rho_air_zw(k) * ( ibit 24 + ibit25 + ibit26)&3119 + ( w_comp(k) * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )& 2890 3120 - w_comp(k-1) * rho_air_zw(k-1) & 2891 3121 * ( & 2892 REAL( IBITS(advc_flags_ 1(k-1,j,i),24,1), KIND = wp ) &2893 + REAL( IBITS(advc_flags_ 1(k-1,j,i),25,1), KIND = wp ) &2894 + REAL( IBITS(advc_flags_ 1(k-1,j,i),26,1), KIND = wp ) &3122 REAL( IBITS(advc_flags_m(k-1,j,i),15,1), KIND = wp ) & 3123 + REAL( IBITS(advc_flags_m(k-1,j,i),16,1), KIND = wp ) & 3124 + REAL( IBITS(advc_flags_m(k-1,j,i),17,1), KIND = wp ) & 2895 3125 ) & 2896 3126 ) * drho_air(k) * ddzw(k) & … … 3011 3241 INTEGER(iwp) :: tn !< number of OpenMP thread 3012 3242 3013 REAL(wp) :: ibit 27!< flag indicating 1st-order scheme along x-direction3014 REAL(wp) :: ibit 28!< flag indicating 3rd-order scheme along x-direction3015 REAL(wp) :: ibit2 9!< flag indicating 5th-order scheme along x-direction3016 REAL(wp) :: ibit 30!< flag indicating 1st-order scheme along y-direction3017 REAL(wp) :: ibit 31!< flag indicating 3rd-order scheme along y-direction3018 REAL(wp) :: ibit 32!< flag indicating 5th-order scheme along y-direction3019 REAL(wp) :: ibit 33!< flag indicating 1st-order scheme along z-direction3020 REAL(wp) :: ibit 34!< flag indicating 3rd-order scheme along z-direction3021 REAL(wp) :: ibit 35!< flag indicating 5th-order scheme along z-direction3243 REAL(wp) :: ibit18 !< flag indicating 1st-order scheme along x-direction 3244 REAL(wp) :: ibit19 !< flag indicating 3rd-order scheme along x-direction 3245 REAL(wp) :: ibit20 !< flag indicating 5th-order scheme along x-direction 3246 REAL(wp) :: ibit21 !< flag indicating 1st-order scheme along y-direction 3247 REAL(wp) :: ibit22 !< flag indicating 3rd-order scheme along y-direction 3248 REAL(wp) :: ibit23 !< flag indicating 5th-order scheme along y-direction 3249 REAL(wp) :: ibit24 !< flag indicating 1st-order scheme along z-direction 3250 REAL(wp) :: ibit25 !< flag indicating 3rd-order scheme along z-direction 3251 REAL(wp) :: ibit26 !< flag indicating 5th-order scheme along z-direction 3022 3252 REAL(wp) :: diss_d !< discretized artificial dissipation at top of the grid box 3023 3253 REAL(wp) :: div !< divergence on w-grid … … 3056 3286 3057 3287 DO k = nzb+1, nzb_max_l 3058 ibit 32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1),KIND = wp )3059 ibit 31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp )3060 ibit 30 = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp )3288 ibit23 = REAL( IBITS(advc_flags_m(k,j-1,i),23,1), KIND = wp ) 3289 ibit22 = REAL( IBITS(advc_flags_m(k,j-1,i),22,1), KIND = wp ) 3290 ibit21 = REAL( IBITS(advc_flags_m(k,j-1,i),21,1), KIND = wp ) 3061 3291 3062 3292 v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv 3063 3293 flux_s_w(k,tn) = v_comp(k) * ( & 3064 ( 37.0_wp * ibit 32* adv_mom_5 &3065 + 7.0_wp * ibit 31* adv_mom_3 &3066 + ibit 30* adv_mom_1 &3294 ( 37.0_wp * ibit23 * adv_mom_5 & 3295 + 7.0_wp * ibit22 * adv_mom_3 & 3296 + ibit21 * adv_mom_1 & 3067 3297 ) * & 3068 3298 ( w(k,j,i) + w(k,j-1,i) ) & 3069 - ( 8.0_wp * ibit 32* adv_mom_5 &3070 + ibit 31* adv_mom_3 &3299 - ( 8.0_wp * ibit23 * adv_mom_5 & 3300 + ibit22 * adv_mom_3 & 3071 3301 ) * & 3072 3302 ( w(k,j+1,i) + w(k,j-2,i) ) & 3073 + ( ibit 32* adv_mom_5 &3303 + ( ibit23 * adv_mom_5 & 3074 3304 ) * & 3075 3305 ( w(k,j+2,i) + w(k,j-3,i) ) & … … 3077 3307 3078 3308 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( & 3079 ( 10.0_wp * ibit 32* adv_mom_5 &3080 + 3.0_wp * ibit 31* adv_mom_3 &3081 + ibit 30* adv_mom_1 &3309 ( 10.0_wp * ibit23 * adv_mom_5 & 3310 + 3.0_wp * ibit22 * adv_mom_3 & 3311 + ibit21 * adv_mom_1 & 3082 3312 ) * & 3083 3313 ( w(k,j,i) - w(k,j-1,i) ) & 3084 - ( 5.0_wp * ibit 32* adv_mom_5 &3085 + ibit 31* adv_mom_3 &3314 - ( 5.0_wp * ibit23 * adv_mom_5 & 3315 + ibit22 * adv_mom_3 & 3086 3316 ) * & 3087 3317 ( w(k,j+1,i) - w(k,j-2,i) ) & 3088 + ( ibit 32* adv_mom_5 &3318 + ( ibit23 * adv_mom_5 & 3089 3319 ) * & 3090 3320 ( w(k,j+2,i) - w(k,j-3,i) ) & … … 3097 3327 v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv 3098 3328 flux_s_w(k,tn) = v_comp(k) * ( & 3099 37.0_wp * ( w(k,j,i) + w(k,j-1,i)) &3100 - 8.0_wp * ( w(k,j+1,i) + w(k,j-2,i)) &3329 37.0_wp * ( w(k,j,i) + w(k,j-1,i) ) & 3330 - 8.0_wp * ( w(k,j+1,i) + w(k,j-2,i) ) & 3101 3331 + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 3102 3332 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( & 3103 10.0_wp * ( w(k,j,i) - w(k,j-1,i)) &3333 10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) & 3104 3334 - 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) & 3105 3335 + ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5 … … 3114 3344 DO k = nzb+1, nzb_max_l 3115 3345 3116 ibit2 9 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp )3117 ibit 28 = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp )3118 ibit 27 = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp )3346 ibit20 = REAL( IBITS(advc_flags_m(k,j,i-1),20,1), KIND = wp ) 3347 ibit19 = REAL( IBITS(advc_flags_m(k,j,i-1),19,1), KIND = wp ) 3348 ibit18 = REAL( IBITS(advc_flags_m(k,j,i-1),18,1), KIND = wp ) 3119 3349 3120 3350 u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu 3121 3351 flux_l_w(k,j,tn) = u_comp(k) * ( & 3122 ( 37.0_wp * ibit2 9* adv_mom_5 &3123 + 7.0_wp * ibit 28* adv_mom_3 &3124 + ibit 27* adv_mom_1 &3352 ( 37.0_wp * ibit20 * adv_mom_5 & 3353 + 7.0_wp * ibit19 * adv_mom_3 & 3354 + ibit18 * adv_mom_1 & 3125 3355 ) * & 3126 3356 ( w(k,j,i) + w(k,j,i-1) ) & 3127 - ( 8.0_wp * ibit2 9* adv_mom_5 &3128 + ibit 28* adv_mom_3 &3357 - ( 8.0_wp * ibit20 * adv_mom_5 & 3358 + ibit19 * adv_mom_3 & 3129 3359 ) * & 3130 3360 ( w(k,j,i+1) + w(k,j,i-2) ) & 3131 + ( ibit2 9* adv_mom_5 &3361 + ( ibit20 * adv_mom_5 & 3132 3362 ) * & 3133 3363 ( w(k,j,i+2) + w(k,j,i-3) ) & … … 3135 3365 3136 3366 diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( & 3137 ( 10.0_wp * ibit2 9* adv_mom_5 &3138 + 3.0_wp * ibit 28* adv_mom_3 &3139 + ibit 27* adv_mom_1 &3367 ( 10.0_wp * ibit20 * adv_mom_5 & 3368 + 3.0_wp * ibit19 * adv_mom_3 & 3369 + ibit18 * adv_mom_1 & 3140 3370 ) * & 3141 3371 ( w(k,j,i) - w(k,j,i-1) ) & 3142 - ( 5.0_wp * ibit2 9* adv_mom_5 &3143 + ibit 28* adv_mom_3 &3372 - ( 5.0_wp * ibit20 * adv_mom_5 & 3373 + ibit19 * adv_mom_3 & 3144 3374 ) * & 3145 3375 ( w(k,j,i+1) - w(k,j,i-2) ) & 3146 + ( ibit2 9* adv_mom_5 &3376 + ( ibit20 * adv_mom_5 & 3147 3377 ) * & 3148 3378 ( w(k,j,i+2) - w(k,j,i-3) ) & … … 3155 3385 u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu 3156 3386 flux_l_w(k,j,tn) = u_comp(k) * ( & 3157 37.0_wp * ( w(k,j,i) + w(k,j,i-1)) &3387 37.0_wp * ( w(k,j,i) + w(k,j,i-1) ) & 3158 3388 - 8.0_wp * ( w(k,j,i+1) + w(k,j,i-2) ) & 3159 3389 + ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5 3160 3390 diss_l_w(k,j,tn) = - ABS( u_comp(k) ) * ( & 3161 10.0_wp * ( w(k,j,i) - w(k,j,i-1)) &3391 10.0_wp * ( w(k,j,i) - w(k,j,i-1) ) & 3162 3392 - 5.0_wp * ( w(k,j,i+1) - w(k,j,i-2) ) & 3163 3393 + ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5 … … 3171 3401 DO k = nzb+1, nzb_max_l 3172 3402 3173 ibit2 9 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )3174 ibit 28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp )3175 ibit 27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )3403 ibit20 = REAL( IBITS(advc_flags_m(k,j,i),20,1), KIND = wp ) 3404 ibit19 = REAL( IBITS(advc_flags_m(k,j,i),19,1), KIND = wp ) 3405 ibit18 = REAL( IBITS(advc_flags_m(k,j,i),18,1), KIND = wp ) 3176 3406 3177 3407 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu 3178 3408 flux_r(k) = u_comp(k) * ( & 3179 ( 37.0_wp * ibit2 9* adv_mom_5 &3180 + 7.0_wp * ibit 28* adv_mom_3 &3181 + ibit 27* adv_mom_1 &3409 ( 37.0_wp * ibit20 * adv_mom_5 & 3410 + 7.0_wp * ibit19 * adv_mom_3 & 3411 + ibit18 * adv_mom_1 & 3182 3412 ) * & 3183 3413 ( w(k,j,i+1) + w(k,j,i) ) & 3184 - ( 8.0_wp * ibit2 9* adv_mom_5 &3185 + ibit 28* adv_mom_3 &3414 - ( 8.0_wp * ibit20 * adv_mom_5 & 3415 + ibit19 * adv_mom_3 & 3186 3416 ) * & 3187 3417 ( w(k,j,i+2) + w(k,j,i-1) ) & 3188 + ( ibit2 9* adv_mom_5 &3418 + ( ibit20 * adv_mom_5 & 3189 3419 ) * & 3190 3420 ( w(k,j,i+3) + w(k,j,i-2) ) & … … 3192 3422 3193 3423 diss_r(k) = - ABS( u_comp(k) ) * ( & 3194 ( 10.0_wp * ibit2 9* adv_mom_5 &3195 + 3.0_wp * ibit 28* adv_mom_3 &3196 + ibit 27* adv_mom_1 &3424 ( 10.0_wp * ibit20 * adv_mom_5 & 3425 + 3.0_wp * ibit19 * adv_mom_3 & 3426 + ibit18 * adv_mom_1 & 3197 3427 ) * & 3198 3428 ( w(k,j,i+1) - w(k,j,i) ) & 3199 - ( 5.0_wp * ibit2 9* adv_mom_5 &3200 + ibit 28* adv_mom_3 &3429 - ( 5.0_wp * ibit20 * adv_mom_5 & 3430 + ibit19 * adv_mom_3 & 3201 3431 ) * & 3202 3432 ( w(k,j,i+2) - w(k,j,i-1) ) & 3203 + ( ibit2 9* adv_mom_5 &3433 + ( ibit20 * adv_mom_5 & 3204 3434 ) * & 3205 3435 ( w(k,j,i+3) - w(k,j,i-2) ) & 3206 3436 ) 3207 3437 3208 ibit 32 = REAL( IBITS(advc_flags_2(k,j,i),0,1),KIND = wp )3209 ibit 31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )3210 ibit 30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )3438 ibit23 = REAL( IBITS(advc_flags_m(k,j,i),23,1), KIND = wp ) 3439 ibit22 = REAL( IBITS(advc_flags_m(k,j,i),22,1), KIND = wp ) 3440 ibit21 = REAL( IBITS(advc_flags_m(k,j,i),21,1), KIND = wp ) 3211 3441 3212 3442 v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv 3213 3443 flux_n(k) = v_comp(k) * ( & 3214 ( 37.0_wp * ibit 32* adv_mom_5 &3215 + 7.0_wp * ibit 31* adv_mom_3 &3216 + ibit 30* adv_mom_1 &3444 ( 37.0_wp * ibit23 * adv_mom_5 & 3445 + 7.0_wp * ibit22 * adv_mom_3 & 3446 + ibit21 * adv_mom_1 & 3217 3447 ) * & 3218 3448 ( w(k,j+1,i) + w(k,j,i) ) & 3219 - ( 8.0_wp * ibit 32* adv_mom_5 &3220 + ibit 31* adv_mom_3 &3449 - ( 8.0_wp * ibit23 * adv_mom_5 & 3450 + ibit22 * adv_mom_3 & 3221 3451 ) * & 3222 3452 ( w(k,j+2,i) + w(k,j-1,i) ) & 3223 + ( ibit 32* adv_mom_5 &3453 + ( ibit23 * adv_mom_5 & 3224 3454 ) * & 3225 3455 ( w(k,j+3,i) + w(k,j-2,i) ) & … … 3227 3457 3228 3458 diss_n(k) = - ABS( v_comp(k) ) * ( & 3229 ( 10.0_wp * ibit 32* adv_mom_5 &3230 + 3.0_wp * ibit 31* adv_mom_3 &3231 + ibit 30* adv_mom_1 &3232 ) * & 3233 ( w(k,j+1,i) - w(k,j,i) )&3234 - ( 5.0_wp * ibit 32* adv_mom_5 &3235 + ibit 31* adv_mom_3 &3236 ) * & 3237 ( w(k,j+2,i) - w(k,j-1,i) )&3238 + ( ibit 32* adv_mom_5 &3239 ) * & 3240 ( w(k,j+3,i) - w(k,j-2,i) )&3459 ( 10.0_wp * ibit23 * adv_mom_5 & 3460 + 3.0_wp * ibit22 * adv_mom_3 & 3461 + ibit21 * adv_mom_1 & 3462 ) * & 3463 ( w(k,j+1,i) - w(k,j,i) ) & 3464 - ( 5.0_wp * ibit23 * adv_mom_5 & 3465 + ibit22 * adv_mom_3 & 3466 ) * & 3467 ( w(k,j+2,i) - w(k,j-1,i) ) & 3468 + ( ibit23 * adv_mom_5 & 3469 ) * & 3470 ( w(k,j+3,i) - w(k,j-2,i) ) & 3241 3471 ) 3242 3472 ENDDO … … 3275 3505 !-- calculated explicetely for the tendency at 3276 3506 !-- the first w-level. For topography wall this is done implicitely by 3277 !-- advc_flags_ 1.3507 !-- advc_flags_m. 3278 3508 k = nzb + 1 3279 3509 w_comp(k) = w(k,j,i) + w(k-1,j,i) … … 3285 3515 !-- k index has to be modified near bottom and top, else array 3286 3516 !-- subscripts will be exceeded. 3287 ibit 35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )3288 ibit 34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )3289 ibit 33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )3290 3291 k_ppp = k + 3 * ibit 353292 k_pp = k + 2 * ( 1 - ibit 33)3293 k_mm = k - 2 * ibit 353517 ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp ) 3518 ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp ) 3519 ibit24 = REAL( IBITS(advc_flags_m(k,j,i),24,1), KIND = wp ) 3520 3521 k_ppp = k + 3 * ibit26 3522 k_pp = k + 2 * ( 1 - ibit24 ) 3523 k_mm = k - 2 * ibit26 3294 3524 3295 3525 w_comp(k) = w(k+1,j,i) + w(k,j,i) 3296 3526 flux_t(k) = w_comp(k) * rho_air(k+1) * ( & 3297 ( 37.0_wp * ibit 35* adv_mom_5 &3298 + 7.0_wp * ibit 34* adv_mom_3 &3299 + ibit 33* adv_mom_1 &3300 ) * & 3301 ( w(k+1,j,i) + w(k,j,i)) &3302 - ( 8.0_wp * ibit 35* adv_mom_5 &3303 + ibit 34* adv_mom_3 &3527 ( 37.0_wp * ibit26 * adv_mom_5 & 3528 + 7.0_wp * ibit25 * adv_mom_3 & 3529 + ibit24 * adv_mom_1 & 3530 ) * & 3531 ( w(k+1,j,i) + w(k,j,i) ) & 3532 - ( 8.0_wp * ibit26 * adv_mom_5 & 3533 + ibit25 * adv_mom_3 & 3304 3534 ) * & 3305 3535 ( w(k_pp,j,i) + w(k-1,j,i) ) & 3306 + ( ibit 35* adv_mom_5 &3536 + ( ibit26 * adv_mom_5 & 3307 3537 ) * & 3308 3538 ( w(k_ppp,j,i) + w(k_mm,j,i) ) & … … 3310 3540 3311 3541 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( & 3312 ( 10.0_wp * ibit 35* adv_mom_5 &3313 + 3.0_wp * ibit 34* adv_mom_3 &3314 + ibit 33* adv_mom_1 &3542 ( 10.0_wp * ibit26 * adv_mom_5 & 3543 + 3.0_wp * ibit25 * adv_mom_3 & 3544 + ibit24 * adv_mom_1 & 3315 3545 ) * & 3316 3546 ( w(k+1,j,i) - w(k,j,i) ) & 3317 - ( 5.0_wp * ibit 35* adv_mom_5 &3318 + ibit 34* adv_mom_3 &3547 - ( 5.0_wp * ibit26 * adv_mom_5 & 3548 + ibit25 * adv_mom_3 & 3319 3549 ) * & 3320 3550 ( w(k_pp,j,i) - w(k-1,j,i) ) & 3321 + ( ibit 35* adv_mom_5 &3551 + ( ibit26 * adv_mom_5 & 3322 3552 ) * & 3323 3553 ( w(k_ppp,j,i) - w(k_mm,j,i) ) & … … 3327 3557 DO k = nzb+3, nzt-2 3328 3558 3329 ibit 35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )3330 ibit 34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )3331 ibit 33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )3559 ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp ) 3560 ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp ) 3561 ibit24 = REAL( IBITS(advc_flags_m(k,j,i),24,1), KIND = wp ) 3332 3562 3333 3563 w_comp(k) = w(k+1,j,i) + w(k,j,i) 3334 3564 flux_t(k) = w_comp(k) * rho_air(k+1) * ( & 3335 ( 37.0_wp * ibit 35* adv_mom_5 &3336 + 7.0_wp * ibit 34* adv_mom_3 &3337 + ibit 33* adv_mom_1 &3338 ) * & 3339 ( w(k+1,j,i) + w(k,j,i) )&3340 - ( 8.0_wp * ibit 35* adv_mom_5 &3341 + ibit 34* adv_mom_3 &3342 ) * & 3343 ( w(k+2,j,i) + w(k-1,j,i) )&3344 + ( ibit 35* adv_mom_5 &3345 ) * & 3346 ( w(k+3,j,i) + w(k-2,j,i) )&3565 ( 37.0_wp * ibit26 * adv_mom_5 & 3566 + 7.0_wp * ibit25 * adv_mom_3 & 3567 + ibit24 * adv_mom_1 & 3568 ) * & 3569 ( w(k+1,j,i) + w(k,j,i) ) & 3570 - ( 8.0_wp * ibit26 * adv_mom_5 & 3571 + ibit25 * adv_mom_3 & 3572 ) * & 3573 ( w(k+2,j,i) + w(k-1,j,i) ) & 3574 + ( ibit26 * adv_mom_5 & 3575 ) * & 3576 ( w(k+3,j,i) + w(k-2,j,i) ) & 3347 3577 ) 3348 3578 3349 3579 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( & 3350 ( 10.0_wp * ibit 35* adv_mom_5 &3351 + 3.0_wp * ibit 34* adv_mom_3 &3352 + ibit 33* adv_mom_1 &3353 ) * & 3354 ( w(k+1,j,i) - w(k,j,i) )&3355 - ( 5.0_wp * ibit 35* adv_mom_5 &3356 + ibit 34* adv_mom_3 &3357 ) * & 3358 ( w(k+2,j,i) - w(k-1,j,i) )&3359 + ( ibit 35* adv_mom_5 &3360 ) * & 3361 ( w(k+3,j,i) - w(k-2,j,i) )&3580 ( 10.0_wp * ibit26 * adv_mom_5 & 3581 + 3.0_wp * ibit25 * adv_mom_3 & 3582 + ibit24 * adv_mom_1 & 3583 ) * & 3584 ( w(k+1,j,i) - w(k,j,i) ) & 3585 - ( 5.0_wp * ibit26 * adv_mom_5 & 3586 + ibit25 * adv_mom_3 & 3587 ) * & 3588 ( w(k+2,j,i) - w(k-1,j,i) ) & 3589 + ( ibit26 * adv_mom_5 & 3590 ) * & 3591 ( w(k+3,j,i) - w(k-2,j,i) ) & 3362 3592 ) 3363 3593 ENDDO … … 3367 3597 !-- k index has to be modified near bottom and top, else array 3368 3598 !-- subscripts will be exceeded. 3369 ibit 35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )3370 ibit 34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )3371 ibit 33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )3372 3373 k_ppp = k + 3 * ibit 353374 k_pp = k + 2 * ( 1 - ibit 33)3375 k_mm = k - 2 * ibit 353599 ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp ) 3600 ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp ) 3601 ibit24 = REAL( IBITS(advc_flags_m(k,j,i),24,1), KIND = wp ) 3602 3603 k_ppp = k + 3 * ibit26 3604 k_pp = k + 2 * ( 1 - ibit24 ) 3605 k_mm = k - 2 * ibit26 3376 3606 3377 3607 w_comp(k) = w(k+1,j,i) + w(k,j,i) 3378 3608 flux_t(k) = w_comp(k) * rho_air(k+1) * ( & 3379 ( 37.0_wp * ibit 35* adv_mom_5 &3380 + 7.0_wp * ibit 34* adv_mom_3 &3381 + ibit 33* adv_mom_1 &3382 ) * & 3383 ( w(k+1,j,i) + w(k,j,i)) &3384 - ( 8.0_wp * ibit 35* adv_mom_5 &3385 + ibit 34* adv_mom_3 &3609 ( 37.0_wp * ibit26 * adv_mom_5 & 3610 + 7.0_wp * ibit25 * adv_mom_3 & 3611 + ibit24 * adv_mom_1 & 3612 ) * & 3613 ( w(k+1,j,i) + w(k,j,i) ) & 3614 - ( 8.0_wp * ibit26 * adv_mom_5 & 3615 + ibit25 * adv_mom_3 & 3386 3616 ) * & 3387 3617 ( w(k_pp,j,i) + w(k-1,j,i) ) & 3388 + ( ibit 35* adv_mom_5 &3618 + ( ibit26 * adv_mom_5 & 3389 3619 ) * & 3390 3620 ( w(k_ppp,j,i) + w(k_mm,j,i) ) & … … 3392 3622 3393 3623 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( & 3394 ( 10.0_wp * ibit 35* adv_mom_5 &3395 + 3.0_wp * ibit 34* adv_mom_3 &3396 + ibit 33* adv_mom_1 &3624 ( 10.0_wp * ibit26 * adv_mom_5 & 3625 + 3.0_wp * ibit25 * adv_mom_3 & 3626 + ibit24 * adv_mom_1 & 3397 3627 ) * & 3398 3628 ( w(k+1,j,i) - w(k,j,i) ) & 3399 - ( 5.0_wp * ibit 35* adv_mom_5 &3400 + ibit 34* adv_mom_3 &3629 - ( 5.0_wp * ibit26 * adv_mom_5 & 3630 + ibit25 * adv_mom_3 & 3401 3631 ) * & 3402 3632 ( w(k_pp,j,i) - w(k-1,j,i) ) & 3403 + ( ibit 35* adv_mom_5 &3633 + ( ibit26 * adv_mom_5 & 3404 3634 ) * & 3405 3635 ( w(k_ppp,j,i) - w(k_mm,j,i) ) & … … 3412 3642 diss_d = diss_t(k-1) 3413 3643 3414 ibit2 9 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp )3415 ibit 28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp )3416 ibit 27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp )3644 ibit20 = REAL( IBITS(advc_flags_m(k,j,i),20,1), KIND = wp ) 3645 ibit19 = REAL( IBITS(advc_flags_m(k,j,i),19,1), KIND = wp ) 3646 ibit18 = REAL( IBITS(advc_flags_m(k,j,i),18,1), KIND = wp ) 3417 3647 3418 ibit 32 = REAL( IBITS(advc_flags_2(k,j,i),0,1),KIND = wp )3419 ibit 31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp )3420 ibit 30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp )3648 ibit23 = REAL( IBITS(advc_flags_m(k,j,i),23,1), KIND = wp ) 3649 ibit22 = REAL( IBITS(advc_flags_m(k,j,i),22,1), KIND = wp ) 3650 ibit21 = REAL( IBITS(advc_flags_m(k,j,i),21,1), KIND = wp ) 3421 3651 3422 ibit 35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp )3423 ibit 34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp )3424 ibit 33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp )3652 ibit26 = REAL( IBITS(advc_flags_m(k,j,i),26,1), KIND = wp ) 3653 ibit25 = REAL( IBITS(advc_flags_m(k,j,i),25,1), KIND = wp ) 3654 ibit24 = REAL( IBITS(advc_flags_m(k,j,i),24,1), KIND = wp ) 3425 3655 ! 3426 3656 !-- Calculate the divergence of the velocity field. A respective 3427 3657 !-- correction is needed to overcome numerical instabilities introduced 3428 3658 !-- by a not sufficient reduction of divergences near topography. 3429 div = ( ( ( u_comp(k) + gu ) * ( ibit 27 + ibit28 + ibit29) &3659 div = ( ( ( u_comp(k) + gu ) * ( ibit18 + ibit19 + ibit20 ) & 3430 3660 - ( u(k+1,j,i) + u(k,j,i) ) & 3431 3661 * ( & 3432 REAL( IBITS(advc_flags_ 1(k,j,i-1),27,1), KIND = wp ) &3433 + REAL( IBITS(advc_flags_ 1(k,j,i-1),28,1), KIND = wp ) &3434 + REAL( IBITS(advc_flags_ 1(k,j,i-1),29,1), KIND = wp ) &3662 REAL( IBITS(advc_flags_m(k,j,i-1),18,1), KIND = wp ) & 3663 + REAL( IBITS(advc_flags_m(k,j,i-1),19,1), KIND = wp ) & 3664 + REAL( IBITS(advc_flags_m(k,j,i-1),20,1), KIND = wp ) & 3435 3665 ) & 3436 3666 ) * ddx & 3437 + ( ( v_comp(k) + gv ) * ( ibit 30 + ibit31 + ibit32) &3667 + ( ( v_comp(k) + gv ) * ( ibit21 + ibit22 + ibit23 ) & 3438 3668 - ( v(k+1,j,i) + v(k,j,i) ) & 3439 3669 * ( & 3440 REAL( IBITS(advc_flags_ 1(k,j-1,i),30,1), KIND = wp ) &3441 + REAL( IBITS(advc_flags_ 1(k,j-1,i),31,1), KIND = wp ) &3442 + REAL( IBITS(advc_flags_ 2(k,j-1,i),0,1),KIND = wp ) &3670 REAL( IBITS(advc_flags_m(k,j-1,i),21,1), KIND = wp ) & 3671 + REAL( IBITS(advc_flags_m(k,j-1,i),22,1), KIND = wp ) & 3672 + REAL( IBITS(advc_flags_m(k,j-1,i),23,1), KIND = wp ) & 3443 3673 ) & 3444 3674 ) * ddy & 3445 3675 + ( w_comp(k) * rho_air(k+1) & 3446 * ( ibit 33 + ibit34 + ibit35) &3676 * ( ibit24 + ibit25 + ibit26 ) & 3447 3677 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 3448 3678 * ( & 3449 REAL( IBITS(advc_flags_ 2(k-1,j,i),1,1), KIND = wp )&3450 + REAL( IBITS(advc_flags_ 2(k-1,j,i),2,1), KIND = wp )&3451 + REAL( IBITS(advc_flags_ 2(k-1,j,i),3,1), KIND = wp )&3679 REAL( IBITS(advc_flags_m(k-1,j,i),24,1), KIND = wp ) & 3680 + REAL( IBITS(advc_flags_m(k-1,j,i),25,1), KIND = wp ) & 3681 + REAL( IBITS(advc_flags_m(k-1,j,i),26,1), KIND = wp ) & 3452 3682 ) & 3453 3683 ) * drho_air_zw(k) * ddzu(k+1) & … … 3531 3761 !> Scalar advection - Call for all grid points 3532 3762 !------------------------------------------------------------------------------! 3533 SUBROUTINE advec_s_ws( sk, sk_char ) 3763 SUBROUTINE advec_s_ws( advc_flag, sk, sk_char, & 3764 non_cyclic_l, non_cyclic_n, & 3765 non_cyclic_r, non_cyclic_s ) 3534 3766 3535 3767 3536 3768 CHARACTER (LEN = *), INTENT(IN) :: sk_char !< string identifier, used for assign fluxes to the correct dimension in the analysis array 3537 INTEGER(iwp) :: sk_num !< integer identifier, used for assign fluxes to the correct dimension in the analysis array3769 INTEGER(iwp) :: sk_num = 0 !< integer identifier, used for assign fluxes to the correct dimension in the analysis array 3538 3770 3539 INTEGER(iwp) :: i !< grid index along x-direction3540 INTEGER(iwp) :: j !< grid index along y-direction3541 INTEGER(iwp) :: k !< grid index along z-direction3542 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults3543 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults3544 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults3545 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms3546 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread3771 INTEGER(iwp) :: i !< grid index along x-direction 3772 INTEGER(iwp) :: j !< grid index along y-direction 3773 INTEGER(iwp) :: k !< grid index along z-direction 3774 INTEGER(iwp) :: k_mm !< k-2 index in disretization, can be modified to avoid segmentation faults 3775 INTEGER(iwp) :: k_pp !< k+2 index in disretization, can be modified to avoid segmentation faults 3776 INTEGER(iwp) :: k_ppp !< k+3 index in disretization, can be modified to avoid segmentation faults 3777 INTEGER(iwp) :: nzb_max_l !< index indicating upper bound for order degradation of horizontal advection terms 3778 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 3547 3779 3780 INTEGER(iwp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: & 3781 advc_flag !< flag array to control order of scalar advection 3782 3783 LOGICAL :: non_cyclic_l !< flag that indicates non-cyclic boundary on the left 3784 LOGICAL :: non_cyclic_n !< flag that indicates non-cyclic boundary on the north 3785 LOGICAL :: non_cyclic_r !< flag that indicates non-cyclic boundary on the right 3786 LOGICAL :: non_cyclic_s !< flag that indicates non-cyclic boundary on the south 3548 3787 ! 3549 3788 !-- sk is an array from parameter list. It should not be a pointer, because … … 3552 3791 !-- even worse, because the compiler cannot assume strided one in the 3553 3792 !-- caller side. 3554 REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar3793 REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar 3555 3794 3556 3795 REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction … … 3612 3851 !-- entire subdomain, in order to avoid unsymmetric loops which might be 3613 3852 !-- an issue for GPUs. 3614 IF( bc_dirichlet_l .OR. bc_radiation_l .OR. & 3615 bc_dirichlet_r .OR. bc_radiation_r .OR. & 3616 bc_dirichlet_s .OR. bc_radiation_s .OR. & 3617 bc_dirichlet_n .OR. bc_radiation_n ) THEN 3853 IF( non_cyclic_l .OR. non_cyclic_r .OR. & 3854 non_cyclic_s .OR. non_cyclic_n ) THEN 3618 3855 nzb_max_l = nzt 3619 3856 ELSE … … 3652 3889 DO k = nzb+1, nzb_max_l 3653 3890 3654 ibit2 = REAL( IBITS(advc_flag s_1(k,j,i-1),2,1), KIND = wp )3655 ibit1 = REAL( IBITS(advc_flag s_1(k,j,i-1),1,1), KIND = wp )3656 ibit0 = REAL( IBITS(advc_flag s_1(k,j,i-1),0,1), KIND = wp )3891 ibit2 = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) 3892 ibit1 = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) 3893 ibit0 = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) 3657 3894 3658 3895 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) … … 3719 3956 !$ACC PRIVATE(flux_t, diss_t, flux_d, diss_d) & 3720 3957 !$ACC PRIVATE(div, u_comp, u_comp_l, v_comp, v_comp_s) & 3721 !$ACC PRESENT(advc_flag s_1) &3958 !$ACC PRESENT(advc_flag) & 3722 3959 !$ACC PRESENT(sk, u, v, w, u_stokes_zu, v_stokes_zu) & 3723 3960 !$ACC PRESENT(drho_air, rho_air_zw, ddzw) & … … 3736 3973 DO k = nzb+1, nzb_max_l 3737 3974 3738 ibit5 = REAL( IBITS(advc_flag s_1(k,j-1,i),5,1), KIND = wp )3739 ibit4 = REAL( IBITS(advc_flag s_1(k,j-1,i),4,1), KIND = wp )3740 ibit3 = REAL( IBITS(advc_flag s_1(k,j-1,i),3,1), KIND = wp )3975 ibit5 = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) 3976 ibit4 = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) 3977 ibit3 = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) 3741 3978 3742 3979 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) … … 3798 4035 DO k = nzb+1, nzb_max_l 3799 4036 3800 ibit2 = REAL( IBITS(advc_flag s_1(k,j,i),2,1), KIND = wp )3801 ibit1 = REAL( IBITS(advc_flag s_1(k,j,i),1,1), KIND = wp )3802 ibit0 = REAL( IBITS(advc_flag s_1(k,j,i),0,1), KIND = wp )4037 ibit2 = REAL( IBITS(advc_flag(k,j,i),2,1), KIND = wp ) 4038 ibit1 = REAL( IBITS(advc_flag(k,j,i),1,1), KIND = wp ) 4039 ibit0 = REAL( IBITS(advc_flag(k,j,i),0,1), KIND = wp ) 3803 4040 3804 4041 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) … … 3836 4073 ! 3837 4074 !-- Recompute the left fluxes. 3838 ibit2_l = REAL( IBITS(advc_flag s_1(k,j,i-1),2,1), KIND = wp )3839 ibit1_l = REAL( IBITS(advc_flag s_1(k,j,i-1),1,1), KIND = wp )3840 ibit0_l = REAL( IBITS(advc_flag s_1(k,j,i-1),0,1), KIND = wp )4075 ibit2_l = REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) 4076 ibit1_l = REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) 4077 ibit0_l = REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) 3841 4078 3842 4079 u_comp_l = u(k,j,i) - u_gtrans + u_stokes_zu(k) … … 3875 4112 #endif 3876 4113 3877 ibit5 = REAL( IBITS(advc_flag s_1(k,j,i),5,1), KIND = wp )3878 ibit4 = REAL( IBITS(advc_flag s_1(k,j,i),4,1), KIND = wp )3879 ibit3 = REAL( IBITS(advc_flag s_1(k,j,i),3,1), KIND = wp )4114 ibit5 = REAL( IBITS(advc_flag(k,j,i),5,1), KIND = wp ) 4115 ibit4 = REAL( IBITS(advc_flag(k,j,i),4,1), KIND = wp ) 4116 ibit3 = REAL( IBITS(advc_flag(k,j,i),3,1), KIND = wp ) 3880 4117 3881 4118 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) … … 3913 4150 ! 3914 4151 !-- Recompute the south fluxes. 3915 ibit5_s = REAL( IBITS(advc_flag s_1(k,j-1,i),5,1), KIND = wp )3916 ibit4_s = REAL( IBITS(advc_flag s_1(k,j-1,i),4,1), KIND = wp )3917 ibit3_s = REAL( IBITS(advc_flag s_1(k,j-1,i),3,1), KIND = wp )4152 ibit5_s = REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) 4153 ibit4_s = REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) 4154 ibit3_s = REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) 3918 4155 3919 4156 v_comp_s = v(k,j,i) - v_gtrans + v_stokes_zu(k) … … 3955 4192 !-- k index has to be modified near bottom and top, else array 3956 4193 !-- subscripts will be exceeded. 3957 ibit8 = REAL( IBITS(advc_flag s_1(k,j,i),8,1), KIND = wp )3958 ibit7 = REAL( IBITS(advc_flag s_1(k,j,i),7,1), KIND = wp )3959 ibit6 = REAL( IBITS(advc_flag s_1(k,j,i),6,1), KIND = wp )4194 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) 4195 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) 4196 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) 3960 4197 3961 4198 k_ppp = k + 3 * ibit8 … … 3998 4235 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & 3999 4236 - u(k,j,i) * ( & 4000 REAL( IBITS(advc_flag s_1(k,j,i-1),0,1), KIND = wp )&4001 + REAL( IBITS(advc_flag s_1(k,j,i-1),1,1), KIND = wp )&4002 + REAL( IBITS(advc_flag s_1(k,j,i-1),2,1), KIND = wp )&4237 REAL( IBITS(advc_flag(k,j,i-1),0,1), KIND = wp ) & 4238 + REAL( IBITS(advc_flag(k,j,i-1),1,1), KIND = wp ) & 4239 + REAL( IBITS(advc_flag(k,j,i-1),2,1), KIND = wp ) & 4003 4240 ) & 4004 4241 ) * ddx & 4005 4242 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 4006 4243 - v(k,j,i) * ( & 4007 REAL( IBITS(advc_flag s_1(k,j-1,i),3,1), KIND = wp )&4008 + REAL( IBITS(advc_flag s_1(k,j-1,i),4,1), KIND = wp )&4009 + REAL( IBITS(advc_flag s_1(k,j-1,i),5,1), KIND = wp )&4244 REAL( IBITS(advc_flag(k,j-1,i),3,1), KIND = wp ) & 4245 + REAL( IBITS(advc_flag(k,j-1,i),4,1), KIND = wp ) & 4246 + REAL( IBITS(advc_flag(k,j-1,i),5,1), KIND = wp ) & 4010 4247 ) & 4011 4248 ) * ddy & … … 4014 4251 - w(k-1,j,i) * rho_air_zw(k-1) * & 4015 4252 ( & 4016 REAL( IBITS(advc_flag s_1(k-1,j,i),6,1), KIND = wp )&4017 + REAL( IBITS(advc_flag s_1(k-1,j,i),7,1), KIND = wp )&4018 + REAL( IBITS(advc_flag s_1(k-1,j,i),8,1), KIND = wp )&4253 REAL( IBITS(advc_flag(k-1,j,i),6,1), KIND = wp ) & 4254 + REAL( IBITS(advc_flag(k-1,j,i),7,1), KIND = wp ) & 4255 + REAL( IBITS(advc_flag(k-1,j,i),8,1), KIND = wp ) & 4019 4256 ) & 4020 4257 ) * drho_air(k) * ddzw(k) … … 4202 4439 !-- k index has to be modified near bottom and top, else array 4203 4440 !-- subscripts will be exceeded. 4204 ibit8 = REAL( IBITS(advc_flag s_1(k,j,i),8,1), KIND = wp )4205 ibit7 = REAL( IBITS(advc_flag s_1(k,j,i),7,1), KIND = wp )4206 ibit6 = REAL( IBITS(advc_flag s_1(k,j,i),6,1), KIND = wp )4441 ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp ) 4442 ibit7 = REAL( IBITS(advc_flag(k,j,i),7,1), KIND = wp ) 4443 ibit6 = REAL( IBITS(advc_flag(k,j,i),6,1), KIND = wp ) 4207 4444 4208 4445 k_ppp = k + 3 * ibit8 … … 4389 4626 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 4390 4627 4391 REAL(wp) :: ibit 9!< flag indicating 1st-order scheme along x-direction4392 REAL(wp) :: ibit1 0!< flag indicating 3rd-order scheme along x-direction4393 REAL(wp) :: ibit 11!< flag indicating 5th-order scheme along x-direction4628 REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction 4629 REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction 4630 REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction 4394 4631 #ifdef _OPENACC 4395 REAL(wp) :: ibit 9_l!< flag indicating 1st-order scheme along x-direction4396 REAL(wp) :: ibit1 0_l !< flag indicating 3rd-order scheme along x-direction4397 REAL(wp) :: ibit 11_l !< flag indicating 5th-order scheme along x-direction4632 REAL(wp) :: ibit0_l !< flag indicating 1st-order scheme along x-direction 4633 REAL(wp) :: ibit1_l !< flag indicating 3rd-order scheme along x-direction 4634 REAL(wp) :: ibit2_l !< flag indicating 5th-order scheme along x-direction 4398 4635 #endif 4399 REAL(wp) :: ibit 12!< flag indicating 1st-order scheme along y-direction4400 REAL(wp) :: ibit 13!< flag indicating 3rd-order scheme along y-direction4401 REAL(wp) :: ibit 14!< flag indicating 5th-order scheme along y-direction4636 REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction 4637 REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction 4638 REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction 4402 4639 #ifdef _OPENACC 4403 REAL(wp) :: ibit 12_s !< flag indicating 1st-order scheme along y-direction4404 REAL(wp) :: ibit 13_s !< flag indicating 3rd-order scheme along y-direction4405 REAL(wp) :: ibit 14_s !< flag indicating 5th-order scheme along y-direction4640 REAL(wp) :: ibit3_s !< flag indicating 1st-order scheme along y-direction 4641 REAL(wp) :: ibit4_s !< flag indicating 3rd-order scheme along y-direction 4642 REAL(wp) :: ibit5_s !< flag indicating 5th-order scheme along y-direction 4406 4643 #endif 4407 REAL(wp) :: ibit 15!< flag indicating 1st-order scheme along z-direction4408 REAL(wp) :: ibit 16!< flag indicating 3rd-order scheme along z-direction4409 REAL(wp) :: ibit 17!< flag indicating 5th-order scheme along z-direction4644 REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction 4645 REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction 4646 REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction 4410 4647 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 4411 4648 REAL(wp) :: div !< diverence on u-grid … … 4468 4705 DO k = nzb+1, nzb_max_l 4469 4706 4470 ibit 11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp )4471 ibit1 0 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp )4472 ibit 9 = REAL( IBITS(advc_flags_1(k,j,i-1),9,1),KIND = wp )4707 ibit2 = REAL( IBITS(advc_flags_m(k,j,i-1),2,1), KIND = wp ) 4708 ibit1 = REAL( IBITS(advc_flags_m(k,j,i-1),1,1), KIND = wp ) 4709 ibit0 = REAL( IBITS(advc_flags_m(k,j,i-1),0,1), KIND = wp ) 4473 4710 4474 4711 u_comp = u(k,j,i) + u(k,j,i-1) - gu 4475 4712 swap_flux_x_local_u(k,j) = u_comp * ( & 4476 ( 37.0_wp * ibit 11 * adv_mom_5&4477 + 7.0_wp * ibit1 0 * adv_mom_3&4478 + ibit 9* adv_mom_1 &4713 ( 37.0_wp * ibit2 * adv_mom_5 & 4714 + 7.0_wp * ibit1 * adv_mom_3 & 4715 + ibit0 * adv_mom_1 & 4479 4716 ) * & 4480 4717 ( u(k,j,i) + u(k,j,i-1) ) & 4481 - ( 8.0_wp * ibit 11 * adv_mom_5&4482 + ibit1 0 * adv_mom_3&4718 - ( 8.0_wp * ibit2 * adv_mom_5 & 4719 + ibit1 * adv_mom_3 &