Changeset 4868 for palm/trunk/SOURCE
- Timestamp:
- Feb 8, 2021 10:07:43 AM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_grid.f90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_grid.f90
r4828 r4868 24 24 ! ----------------- 25 25 ! $Id$ 26 ! Height of level k=0 for the u,v-grid is always set 0.0, 27 ! small changes to follow coding standard 28 ! 29 ! 4828 2021-01-05 11:21:41Z Giersch 26 30 ! Bugfix in building mapping when static driver is available but no actual building is present 27 31 ! within the model domain … … 161 165 dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end, & 162 166 dz_stretch_level_end_index, dz_stretch_level_start_index, & 163 dz_stretch_level_start, ibc_uv_b, message_string,&167 dz_stretch_level_start, message_string, & 164 168 number_stretch_level_end, & 165 169 number_stretch_level_start, & … … 386 390 387 391 ! 388 !-- Grid for atmosphere with surface at z=0 (k=0, w-grid). 389 !-- First compute the u- and v-levels. In case of dirichlet bc for u and v the first u/v- and 390 !-- w-level (k=0) are defined at same height (z=0). 391 !-- The second u-level (k=1) corresponds to the top of the surface layer. In case of symmetric 392 !-- boundaries (closed channel flow), the first grid point is always at z=0. 393 IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 .OR. topography == 'closed_channel' ) THEN 394 zu(0) = 0.0_wp 395 ELSE 396 zu(0) = - dz(1) * 0.5_wp 397 ENDIF 398 399 zu(1) = dz(1) * 0.5_wp 392 !-- Grid for atmosphere with surface at z=0. This corresponds to k=0 on the w-grid AND 393 !-- the u-grid. 394 !-- First compute the u- and v-levels. 395 !-- The second u-level (k=1) corresponds to the top of the surface layer. 396 zu(0) = 0.0_wp 397 zu(1) = dz(1) * 0.5_wp 400 398 401 399 ! … … 430 428 dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) ) 431 429 432 IF ( dz_level_end < dz(n+1)/3.0 ) THEN430 IF ( dz_level_end < dz(n+1) / 3.0 ) THEN 433 431 zu(k) = dz_stretch_level_end(n) 434 432 dz_stretched = dz(n+1) … … 572 570 ! 573 571 !-- Compute the w-levels. They are always staggered half-way between the corresponding u-levels, 574 !-- except in case of dirichlet bcfor u and v at the ground. In this case the first u- and572 !-- except for u and v at the ground. In this case the first u- and 575 573 !-- w-level are defined at same height. The top w-level (nzt+1) is not used but set for 576 !-- consistency, since w and all scalar variables are defined up t pnzt+1.574 !-- consistency, since w and all scalar variables are defined up to nzt+1. 577 575 zw(nzt+1) = dz(1) 578 576 zw(nzt) = 0.0_wp … … 582 580 583 581 ! 584 !-- In case of dirichlet bc for u and v the first u- and w-level are defined at same height. 585 IF ( ibc_uv_b == 0 ) THEN 586 zu(0) = zw(0) 587 ENDIF 582 !-- Like for the atmosphere, the first u-grid and w-grid-level are defined at same height. 583 zu(0) = zw(0) 588 584 589 585 ENDIF !End of defining the vertical grid levels … … 822 818 !-- Second branch for stretching from fine to rough 823 819 ELSEIF ( dz(n) < dz(n+1) ) THEN 820 824 821 DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit ) 825 822 … … 853 850 854 851 ELSE 852 855 853 message_string= 'Two adjacent values of dz must be different' 856 854 CALL message( 'init_grid', 'PA0228', 1, 2, 0, 6, 0 ) … … 872 870 873 871 ENDIF 872 874 873 ENDDO 875 874 … … 1280 1279 !-- Topography input via ASCII format. 1281 1280 ELSE 1282 ocean_offset = MERGE( zw(0), 0.0_wp, ocean_mode ) 1281 1282 ocean_offset = MERGE( zw(0), 0.0_wp, ocean_mode ) 1283 1283 ! 1284 1284 !-- Initialize topography bit 0 (indicates obstacle) everywhere to zero and clear all grid points … … 1305 1305 ENDDO 1306 1306 ENDDO 1307 1307 1308 ENDIF 1308 1309 … … 1353 1354 INTEGER(iwp) :: num_wall !< number of surrounding vertical walls for a single grid point 1354 1355 1355 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo_tmp 1356 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo_3d 1357 1356 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo_tmp !< temporary 3D-topography used to fill holes 1357 INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: topo_3d !< 3D-topography array merging buildings and 1358 !< orography 1358 1359 1359 1360 LOGICAL :: filled = .FALSE. !< flag indicating if holes were filled … … 1790 1791 1791 1792 CASE ( 'tunnel' ) 1792 1793 1793 ! 1794 1794 !-- Tunnel height … … 2075 2075 ENDIF 2076 2076 2077 2078 2077 END SUBROUTINE init_topo 2078 2079 2079 2080 2080 SUBROUTINE set_topo_flags(topo) … … 2126 2126 ENDDO 2127 2127 2128 DO k = nzb, nzt2128 DO k = nzb, nzt 2129 2129 ! 2130 2130 !-- w grid … … 2145 2145 !-- Set outer array for scalars to mask near-surface grid points. Note, on basis of flag 24 futher 2146 2146 !-- flags will be derived which are used to control production of subgrid TKE production near walls. 2147 2148 2147 ALLOCATE( wall_flags_total_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 2149 2148 wall_flags_total_0 = 0 2150 2149 2151 DO i = nxl, nxr2152 DO j = nys, nyn2153 DO k = nzb, nzt+12150 DO i = nxl, nxr 2151 DO j = nys, nyn 2152 DO k = nzb, nzt+1 2154 2153 wall_flags_total_0(k,j,i) = IOR( wall_flags_total_0(k,j,i), wall_flags_static_0(k,j,i) ) 2155 2154 ENDDO … … 2159 2158 CALL exchange_horiz_int( wall_flags_total_0, nys, nyn, nxl, nxr, nzt, nbgp ) 2160 2159 2161 DO i = nxl, nxr2162 DO j = nys, nyn2163 DO k = nzb, nzt+12160 DO i = nxl, nxr 2161 DO j = nys, nyn 2162 DO k = nzb, nzt+1 2164 2163 IF ( BTEST( wall_flags_total_0(k,j-1,i), 0 ) .AND. & 2165 2164 BTEST( wall_flags_total_0(k,j+1,i), 0 ) .AND. & … … 2176 2175 ! 2177 2176 !-- Set further special flags 2178 DO i = nxl, nxr2179 DO j = nys, nyn2180 DO k = nzb, nzt+12177 DO i = nxl, nxr 2178 DO j = nys, nyn 2179 DO k = nzb, nzt+1 2181 2180 ! 2182 2181 !-- scalar grid, former nzb_diff_s_inner. … … 2205 2204 2206 2205 2207 DO k = nzb+1, nzt2206 DO k = nzb+1, nzt 2208 2207 ! 2209 2208 !-- Special flag on u grid, former nzb_u_inner + 1, required for disturb_field and … … 2254 2253 ! 2255 2254 !-- Flags indicating downward facing walls 2256 DO k = nzb+1, nzt+12255 DO k = nzb+1, nzt+1 2257 2256 ! 2258 2257 !-- Scalar grid … … 2278 2277 ! 2279 2278 !-- Flags indicating upward facing walls 2280 DO k = nzb, nzt2279 DO k = nzb, nzt 2281 2280 ! 2282 2281 !-- Upward facing wall on scalar grid … … 2331 2330 !-- species or aerosols. 2332 2331 IF ( scalar_advec == 'ws-scheme' ) THEN 2333 DO k = nzb, nzt2332 DO k = nzb, nzt 2334 2333 IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) .AND. ( & 2335 2334 ANY( .NOT. BTEST( wall_flags_total_0(k,j-3:j+3,i-1), 0 ) ) .OR. & … … 2364 2363 wall_flags_static_0(nzb,:,:) = IBSET( wall_flags_static_0(nzb,:,:), 5 ) 2365 2364 ELSE 2366 DO i = nxl, nxr2367 DO j = nys, nyn2368 DO k = nzb, nzt+12365 DO i = nxl, nxr 2366 DO j = nys, nyn 2367 DO k = nzb, nzt+1 2369 2368 ! 2370 2369 !-- Natural terrain grid point … … 2377 2376 ! 2378 2377 !-- Building grid points. 2379 DO i = nxl, nxr2380 DO j = nys, nyn2381 DO k = nzb, nzt+12378 DO i = nxl, nxr 2379 DO j = nys, nyn 2380 DO k = nzb, nzt+1 2382 2381 IF ( BTEST( topo(k,j,i), 2 ) ) & 2383 2382 wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 6 ) … … 2387 2386 ! 2388 2387 !-- Set flag 4, indicating new topography grid points due to filtering. 2389 DO i = nxl, nxr2390 DO j = nys, nyn2391 DO k = nzb, nzt+12388 DO i = nxl, nxr 2389 DO j = nys, nyn 2390 DO k = nzb, nzt+1 2392 2391 IF ( BTEST( topo(k,j,i), 4 ) ) & 2393 2392 wall_flags_static_0(k,j,i) = IBSET( wall_flags_static_0(k,j,i), 4 ) … … 2398 2397 CALL exchange_horiz_int( wall_flags_static_0, nys, nyn, nxl, nxr, nzt, nbgp ) 2399 2398 2400 DO i = nxl, nxr2401 DO j = nys, nyn2402 DO k = nzb, nzt+12399 DO i = nxl, nxr 2400 DO j = nys, nyn 2401 DO k = nzb, nzt+1 2403 2402 wall_flags_total_0(k,j,i) = IOR( wall_flags_total_0(k,j,i), wall_flags_static_0(k,j,i) ) 2404 2403 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
