- Timestamp:
- Jan 27, 2020 3:07:30 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_grid.f90 (modified) (36 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_grid.f90
r4360 r4386 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Allocation statements, comments, naming of variables revised and _wp added to 28 ! real type values 29 ! 30 ! 4360 2020-01-07 11:25:50Z suehring 27 31 ! Revise error messages for generic tunnel setup. 28 32 ! … … 199 203 ! 200 204 !-- Allocate grid arrays 201 ALLOCATE( x(0:nx), xu(0:nx) ) 205 ALLOCATE( x(0:nx) ) 206 ALLOCATE( xu(0:nx) ) 207 202 208 DO i = 0, nx 203 209 xu(i) = i * dx … … 205 211 ENDDO 206 212 207 ALLOCATE( y(0:ny), yv(0:ny) ) 213 ALLOCATE( y(0:ny) ) 214 ALLOCATE( yv(0:ny) ) 215 208 216 DO j = 0, ny 209 217 yv(j) = j * dy … … 211 219 ENDDO 212 220 213 ! 214 !-- Allocate grid arrays 215 ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), & 216 dzw(1:nzt+1), zu(nzb:nzt+1), zw(nzb:nzt+1) ) 217 218 ! 219 !-- Compute height of u-levels from constant grid length and dz stretch factors 221 ALLOCATE( ddzu(1:nzt+1) ) 222 ALLOCATE( ddzw(1:nzt+1) ) 223 ALLOCATE( dd2zu(1:nzt) ) 224 ALLOCATE( dzu(1:nzt+1) ) 225 ALLOCATE( dzw(1:nzt+1) ) 226 ALLOCATE( zu(nzb:nzt+1) ) 227 ALLOCATE( zw(nzb:nzt+1) ) 228 229 ! 230 !-- For constructing an appropriate grid, the vertical grid spacing dz has to 231 !-- be specified with a non-negative value in the parameter file 220 232 IF ( dz(1) == -1.0_wp ) THEN 221 233 message_string = 'missing dz' … … 248 260 249 261 ! 250 !-- The number of specified end levels +1 has to be the same thanthe number262 !-- The number of specified end levels +1 has to be the same as the number 251 263 !-- of specified dz values 252 264 IF ( number_dz /= number_stretch_level_end + 1 ) THEN 253 WRITE( message_string, * ) 'The number of values for dz = ', &254 number_dz, 'has to be the same than& ',&255 'the number of values for ', &256 'dz_stretch_level_end + 1 = ', &265 WRITE( message_string, * ) 'The number of values for dz = ', & 266 number_dz, 'has to be the same as& ', & 267 'the number of values for ', & 268 'dz_stretch_level_end + 1 = ', & 257 269 number_stretch_level_end+1 258 270 CALL message( 'init_grid', 'PA0156', 1, 2, 0, 6, 0 ) … … 260 272 261 273 ! 262 !-- 263 !-- 264 IF ( number_dz /= number_stretch_level_start + 1 .AND. &274 !-- The number of specified start levels has to be the same or one less than 275 !-- the number of specified dz values 276 IF ( number_dz /= number_stretch_level_start + 1 .AND. & 265 277 number_dz /= number_stretch_level_start ) THEN 266 WRITE( message_string, * ) 'The number of values for dz = ', &267 number_dz, 'has to be the same or one ', &268 'more than& the number of values for ', &269 'dz_stretch_level_start = ', &278 WRITE( message_string, * ) 'The number of values for dz = ', & 279 number_dz, 'has to be the same as or one ', & 280 'more than& the number of values for ', & 281 'dz_stretch_level_start = ', & 270 282 number_stretch_level_start 271 283 CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 ) 272 284 ENDIF 273 285 274 !-- 275 !-- 276 IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND. &286 !-- The number of specified start levels has to be the same or one more than 287 !-- the number of specified end levels 288 IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND. & 277 289 number_stretch_level_start /= number_stretch_level_end ) THEN 278 WRITE( message_string, * ) 'The number of values for ', &279 'dz_stretch_level_start = ', &280 dz_stretch_level_start, 'has to be the ', &281 'same or one more than& the number of ', &282 'values for dz_stretch_level_end = ', &290 WRITE( message_string, * ) 'The number of values for ', & 291 'dz_stretch_level_start = ', & 292 dz_stretch_level_start, 'has to be the ', & 293 'same or one more than& the number of ', & 294 'values for dz_stretch_level_end = ', & 283 295 number_stretch_level_end 284 296 CALL message( 'init_grid', 'PA0216', 1, 2, 0, 6, 0 ) … … 306 318 307 319 ! 308 !-- Define the vertical grid levels 320 !-- Define the vertical grid levels. Start with atmosphere branch 309 321 IF ( .NOT. ocean_mode ) THEN 310 322 311 323 ! 312 324 !-- The stretching region has to be large enough to allow for a smooth 313 !-- transition between two different grid spacings 325 !-- transition between two different grid spacings. The number 4 is an 326 !-- empirical value 314 327 DO n = 1, number_stretch_level_start 315 328 min_dz_stretch_level_end(n) = dz_stretch_level_start(n) + & … … 319 332 IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) > & 320 333 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN 321 message_string= 'E each dz_stretch_level_end has to be larger ' //&334 message_string= 'Each dz_stretch_level_end has to be larger ' // & 322 335 'than its corresponding value for &' // & 323 336 'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//& … … 327 340 328 341 ! 329 !-- Stretching must not be applied within the prandtl_layer342 !-- Stretching must not be applied within the surface layer 330 343 !-- (first two grid points). For the default case dz_stretch_level_start 331 344 !-- is negative. Therefore the absolut value is checked here. 332 IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_wp ) ) THEN333 WRITE( message_string, * ) 'E each dz_stretch_level_start has to be ',&345 IF ( ANY( ABS( dz_stretch_level_start ) <= dz(1) * 1.5_wp ) ) THEN 346 WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',& 334 347 'larger than ', dz(1) * 1.5 335 348 CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 ) … … 338 351 ! 339 352 !-- The stretching has to start and end on a grid level. Therefore 340 !-- user-specified values have to ''interpolate'' to the next lowest level 353 !-- user-specified values are mapped to the next lowest level. The 354 !-- calculation of the first level is realized differently just because of 355 !-- historical reasons (the advanced/new stretching mechanism was realized 356 !-- in a way that results don't change if the old parameters 357 !-- dz_stretch_level, dz_stretch_factor and dz_max are used) 341 358 IF ( number_stretch_level_start /= 0 ) THEN 342 359 dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) - & … … 358 375 ENDDO 359 376 ENDIF 360 377 361 378 ! 362 379 !-- Determine stretching factor if necessary … … 370 387 !-- the first u/v- and w-level (k=0) are defined at same height (z=0). 371 388 !-- The second u-level (k=1) corresponds to the top of the 372 !-- Prandtl-layer. In case of symmetric boundaries (closed channel flow),389 !-- surface layer. In case of symmetric boundaries (closed channel flow), 373 390 !-- the first grid point is always at z=0. 374 391 IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 .OR. & … … 390 407 391 408 !-- The default value of dz_stretch_level_start is negative, thus the first 392 !-- condition is always true. Hence, the second condition is necessary. 409 !-- condition is true even if no stretching shall be applied. Hence, the 410 !-- second condition is also necessary. 393 411 DO k = 2, nzt+1-symmetry_flag 394 412 IF ( dz_stretch_level_start(n) <= zu(k-1) .AND. & … … 452 470 ENDIF 453 471 454 ELSE 472 ELSE !ocean branch 455 473 456 474 ! 457 475 !-- The stretching region has to be large enough to allow for a smooth 458 !-- transition between two different grid spacings 476 !-- transition between two different grid spacings. The number 4 is an 477 !-- empirical value 459 478 DO n = 1, number_stretch_level_start 460 479 min_dz_stretch_level_end(n) = dz_stretch_level_start(n) - & … … 464 483 IF ( ANY( min_dz_stretch_level_end (1:number_stretch_level_start) < & 465 484 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN 466 message_string= 'E each dz_stretch_level_end has to be less ' // &485 message_string= 'Each dz_stretch_level_end has to be less ' // & 467 486 'than its corresponding value for &' // & 468 487 'dz_stretch_level_start - 4*MAX(dz(n),dz(n+1)) '//& … … 474 493 !-- Stretching must not be applied close to the surface (last two grid 475 494 !-- points). For the default case dz_stretch_level_start is negative. 476 IF ( ANY( dz_stretch_level_start > - dz(1) * 1.5_wp ) ) THEN477 WRITE( message_string, * ) 'E each dz_stretch_level_start has to be ',&478 'less than ', dz(1) * 1.5495 IF ( ANY( dz_stretch_level_start >= - dz(1) * 1.5_wp ) ) THEN 496 WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',& 497 'less than ', -dz(1) * 1.5 479 498 CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 ) 480 499 ENDIF … … 482 501 ! 483 502 !-- The stretching has to start and end on a grid level. Therefore 484 !-- user-specified values have to ''interpolate'' to the next highest level 503 !-- user-specified values are mapped to the next highest level. The 504 !-- calculation of the first level is realized differently just because of 505 !-- historical reasons (the advanced/new stretching mechanism was realized 506 !-- in a way that results don't change if the old parameters 507 !-- dz_stretch_level, dz_stretch_factor and dz_max are used) 485 508 IF ( number_stretch_level_start /= 0 ) THEN 486 509 dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) + & … … 514 537 !-- below the first w-level (k=0). In case of dirichlet bc the first u- and 515 538 !-- w-level are defined at same height, but staggered from the second level. 516 !-- The second u-level (k=1) corresponds to the top of the Prandtl-layer.539 !-- The second u-level (k=1) corresponds to the top of the surface layer. 517 540 !-- z values are negative starting from z=0 (surface) 518 541 zu(nzt+1) = dz(1) * 0.5_wp … … 576 599 ENDIF 577 600 578 ENDIF 601 ENDIF !End of defining the vertical grid levels 579 602 580 603 ! … … 617 640 topo = 0 618 641 ! 619 !-- Initialize topography by generic topography or read fromtopography from file.642 !-- Initialize topography by generic topography or read topography from file. 620 643 CALL init_topo( topo ) 621 644 ! … … 691 714 !-- Allocate and set the arrays containing the topography height 692 715 IF ( nxr == nx .AND. nyn /= ny ) THEN 693 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn) , &694 716 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn) ) 717 ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn) ) 695 718 ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN 696 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1) , &697 719 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1) ) 720 ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn+1) ) 698 721 ELSEIF ( nxr == nx .AND. nyn == ny ) THEN 699 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1) , &700 722 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1) ) 723 ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn+1) ) 701 724 ELSE 702 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn) , &703 725 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn) ) 726 ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn) ) 704 727 ENDIF 705 728 … … 731 754 !-- Allocate outer and inner index arrays for topography and set 732 755 !-- defaults. 733 ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg) , &734 nzb_s_outer(nysg:nyng,nxlg:nxrg), &735 nzb_u_inner(nysg:nyng,nxlg:nxrg), &736 nzb_u_outer(nysg:nyng,nxlg:nxrg), &737 nzb_v_inner(nysg:nyng,nxlg:nxrg), &738 nzb_v_outer(nysg:nyng,nxlg:nxrg), &739 nzb_w_inner(nysg:nyng,nxlg:nxrg), &740 nzb_w_outer(nysg:nyng,nxlg:nxrg), &741 nzb_diff_s_inner(nysg:nyng,nxlg:nxrg), &742 nzb_diff_s_outer(nysg:nyng,nxlg:nxrg), &743 nzb_local(nysg:nyng,nxlg:nxrg), &744 756 ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg) ) 757 ALLOCATE( nzb_s_outer(nysg:nyng,nxlg:nxrg) ) 758 ALLOCATE( nzb_u_inner(nysg:nyng,nxlg:nxrg) ) 759 ALLOCATE( nzb_u_outer(nysg:nyng,nxlg:nxrg) ) 760 ALLOCATE( nzb_v_inner(nysg:nyng,nxlg:nxrg) ) 761 ALLOCATE( nzb_v_outer(nysg:nyng,nxlg:nxrg) ) 762 ALLOCATE( nzb_w_inner(nysg:nyng,nxlg:nxrg) ) 763 ALLOCATE( nzb_w_outer(nysg:nyng,nxlg:nxrg) ) 764 ALLOCATE( nzb_diff_s_inner(nysg:nyng,nxlg:nxrg) ) 765 ALLOCATE( nzb_diff_s_outer(nysg:nyng,nxlg:nxrg) ) 766 ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) ) 767 ALLOCATE( nzb_tmp(nysg:nyng,nxlg:nxrg) ) 745 768 ! 746 769 !-- Initialize 2D-index arrays. Note, these will be removed soon! … … 996 1019 REAL(wp) :: dz_stretch_factor_array_2(9) = 1.08_wp !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1 997 1020 998 REAL(wp), PARAMETER :: stretch_factor_interval = 1.0E-06 !< interval for sampling possible stretching factors999 REAL(wp), PARAMETER :: stretch_factor_lower_limit = 0.88 !< lowest possible stretching factor1000 REAL(wp), PARAMETER :: stretch_factor_upper_limit = 1.12 !< highest possible stretching factor1021 REAL(wp), PARAMETER :: stretch_factor_interval = 1.0E-06_wp !< interval for sampling possible stretching factors 1022 REAL(wp), PARAMETER :: stretch_factor_lower_limit = 0.88_wp !< lowest possible stretching factor 1023 REAL(wp), PARAMETER :: stretch_factor_upper_limit = 1.12_wp !< highest possible stretching factor 1001 1024 1002 1025 … … 1005 1028 1006 1029 iterations = 1 1007 stretch_factor_1 = 1.0 1008 stretch_factor_2 = 1.0 1009 delta_total_old = 1.0 1030 stretch_factor_1 = 1.0_wp 1031 stretch_factor_2 = 1.0_wp 1032 delta_total_old = 1.0_wp 1010 1033 1034 ! 1035 !-- First branch for stretching from rough to fine 1011 1036 IF ( dz(n) > dz(n+1) ) THEN 1012 1037 DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) 1013 1038 1014 stretch_factor_1 = 1.0 - iterations * stretch_factor_interval1015 distance = ABS( dz_stretch_level_end(n) - &1016 dz_stretch_level_start(n) ) 1017 numerator = distance*stretch_factor_1/dz(n) + &1039 stretch_factor_1 = 1.0_wp - iterations * stretch_factor_interval 1040 distance = ABS( dz_stretch_level_end(n) - & 1041 dz_stretch_level_start(n) ) 1042 numerator = distance*stretch_factor_1/dz(n) + & 1018 1043 stretch_factor_1 - distance/dz(n) 1019 1044 1020 IF ( numerator > 0.0 ) THEN1021 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0 1045 IF ( numerator > 0.0_wp ) THEN 1046 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp 1022 1047 l_rounded = NINT( l ) 1023 1048 delta_l = ABS( l_rounded - l ) / l … … 1026 1051 stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) ) 1027 1052 1028 delta_stretch_factor = ABS( stretch_factor_1 - &1029 stretch_factor_2 ) / &1053 delta_stretch_factor = ABS( stretch_factor_1 - & 1054 stretch_factor_2 ) / & 1030 1055 stretch_factor_2 1031 1056 … … 1033 1058 1034 1059 ! 1035 !-- 1036 !-- 1037 !-- 1038 !-- 1060 !-- stretch_factor_1 is taken to guarantee that the stretching 1061 !-- procedure ends as close as possible to dz_stretch_level_end. 1062 !-- stretch_factor_2 would guarantee that the stretched dz(n) is 1063 !-- equal to dz(n+1) after l_rounded grid levels. 1039 1064 IF (delta_total_new < delta_total_old) THEN 1040 1065 dz_stretch_factor_array(n) = stretch_factor_1 … … 1046 1071 1047 1072 ENDDO 1048 1073 1074 ! 1075 !-- Second branch for stretching from fine to rough 1049 1076 ELSEIF ( dz(n) < dz(n+1) ) THEN 1050 1077 DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit ) 1051 1078 1052 stretch_factor_1 = 1.0 + iterations * stretch_factor_interval1053 distance = ABS( dz_stretch_level_end(n) - &1079 stretch_factor_1 = 1.0_wp + iterations * stretch_factor_interval 1080 distance = ABS( dz_stretch_level_end(n) - & 1054 1081 dz_stretch_level_start(n) ) 1055 numerator = distance*stretch_factor_1/dz(n) + &1082 numerator = distance*stretch_factor_1/dz(n) + & 1056 1083 stretch_factor_1 - distance/dz(n) 1057 1084 1058 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0 1085 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp 1059 1086 l_rounded = NINT( l ) 1060 1087 delta_l = ABS( l_rounded - l ) / l … … 1062 1089 stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) ) 1063 1090 1064 delta_stretch_factor = ABS( stretch_factor_1 - &1065 stretch_factor_2 ) / &1091 delta_stretch_factor = ABS( stretch_factor_1 - & 1092 stretch_factor_2 ) / & 1066 1093 stretch_factor_2 1067 1094 … … 1069 1096 1070 1097 ! 1071 !-- 1072 !-- 1073 !-- 1074 !-- 1098 !-- stretch_factor_1 is taken to guarantee that the stretching 1099 !-- procedure ends as close as possible to dz_stretch_level_end. 1100 !-- stretch_factor_2 would guarantee that the stretched dz(n) is 1101 !-- equal to dz(n+1) after l_rounded grid levels. 1075 1102 IF (delta_total_new < delta_total_old) THEN 1076 1103 dz_stretch_factor_array(n) = stretch_factor_1 … … 1794 1821 IMPLICIT NONE 1795 1822 1796 INTEGER(iwp) :: bh !< temporary vertical index of building height1797 INTEGER(iwp) :: blx!< grid point number of building size along x1798 INTEGER(iwp) :: bly!< grid point number of building size along y1799 INTEGER(iwp) :: bxl!< index for left building wall1800 INTEGER(iwp) :: bxr!< index for right building wall1801 INTEGER(iwp) :: byn!< index for north building wall1802 INTEGER(iwp) :: bys!< index for south building wall1803 INTEGER(iwp) :: ch !< temporary vertical index for canyon height1804 INTEGER(iwp) :: cwx!< grid point number of canyon size along x1805 INTEGER(iwp) :: cwy!< grid point number of canyon size along y1806 INTEGER(iwp) :: cxl!< index for left canyon wall1807 INTEGER(iwp) :: cxr!< index for right canyon wall1808 INTEGER(iwp) :: cyn!< index for north canyon wall1809 INTEGER(iwp) :: cys!< index for south canyon wall1810 INTEGER(iwp) :: i !< index variable along x1811 INTEGER(iwp) :: j !< index variable along y1812 INTEGER(iwp) :: k !< index variable along z1813 INTEGER(iwp) :: hv_in !< heavyside function to model inner tunnel surface1814 INTEGER(iwp) :: hv_out !< heavyside function to model outer tunnel surface1815 INTEGER(iwp) :: txe_out !< end position of outer tunnel wall in x1816 INTEGER(iwp) :: txs_out !< start position of outer tunnel wall in x1817 INTEGER(iwp) :: tye_out !< end position of outer tunnel wall in y1818 INTEGER(iwp) :: tys_out !< start position of outer tunnel wall in y1819 INTEGER(iwp) :: txe_in !< end position of inner tunnel wall in x1820 INTEGER(iwp) :: txs_in !< start position of inner tunnel wall in x1821 INTEGER(iwp) :: tye_in !< end position of inner tunnel wall in y1822 INTEGER(iwp) :: tys_in !< start position of inner tunnel wall in y1823 INTEGER(iwp) :: td !< tunnel wall depth1824 INTEGER(iwp) :: th !< height of outer tunnel wall1823 INTEGER(iwp) :: bh !< temporary vertical index of building height 1824 INTEGER(iwp) :: ngp_bx !< grid point number of building size along x 1825 INTEGER(iwp) :: ngp_by !< grid point number of building size along y 1826 INTEGER(iwp) :: index_left_bwall !< index for left building wall 1827 INTEGER(iwp) :: index_right_bwall !< index for right building wall 1828 INTEGER(iwp) :: index_north_bwall !< index for north building wall 1829 INTEGER(iwp) :: index_south_bwall !< index for south building wall 1830 INTEGER(iwp) :: ch !< temporary vertical index for canyon height 1831 INTEGER(iwp) :: ngp_cx !< grid point number of canyon size along x 1832 INTEGER(iwp) :: ngp_cy !< grid point number of canyon size along y 1833 INTEGER(iwp) :: index_left_cwall !< index for left canyon wall 1834 INTEGER(iwp) :: index_right_cwall !< index for right canyon wall 1835 INTEGER(iwp) :: index_north_cwall !< index for north canyon wall 1836 INTEGER(iwp) :: index_south_cwall !< index for south canyon wall 1837 INTEGER(iwp) :: i !< index variable along x 1838 INTEGER(iwp) :: j !< index variable along y 1839 INTEGER(iwp) :: k !< index variable along z 1840 INTEGER(iwp) :: hv_in !< heavyside function to model inner tunnel surface 1841 INTEGER(iwp) :: hv_out !< heavyside function to model outer tunnel surface 1842 INTEGER(iwp) :: txe_out !< end position of outer tunnel wall in x 1843 INTEGER(iwp) :: txs_out !< start position of outer tunnel wall in x 1844 INTEGER(iwp) :: tye_out !< end position of outer tunnel wall in y 1845 INTEGER(iwp) :: tys_out !< start position of outer tunnel wall in y 1846 INTEGER(iwp) :: txe_in !< end position of inner tunnel wall in x 1847 INTEGER(iwp) :: txs_in !< start position of inner tunnel wall in x 1848 INTEGER(iwp) :: tye_in !< end position of inner tunnel wall in y 1849 INTEGER(iwp) :: tys_in !< start position of inner tunnel wall in y 1850 INTEGER(iwp) :: td !< tunnel wall depth 1851 INTEGER(iwp) :: th !< height of outer tunnel wall 1825 1852 1826 1853 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography top at cell-center … … 1829 1856 !-- Check for correct setting of the namelist parameter topography. If 1830 1857 !-- topography information is read from file but topography = 'flat', 1831 !-- initialization does not properly.1858 !-- initialization does not work properly. 1832 1859 IF ( ( buildings_f%from_file .OR. terrain_height_f%from_file ) .AND. & 1833 1860 TRIM( topography ) /= 'read_from_file' ) THEN … … 1860 1887 !-- Single rectangular building, by default centered in the middle of the 1861 1888 !-- total domain 1862 blx = NINT( building_length_x / dx )1863 bly = NINT( building_length_y / dy )1889 ngp_bx = NINT( building_length_x / dx ) 1890 ngp_by = NINT( building_length_y / dy ) 1864 1891 bh = MINLOC( ABS( zw - building_height ), 1 ) - 1 1865 1892 IF ( ABS( zw(bh) - building_height ) == & 1866 1893 ABS( zw(bh+1) - building_height ) ) bh = bh + 1 1867 1894 IF ( building_wall_left == 9999999.9_wp ) THEN 1868 building_wall_left = ( nx + 1 - blx ) / 2 * dx1869 ENDIF 1870 bxl = NINT( building_wall_left / dx )1871 bxr = bxl + blx1895 building_wall_left = ( nx + 1 - ngp_bx ) / 2 * dx 1896 ENDIF 1897 index_left_bwall = NINT( building_wall_left / dx ) 1898 index_right_bwall = index_left_bwall + ngp_bx 1872 1899 1873 1900 IF ( building_wall_south == 9999999.9_wp ) THEN 1874 building_wall_south = ( ny + 1 - bly ) / 2 * dy1875 ENDIF 1876 bys= NINT( building_wall_south / dy )1877 byn = bys + bly1901 building_wall_south = ( ny + 1 - ngp_by ) / 2 * dy 1902 ENDIF 1903 index_south_bwall = NINT( building_wall_south / dy ) 1904 index_north_bwall = index_south_bwall + ngp_by 1878 1905 1879 1906 ! 1880 1907 !-- Building size has to meet some requirements 1881 IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR. & 1882 ( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) ) THEN 1908 IF ( ( index_left_bwall < 1 ) .OR. ( index_right_bwall > nx-1 ) .OR. & 1909 ( index_right_bwall < index_left_bwall+3 ) .OR. & 1910 ( index_south_bwall < 1 ) .OR. ( index_north_bwall > ny-1 ) .OR.& 1911 ( index_north_bwall < index_south_bwall+3 ) ) THEN 1883 1912 WRITE( message_string, * ) 'inconsistent building parameters:', & 1884 '&bxl=', bxl, 'bxr=', bxr, 'bys=', bys, & 1885 'byn=', byn, 'nx=', nx, 'ny=', ny 1913 '&index_left_bwall=', index_left_bwall, & 1914 'index_right_bwall=', index_right_bwall, & 1915 'index_south_bwall=', index_south_bwall, & 1916 'index_north_bwall=', index_north_bwall, & 1917 'nx=', nx, 'ny=', ny 1886 1918 CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 ) 1887 1919 ENDIF … … 1891 1923 ! 1892 1924 !-- Define the building. 1893 IF ( bxl <= nxr .AND. bxr >= nxl .AND. & 1894 bys <= nyn .AND. byn >= nys ) & 1895 nzb_local(MAX(nys,bys):MIN(nyn,byn),MAX(nxl,bxl):MIN(nxr,bxr)) = bh 1925 IF ( index_left_bwall <= nxr .AND. index_right_bwall >= nxl .AND. & 1926 index_south_bwall <= nyn .AND. index_north_bwall >= nys ) & 1927 nzb_local(MAX(nys,index_south_bwall):MIN(nyn,index_north_bwall), & 1928 MAX(nxl,index_left_bwall):MIN(nxr,index_right_bwall)) = bh 1896 1929 ! 1897 1930 !-- Set bit array on basis of nzb_local … … 1941 1974 ! 1942 1975 !-- Street canyon in y direction 1943 cwx = NINT( canyon_width_x / dx )1976 ngp_cx = NINT( canyon_width_x / dx ) 1944 1977 IF ( canyon_wall_left == 9999999.9_wp ) THEN 1945 canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx1946 ENDIF 1947 cxl= NINT( canyon_wall_left / dx )1948 cxr = cxl + cwx1978 canyon_wall_left = ( nx + 1 - ngp_cx ) / 2 * dx 1979 ENDIF 1980 index_left_cwall= NINT( canyon_wall_left / dx ) 1981 index_right_cwall= index_left_cwall+ ngp_cx 1949 1982 ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN 1950 1983 ! 1951 1984 !-- Street canyon in x direction 1952 cwy = NINT( canyon_width_y / dy )1985 ngp_cy = NINT( canyon_width_y / dy ) 1953 1986 IF ( canyon_wall_south == 9999999.9_wp ) THEN 1954 canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy1955 ENDIF 1956 cys= NINT( canyon_wall_south / dy )1957 cyn = cys + cwy1987 canyon_wall_south = ( ny + 1 - ngp_cy ) / 2 * dy 1988 ENDIF 1989 index_south_cwall = NINT( canyon_wall_south / dy ) 1990 index_north_cwall = index_south_cwall + ngp_cy 1958 1991 1959 1992 ELSE … … 1971 2004 !-- Street canyon size has to meet some requirements 1972 2005 IF ( canyon_width_x /= 9999999.9_wp ) THEN 1973 IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR. & 2006 IF ( ( index_left_cwall< 1 ) .OR. ( index_right_cwall> nx-1 ) .OR.& 2007 ( ngp_cx < 3 ) .OR. ( ch < 3 ) ) THEN 2008 WRITE( message_string, * ) 'inconsistent canyon parameters:', & 2009 '&index_left_cwall=', index_left_cwall, & 2010 ' index_right_cwall=', index_right_cwall, & 2011 ' ngp_cx=', ngp_cx, & 2012 ' ch=', ch, ' nx=', nx, ' ny=', ny 2013 CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 ) 2014 ENDIF 2015 ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN 2016 IF ( ( index_south_cwall < 1 ) .OR. & 2017 ( index_north_cwall > ny-1 ) .OR. ( ngp_cy < 3 ) .OR. & 1974 2018 ( ch < 3 ) ) THEN 1975 2019 WRITE( message_string, * ) 'inconsistent canyon parameters:', & 1976 '&cxl=', cxl, ' cxr=', cxr, & 1977 ' cwx=', cwx, & 1978 ' ch=', ch, ' nx=', nx, ' ny=', ny 1979 CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 ) 1980 ENDIF 1981 ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN 1982 IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR. & 1983 ( ch < 3 ) ) THEN 1984 WRITE( message_string, * ) 'inconsistent canyon parameters:', & 1985 '&cys=', cys, ' cyn=', cyn, & 1986 ' cwy=', cwy, & 2020 '&index_south_cwall=', index_south_cwall, & 2021 ' index_north_cwall=', index_north_cwall, & 2022 ' ngp_cy=', ngp_cy, & 1987 2023 ' ch=', ch, ' nx=', nx, ' ny=', ny 1988 2024 CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 ) … … 2000 2036 nzb_local = ch 2001 2037 IF ( canyon_width_x /= 9999999.9_wp ) THEN 2002 IF ( cxl <= nxr .AND. cxr >= nxl )&2003 nzb_local(:,MAX(nxl, cxl+1):MIN(nxr,cxr-1)) = 02038 IF ( index_left_cwall<= nxr .AND. index_right_cwall>= nxl ) & 2039 nzb_local(:,MAX(nxl,index_left_cwall+1):MIN(nxr,index_right_cwall-1)) = 0 2004 2040 ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN 2005 IF ( cys <= nyn .AND. cyn >= nys )&2006 nzb_local(MAX(nys, cys+1):MIN(nyn,cyn-1),:) = 02041 IF ( index_south_cwall <= nyn .AND. index_north_cwall >= nys ) & 2042 nzb_local(MAX(nys,index_south_cwall+1):MIN(nyn,index_north_cwall-1),:) = 0 2007 2043 ENDIF 2008 2044 !
Note: See TracChangeset
for help on using the changeset viewer.
