Changeset 2232 for palm/trunk/SOURCE/init_grid.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_grid.f90 (modified) (32 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_grid.f90
r2201 r2232 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! - Adjustments according to new topography representation 23 ! - Bugfix: Move determination of nzb_max behind topography modification in 24 ! cell-edge case 25 ! - Get rid off global arrays required for topography output 26 ! - Enable topography input via netcdf 27 ! - Generic tunnel set-up added 23 28 ! 24 29 ! Former revisions: … … 212 217 ! ------------ 213 218 !> Creating grid depending constants 219 !> To Do: Setting topo flags only based on topo_3d array - flags for former 220 !> nzb_outer arrays are still not set properly. 221 !> To Do: Rearrange topo flag list 214 222 !------------------------------------------------------------------------------! 215 223 SUBROUTINE init_grid … … 230 238 dz_stretch_level, dz_stretch_level_index, grid_level, ibc_uv_b, & 231 239 io_blocks, io_group, inflow_l, inflow_n, inflow_r, inflow_s, & 232 masking_method, maximum_grid_level, message_string,&240 lod, masking_method, maximum_grid_level, message_string, & 233 241 momentum_advec, nest_domain, nest_bound_l, nest_bound_n, & 234 242 nest_bound_r, nest_bound_s, ocean, outflow_l, outflow_n, & 235 243 outflow_r, outflow_s, psolver, scalar_advec, topography, & 236 topography_grid_convention, use_surface_fluxes, use_top_fluxes, & 237 wall_adjustment_factor 244 topography_grid_convention, tunnel_height, tunnel_length, & 245 tunnel_width_x, tunnel_width_y, tunnel_wall_depth, & 246 use_surface_fluxes, use_top_fluxes, wall_adjustment_factor 238 247 239 248 USE grid_variables, & 240 ONLY: ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, fwxm, & 241 fwxp, fwym, fwyp, fxm, fxp, fym, fyp, wall_e_x, wall_e_y, & 242 wall_u, wall_v, wall_w_x, wall_w_y, zu_s_inner, zw_w_inner 249 ONLY: ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, zu_s_inner, zw_w_inner 243 250 244 251 USE indices, & 245 ONLY: flags, nbgp, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg,&246 n y, nyn, nyng, nyn_mg, nys, nys_mg, nysg, nz, nzb,&247 nzb _diff, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u,&248 nzb_ diff_v, nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner,&252 ONLY: advc_flags_1, advc_flags_2, flags, nbgp, nx, nxl, nxlg, nxl_mg, & 253 nxr, nxrg, nxr_mg, ny, nyn, nyng, nyn_mg, nys, nys_mg, nysg, nz,& 254 nzb, nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer, & 255 nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner, & 249 256 nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner, & 250 nzb_w_outer, nzt, nzt_diff, nzt_mg, rflags_invers, & 251 rflags_s_inner, wall_flags_0, wall_flags_00, wall_flags_1, & 257 nzb_w_outer, nzt, nzt_mg, wall_flags_0, wall_flags_1, & 252 258 wall_flags_10, wall_flags_2, wall_flags_3, wall_flags_4, & 253 259 wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8, & … … 255 261 256 262 USE kinds 257 263 #if defined ( __netcdf ) 264 USE netcdf_interface, & 265 ONLY: netcdf_close_file, netcdf_open_read_file, netcdf_get_attribute, & 266 netcdf_get_variable 267 #endif 258 268 USE pegrid 269 270 USE surface_mod, & 271 ONLY: init_bc 259 272 260 273 IMPLICIT NONE … … 275 288 INTEGER(iwp) :: cys !< index for south canyon wall 276 289 INTEGER(iwp) :: i !< index variable along x 290 INTEGER(iwp) :: id_topo !< NetCDF id of topograhy input file 277 291 INTEGER(iwp) :: ii !< loop variable for reading topography file 278 292 INTEGER(iwp) :: inc !< incremental parameter for coarsening grid level 279 293 INTEGER(iwp) :: j !< index variable along y 280 294 INTEGER(iwp) :: k !< index variable along z 295 INTEGER(iwp) :: k_top !< topography top index on local PE 281 296 INTEGER(iwp) :: l !< loop variable 282 297 INTEGER(iwp) :: nxl_l !< index of left PE boundary for multigrid level … … 286 301 INTEGER(iwp) :: nzb_local_max !< vertical grid index of maximum topography height 287 302 INTEGER(iwp) :: nzb_local_min !< vertical grid index of minimum topography height 288 INTEGER(iwp) :: nzb_si !< dummy index for local nzb_s_inner289 303 INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level 290 304 INTEGER(iwp) :: num_hole !< number of holes (in topography) resolved by only one grid point … … 292 306 INTEGER(iwp) :: num_wall !< number of surrounding vertical walls for a single grid point 293 307 INTEGER(iwp) :: skip_n_rows !< counting variable to skip rows while reading topography file 294 INTEGER(iwp) :: vi !< dummy for vertical influence 295 296 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: & 297 vertical_influence !< number of vertical grid points above obstacle where adjustment of near-wall mixing length is required 298 299 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_nl !< index of north-left corner location to limit near-wall mixing length 300 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_nr !< north-right 301 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_sl !< south-left 302 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_sr !< south-right 308 INTEGER(iwp) :: hv_in !< heavyside function to model inner tunnel surface 309 INTEGER(iwp) :: hv_out !< heavyside function to model outer tunnel surface 310 INTEGER(iwp) :: txe_out !< end position of outer tunnel wall in x 311 INTEGER(iwp) :: txs_out !< start position of outer tunnel wall in x 312 INTEGER(iwp) :: tye_out !< end position of outer tunnel wall in y 313 INTEGER(iwp) :: tys_out !< start position of outer tunnel wall in y 314 INTEGER(iwp) :: txe_in !< end position of inner tunnel wall in x 315 INTEGER(iwp) :: txs_in !< start position of inner tunnel wall in x 316 INTEGER(iwp) :: tye_in !< end position of inner tunnel wall in y 317 INTEGER(iwp) :: tys_in !< start position of inner tunnel wall in y 318 INTEGER(iwp) :: td !< tunnel wall depth 319 INTEGER(iwp) :: th !< height of outer tunnel wall 320 303 321 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography top at cell-center 304 322 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp !< dummy to calculate topography indices on u- and v-grid 305 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_l !< distance to adjacent left-facing wall 306 INTEGER(iwp), DIMENSION(:,: ), ALLOCATABLE :: wall_n !< north-facing307 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_r !< right-facing 308 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_s !< right-facing323 324 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo_3d !< input array for 3D topography and dummy array for setting "outer"-flags 325 326 LOGICAL :: netcdf_extend = .FALSE. !< Flag indicating wether netcdf topography input file or not 309 327 310 328 REAL(wp) :: dum !< dummy variable to skip columns while reading topography file 311 329 REAL(wp) :: dz_stretched !< stretched vertical grid spacing 312 330 331 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: oro_height !< input variable for terrain height 313 332 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: topo_height !< input variable for topography height 314 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zu_s_inner_l !< dummy array on global scale to write topography output array 315 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw_w_inner_l !< dummy array on global scale to write topography output array 316 317 333 318 334 ! 319 335 !-- Calculation of horizontal array bounds including ghost layers … … 464 480 ! 465 481 !-- Allocate outer and inner index arrays for topography and set 466 !-- defaults. 467 468 ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), & 469 corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), & 470 wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr), & 471 wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) ) 472 473 ALLOCATE( fwxm(nysg:nyng,nxlg:nxrg), fwxp(nysg:nyng,nxlg:nxrg), & 474 fwym(nysg:nyng,nxlg:nxrg), fwyp(nysg:nyng,nxlg:nxrg), & 475 fxm(nysg:nyng,nxlg:nxrg), fxp(nysg:nyng,nxlg:nxrg), & 476 fym(nysg:nyng,nxlg:nxrg), fyp(nysg:nyng,nxlg:nxrg), & 477 nzb_s_inner(nysg:nyng,nxlg:nxrg), & 482 !-- defaults. 483 ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg), & 478 484 nzb_s_outer(nysg:nyng,nxlg:nxrg), & 479 485 nzb_u_inner(nysg:nyng,nxlg:nxrg), & … … 485 491 nzb_diff_s_inner(nysg:nyng,nxlg:nxrg), & 486 492 nzb_diff_s_outer(nysg:nyng,nxlg:nxrg), & 487 nzb_diff_u(nysg:nyng,nxlg:nxrg), &488 nzb_diff_v(nysg:nyng,nxlg:nxrg), &489 493 nzb_local(nysg:nyng,nxlg:nxrg), & 490 494 nzb_tmp(nysg:nyng,nxlg:nxrg), & 491 rflags_s_inner(nzb:nzt+2,nysg:nyng,nxlg:nxrg), & 492 rflags_invers(nysg:nyng,nxlg:nxrg,nzb:nzt+2), & 493 wall_e_x(nysg:nyng,nxlg:nxrg), & 494 wall_e_y(nysg:nyng,nxlg:nxrg), & 495 wall_u(nysg:nyng,nxlg:nxrg), & 496 wall_v(nysg:nyng,nxlg:nxrg), & 497 wall_w_x(nysg:nyng,nxlg:nxrg), & 498 wall_w_y(nysg:nyng,nxlg:nxrg) ) 499 500 495 wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 496 497 ALLOCATE( topo_3d(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 498 topo_3d = 0 501 499 502 500 ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 503 504 501 505 502 nzb_s_inner = nzb; nzb_s_outer = nzb … … 507 504 nzb_v_inner = nzb; nzb_v_outer = nzb 508 505 nzb_w_inner = nzb; nzb_w_outer = nzb 509 510 rflags_s_inner = 1.0_wp511 rflags_invers = 1.0_wp512 506 513 507 ! … … 519 513 nzb_diff = nzb + 1 520 514 ENDIF 521 IF ( use_top_fluxes ) THEN522 nzt_diff = nzt - 1523 ELSE524 nzt_diff = nzt525 ENDIF526 515 527 516 nzb_diff_s_inner = nzb_diff; nzb_diff_s_outer = nzb_diff 528 nzb_diff_u = nzb_diff; nzb_diff_v = nzb_diff529 530 wall_e_x = 0.0_wp; wall_e_y = 0.0_wp; wall_u = 0.0_wp; wall_v = 0.0_wp531 wall_w_x = 0.0_wp; wall_w_y = 0.0_wp532 fwxp = 1.0_wp; fwxm = 1.0_wp; fwyp = 1.0_wp; fwym = 1.0_wp533 fxp = 1.0_wp; fxm = 1.0_wp; fyp = 1.0_wp; fym = 1.0_wp534 517 535 518 ! … … 542 525 ENDDO 543 526 l_wall(nzt+1,:,:) = l_grid(nzt) 544 545 ALLOCATE ( vertical_influence(nzb:nzt) )546 DO k = 1, nzt547 vertical_influence(k) = MIN ( INT( l_grid(k) / &548 ( wall_adjustment_factor * dzw(k) ) + 0.5_wp ), nzt - k )549 ENDDO550 527 551 528 DO k = 1, nzt … … 561 538 ENDIF 562 539 ENDDO 563 vertical_influence(0) = vertical_influence(1)564 565 DO k = nzb + 1, nzb + vertical_influence(nzb)566 l_wall(k,:,:) = zu(k) - zw(nzb)567 ENDDO568 569 540 ! 570 541 !-- Set outer and inner index arrays for non-flat topography. … … 580 551 !-- nzb_local is required for the multigrid solver 581 552 nzb_local = 0 553 ! 554 !-- Initialilize 3D topography array, used later for initializing flags 555 topo_3d(nzb+1:nzt+1,:,:) = IBSET( topo_3d(nzb+1:nzt+1,:,:), 0 ) 556 ! 557 !-- level of detail is required for output routines 558 lod = 1 582 559 583 560 CASE ( 'single_building' ) … … 587 564 blx = NINT( building_length_x / dx ) 588 565 bly = NINT( building_length_y / dy ) 589 IF ( .NOT. ocean ) THEN 590 bh = MINLOC( ABS( zw - building_height ), 1 ) - 1 591 ELSE 592 bh = MINLOC( ABS( zw - zw(0) - building_height ), 1 ) - 1 593 ENDIF 594 595 IF ( ABS( zw(bh ) - building_height ) == & 566 bh = MINLOC( ABS( zw - building_height ), 1 ) - 1 567 IF ( ABS( zw(bh) - building_height ) == & 596 568 ABS( zw(bh+1) - building_height ) ) bh = bh + 1 597 569 … … 626 598 627 599 CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp ) 600 ! 601 !-- Set bit array to mask topography 602 DO i = nxlg, nxrg 603 DO j = nysg, nyng 604 605 topo_3d(nzb_local(j,i)+1:nzt+1,j,i) = & 606 IBSET( topo_3d(nzb_local(j,i)+1:nzt+1,j,i), 0 ) 607 ENDDO 608 ENDDO 609 ! 610 !-- level of detail is required for output routines. Here, 2D topography. 611 lod = 1 612 613 CALL exchange_horiz_int( topo_3d, nbgp ) 628 614 629 615 CASE ( 'single_street_canyon' ) … … 658 644 ENDIF 659 645 660 IF ( .NOT. ocean ) THEN 661 ch = MINLOC( ABS( zw - canyon_height ), 1 ) - 1 662 ELSE 663 ch = MINLOC( ABS( zw - zw(0) - canyon_height ), 1 ) - 1 664 ENDIF 665 666 IF ( ABS( zw(ch ) - canyon_height ) == & 646 ch = MINLOC( ABS( zw - canyon_height ), 1 ) - 1 647 IF ( ABS( zw(ch) - canyon_height ) == & 667 648 ABS( zw(ch+1) - canyon_height ) ) ch = ch + 1 668 649 … … 707 688 708 689 CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp ) 690 ! 691 !-- Set bit array to mask topography 692 DO i = nxlg, nxrg 693 DO j = nysg, nyng 694 topo_3d(nzb_local(j,i)+1:nzt+1,j,i) = & 695 IBSET( topo_3d(nzb_local(j,i)+1:nzt+1,j,i), 0 ) 696 ENDDO 697 ENDDO 698 ! 699 !-- level of detail is required for output routines. Here, 2D topography. 700 lod = 1 701 702 CALL exchange_horiz_int( topo_3d, nbgp ) 703 704 CASE ( 'tunnel' ) 705 706 ! 707 !-- Tunnel height 708 IF ( tunnel_height == 9999999.9_wp ) THEN 709 th = zw( INT( 0.2 * nz) ) 710 ELSE 711 th = tunnel_height 712 ENDIF 713 ! 714 !-- Tunnel-wall depth 715 IF ( tunnel_wall_depth == 9999999.9_wp ) THEN 716 td = MAX ( dx, dy, dz ) 717 ELSE 718 td = tunnel_wall_depth 719 ENDIF 720 ! 721 !-- Check for tunnel width 722 IF ( tunnel_width_x == 9999999.9_wp .AND. & 723 tunnel_width_y == 9999999.9_wp ) THEN 724 message_string = 'No tunnel width is given. ' 725 CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 ) 726 ENDIF 727 IF ( tunnel_width_x /= 9999999.9_wp .AND. & 728 tunnel_width_y /= 9999999.9_wp ) THEN 729 message_string = 'Inconsistent tunnel parameters:' // & 730 'tunnel can only be oriented' // & 731 'either in x- or in y-direction.' 732 CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 ) 733 ENDIF 734 ! 735 !-- Tunnel axis along y 736 IF ( tunnel_width_x /= 9999999.9_wp ) THEN 737 IF ( tunnel_width_x > ( nx + 1 ) * dx ) THEN 738 message_string = 'Tunnel width too large' 739 CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 ) 740 ENDIF 741 742 txs_out = INT( ( nx + 1 ) * 0.5_wp * dx - tunnel_width_x * 0.5_wp ) 743 txe_out = INT( ( nx + 1 ) * 0.5_wp * dx + tunnel_width_x * 0.5_wp ) 744 txs_in = INT( ( nx + 1 ) * 0.5_wp * dx - & 745 ( tunnel_width_x * 0.5_wp - td ) ) 746 txe_in = INT( ( nx + 1 ) * 0.5_wp * dx + & 747 ( tunnel_width_x * 0.5_wp - td ) ) 748 749 tys_out = INT( ( ny + 1 ) * 0.5_wp * dy - tunnel_length * 0.5_wp ) 750 tye_out = INT( ( ny + 1 ) * 0.5_wp * dy + tunnel_length * 0.5_wp ) 751 tys_in = tys_out 752 tye_in = tye_out 753 ENDIF 754 IF ( tunnel_width_x /= 9999999.9_wp .AND. & 755 tunnel_width_x - 2.0_wp * tunnel_wall_depth <= 2.0_wp * dx ) THEN 756 message_string = 'Tunnel width too small' 757 CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 ) 758 ENDIF 759 IF ( tunnel_width_y /= 9999999.9_wp .AND. & 760 tunnel_width_y - 2.0_wp * tunnel_wall_depth <= 2.0_wp * dy ) THEN 761 message_string = 'Tunnel width too small' 762 CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 ) 763 ENDIF 764 ! 765 !-- Tunnel axis along x 766 IF ( tunnel_width_y /= 9999999.9_wp ) THEN 767 IF ( tunnel_width_y > ( ny + 1 ) * dy ) THEN 768 message_string = 'Tunnel width too large' 769 CALL message( 'init_grid', 'PA0207', 1, 2, 0, 6, 0 ) 770 ENDIF 771 772 txs_out = INT( ( nx + 1 ) * 0.5_wp * dx - tunnel_length * 0.5_wp ) 773 txe_out = INT( ( nx + 1 ) * 0.5_wp * dx + tunnel_length * 0.5_wp ) 774 txs_in = txs_out 775 txe_in = txe_out 776 777 tys_out = INT( ( ny + 1 ) * 0.5_wp * dy - tunnel_width_y * 0.5_wp ) 778 tye_out = INT( ( ny + 1 ) * 0.5_wp * dy + tunnel_width_y * 0.5_wp ) 779 tys_in = INT( ( ny + 1 ) * 0.5_wp * dy - & 780 ( tunnel_width_y * 0.5_wp - td ) ) 781 tye_in = INT( ( ny + 1 ) * 0.5_wp * dy + & 782 ( tunnel_width_y * 0.5_wp - td ) ) 783 ENDIF 784 785 topo_3d = 0 786 DO i = nxl, nxr 787 DO j = nys, nyn 788 ! 789 !-- Use heaviside function to model outer tunnel surface 790 hv_out = th * 0.5_wp * & 791 ( ( SIGN( 1.0_wp, i * dx - txs_out ) + 1.0_wp ) & 792 - ( SIGN( 1.0_wp, i * dx - txe_out ) + 1.0_wp ) ) 793 794 hv_out = hv_out * 0.5_wp * & 795 ( ( SIGN( 1.0_wp, j * dy - tys_out ) + 1.0_wp ) & 796 - ( SIGN( 1.0_wp, j * dy - tye_out ) + 1.0_wp ) ) 797 ! 798 !-- Use heaviside function to model inner tunnel surface 799 hv_in = ( th - td ) * 0.5_wp * & 800 ( ( SIGN( 1.0_wp, i * dx - txs_in ) + 1.0_wp ) & 801 - ( SIGN( 1.0_wp, i * dx - txe_in ) + 1.0_wp ) ) 802 803 hv_in = hv_in * 0.5_wp * & 804 ( ( SIGN( 1.0_wp, j * dy - tys_in ) + 1.0_wp ) & 805 - ( SIGN( 1.0_wp, j * dy - tye_in ) + 1.0_wp ) ) 806 ! 807 !-- Set flags at x-y-positions without any tunnel surface 808 IF ( hv_out - hv_in == 0.0_wp ) THEN 809 topo_3d(nzb+1:nzt+1,j,i) = IBSET( topo_3d(nzb+1:nzt+1,j,i), 0 ) 810 ! 811 !-- Set flags at x-y-positions with tunnel surfaces 812 ELSE 813 DO k = nzb + 1, nzt + 1 814 ! 815 !-- Inner tunnel 816 IF ( hv_out - hv_in == th ) THEN 817 IF ( zw(k) <= hv_out ) THEN 818 topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 ) 819 ELSE 820 topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 ) 821 ENDIF 822 ENDIF 823 ! 824 !-- Lateral tunnel walls 825 IF ( hv_out - hv_in == td ) THEN 826 IF ( zw(k) <= hv_in ) THEN 827 topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 ) 828 ELSEIF ( zw(k) > hv_in .AND. zw(k) <= hv_out ) THEN 829 topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 ) 830 ELSEIF ( zw(k) > hv_out ) THEN 831 topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 ) 832 ENDIF 833 ENDIF 834 ENDDO 835 ENDIF 836 ENDDO 837 ENDDO 838 839 nzb_local = 0 840 ! 841 !-- level of detail is required for output routines. Here, 3D topography. 842 lod = 2 843 844 CALL exchange_horiz_int( topo_3d, nbgp ) 709 845 710 846 CASE ( 'read_from_file' ) 711 847 848 ALLOCATE ( oro_height(nys:nyn,nxl:nxr) ) 712 849 ALLOCATE ( topo_height(nys:nyn,nxl:nxr) ) 850 oro_height = 0.0_wp 851 topo_height = 0.0_wp 713 852 714 853 DO ii = 0, io_blocks-1 … … 717 856 ! 718 857 !-- Arbitrary irregular topography data in PALM format (exactly 719 !-- matching the grid size and total domain size) 720 OPEN( 90, FILE='TOPOGRAPHY_DATA'//TRIM( coupling_char ), & 721 STATUS='OLD', FORM='FORMATTED', ERR=10 ) 722 ! 723 !-- Read topography PE-wise. Rows are read from nyn to nys, columns 724 !-- are read from nxl to nxr. At first, ny-nyn rows need to be skipped. 725 skip_n_rows = 0 726 DO WHILE ( skip_n_rows < ny - nyn ) 727 READ( 90, * ) 728 skip_n_rows = skip_n_rows + 1 729 ENDDO 730 ! 731 !-- Read data from nyn to nys and nxl to nxr. Therefore, skip 732 !-- column until nxl-1 is reached 733 DO j = nyn, nys, -1 734 READ( 90, *, ERR=11, END=11 ) & 858 !-- matching the grid size and total domain size). 859 !-- First, check if NetCDF file for topography exist or not. 860 !-- This case, read topography from NetCDF, else read it from 861 !-- ASCII file. 862 #if defined ( __netcdf ) 863 INQUIRE( FILE='TOPOGRAPHY_DATA_NC'//TRIM( coupling_char ), & 864 EXIST=netcdf_extend ) 865 ! 866 !-- NetCDF branch 867 IF ( netcdf_extend ) THEN 868 ! 869 !-- Open file in read-only mode 870 CALL netcdf_open_read_file( 'TOPOGRAPHY_DATA_NC', & 871 id_topo, 20 ) !Error number still need to be set properly 872 873 ! 874 !-- Read terrain height. Reading is done PE-wise, i.e. each 875 !-- processor reads its own domain. Reading is realized 876 !-- via looping over x-dimension, i.e. calling 877 !-- netcdf_get_variable reads topography along y for given x. 878 !-- Orography is 2D. 879 DO i = nxl, nxr 880 CALL netcdf_get_variable( id_topo, 'orography_0', & 881 i, oro_height(:,i), 20 ) !Error number still need to be set properly 882 ENDDO 883 ! 884 !-- Read attribute lod (level of detail), required for variable 885 !-- buildings_0 886 CALL netcdf_get_attribute( id_topo, "lod", lod, .FALSE., & 887 20, 'buildings_0' ) !Error number still need to be set properly 888 ! 889 !-- Read building height 890 !-- 2D for lod = 1, 3D for lod = 2 891 IF ( lod == 1 ) THEN 892 DO i = nxl, nxr 893 CALL netcdf_get_variable( id_topo, 'buildings_0', & 894 i, topo_height(:,i), 20 ) !Error number still need to be set properly 895 ENDDO 896 897 ELSEIF ( lod == 2 ) THEN 898 ! 899 !-- Read data PE-wise. Read yz-slices. 900 DO i = nxl, nxr 901 DO j = nys, nyn 902 CALL netcdf_get_variable( id_topo, 'buildings_0', & 903 i, j, topo_3d(:,j,i), 20 ) !Error number still need to be set properly 904 ENDDO 905 ENDDO 906 ELSE 907 message_string = 'NetCDF attribute lod ' // & 908 '(level of detail) is not set properly.' 909 CALL message( 'init_grid', 'PA0208', 1, 2, 0, 6, 0 ) !Error number still need to be set properly 910 ENDIF 911 ! 912 !-- Close topography input file 913 CALL netcdf_close_file( id_topo, 20 ) 914 #endif 915 916 ! 917 !-- ASCII branch. Please note, reading of 3D topography is not 918 !-- supported in ASCII format. Further, no distinction is made 919 !-- between orography and buildings 920 ELSE 921 922 OPEN( 90, FILE='TOPOGRAPHY_DATA'//TRIM( coupling_char ), & 923 STATUS='OLD', FORM='FORMATTED', ERR=10 ) 924 ! 925 !-- Read topography PE-wise. Rows are read from nyn to nys, columns 926 !-- are read from nxl to nxr. At first, ny-nyn rows need to be skipped. 927 skip_n_rows = 0 928 DO WHILE ( skip_n_rows < ny - nyn ) 929 READ( 90, * ) 930 skip_n_rows = skip_n_rows + 1 931 ENDDO 932 ! 933 !-- Read data from nyn to nys and nxl to nxr. Therefore, skip 934 !-- column until nxl-1 is reached 935 DO j = nyn, nys, -1 936 READ( 90, *, ERR=11, END=11 ) & 735 937 ( dum, i = 0, nxl-1 ), & 736 938 ( topo_height(j,i), i = nxl, nxr ) 737 ENDDO738 739 GOTO 12939 ENDDO 940 941 GOTO 12 740 942 741 10 message_string = 'file TOPOGRAPHY'//TRIM( coupling_char )// & 742 ' does not exist' 743 CALL message( 'init_grid', 'PA0208', 1, 2, 0, 6, 0 ) 744 745 11 message_string = 'errors in file TOPOGRAPHY_DATA'// & 746 TRIM( coupling_char ) 747 CALL message( 'init_grid', 'PA0209', 1, 2, 0, 6, 0 ) 748 749 12 CLOSE( 90 ) 943 10 message_string = 'file TOPOGRAPHY'//TRIM( coupling_char )// & 944 ' does not exist' 945 CALL message( 'init_grid', 'PA0208', 1, 2, 0, 6, 0 ) 946 947 11 message_string = 'errors in file TOPOGRAPHY_DATA'// & 948 TRIM( coupling_char ) 949 CALL message( 'init_grid', 'PA0209', 1, 2, 0, 6, 0 ) 950 951 12 CLOSE( 90 ) 952 953 ENDIF 750 954 751 955 ENDIF … … 761 965 DO j = nys, nyn 762 966 IF ( .NOT. ocean ) THEN 763 nzb_local(j,i) = MINLOC( ABS( zw - topo_height(j,i) ), 1 ) - 1 764 IF ( ABS( zw(nzb_local(j,i) ) - topo_height(j,i) ) == & 765 ABS( zw(nzb_local(j,i)+1) - topo_height(j,i) ) ) & 967 nzb_local(j,i) = MINLOC( ABS( zw - topo_height(j,i) & 968 - oro_height(j,i) ), 1 ) - 1 969 IF ( ABS( zw(nzb_local(j,i) ) - topo_height(j,i) & 970 - oro_height(j,i) ) == & 971 ABS( zw(nzb_local(j,i)+1) - topo_height(j,i) & 972 - oro_height(j,i) ) ) & 766 973 nzb_local(j,i) = nzb_local(j,i) + 1 767 974 ELSE 768 975 nzb_local(j,i) = MINLOC( ABS( zw - zw(0) & 769 - topo_height(j,i) ), 1 ) - 1 976 - topo_height(j,i) & 977 - oro_height(j,i) ), 1 ) - 1 770 978 IF ( ABS( zw(nzb_local(j,i) ) - zw(0) & 771 - topo_height(j,i) ) == & 979 - topo_height(j,i) & 980 - oro_height(j,i) ) == & 772 981 ABS( zw(nzb_local(j,i)+1) - zw(0) & 773 - topo_height(j,i) ) ) & 982 - topo_height(j,i) & 983 - oro_height(j,i) ) ) & 774 984 nzb_local(j,i) = nzb_local(j,i) + 1 775 985 ENDIF … … 778 988 ENDDO 779 989 990 DEALLOCATE ( oro_height ) 780 991 DEALLOCATE ( topo_height ) 781 992 ! … … 836 1047 CALL message( 'init_grid', 'PA0430', 0, 0, 0, 6, 0 ) 837 1048 ENDIF 1049 1050 ! 1051 !-- Set bit array to mask topography. Only required for lod = 1 1052 IF ( lod == 1 ) THEN 1053 DO i = nxlg, nxrg 1054 DO j = nysg, nyng 1055 topo_3d(nzb_local(j,i)+1:nzt+1,j,i) = & 1056 IBSET( topo_3d(nzb_local(j,i)+1:nzt+1,j,i), 0 ) 1057 ENDDO 1058 ENDDO 1059 ENDIF 838 1060 ! 839 1061 !-- Exchange ghost-points, as well as add cyclic or Neumann boundary 840 1062 !-- conditions. 1063 CALL exchange_horiz_int( topo_3d, nbgp ) 841 1064 CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp ) 842 1065 843 1066 IF ( .NOT. bc_ns_cyc ) THEN 1067 IF ( nys == 0 ) topo_3d(:,-1,:) = topo_3d(:,0,:) 1068 IF ( nyn == ny ) topo_3d(:,ny+1,:) = topo_3d(:,ny,:) 1069 844 1070 IF ( nys == 0 ) nzb_local(-1,:) = nzb_local(0,:) 845 1071 IF ( nyn == ny ) nzb_local(ny+1,:) = nzb_local(ny,:) … … 847 1073 848 1074 IF ( .NOT. bc_lr_cyc ) THEN 1075 IF ( nxl == 0 ) topo_3d(:,:,-1) = topo_3d(:,:,0) 1076 IF ( nxr == nx ) topo_3d(:,:,nx+1) = topo_3d(:,:,nx) 1077 849 1078 IF ( nxl == 0 ) nzb_local(:,-1) = nzb_local(:,0) 850 1079 IF ( nxr == nx ) nzb_local(:,nx+1) = nzb_local(:,nx) 851 1080 ENDIF 1081 852 1082 853 1083 CASE DEFAULT … … 857 1087 !-- case in the user interface. There, the subroutine user_init_grid 858 1088 !-- checks which of these two conditions applies. 859 CALL user_init_grid( nzb_local)1089 CALL user_init_grid( topo_3d ) 860 1090 861 1091 END SELECT 862 !863 !-- Determine the maximum level of topography. Furthermore it is used for864 !-- steering the degradation of order of the applied advection scheme.865 !-- In case of non-cyclic lateral boundaries, the order of the advection866 !-- scheme has to be reduced up to nzt (required at the lateral boundaries).867 #if defined( __parallel )868 CALL MPI_ALLREDUCE( MAXVAL( nzb_local ) + 1, nzb_max, 1, MPI_INTEGER, &869 MPI_MAX, comm2d, ierr )870 #else871 nzb_max = MAXVAL( nzb_local ) + 1872 #endif873 IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR. &874 inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s .OR. &875 nest_domain ) &876 THEN877 nzb_max = nzt878 ENDIF879 880 1092 ! 881 1093 !-- Consistency checks and index array initialization are only required for … … 901 1113 '&MAXVAL( nzb_local ) = ', nzb_local_max 902 1114 CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 ) 1115 ENDIF 1116 ! 1117 !-- In case of non-flat topography, check whether the convention how to 1118 !-- define the topography grid has been set correctly, or whether the default 1119 !-- is applicable. If this is not possible, abort. 1120 IF ( TRIM( topography_grid_convention ) == ' ' ) THEN 1121 IF ( TRIM( topography ) /= 'single_building' .AND. & 1122 TRIM( topography ) /= 'single_street_canyon' .AND. & 1123 TRIM( topography ) /= 'tunnel' .AND. & 1124 TRIM( topography ) /= 'read_from_file') THEN 1125 !-- The default value is not applicable here, because it is only valid 1126 !-- for the two standard cases 'single_building' and 'read_from_file' 1127 !-- defined in init_grid. 1128 WRITE( message_string, * ) & 1129 'The value for "topography_grid_convention" ', & 1130 'is not set. Its default value is & only valid for ', & 1131 '"topography" = ''single_building'', ', & 1132 '''single_street_canyon'' & or ''read_from_file''.', & 1133 ' & Choose ''cell_edge'' or ''cell_center''.' 1134 CALL message( 'init_grid', 'PA0239', 1, 2, 0, 6, 0 ) 1135 ELSE 1136 !-- The default value is applicable here. 1137 !-- Set convention according to topography. 1138 IF ( TRIM( topography ) == 'single_building' .OR. & 1139 TRIM( topography ) == 'single_street_canyon' ) THEN 1140 topography_grid_convention = 'cell_edge' 1141 ELSEIF ( TRIM( topography ) == 'read_from_file' .OR. & 1142 TRIM( topography ) == 'tunnel') THEN 1143 topography_grid_convention = 'cell_center' 1144 ENDIF 1145 ENDIF 1146 ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND. & 1147 TRIM( topography_grid_convention ) /= 'cell_center' ) THEN 1148 WRITE( message_string, * ) & 1149 'The value for "topography_grid_convention" is ', & 1150 'not recognized. & Choose ''cell_edge'' or ''cell_center''.' 1151 CALL message( 'init_grid', 'PA0240', 1, 2, 0, 6, 0 ) 903 1152 ENDIF 904 1153 … … 964 1213 nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j+1,i) ) 965 1214 ENDDO 966 ENDDO 1215 ENDDO 967 1216 ! 968 1217 !-- Exchange ghost points 969 1218 CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp ) 1219 ! 1220 !-- Apply cell-edge convention also for 3D topo array. The former setting 1221 !-- of nzb_local will be removed later. 1222 DO j = nys+1, nyn+1 1223 DO i = nxl-1, nxr 1224 DO k = nzb, nzt+1 1225 IF ( BTEST( topo_3d(k,j,i), 0 ) .OR. & 1226 BTEST( topo_3d(k,j,i+1), 0 ) ) & 1227 topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 ) 1228 ENDDO 1229 ENDDO 1230 ENDDO 1231 CALL exchange_horiz_int( topo_3d, nbgp ) 1232 1233 DO i = nxl, nxr+1 1234 DO j = nys-1, nyn 1235 DO k = nzb, nzt+1 1236 IF ( BTEST( topo_3d(k,j,i), 0 ) .OR. & 1237 BTEST( topo_3d(k,j+1,i), 0 ) ) & 1238 topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 0 ) 1239 ENDDO 1240 ENDDO 1241 ENDDO 1242 CALL exchange_horiz_int( topo_3d, nbgp ) 1243 970 1244 ENDIF 1245 971 1246 ! 972 1247 !-- Initialize index arrays nzb_s_inner and nzb_w_inner 1248 973 1249 nzb_s_inner = nzb_local 974 1250 nzb_w_inner = nzb_local … … 1103 1379 CALL exchange_horiz_2d_int( nzb_s_outer, nys, nyn, nxl, nxr, nbgp ) 1104 1380 1105 !1106 !-- Allocate and set the arrays containing the topography height1107 ALLOCATE( zu_s_inner(0:nx+1,0:ny+1), zw_w_inner(0:nx+1,0:ny+1), &1108 zu_s_inner_l(0:nx+1,0:ny+1), zw_w_inner_l(0:nx+1,0:ny+1) )1109 1110 zu_s_inner = 0.0_wp1111 zw_w_inner = 0.0_wp1112 zu_s_inner_l = 0.0_wp1113 zw_w_inner_l = 0.0_wp1114 1115 DO i = nxl, nxr1116 DO j = nys, nyn1117 zu_s_inner_l(i,j) = zu(nzb_local(j,i))1118 zw_w_inner_l(i,j) = zw(nzb_local(j,i))1119 ENDDO1120 ENDDO1121 1122 #if defined( __parallel )1123 CALL MPI_REDUCE( zu_s_inner_l, zu_s_inner, (nx+2)*(ny+2), &1124 MPI_REAL, MPI_SUM, 0, comm2d, ierr )1125 CALL MPI_REDUCE( zw_w_inner_l, zw_w_inner, (nx+2)*(ny+2), &1126 MPI_REAL, MPI_SUM, 0, comm2d, ierr )1127 #else1128 zu_s_inner = zu_s_inner_l1129 zw_w_inner = zw_w_inner_l1130 #endif1131 1132 DEALLOCATE( zu_s_inner_l, zw_w_inner_l )1133 IF ( myid /= 0 ) DEALLOCATE( zu_s_inner, zw_w_inner )1134 !1135 !-- Set south and left ghost points, required for netcdf output1136 IF ( myid == 0 ) THEN1137 IF( bc_lr_cyc ) THEN1138 zu_s_inner(nx+1,:) = zu_s_inner(0,:)1139 zw_w_inner(nx+1,:) = zw_w_inner(0,:)1140 ELSE1141 zu_s_inner(nx+1,:) = zu_s_inner(nx,:)1142 zw_w_inner(nx+1,:) = zw_w_inner(nx,:)1143 ENDIF1144 IF( bc_ns_cyc ) THEN1145 zu_s_inner(:,ny+1) = zu_s_inner(:,0)1146 zw_w_inner(:,ny+1) = zw_w_inner(:,0)1147 ELSE1148 zu_s_inner(:,ny+1) = zu_s_inner(:,ny)1149 zw_w_inner(:,ny+1) = zw_w_inner(:,ny)1150 ENDIF1151 ENDIF1152 !1153 !-- Set flag arrays to be used for masking of grid points1154 DO i = nxlg, nxrg1155 DO j = nysg, nyng1156 DO k = nzb, nzt+11157 IF ( k <= nzb_s_inner(j,i) ) rflags_s_inner(k,j,i) = 0.0_wp1158 IF ( k <= nzb_s_inner(j,i) ) rflags_invers(j,i,k) = 0.0_wp1159 ENDDO1160 ENDDO1161 ENDDO1162 1163 1381 ENDIF 1164 1382 ! … … 1166 1384 !-- grid-levels further below. 1167 1385 DEALLOCATE( nzb_tmp ) 1386 1387 ! 1388 !-- Determine the maximum level of topography. It is used for 1389 !-- steering the degradation of order of the applied advection scheme. 1390 !-- In case of non-cyclic lateral boundaries, the order of the advection 1391 !-- scheme has to be reduced up to nzt (required at the lateral boundaries). 1392 k_top = 0 1393 DO i = nxl, nxr 1394 DO j = nys, nyn 1395 DO k = nzb, nzt + 1 1396 k_top = MAX( k_top, MERGE( k, 0, & 1397 .NOT. BTEST( topo_3d(k,j,i), 0 ) ) ) 1398 ENDDO 1399 ENDDO 1400 ENDDO 1401 #if defined( __parallel ) 1402 CALL MPI_ALLREDUCE( k_top + 1, nzb_max, 1, MPI_INTEGER, & !is +1 really necessary here? 1403 MPI_MAX, comm2d, ierr ) 1404 #else 1405 nzb_max = k_top + 1 1406 #endif 1407 IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR. & 1408 inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s .OR. & 1409 nest_domain ) & 1410 THEN 1411 nzb_max = nzt 1412 ENDIF 1413 ! 1414 !-- Finally, if topography extents up to the model top, limit nzb_max to nzt. 1415 nzb_max = MIN( nzb_max, nzt ) 1168 1416 1169 1417 ! … … 1172 1420 !-- applied 1173 1421 IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN 1174 nzb_diff_u = nzb_u_inner + 21175 nzb_diff_v = nzb_v_inner + 21176 1422 nzb_diff_s_inner = nzb_s_inner + 2 1177 1423 nzb_diff_s_outer = nzb_s_outer + 2 1178 1424 ELSE 1179 nzb_diff_u = nzb_u_inner + 11180 nzb_diff_v = nzb_v_inner + 11181 1425 nzb_diff_s_inner = nzb_s_inner + 1 1182 1426 nzb_diff_s_outer = nzb_s_outer + 1 … … 1184 1428 1185 1429 ! 1186 !-- Calculation of wall switches and factors required by diffusion_u/v.f90 and 1187 !-- for limitation of near-wall mixing length l_wall further below 1188 corner_nl = 0 1189 corner_nr = 0 1190 corner_sl = 0 1191 corner_sr = 0 1192 wall_l = 0 1193 wall_n = 0 1194 wall_r = 0 1195 wall_s = 0 1196 1197 DO i = nxl, nxr 1198 DO j = nys, nyn 1199 ! 1200 !-- u-component 1201 IF ( nzb_u_outer(j,i) > nzb_u_outer(j+1,i) ) THEN 1202 wall_u(j,i) = 1.0_wp ! north wall (location of adjacent fluid) 1203 fym(j,i) = 0.0_wp 1204 fyp(j,i) = 1.0_wp 1205 ELSEIF ( nzb_u_outer(j,i) > nzb_u_outer(j-1,i) ) THEN 1206 wall_u(j,i) = 1.0_wp ! south wall (location of adjacent fluid) 1207 fym(j,i) = 1.0_wp 1208 fyp(j,i) = 0.0_wp 1209 ENDIF 1210 ! 1211 !-- v-component 1212 IF ( nzb_v_outer(j,i) > nzb_v_outer(j,i+1) ) THEN 1213 wall_v(j,i) = 1.0_wp ! rigth wall (location of adjacent fluid) 1214 fxm(j,i) = 0.0_wp 1215 fxp(j,i) = 1.0_wp 1216 ELSEIF ( nzb_v_outer(j,i) > nzb_v_outer(j,i-1) ) THEN 1217 wall_v(j,i) = 1.0_wp ! left wall (location of adjacent fluid) 1218 fxm(j,i) = 1.0_wp 1219 fxp(j,i) = 0.0_wp 1220 ENDIF 1221 ! 1222 !-- w-component, also used for scalars, separate arrays for shear 1223 !-- production of tke 1224 IF ( nzb_w_outer(j,i) > nzb_w_outer(j+1,i) ) THEN 1225 wall_e_y(j,i) = 1.0_wp ! north wall (location of adjacent fluid) 1226 wall_w_y(j,i) = 1.0_wp 1227 fwym(j,i) = 0.0_wp 1228 fwyp(j,i) = 1.0_wp 1229 ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j-1,i) ) THEN 1230 wall_e_y(j,i) = -1.0_wp ! south wall (location of adjacent fluid) 1231 wall_w_y(j,i) = 1.0_wp 1232 fwym(j,i) = 1.0_wp 1233 fwyp(j,i) = 0.0_wp 1234 ENDIF 1235 IF ( nzb_w_outer(j,i) > nzb_w_outer(j,i+1) ) THEN 1236 wall_e_x(j,i) = 1.0_wp ! right wall (location of adjacent fluid) 1237 wall_w_x(j,i) = 1.0_wp 1238 fwxm(j,i) = 0.0_wp 1239 fwxp(j,i) = 1.0_wp 1240 ELSEIF ( nzb_w_outer(j,i) > nzb_w_outer(j,i-1) ) THEN 1241 wall_e_x(j,i) = -1.0_wp ! left wall (location of adjacent fluid) 1242 wall_w_x(j,i) = 1.0_wp 1243 fwxm(j,i) = 1.0_wp 1244 fwxp(j,i) = 0.0_wp 1245 ENDIF 1246 ! 1247 !-- Wall and corner locations inside buildings for limitation of 1248 !-- near-wall mixing length l_wall 1249 IF ( nzb_s_inner(j,i) > nzb_s_inner(j+1,i) ) THEN 1250 1251 wall_n(j,i) = nzb_s_inner(j+1,i) + 1 ! North wall 1252 1253 IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN 1254 corner_nl(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northleft corner 1255 nzb_s_inner(j,i-1) ) + 1 1256 ENDIF 1257 1258 IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN 1259 corner_nr(j,i) = MAX( nzb_s_inner(j+1,i), & ! Northright corner 1260 nzb_s_inner(j,i+1) ) + 1 1261 ENDIF 1262 1263 ENDIF 1264 1265 IF ( nzb_s_inner(j,i) > nzb_s_inner(j-1,i) ) THEN 1266 1267 wall_s(j,i) = nzb_s_inner(j-1,i) + 1 ! South wall 1268 IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN 1269 corner_sl(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southleft corner 1270 nzb_s_inner(j,i-1) ) + 1 1271 ENDIF 1272 1273 IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN 1274 corner_sr(j,i) = MAX( nzb_s_inner(j-1,i), & ! Southright corner 1275 nzb_s_inner(j,i+1) ) + 1 1276 ENDIF 1277 1278 ENDIF 1279 1280 IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i-1) ) THEN 1281 wall_l(j,i) = nzb_s_inner(j,i-1) + 1 ! Left wall 1282 ENDIF 1283 1284 IF ( nzb_s_inner(j,i) > nzb_s_inner(j,i+1) ) THEN 1285 wall_r(j,i) = nzb_s_inner(j,i+1) + 1 ! Right wall 1286 ENDIF 1430 !-- Set-up topography flags. First, set flags only for s, u, v and w-grid. 1431 !-- Further special flags will be set in following loops. 1432 wall_flags_0 = 0 1433 DO j = nys, nyn 1434 DO i = nxl, nxr 1435 DO k = nzb, nzt+1 1436 ! 1437 !-- scalar grid 1438 IF ( BTEST( topo_3d(k,j,i), 0 ) ) & 1439 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 ) 1440 ! 1441 !-- v grid 1442 IF ( BTEST( topo_3d(k,j,i), 0 ) .AND. & 1443 BTEST( topo_3d(k,j-1,i), 0 ) ) & 1444 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 ) 1445 ! 1446 !-- To do: set outer arrays on basis of topo_3d array, adjust for downward-facing walls 1447 !-- s grid outer array 1448 IF ( k >= nzb_s_outer(j,i) ) & 1449 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 ) 1450 ! 1451 !-- s grid outer array 1452 IF ( k >= nzb_u_outer(j,i) ) & 1453 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 ) 1454 ! 1455 !-- s grid outer array 1456 IF ( k >= nzb_v_outer(j,i) ) & 1457 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 ) 1458 ! 1459 !-- w grid outer array 1460 IF ( k >= nzb_w_outer(j,i) ) & 1461 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 ) 1462 ENDDO 1463 1464 DO k = nzb, nzt 1465 ! 1466 !-- w grid 1467 IF ( BTEST( topo_3d(k,j,i), 0 ) .AND. & 1468 BTEST( topo_3d(k+1,j,i), 0 ) ) & 1469 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 ) 1470 ENDDO 1471 wall_flags_0(nzt+1,j,i) = IBSET( wall_flags_0(nzt+1,j,i), 3 ) 1287 1472 1288 1473 ENDDO 1289 1474 ENDDO 1475 ! 1476 !-- u grid. Note, reverse 1477 !-- memory access is required for setting flag on u-grid 1478 DO j = nys, nyn 1479 DO i = nxl, nxr 1480 DO k = nzb, nzt+1 1481 IF ( BTEST( topo_3d(k,j,i), 0 ) .AND. & 1482 BTEST( topo_3d(k,j,i-1), 0 ) ) & 1483 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 ) 1484 ENDDO 1485 ENDDO 1486 ENDDO 1487 ! 1488 !-- Set further special flags 1489 DO i = nxl, nxr 1490 DO j = nys, nyn 1491 DO k = nzb, nzt+1 1492 ! 1493 !-- scalar grid, former nzb_diff_s_inner. 1494 !-- Note, use this flag also to mask topography in diffusion_u and 1495 !-- diffusion_v along the vertical direction. In case of 1496 !-- use_surface_fluxes, fluxes are calculated via MOST, else, simple 1497 !-- gradient approach is applied. Please note, in case of u- and v- 1498 !-- diffuison, a small error is made at edges (on the east side for u, 1499 !-- at the north side for v), since topography on scalar grid point 1500 !-- is used instead of topography on u/v-grid. As number of topography grid 1501 !-- points on uv-grid is different than s-grid, different number of 1502 !-- surface elements would be required. In order to avoid this, 1503 !-- treat edges (u(k,j,i+1)) simply by a gradient approach, i.e. these 1504 !-- points are not masked within diffusion_u. Tests had shown that the 1505 !-- effect on the flow is negligible. 1506 IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN 1507 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) & 1508 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 ) 1509 ELSE 1510 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 ) 1511 ENDIF 1512 1513 ENDDO 1514 ! 1515 !-- Special flag to control vertical diffusion at model top - former 1516 !-- nzt_diff 1517 wall_flags_0(:,j,i) = IBSET( wall_flags_0(:,j,i), 9 ) 1518 IF ( use_top_fluxes ) & 1519 wall_flags_0(nzt:nzt+1,j,i) = IBCLR( wall_flags_0(nzt:nzt+1,j,i), 9 ) 1520 1521 DO k = nzb+1, nzt 1522 ! 1523 !-- Special flag on u grid, former nzb_u_inner + 1, required 1524 !-- for disturb_field and initialization. Do not disturb directly at 1525 !-- topography, as well as initialize u with zero one grid point outside 1526 !-- of topography. 1527 IF ( BTEST( wall_flags_0(k-1,j,i), 1 ) .AND. & 1528 BTEST( wall_flags_0(k,j,i), 1 ) .AND. & 1529 BTEST( wall_flags_0(k+1,j,i), 1 ) ) & 1530 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 ) 1531 ! 1532 !-- Special flag on v grid, former nzb_v_inner + 1, required 1533 !-- for disturb_field and initialization. Do not disturb directly at 1534 !-- topography, as well as initialize v with zero one grid point outside 1535 !-- of topography. 1536 IF ( BTEST( wall_flags_0(k-1,j,i), 2 ) .AND. & 1537 BTEST( wall_flags_0(k,j,i), 2 ) .AND. & 1538 BTEST( wall_flags_0(k+1,j,i), 2 ) ) & 1539 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 ) 1540 ! 1541 !-- Special flag on scalar grid, former nzb_s_inner+1. Used for 1542 !-- lpm_sgs_tke 1543 IF ( BTEST( wall_flags_0(k,j,i), 0 ) .AND. & 1544 BTEST( wall_flags_0(k-1,j,i), 0 ) .AND. & 1545 BTEST( wall_flags_0(k+1,j,i), 0 ) ) & 1546 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 ) 1547 ! 1548 !-- Special flag on scalar grid, nzb_diff_s_outer - 1, required in 1549 !-- in production_e 1550 IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN 1551 IF ( BTEST( wall_flags_0(k,j,i), 24 ) .AND. & 1552 BTEST( wall_flags_0(k-1,j,i), 24 ) .AND. & 1553 BTEST( wall_flags_0(k+1,j,i), 0 ) ) & 1554 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 29 ) 1555 ELSE 1556 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) & 1557 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 29 ) 1558 ENDIF 1559 ! 1560 !-- Special flag on scalar grid, nzb_diff_s_outer - 1, required in 1561 !-- in production_e 1562 IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN 1563 IF ( BTEST( wall_flags_0(k,j,i), 0 ) .AND. & 1564 BTEST( wall_flags_0(k-1,j,i), 0 ) .AND. & 1565 BTEST( wall_flags_0(k+1,j,i), 0 ) ) & 1566 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 ) 1567 ELSE 1568 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) & 1569 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 ) 1570 ENDIF 1571 ENDDO 1572 ! 1573 !-- Flags indicating downward facing walls 1574 DO k = nzb+1, nzt 1575 ! 1576 !-- Scalar grid 1577 IF ( BTEST( wall_flags_0(k-1,j,i), 0 ) .AND. & 1578 .NOT. BTEST( wall_flags_0(k,j,i), 0 ) ) & 1579 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 ) 1580 ! 1581 !-- Downward facing wall on u grid 1582 IF ( BTEST( wall_flags_0(k-1,j,i), 1 ) .AND. & 1583 .NOT. BTEST( wall_flags_0(k,j,i), 1 ) ) & 1584 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 ) 1585 ! 1586 !-- Downward facing wall on v grid 1587 IF ( BTEST( wall_flags_0(k-1,j,i), 2 ) .AND. & 1588 .NOT. BTEST( wall_flags_0(k,j,i), 2 ) ) & 1589 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 ) 1590 ! 1591 !-- Downward facing wall on w grid 1592 IF ( BTEST( wall_flags_0(k-1,j,i), 3 ) .AND. & 1593 .NOT. BTEST( wall_flags_0(k,j,i), 3 ) ) & 1594 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 ) 1595 ENDDO 1596 ! 1597 !-- Flags indicating upward facing walls 1598 DO k = nzb, nzt 1599 ! 1600 !-- Upward facing wall on scalar grid 1601 IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) .AND. & 1602 BTEST( wall_flags_0(k+1,j,i), 0 ) ) & 1603 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 ) 1604 ! 1605 !-- Upward facing wall on u grid 1606 IF ( .NOT. BTEST( wall_flags_0(k,j,i), 1 ) .AND. & 1607 BTEST( wall_flags_0(k+1,j,i), 1 ) ) & 1608 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 ) 1609 1610 ! 1611 !-- Upward facing wall on v grid 1612 IF ( .NOT. BTEST( wall_flags_0(k,j,i), 2 ) .AND. & 1613 BTEST( wall_flags_0(k+1,j,i), 2 ) ) & 1614 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 ) 1615 1616 ! 1617 !-- Upward facing wall on w grid 1618 IF ( .NOT. BTEST( wall_flags_0(k,j,i), 3 ) .AND. & 1619 BTEST( wall_flags_0(k+1,j,i), 3 ) ) & 1620 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 ) 1621 ! 1622 !-- Special flag on scalar grid, former nzb_s_inner 1623 IF ( BTEST( wall_flags_0(k,j,i), 0 ) .OR. & 1624 BTEST( wall_flags_0(k,j,i), 12 ) .OR. & 1625 BTEST( wall_flags_0(k,j,i), 13 ) ) & 1626 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 ) 1627 ! 1628 !-- Special flag on scalar grid, nzb_diff_s_inner - 1, required for 1629 !-- flow_statistics 1630 IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN 1631 IF ( BTEST( wall_flags_0(k,j,i), 0 ) .AND. & 1632 BTEST( wall_flags_0(k+1,j,i), 0 ) ) & 1633 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 ) 1634 ELSE 1635 IF ( BTEST( wall_flags_0(k,j,i), 22 ) ) & 1636 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 ) 1637 ENDIF 1638 1639 1640 ENDDO 1641 wall_flags_0(nzt+1,j,i) = IBSET( wall_flags_0(nzt+1,j,i), 22 ) 1642 wall_flags_0(nzt+1,j,i) = IBSET( wall_flags_0(nzt+1,j,i), 23 ) 1643 ENDDO 1644 ENDDO 1645 ! 1646 !-- Exchange ghost points for wall flags 1647 CALL exchange_horiz_int( wall_flags_0, nbgp ) 1648 ! 1649 !-- Set boundary conditions also for flags. Can be interpreted as Neumann 1650 !-- boundary conditions for topography. 1651 IF ( .NOT. bc_ns_cyc ) THEN 1652 IF ( nys == 0 ) wall_flags_0(:,-1,:) = wall_flags_0(:,0,:) 1653 IF ( nyn == ny ) wall_flags_0(:,ny+1,:) = wall_flags_0(:,ny,:) 1654 ENDIF 1655 IF ( .NOT. bc_lr_cyc ) THEN 1656 IF ( nxl == 0 ) wall_flags_0(:,:,-1) = wall_flags_0(:,:,0) 1657 IF ( nxr == nx ) wall_flags_0(:,:,nx+1) = wall_flags_0(:,:,nx) 1658 ENDIF 1659 1660 ! 1661 !-- Initialize boundary conditions via surface type 1662 CALL init_bc 1663 ! 1664 !-- Allocate and set topography height arrays required for data output 1665 IF ( TRIM( topography ) /= 'flat' ) THEN 1666 ! 1667 !-- Allocate and set the arrays containing the topography height 1668 1669 IF ( nxr == nx .AND. nyn /= ny ) THEN 1670 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn), & 1671 zw_w_inner(nxl:nxr+1,nys:nyn) ) 1672 ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN 1673 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1), & 1674 zw_w_inner(nxl:nxr,nys:nyn+1) ) 1675 ELSEIF ( nxr == nx .AND. nyn == ny ) THEN 1676 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1), & 1677 zw_w_inner(nxl:nxr+1,nys:nyn+1) ) 1678 ELSE 1679 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn), & 1680 zw_w_inner(nxl:nxr,nys:nyn) ) 1681 ENDIF 1682 1683 zu_s_inner = 0.0_wp 1684 zw_w_inner = 0.0_wp 1685 ! 1686 !-- Determine local topography height on scalar and w-grid. Note, setting 1687 !-- lateral boundary values is not necessary, realized via wall_flags_0 1688 !-- array. Further, please note that loop bounds are different from 1689 !-- nxl to nxr and nys to nyn on south and right model boundary, hence, 1690 !-- use intrinsic lbound and ubound functions to infer array bounds. 1691 DO i = lbound(zu_s_inner, 1), ubound(zu_s_inner, 1) 1692 DO j = lbound(zu_s_inner, 2), ubound(zu_s_inner, 2) 1693 ! 1694 !-- Topography height on scalar grid. Therefore, determine index of 1695 !-- upward-facing surface element on scalar grid (bit 12). 1696 zu_s_inner(i,j) = zu( MAXLOC( MERGE( & 1697 1, 0, BTEST( wall_flags_0(:,j,i), 12 )& 1698 ), DIM = 1 & 1699 ) - 1 & 1700 ) 1701 ! 1702 !-- Topography height on w grid. Therefore, determine index of 1703 !-- upward-facing surface element on w grid (bit 18). 1704 zw_w_inner(i,j) = zw( MAXLOC( MERGE( & 1705 1, 0, BTEST( wall_flags_0(:,j,i), 18 )& 1706 ), DIM = 1 & 1707 ) - 1 & 1708 ) 1709 ENDDO 1710 ENDDO 1711 1712 1713 ENDIF 1714 1290 1715 ! 1291 1716 !-- Calculate wall flag arrays for the multigrid method. … … 1450 1875 !-- required in the ws-scheme, the arrays need to be allocated here as they are 1451 1876 !-- used in OpenACC directives. 1452 ALLOCATE( wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &1453 wall_flags_00(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )1454 wall_flags_0= 01455 wall_flags_00= 01877 ALLOCATE( advc_flags_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 1878 advc_flags_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1879 advc_flags_1 = 0 1880 advc_flags_2 = 0 1456 1881 ! 1457 1882 !-- Init flags for ws-scheme to degrade order of the numerics near walls, i.e. … … 1466 1891 !-- Go through all points of the subdomain one by one and look for the closest 1467 1892 !-- surface 1468 IF ( TRIM(topography) /= 'flat' ) THEN 1469 DO i = nxl, nxr 1470 DO j = nys, nyn 1471 1472 nzb_si = nzb_s_inner(j,i) 1473 vi = vertical_influence(nzb_si) 1474 1475 IF ( wall_n(j,i) > 0 ) THEN 1476 ! 1477 !-- North wall (y distance) 1478 DO k = wall_n(j,i), nzb_si 1479 l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), 0.5_wp * dy ) 1480 ENDDO 1481 ! 1482 !-- Above North wall (yz distance) 1483 DO k = nzb_si + 1, nzb_si + vi 1484 l_wall(k,j+1,i) = MIN( l_wall(k,j+1,i), & 1485 SQRT( 0.25_wp * dy**2 + & 1486 ( zu(k) - zw(nzb_si) )**2 ) ) 1487 ENDDO 1488 ! 1489 !-- Northleft corner (xy distance) 1490 IF ( corner_nl(j,i) > 0 ) THEN 1491 DO k = corner_nl(j,i), nzb_si 1492 l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), & 1493 0.5_wp * SQRT( dx**2 + dy**2 ) ) 1494 ENDDO 1495 ! 1496 !-- Above Northleft corner (xyz distance) 1497 DO k = nzb_si + 1, nzb_si + vi 1498 l_wall(k,j+1,i-1) = MIN( l_wall(k,j+1,i-1), & 1499 SQRT( 0.25_wp * (dx**2 + dy**2) + & 1500 ( zu(k) - zw(nzb_si) )**2 ) ) 1501 ENDDO 1502 ENDIF 1503 ! 1504 !-- Northright corner (xy distance) 1505 IF ( corner_nr(j,i) > 0 ) THEN 1506 DO k = corner_nr(j,i), nzb_si 1507 l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), & 1508 0.5_wp * SQRT( dx**2 + dy**2 ) ) 1509 ENDDO 1510 ! 1511 !-- Above northright corner (xyz distance) 1512 DO k = nzb_si + 1, nzb_si + vi 1513 l_wall(k,j+1,i+1) = MIN( l_wall(k,j+1,i+1), & 1514 SQRT( 0.25_wp * (dx**2 + dy**2) + & 1515 ( zu(k) - zw(nzb_si) )**2 ) ) 1516 ENDDO 1517 ENDIF 1893 DO i = nxl, nxr 1894 DO j = nys, nyn 1895 DO k = nzb+1, nzt 1896 ! 1897 !-- Check if current gridpoint belongs to the atmosphere 1898 IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN 1899 ! 1900 !-- Check for neighbouring grid-points. 1901 !-- Vertical distance, down 1902 IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) & 1903 l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) ) 1904 ! 1905 !-- Vertical distance, up 1906 IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) & 1907 l_wall(k,j,i) = MIN( l_grid(k), zw(k) - zu(k) ) 1908 ! 1909 !-- y-distance 1910 IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) .OR. & 1911 .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) & 1912 l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy ) 1913 ! 1914 !-- x-distance 1915 IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) .OR. & 1916 .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) & 1917 l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx ) 1918 ! 1919 !-- yz-distance (vertical edges, down) 1920 IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i), 0 ) .OR. & 1921 .NOT. BTEST( wall_flags_0(k-1,j+1,i), 0 ) ) & 1922 l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 1923 SQRT( 0.25_wp * dy**2 + & 1924 ( zu(k) - zw(k-1) )**2 ) ) 1925 ! 1926 !-- yz-distance (vertical edges, up) 1927 IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i), 0 ) .OR. & 1928 .NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 ) ) & 1929 l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 1930 SQRT( 0.25_wp * dy**2 + & 1931 ( zw(k) - zu(k) )**2 ) ) 1932 ! 1933 !-- xz-distance (vertical edges, down) 1934 IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 ) .OR. & 1935 .NOT. BTEST( wall_flags_0(k-1,j,i+1), 0 ) ) & 1936 l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 1937 SQRT( 0.25_wp * dx**2 + & 1938 ( zu(k) - zw(k-1) )**2 ) ) 1939 ! 1940 !-- xz-distance (vertical edges, up) 1941 IF ( .NOT. BTEST( wall_flags_0(k+1,j,i-1), 0 ) .OR. & 1942 .NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 ) ) & 1943 l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 1944 SQRT( 0.25_wp * dx**2 + & 1945 ( zw(k) - zu(k) )**2 ) ) 1946 ! 1947 !-- xy-distance (horizontal edges) 1948 IF ( .NOT. BTEST( wall_flags_0(k,j-1,i-1), 0 ) .OR. & 1949 .NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 ) .OR. & 1950 .NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 ) .OR. & 1951 .NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) ) & 1952 l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 1953 SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) ) 1954 ! 1955 !-- xyz distance (vertical and horizontal edges, down) 1956 IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i-1), 0 ) .OR. & 1957 .NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 ) .OR. & 1958 .NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 ) .OR. & 1959 .NOT. BTEST( wall_flags_0(k-1,j+1,i+1), 0 ) ) & 1960 l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 1961 SQRT( 0.25_wp * ( dx**2 + dy**2 ) & 1962 + ( zu(k) - zw(k-1) )**2 ) ) 1963 ! 1964 !-- xyz distance (vertical and horizontal edges, up) 1965 IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i-1), 0 ) .OR. & 1966 .NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 ) .OR. & 1967 .NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 ) .OR. & 1968 .NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) ) & 1969 l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 1970 SQRT( 0.25_wp * ( dx**2 + dy**2 ) & 1971 + ( zw(k) - zu(k) )**2 ) ) 1972 1518 1973 ENDIF 1519 1520 IF ( wall_s(j,i) > 0 ) THEN1521 !1522 !-- South wall (y distance)1523 DO k = wall_s(j,i), nzb_si1524 l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), 0.5_wp * dy )1525 ENDDO1526 !1527 !-- Above south wall (yz distance)1528 DO k = nzb_si + 1, nzb_si + vi1529 l_wall(k,j-1,i) = MIN( l_wall(k,j-1,i), &1530 SQRT( 0.25_wp * dy**2 + &1531 ( zu(k) - zw(nzb_si) )**2 ) )1532 ENDDO1533 !1534 !-- Southleft corner (xy distance)1535 IF ( corner_sl(j,i) > 0 ) THEN1536 DO k = corner_sl(j,i), nzb_si1537 l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &1538 0.5_wp * SQRT( dx**2 + dy**2 ) )1539 ENDDO1540 !1541 !-- Above southleft corner (xyz distance)1542 DO k = nzb_si + 1, nzb_si + vi1543 l_wall(k,j-1,i-1) = MIN( l_wall(k,j-1,i-1), &1544 SQRT( 0.25_wp * (dx**2 + dy**2) + &1545 ( zu(k) - zw(nzb_si) )**2 ) )1546 ENDDO1547 ENDIF1548 !1549 !-- Southright corner (xy distance)1550 IF ( corner_sr(j,i) > 0 ) THEN1551 DO k = corner_sr(j,i), nzb_si1552 l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &1553 0.5_wp * SQRT( dx**2 + dy**2 ) )1554 ENDDO1555 !1556 !-- Above southright corner (xyz distance)1557 DO k = nzb_si + 1, nzb_si + vi1558 l_wall(k,j-1,i+1) = MIN( l_wall(k,j-1,i+1), &1559 SQRT( 0.25_wp * (dx**2 + dy**2) + &1560 ( zu(k) - zw(nzb_si) )**2 ) )1561 ENDDO1562 ENDIF1563 1564 ENDIF1565 1566 IF ( wall_l(j,i) > 0 ) THEN1567 !1568 !-- Left wall (x distance)1569 DO k = wall_l(j,i), nzb_si1570 l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), 0.5_wp * dx )1571 ENDDO1572 !1573 !-- Above left wall (xz distance)1574 DO k = nzb_si + 1, nzb_si + vi1575 l_wall(k,j,i-1) = MIN( l_wall(k,j,i-1), &1576 SQRT( 0.25_wp * dx**2 + &1577 ( zu(k) - zw(nzb_si) )**2 ) )1578 ENDDO1579 ENDIF1580 1581 IF ( wall_r(j,i) > 0 ) THEN1582 !1583 !-- Right wall (x distance)1584 DO k = wall_r(j,i), nzb_si1585 l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), 0.5_wp * dx )1586 ENDDO1587 !1588 !-- Above right wall (xz distance)1589 DO k = nzb_si + 1, nzb_si + vi1590 l_wall(k,j,i+1) = MIN( l_wall(k,j,i+1), &1591 SQRT( 0.25_wp * dx**2 + &1592 ( zu(k) - zw(nzb_si) )**2 ) )1593 ENDDO1594 1595 ENDIF1596 1597 1974 ENDDO 1598 1975 ENDDO 1599 1600 ENDIF 1601 1602 ! 1603 !-- Multiplication with wall_adjustment_factor 1604 l_wall = wall_adjustment_factor * l_wall 1605 1976 ENDDO 1606 1977 ! 1607 1978 !-- Set lateral boundary conditions for l_wall 1608 CALL exchange_horiz( l_wall, nbgp ) 1609 1610 DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, & 1611 vertical_influence, wall_l, wall_n, wall_r, wall_s ) 1979 CALL exchange_horiz( l_wall, nbgp ) 1612 1980 1613 1981
Note: See TracChangeset
for help on using the changeset viewer.
