Changeset 3065 for palm/trunk/SOURCE
- Timestamp:
- Jun 12, 2018 7:03:02 AM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 21 edited
-
check_parameters.f90 (modified) (3 diffs)
-
header.f90 (modified) (3 diffs)
-
init_grid.f90 (modified) (9 diffs)
-
init_masks.f90 (modified) (5 diffs)
-
lpm_advec.f90 (modified) (11 diffs)
-
lpm_exchange_horiz.f90 (modified) (2 diffs)
-
lpm_init.f90 (modified) (4 diffs)
-
lpm_set_attributes.f90 (modified) (5 diffs)
-
modules.f90 (modified) (7 diffs)
-
parin.f90 (modified) (3 diffs)
-
plant_canopy_model_mod.f90 (modified) (5 diffs)
-
pmc_interface_mod.f90 (modified) (3 diffs)
-
radiation_model_mod.f90 (modified) (14 diffs)
-
read_restart_data_mod.f90 (modified) (3 diffs)
-
synthetic_turbulence_generator_mod.f90 (modified) (4 diffs)
-
urban_surface_mod.f90 (modified) (6 diffs)
-
user_init_grid.f90 (modified) (2 diffs)
-
vertical_nesting_mod.f90 (modified) (11 diffs)
-
virtual_flight_mod.f90 (modified) (3 diffs)
-
wind_turbine_model_mod.f90 (modified) (22 diffs)
-
write_restart_data_mod.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/check_parameters.f90
r3049 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dz was replaced by dz(1), error message revised 28 ! 29 ! 3049 2018-05-29 13:52:36Z Giersch 27 30 ! add variable description 28 31 ! … … 4178 4181 IF ( ( constant_flux_layer .OR. & 4179 4182 INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) & 4180 .AND. roughness_length >= 0.5 * dz ) THEN4183 .AND. roughness_length >= 0.5 * dz(1) ) THEN 4181 4184 message_string = 'roughness_length must be smaller than dz/2' 4182 4185 CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 ) … … 4198 4201 !-- Check if vertical grid stretching is switched off in case of complex 4199 4202 !-- terrain simulations 4200 IF ( complex_terrain .AND. dz_stretch_level < 100000.0_wp ) THEN 4203 IF ( complex_terrain .AND. & 4204 ANY( dz_stretch_level_start /= -9999999.9_wp ) ) THEN 4201 4205 message_string = 'Vertical grid stretching is not allowed for ' // & 4202 4206 'complex_terrain = .T.' -
palm/trunk/SOURCE/header.f90
r3045 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Header output concerning stretching revised 28 ! 29 ! 3045 2018-05-28 07:55:41Z Giersch 27 30 ! Error messages revised 28 31 ! … … 800 803 !-- Computational grid 801 804 IF ( .NOT. ocean ) THEN 802 WRITE ( io, 250 ) dx, dy, dz, (nx+1)*dx, (ny+1)*dy, zu(nzt+1) 803 IF ( dz_stretch_level_index < nzt+1 ) THEN 804 WRITE ( io, 252 ) dz_stretch_level, dz_stretch_level_index, & 805 dz_stretch_factor, dz_max 806 ENDIF 805 WRITE ( io, 250 ) dx, dy 806 807 DO i = 1, number_stretch_level_start+1 808 WRITE ( io, 253 ) i, dz(i) 809 ENDDO 810 811 WRITE( io, 251 ) (nx+1)*dx, (ny+1)*dy, zu(nzt+1) 812 813 IF ( ANY( dz_stretch_level_start_index < nzt+1 ) ) THEN 814 WRITE( io, '(A)', advance='no') ' Vertical stretching starts at height:' 815 DO i = 1, number_stretch_level_start 816 WRITE ( io, '(F10.1,A3)', advance='no' ) dz_stretch_level_start(i), ' m,' 817 ENDDO 818 WRITE( io, '(/,A)', advance='no') ' Vertical stretching starts at index: ' 819 DO i = 1, number_stretch_level_start 820 WRITE ( io, '(I12,A1)', advance='no' ) dz_stretch_level_start_index(i), ',' 821 ENDDO 822 WRITE( io, '(/,A)', advance='no') ' Vertical stretching ends at height: ' 823 DO i = 1, number_stretch_level_start 824 WRITE ( io, '(F10.1,A3)', advance='no' ) dz_stretch_level_end(i), ' m,' 825 ENDDO 826 WRITE( io, '(/,A)', advance='no') ' Vertical stretching ends at index: ' 827 DO i = 1, number_stretch_level_start 828 WRITE ( io, '(I12,A1)', advance='no' ) dz_stretch_level_end_index(i), ',' 829 ENDDO 830 WRITE( io, '(/,A)', advance='no') ' Factor used for stretching: ' 831 DO i = 1, number_stretch_level_start 832 WRITE ( io, '(F12.3,A1)', advance='no' ) dz_stretch_factor_array(i), ',' 833 ENDDO 834 ENDIF 835 807 836 ELSE 808 WRITE ( io, 250 ) dx, dy, dz, (nx+1)*dx, (ny+1)*dy, zu(0) 809 IF ( dz_stretch_level_index > 0 ) THEN 810 WRITE ( io, 252 ) dz_stretch_level, dz_stretch_level_index, & 811 dz_stretch_factor, dz_max 812 ENDIF 813 ENDIF 814 WRITE ( io, 254 ) nx, ny, nzt+1, MIN( nnx, nx+1 ), MIN( nny, ny+1 ), & 837 WRITE ( io, 250 ) dx, dy 838 DO i = 1, number_stretch_level_start+1 839 WRITE ( io, 253 ) i, dz(i) 840 ENDDO 841 842 WRITE ( io, 251 ) (nx+1)*dx, (ny+1)*dy, zu(0) 843 844 IF ( ANY( dz_stretch_level_start_index > 0 ) ) THEN 845 WRITE( io, '(A)', advance='no') ' Vertical stretching starts at height:' 846 DO i = 1, number_stretch_level_start 847 WRITE ( io, '(F10.1,A3)', advance='no' ) dz_stretch_level_start(i), ' m,' 848 ENDDO 849 WRITE( io, '(/,A)', advance='no') ' Vertical stretching starts at index: ' 850 DO i = 1, number_stretch_level_start 851 WRITE ( io, '(I12,A1)', advance='no' ) dz_stretch_level_start_index(i), ',' 852 ENDDO 853 WRITE( io, '(/,A)', advance='no') ' Vertical stretching ends at height: ' 854 DO i = 1, number_stretch_level_start 855 WRITE ( io, '(F10.1,A3)', advance='no' ) dz_stretch_level_end(i), ' m,' 856 ENDDO 857 WRITE( io, '(/,A)', advance='no') ' Vertical stretching ends at index: ' 858 DO i = 1, number_stretch_level_start 859 WRITE ( io, '(I12,A1)', advance='no' ) dz_stretch_level_end_index(i), ',' 860 ENDDO 861 WRITE( io, '(/,A)', advance='no') ' Factor used for stretching: ' 862 DO i = 1, number_stretch_level_start 863 WRITE ( io, '(F12.3,A1)', advance='no' ) dz_stretch_factor_array(i), ',' 864 ENDDO 865 ENDIF 866 ENDIF 867 WRITE ( io, 254 ) nx, ny, nzt+1, MIN( nnx, nx+1 ), MIN( nny, ny+1 ), & 815 868 MIN( nnz+2, nzt+2 ) 816 869 IF ( sloping_surface ) WRITE ( io, 260 ) alpha_surface … … 2083 2136 250 FORMAT (//' Computational grid and domain size:'/ & 2084 2137 ' ----------------------------------'// & 2085 ' Grid length: dx = ',F7.3,' m dy = ',F7.3, & 2086 ' m dz = ',F7.3,' m'/ & 2087 ' Domain size: x = ',F10.3,' m y = ',F10.3, & 2138 ' Grid length: dx = ',F8.3,' m dy = ',F8.3, ' m') 2139 251 FORMAT ( /' Domain size: x = ',F10.3,' m y = ',F10.3, & 2088 2140 ' m z(u) = ',F10.3,' m'/) 2089 252 FORMAT (' dz constant up to ',F10.3,' m (k=',I4,'), then stretched by', & 2090 ' factor:',F6.3/ & 2091 ' maximum dz not to be exceeded is dz_max = ',F10.3,' m'/) 2092 254 FORMAT (' Number of gridpoints (x,y,z): (0:',I4,', 0:',I4,', 0:',I4,')'/ & 2141 253 FORMAT ( ' dz(',I1,') = ', F8.3, ' m') 2142 254 FORMAT (//' Number of gridpoints (x,y,z): (0:',I4,', 0:',I4,', 0:',I4,')'/ & 2093 2143 ' Subdomain size (x,y,z): ( ',I4,', ',I4,', ',I4,')'/) 2094 2144 260 FORMAT (/' The model has a slope in x-direction. Inclination angle: ',F6.2,& -
palm/trunk/SOURCE/init_grid.f90
r3051 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! New vertical stretching mechanism introduced 28 ! 29 ! 3051 2018-05-30 17:43:55Z suehring 27 30 ! Minor bugfix concerning mapping 3D buildings on top of terrain 28 31 ! … … 325 328 canyon_height, canyon_wall_left, canyon_wall_south, & 326 329 canyon_width_x, canyon_width_y, constant_flux_layer, & 327 dp_level_ind_b, dz, dz_max, dz_stretch_factor, & 328 dz_stretch_level, dz_stretch_level_index, grid_level, & 330 dp_level_ind_b, dz, dz_max, dz_stretch_factor, & 331 dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end,& 332 dz_stretch_level_end_index, dz_stretch_level_start_index, & 333 dz_stretch_level_start, grid_level, & 329 334 force_bound_l, force_bound_r, force_bound_n, force_bound_s, & 330 335 ibc_uv_b, inflow_l, inflow_n, inflow_r, inflow_s, & 331 336 masking_method, maximum_grid_level, message_string, & 332 momentum_advec, nest_domain, nest_bound_l, nest_bound_n, & 333 nest_bound_r, nest_bound_s, ocean, outflow_l, outflow_n, & 334 outflow_r, outflow_s, psolver, scalar_advec, topography, & 335 topography_grid_convention, tunnel_height, tunnel_length, & 336 tunnel_width_x, tunnel_width_y, tunnel_wall_depth, & 337 use_surface_fluxes, use_top_fluxes, wall_adjustment_factor 337 momentum_advec, nest_domain, nest_bound_l, & 338 nest_bound_n, nest_bound_r, nest_bound_s, & 339 number_stretch_level_end, number_stretch_level_start, ocean, & 340 outflow_l, outflow_n, outflow_r, outflow_s, psolver, & 341 scalar_advec, topography, topography_grid_convention, & 342 tunnel_height, tunnel_length, tunnel_width_x, tunnel_width_y, & 343 tunnel_wall_depth, use_surface_fluxes, use_top_fluxes, & 344 wall_adjustment_factor 338 345 339 346 USE grid_variables, & … … 361 368 IMPLICIT NONE 362 369 363 INTEGER(iwp) :: i !< index variable along x 364 INTEGER(iwp) :: j !< index variable along y 365 INTEGER(iwp) :: k !< index variable along z 366 INTEGER(iwp) :: k_top !< topography top index on local PE 367 INTEGER(iwp) :: l !< loop variable 368 INTEGER(iwp) :: nzb_local_max !< vertical grid index of maximum topography height 369 INTEGER(iwp) :: nzb_local_min !< vertical grid index of minimum topography height 370 INTEGER(iwp) :: i !< index variable along x 371 INTEGER(iwp) :: j !< index variable along y 372 INTEGER(iwp) :: k !< index variable along z 373 INTEGER(iwp) :: k_top !< topography top index on local PE 374 INTEGER(iwp) :: n !< loop variable for stretching 375 INTEGER(iwp) :: number_dz !< number of user-specified dz values 376 INTEGER(iwp) :: nzb_local_max !< vertical grid index of maximum topography height 377 INTEGER(iwp) :: nzb_local_min !< vertical grid index of minimum topography height 370 378 371 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local 372 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp 379 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography top at cell-center 380 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp !< dummy to calculate topography indices on u- and v-grid 373 381 374 382 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo !< input array for 3D topography and dummy array for setting "outer"-flags 375 383 384 REAL(wp) :: dz_level_end !< distance between calculated height level for u/v-grid and user-specified end level for stretching 376 385 REAL(wp) :: dz_stretched !< stretched vertical grid spacing 386 387 REAL(wp), DIMENSION(:), ALLOCATABLE :: min_dz_stretch_level_end !< Array that contains all minimum heights where the stretching can end 377 388 378 389 … … 391 402 ! 392 403 !-- Compute height of u-levels from constant grid length and dz stretch factors 393 IF ( dz == -1.0_wp ) THEN404 IF ( dz(1) == -1.0_wp ) THEN 394 405 message_string = 'missing dz' 395 406 CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 ) 396 ELSEIF ( dz <= 0.0_wp ) THEN397 WRITE( message_string, * ) 'dz=',dz ,' <= 0.0'407 ELSEIF ( dz(1) <= 0.0_wp ) THEN 408 WRITE( message_string, * ) 'dz=',dz(1),' <= 0.0' 398 409 CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 ) 399 410 ENDIF 400 411 401 412 ! 413 !-- Initialize dz_stretch_level_start with the value of dz_stretch_level 414 !-- if it was set by the user 415 IF ( dz_stretch_level /= -9999999.9_wp ) THEN 416 dz_stretch_level_start(1) = dz_stretch_level 417 ENDIF 418 419 ! 420 !-- Determine number of dz values and stretching levels specified by the 421 !-- user to allow right controlling of the stretching mechanism and to 422 !-- perform error checks 423 number_dz = COUNT( dz /= -1.0_wp ) 424 number_stretch_level_start = COUNT( dz_stretch_level_start /= & 425 -9999999.9_wp ) 426 number_stretch_level_end = COUNT( dz_stretch_level_end /= & 427 9999999.9_wp ) 428 429 ! 430 !-- The number of specified end levels +1 has to be the same than the number 431 !-- of specified dz values 432 IF ( number_dz /= number_stretch_level_end + 1 ) THEN 433 WRITE( message_string, * ) 'The number of values for dz = ', & 434 number_dz, 'has to be the same than& ', & 435 'the number of values for ', & 436 'dz_stretch_level_end + 1 = ', & 437 number_stretch_level_end+1 438 CALL message( 'init_grid', 'PA0156', 1, 2, 0, 6, 0 ) 439 ENDIF 440 441 ! 442 !-- The number of specified start levels has to be the same or one less than 443 !-- the number of specified dz values 444 IF ( number_dz /= number_stretch_level_start + 1 .AND. & 445 number_dz /= number_stretch_level_start ) THEN 446 WRITE( message_string, * ) 'The number of values for dz = ', & 447 number_dz, 'has to be the same or one ', & 448 'more than& the number of values for ', & 449 'dz_stretch_level_start = ', & 450 number_stretch_level_start 451 CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 ) 452 ENDIF 453 454 !-- The number of specified start levels has to be the same or one more than 455 !-- the number of specified end levels 456 IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND. & 457 number_stretch_level_start /= number_stretch_level_end ) THEN 458 WRITE( message_string, * ) 'The number of values for ', & 459 'dz_stretch_level_start = ', & 460 dz_stretch_level_start, 'has to be the ',& 461 'same or one more than& the number of ', & 462 'values for dz_stretch_level_end = ', & 463 number_stretch_level_end 464 CALL message( 'init_grid', 'PA0216', 1, 2, 0, 6, 0 ) 465 ENDIF 466 467 ! 468 !-- Initialize dz for the free atmosphere with the value of dz_max 469 IF ( dz(number_stretch_level_start+1) == -1.0_wp .AND. & 470 number_stretch_level_start /= 0 ) THEN 471 dz(number_stretch_level_start+1) = dz_max 472 ENDIF 473 474 ! 475 !-- Initialize the stretching factor if (infinitely) stretching in the free 476 !-- atmosphere is desired (dz_stretch_level_end was not specified for the 477 !-- free atmosphere) 478 IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN 479 dz_stretch_factor_array(number_stretch_level_start) = & 480 dz_stretch_factor 481 ENDIF 482 483 ! 484 !-- Allocation of arrays for stretching 485 ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) ) 486 487 ! 402 488 !-- Define the vertical grid levels 403 489 IF ( .NOT. ocean ) THEN 490 491 ! 492 !-- The stretching region has to be large enough to allow for a smooth 493 !-- transition between two different grid spacings 494 DO n = 1, number_stretch_level_start 495 min_dz_stretch_level_end(n) = dz_stretch_level_start(n) + & 496 4 * MAX( dz(n),dz(n+1) ) 497 ENDDO 498 499 IF ( ANY( min_dz_stretch_level_end > dz_stretch_level_end ) ) THEN 500 message_string= 'Eeach dz_stretch_level_end has to be larger ' // & 501 'than its corresponding value for &' // & 502 'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//& 503 'to allow for smooth grid stretching' 504 CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 ) 505 ENDIF 506 507 ! 508 !-- Stretching must not be applied within the prandtl_layer 509 !-- (first two grid points). For the default case dz_stretch_level_start 510 !-- is negative. Therefore the absolut value is checked here. 511 IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_wp ) ) THEN 512 WRITE( message_string, * ) 'Eeach dz_stretch_level_start has to be ',& 513 'larger than ', dz(1) * 1.5 514 CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 ) 515 ENDIF 516 517 ! 518 !-- The stretching has to start and end on a grid level. Therefore 519 !-- user-specified values have to ''interpolate'' to the next lowest level 520 IF ( number_stretch_level_start /= 0 ) THEN 521 dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) - & 522 dz(1)/2.0) / dz(1) ) & 523 * dz(1) + dz(1)/2.0 524 ENDIF 525 526 IF ( number_stretch_level_start > 1 ) THEN 527 DO n = 2, number_stretch_level_start 528 dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / & 529 dz(n) ) * dz(n) 530 ENDDO 531 ENDIF 532 533 IF ( number_stretch_level_end /= 0 ) THEN 534 DO n = 1, number_stretch_level_end 535 dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / & 536 dz(n+1) ) * dz(n+1) 537 ENDDO 538 ENDIF 539 540 ! 541 !-- Determine stretching factor if necessary 542 IF ( number_stretch_level_end >= 1 ) THEN 543 CALL calculate_stretching_factor( number_stretch_level_end ) 544 ENDIF 545 404 546 ! 405 547 !-- Grid for atmosphere with surface at z=0 (k=0, w-grid). 548 !-- First compute the u- and v-levels. In case of dirichlet bc for u and v 549 !-- the first u/v- and w-level (k=0) are defined at same height (z=0). 406 550 !-- The second u-level (k=1) corresponds to the top of the 407 551 !-- Prandtl-layer. 408 409 552 IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN 410 553 zu(0) = 0.0_wp 411 ! zu(0) = - dz * 0.5_wp412 554 ELSE 413 zu(0) = - dz * 0.5_wp 414 ENDIF 415 zu(1) = dz * 0.5_wp 416 417 dz_stretch_level_index = nzt+1 418 dz_stretched = dz 555 zu(0) = - dz(1) * 0.5_wp 556 ENDIF 557 558 zu(1) = dz(1) * 0.5_wp 559 560 ! 561 !-- Determine u and v height levels considering the possibility of grid 562 !-- stretching in several heights. 563 n = 1 564 dz_stretch_level_start_index = nzt+1 565 dz_stretch_level_end_index = nzt+1 566 dz_stretched = dz(1) 567 568 !-- The default value of dz_stretch_level_start is negative, thus the first 569 !-- condition is always true. Hence, the second condition is necessary. 419 570 DO k = 2, nzt+1 420 IF ( dz_stretch_level <= zu(k-1) .AND. dz_stretched < dz_max ) THEN 421 dz_stretched = dz_stretched * dz_stretch_factor 422 dz_stretched = MIN( dz_stretched, dz_max ) 423 IF ( dz_stretch_level_index == nzt+1 ) dz_stretch_level_index = k-1 424 ENDIF 571 IF ( dz_stretch_level_start(n) <= zu(k-1) .AND. & 572 dz_stretch_level_start(n) /= -9999999.9_wp ) THEN 573 dz_stretched = dz_stretched * dz_stretch_factor_array(n) 574 575 IF ( dz(n) > dz(n+1) ) THEN 576 dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz 577 ELSE 578 dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz 579 ENDIF 580 581 IF ( dz_stretch_level_start_index(n) == nzt+1 ) & 582 dz_stretch_level_start_index(n) = k-1 583 584 ENDIF 585 425 586 zu(k) = zu(k-1) + dz_stretched 587 588 ! 589 !-- Make sure that the stretching ends exactly at dz_stretch_level_end 590 dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) ) 591 592 IF ( dz_level_end < dz(n+1)/3.0 ) THEN 593 zu(k) = dz_stretch_level_end(n) 594 dz_stretched = dz(n+1) 595 dz_stretch_level_end_index(n) = k 596 n = n + 1 597 ENDIF 426 598 ENDDO 427 599 … … 438 610 439 611 ELSE 612 613 ! 614 !-- The stretching region has to be large enough to allow for a smooth 615 !-- transition between two different grid spacings 616 DO n = 1, number_stretch_level_start 617 min_dz_stretch_level_end(n) = dz_stretch_level_start(n) - & 618 4 * MAX( dz(n),dz(n+1) ) 619 ENDDO 620 621 IF ( ANY( min_dz_stretch_level_end < dz_stretch_level_end ) ) THEN 622 message_string= 'Eeach dz_stretch_level_end has to be less ' // & 623 'than its corresponding value for &' // & 624 'dz_stretch_level_start - 4*MAX(dz(n),dz(n+1)) '//& 625 'to allow for smooth grid stretching' 626 CALL message( 'init_grid', 'PA0224', 1, 2, 0, 6, 0 ) 627 ENDIF 628 629 ! 630 !-- Stretching must not be applied within the prandtl_layer 631 !-- (first two grid points). For the default case dz_stretch_level_start 632 !-- is negative. Therefore the absolut value is checked here. 633 IF ( ANY( dz_stretch_level_start > dz(1) * 1.5_wp ) ) THEN 634 WRITE( message_string, * ) 'Eeach dz_stretch_level_start has to be ',& 635 'less than ', dz(1) * 1.5 636 CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 ) 637 ENDIF 638 639 ! 640 !-- The stretching has to start and end on a grid level. Therefore 641 !-- user-specified values have to ''interpolate'' to the next highest level 642 IF ( number_stretch_level_start /= 0 ) THEN 643 dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) + & 644 dz(1)/2.0) / dz(1) ) & 645 * dz(1) - dz(1)/2.0 646 ENDIF 647 648 IF ( number_stretch_level_start > 1 ) THEN 649 DO n = 2, number_stretch_level_start 650 dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / & 651 dz(n) ) * dz(n) 652 ENDDO 653 ENDIF 654 655 IF ( number_stretch_level_end /= 0 ) THEN 656 DO n = 1, number_stretch_level_end 657 dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / & 658 dz(n+1) ) * dz(n+1) 659 ENDDO 660 ENDIF 661 662 ! 663 !-- Determine stretching factor if necessary 664 IF ( number_stretch_level_end >= 1 ) THEN 665 CALL calculate_stretching_factor( number_stretch_level_end ) 666 ENDIF 667 440 668 ! 441 669 !-- Grid for ocean with free water surface is at k=nzt (w-grid). … … 444 672 !-- w-level are defined at same height, but staggered from the second level. 445 673 !-- The second u-level (k=1) corresponds to the top of the Prandtl-layer. 446 zu(nzt+1) = dz * 0.5_wp 447 zu(nzt) = - dz * 0.5_wp 448 449 dz_stretch_level_index = 0 450 dz_stretched = dz 674 !-- z values are negative starting from z=0 (surface) 675 zu(nzt+1) = dz(1) * 0.5_wp 676 zu(nzt) = - dz(1) * 0.5_wp 677 678 ! 679 !-- Determine u and v height levels considering the possibility of grid 680 !-- stretching in several heights. 681 n = 1 682 dz_stretch_level_start_index = 0 683 dz_stretch_level_end_index = 0 684 dz_stretched = dz(1) 685 451 686 DO k = nzt-1, 0, -1 452 ! 453 !-- The default value of dz_stretch_level is positive, thus the first 454 !-- condition is always true. Hence, the second condition is necessary. 455 IF ( dz_stretch_level >= zu(k+1) .AND. dz_stretch_level <= 0.0 & 456 .AND. dz_stretched < dz_max ) THEN 457 dz_stretched = dz_stretched * dz_stretch_factor 458 dz_stretched = MIN( dz_stretched, dz_max ) 459 IF ( dz_stretch_level_index == 0 ) dz_stretch_level_index = k+1 460 ENDIF 687 688 IF ( dz_stretch_level_start(n) >= zu(k+1) ) THEN 689 dz_stretched = dz_stretched * dz_stretch_factor_array(n) 690 691 IF ( dz(n) > dz(n+1) ) THEN 692 dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz 693 ELSE 694 dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz 695 ENDIF 696 697 IF ( dz_stretch_level_start_index(n) == 0 ) & 698 dz_stretch_level_start_index(n) = k+1 699 700 ENDIF 701 461 702 zu(k) = zu(k+1) - dz_stretched 703 704 ! 705 !-- Make sure that the stretching ends exactly at dz_stretch_level_end 706 dz_level_end = ABS( zu(k) - dz_stretch_level_end(n) ) 707 708 IF ( dz_level_end < dz(n+1)/3.0 ) THEN 709 zu(k) = dz_stretch_level_end(n) 710 dz_stretched = dz(n+1) 711 dz_stretch_level_end_index(n) = k 712 n = n + 1 713 ENDIF 462 714 ENDDO 463 715 464 716 ! 465 717 !-- Compute the w-levels. They are always staggered half-way between the … … 468 720 !-- same height. The top w-level (nzt+1) is not used but set for 469 721 !-- consistency, since w and all scalar variables are defined up tp nzt+1. 470 zw(nzt+1) = dz 722 zw(nzt+1) = dz(1) 471 723 zw(nzt) = 0.0_wp 472 724 DO k = 0, nzt … … 855 1107 END SUBROUTINE init_grid 856 1108 1109 1110 ! Description: 1111 ! -----------------------------------------------------------------------------! 1112 !> Calculation of the stretching factor through an iterative method. Ideas were 1113 !> taken from the paper "Regional stretched grid generation and its application 1114 !> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the 1115 !> system of equations has two variables (r,l) but four requirements 1116 !> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which 1117 !> results into an overdetermined system. 1118 !------------------------------------------------------------------------------! 1119 SUBROUTINE calculate_stretching_factor( number_end ) 1120 1121 USE control_parameters, & 1122 ONLY: dz, dz_stretch_factor, dz_stretch_factor_array, & 1123 dz_stretch_level_end, dz_stretch_level_start, message_string 1124 1125 USE kinds 1126 1127 IMPLICIT NONE 1128 1129 INTEGER(iwp) :: iterations !< number of iterations until stretch_factor_lower/upper_limit is reached 1130 INTEGER(iwp) :: l_rounded !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2 1131 INTEGER(iwp) :: n !< loop variable for stretching 1132 1133 INTEGER(iwp), INTENT(IN) :: number_end !< number of user-specified end levels for stretching 1134 1135 REAL(wp) :: delta_l !< absolute difference between l and l_rounded 1136 REAL(wp) :: delta_stretch_factor !< absolute difference between stretch_factor_1 and stretch_factor_2 1137 REAL(wp) :: delta_total_new !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible) 1138 REAL(wp) :: delta_total_old !< sum of delta_l and delta_stretch_factor for the last iteration 1139 REAL(wp) :: distance !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region) 1140 REAL(wp) :: l !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly 1141 REAL(wp) :: numerator !< numerator of the quotient 1142 REAL(wp) :: stretch_factor_1 !< stretching factor that fulfil Eq. (5) togehter with l exactly 1143 REAL(wp) :: stretch_factor_2 !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly 1144 1145 REAL(wp), PARAMETER :: stretch_factor_interval = 1.0E-06 !< interval for sampling possible stretching factors 1146 REAL(wp), PARAMETER :: stretch_factor_lower_limit = 0.88 !< lowest possible stretching factor 1147 REAL(wp), PARAMETER :: stretch_factor_upper_limit = 1.12 !< highest possible stretching factor 1148 1149 1150 l = 0 1151 DO n = 1, number_end 1152 1153 iterations = 1 1154 stretch_factor_1 = 1.0 1155 stretch_factor_2 = 1.0 1156 delta_total_old = 1.0 1157 1158 IF ( dz(n) > dz(n+1) ) THEN 1159 DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) 1160 1161 stretch_factor_1 = 1.0 - iterations * stretch_factor_interval 1162 distance = ABS( dz_stretch_level_end(n) - & 1163 dz_stretch_level_start(n) ) 1164 numerator = distance*stretch_factor_1/dz(n) + & 1165 stretch_factor_1 - distance/dz(n) 1166 1167 IF ( numerator > 0.0 ) THEN 1168 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0 1169 l_rounded = NINT( l ) 1170 delta_l = ABS( l_rounded - l ) / l 1171 ENDIF 1172 1173 stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) ) 1174 1175 delta_stretch_factor = ABS( stretch_factor_1 - & 1176 stretch_factor_2 ) / & 1177 stretch_factor_2 1178 1179 delta_total_new = delta_l + delta_stretch_factor 1180 1181 ! 1182 !-- stretch_factor_1 is taken to guarantee that the stretching 1183 !-- procedure ends as close as possible to dz_stretch_level_end. 1184 !-- stretch_factor_2 would guarantee that the stretched dz(n) is 1185 !-- equal to dz(n+1) after l_rounded grid levels. 1186 IF (delta_total_new < delta_total_old) THEN 1187 dz_stretch_factor_array(n) = stretch_factor_1 1188 delta_total_old = delta_total_new 1189 ENDIF 1190 1191 iterations = iterations + 1 1192 1193 ENDDO 1194 1195 ELSEIF ( dz(n) < dz(n+1) ) THEN 1196 DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit ) 1197 1198 stretch_factor_1 = 1.0 + iterations * stretch_factor_interval 1199 distance = ABS( dz_stretch_level_end(n) - & 1200 dz_stretch_level_start(n) ) 1201 numerator = distance*stretch_factor_1/dz(n) + & 1202 stretch_factor_1 - distance/dz(n) 1203 1204 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0 1205 l_rounded = NINT( l ) 1206 delta_l = ABS( l_rounded - l ) / l 1207 1208 stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) ) 1209 1210 delta_stretch_factor = ABS( stretch_factor_1 - & 1211 stretch_factor_2 ) / & 1212 stretch_factor_2 1213 1214 delta_total_new = delta_l + delta_stretch_factor 1215 1216 ! 1217 !-- stretch_factor_1 is taken to guarantee that the stretching 1218 !-- procedure ends as close as possible to dz_stretch_level_end. 1219 !-- stretch_factor_2 would guarantee that the stretched dz(n) is 1220 !-- equal to dz(n+1) after l_rounded grid levels. 1221 IF (delta_total_new < delta_total_old) THEN 1222 dz_stretch_factor_array(n) = stretch_factor_1 1223 delta_total_old = delta_total_new 1224 ENDIF 1225 1226 iterations = iterations + 1 1227 ENDDO 1228 1229 ELSE 1230 message_string= 'Two adjacent values of dz must be different' 1231 CALL message( 'init_grid', 'PA0228', 1, 2, 0, 6, 0 ) 1232 1233 ENDIF 1234 ENDDO 1235 1236 END SUBROUTINE calculate_stretching_factor 1237 1238 857 1239 ! Description: 858 1240 ! -----------------------------------------------------------------------------! … … 1703 2085 !-- Tunnel-wall depth 1704 2086 IF ( tunnel_wall_depth == 9999999.9_wp ) THEN 1705 td = MAX ( dx, dy, dz )2087 td = MAX ( dx, dy, dz(1) ) 1706 2088 ELSE 1707 2089 td = tunnel_wall_depth -
palm/trunk/SOURCE/init_masks.f90
r3049 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dz_stretch_level was replaced by dz_stretch_level_start 28 ! 29 ! 3049 2018-05-29 13:52:36Z Giersch 27 30 ! Error messages revised 28 31 ! … … 133 136 ONLY: constant_diffusion, cloud_droplets, cloud_physics, & 134 137 data_output_masks, data_output_masks_user, & 135 doav, doav_n, domask, domask_no, dz, dz_stretch_level , humidity,&136 mask, masks, mask_scale, mask_i,&138 doav, doav_n, domask, domask_no, dz, dz_stretch_level_start, & 139 humidity, mask, masks, mask_scale, mask_i, & 137 140 mask_i_global, mask_j, mask_j_global, mask_k, mask_k_global, & 138 141 mask_loop, mask_size, mask_size_l, mask_start_l, mask_x, & … … 141 144 microphysics_morrison, microphysics_seifert, passive_scalar, & 142 145 ocean, varnamelength 146 143 147 144 148 USE grid_variables, & … … 490 494 CALL set_mask_locations( 1, dx, 'dx', nx, 'nx', nxl, nxr ) 491 495 CALL set_mask_locations( 2, dy, 'dy', ny, 'ny', nys, nyn ) 492 CALL set_mask_locations( 3, dz , 'dz', nz, 'nz', nzb, nzt )496 CALL set_mask_locations( 3, dz(1), 'dz', nz, 'nz', nzb, nzt ) 493 497 ! 494 498 !-- Set global masks along all three dimensions (required by … … 727 731 !-- The following line assumes a constant vertical grid spacing within 728 732 !-- the vertical mask range; it fails for vertical grid stretching. 729 !-- Maybe revise later. Issue warning but continue execution. 733 !-- Maybe revise later. Issue warning but continue execution. ABS(...) 734 !-- within the IF statement is necessary because the default value of 735 !-- dz_stretch_level_start is -9999999.9_wp. 730 736 loop_stride = NINT( mask_loop(mid,dim,3) * mask_scale(dim) * ddxyz ) 731 737 732 IF ( mask_loop(mid,dim,2) * mask_scale(dim) > dz_stretch_level )&733 THEN738 IF ( mask_loop(mid,dim,2) * mask_scale(dim) > & 739 ABS( dz_stretch_level_start(1) ) ) THEN 734 740 WRITE ( message_string, '(A,I3,A,I1,A,F9.3,A,F8.2,3A)' ) & 735 741 'mask_loop(',mid,',',dim,',2)=', mask_loop(mid,dim,2), & 736 ' exceeds dz_stretch_level=',dz_stretch_level ,&742 ' exceeds dz_stretch_level=',dz_stretch_level_start(1), & 737 743 '.&Vertical mask locations will not ', & 738 744 'match the desired heights within the stretching ', & -
palm/trunk/SOURCE/lpm_advec.f90
r2969 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dz values were replaced by dzw or dz(1) to allow for right vertical stretching 28 ! 29 ! 2969 2018-04-13 11:55:09Z thiele 27 30 ! Bugfix in Interpolation indices. 28 31 ! … … 158 161 159 162 USE arrays_3d, & 160 ONLY: de_dx, de_dy, de_dz, diss, e, km, u, v, w, zu, zw163 ONLY: de_dx, de_dy, de_dz, diss, dzw, e, km, u, v, w, zu, zw 161 164 162 165 USE cpulog … … 409 412 + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * & 410 413 u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans 411 u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *&414 u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 412 415 ( u_int_u - u_int_l ) 413 416 ENDIF … … 503 506 + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) & 504 507 ) / ( 3.0_wp * gg ) - v_gtrans 505 v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz *&508 v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 506 509 ( v_int_u - v_int_l ) 507 510 ENDIF … … 540 543 ( gg-dd ) * w(k+1,j+1,i+1) & 541 544 ) / ( 3.0_wp * gg ) 542 w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz *&545 w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k) * & 543 546 ( w_int_u - w_int_l ) 544 547 ENDIF … … 611 614 ( gg - dd ) * e(k+1,j+1,i+1) & 612 615 ) / ( 3.0_wp * gg ) 613 e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz *&616 e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 614 617 ( e_int_u - e_int_l ) 615 618 ENDIF … … 637 640 ( gg - dd ) * de_dx(k+1,j+1,i+1) & 638 641 ) / ( 3.0_wp * gg ) 639 de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz *&642 de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 640 643 ( de_dx_int_u - de_dx_int_l ) 641 644 ENDIF … … 655 658 ( gg - dd ) * de_dy(k+1,j+1,i+1) & 656 659 ) / ( 3.0_wp * gg ) 657 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &660 de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 658 661 ( de_dy_int_u - de_dy_int_l ) 659 662 ENDIF … … 661 664 ! 662 665 !-- Interpolate the TKE gradient along z 663 IF ( zv(n) < 0.5_wp * dz ) THEN666 IF ( zv(n) < 0.5_wp * dz(1) ) THEN 664 667 de_dz_int(n) = 0.0_wp 665 668 ELSE … … 678 681 ( gg - dd ) * de_dz(k+1,j+1,i+1) & 679 682 ) / ( 3.0_wp * gg ) 680 de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &683 de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 681 684 ( de_dz_int_u - de_dz_int_l ) 682 685 ENDIF … … 699 702 ( gg - dd ) * diss(k+1,j+1,i+1) & 700 703 ) / ( 3.0_wp * gg ) 701 diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&704 diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 702 705 ( diss_int_u - diss_int_l ) 703 706 ENDIF -
palm/trunk/SOURCE/lpm_exchange_horiz.f90
r3049 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dz was replaced by dz(1) to allow for right vertical stretching 28 ! 29 ! 3049 2018-05-29 13:52:36Z Giersch 27 30 ! Error messages revised 28 31 ! … … 910 913 ip = particle_array(n)%x * ddx 911 914 jp = particle_array(n)%y * ddy 912 kp = particle_array(n)%z / dz + 1 + offset_ocean_nzt915 kp = particle_array(n)%z / dz(1) + 1 + offset_ocean_nzt 913 916 ! 914 917 !-- In case of grid stretching a particle might be above or below the -
palm/trunk/SOURCE/lpm_init.f90
r3049 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dz was replaced by dzw or dz(1) to allow for right vertical stretching 28 ! 29 ! 3049 2018-05-29 13:52:36Z Giersch 27 30 ! Error messages revised 28 31 ! … … 794 797 ip = tmp_particle%x * ddx 795 798 jp = tmp_particle%y * ddy 796 kp = tmp_particle%z / dz + 1 + offset_ocean_nzt799 kp = tmp_particle%z / dz(1) + 1 + offset_ocean_nzt 797 800 DO WHILE( zw(kp) < tmp_particle%z ) 798 801 kp = kp + 1 … … 969 972 pdz(particles(n)%group) 970 973 particles(n)%z = particles(n)%z + & 971 MERGE( rand_contr, SIGN( dz , rand_contr ),&972 ABS( rand_contr ) < dz 974 MERGE( rand_contr, SIGN( dzw(kp), rand_contr ), & 975 ABS( rand_contr ) < dzw(kp) & 973 976 ) 974 977 ENDIF … … 987 990 i = particles(n)%x * ddx 988 991 j = particles(n)%y * ddy 989 k = particles(n)%z / dz + 1 + offset_ocean_nzt992 k = particles(n)%z / dz(1) + 1 + offset_ocean_nzt 990 993 DO WHILE( zw(k) < particles(n)%z ) 991 994 k = k + 1 -
palm/trunk/SOURCE/lpm_set_attributes.f90
r2718 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dz was replaced by dzw to allow for right vertical stretching 28 ! 29 ! 2718 2018-01-02 08:49:38Z maronga 27 30 ! Corrected "Former revisions" section 28 31 ! … … 80 83 81 84 USE arrays_3d, & 82 ONLY: pt, u, v, w, zu, zw85 ONLY: dzw, pt, u, v, w, zu, zw 83 86 84 87 USE control_parameters, & 85 ONLY: u_gtrans, v_gtrans , dz88 ONLY: u_gtrans, v_gtrans 86 89 87 90 USE cpulog, & … … 207 210 ( gg - dd ) * u(k+1,j+1,i+1) & 208 211 ) / ( 3.0_wp * gg ) - u_gtrans 209 u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz *&212 u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 210 213 ( u_int_u - u_int_l ) 211 214 ENDIF … … 240 243 ( gg - dd ) * v(k+1,j+1,i+1) & 241 244 ) / ( 3.0_wp * gg ) - v_gtrans 242 v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz *&245 v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 243 246 ( v_int_u - v_int_l ) 244 247 ENDIF … … 344 347 ) / ( 3.0_wp * gg ) - sums(k,4) 345 348 346 pt_int = pt_int_l + ( zv(n) - zu(k) ) / dz *&349 pt_int = pt_int_l + ( zv(n) - zu(k) ) / dzw(k) * & 347 350 ( pt_int_u - pt_int_l ) 348 351 -
palm/trunk/SOURCE/modules.f90
r3045 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Variables concerning stretching introduced or revised 28 ! 29 ! 3045 2018-05-28 07:55:41Z Giersch 27 30 ! z_max_do2d removed 28 31 ! … … 1148 1151 INTEGER(iwp) :: dp_level_ind_b = 0 !< lowest grid index for external pressure gradient forcing 1149 1152 INTEGER(iwp) :: dvrp_filecount = 0 !< parameter for dvr visualization software 1150 INTEGER(iwp) :: dz_stretch_level_index !< vertical grid level index above which the vertical grid spacing is stretched1151 1153 INTEGER(iwp) :: ensemble_member_nr = 0 !< namelist parameter 1152 1154 INTEGER(iwp) :: gamma_mg !< switch for steering the multigrid cycle: 1: v-cycle, 2: w-cycle … … 1190 1192 INTEGER(iwp) :: num_var_fl !< number of sampling/output variables in virtual flight measurements 1191 1193 INTEGER(iwp) :: num_var_fl_user=0 !< number of user-defined sampling/output variables in virtual flight measurements 1194 INTEGER(iwp) :: number_stretch_level_start !< number of user-specified start levels for stretching 1195 INTEGER(iwp) :: number_stretch_level_end !< number of user-specified end levels for stretching 1192 1196 INTEGER(iwp) :: nz_do3d = -9999 !< namelist parameter 1193 1197 INTEGER(iwp) :: prt_time_count = 0 !< number of output intervals for particle data output … … 1199 1203 INTEGER(iwp) :: timestep_count = 0 !< number of timesteps carried out since the beginning of the initial run 1200 1204 INTEGER(iwp) :: y_shift = 0 !< namelist parameter 1205 1201 1206 INTEGER(iwp) :: dist_nxl(0:1) !< left boundary of disturbance region 1202 1207 INTEGER(iwp) :: dist_nxr(0:1) !< right boundary of disturbance region … … 1211 1216 INTEGER(iwp) :: domask_no(max_masks,0:1) = 0 !< number of masked output quantities 1212 1217 INTEGER(iwp) :: domask_time_count(max_masks,0:1) !< number of output intervals for masked data 1218 INTEGER(iwp) :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched 1219 INTEGER(iwp) :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched 1213 1220 INTEGER(iwp) :: mask_size(max_masks,3) = -1 !< size of mask array per mask and dimension (for netcdf output) 1214 1221 INTEGER(iwp) :: mask_size_l(max_masks,3) = -1 !< subdomain size of mask array per mask and dimension (for netcdf output) … … 1416 1423 REAL(wp) :: dt_spinup = 60.0_wp !< namelist parameter 1417 1424 REAL(wp) :: dt_3d = 1.0_wp !< time step 1418 REAL(wp) :: dz = -1.0_wp !< namelist parameter 1419 REAL(wp) :: dz_max = 9999999.9_wp !< namelist parameter 1425 REAL(wp) :: dz_max = 1000.0_wp !< namelist parameter 1420 1426 REAL(wp) :: dz_stretch_factor = 1.08_wp !< namelist parameter 1421 REAL(wp) :: dz_stretch_level = 100000.0_wp!< namelist parameter1427 REAL(wp) :: dz_stretch_level = -9999999.9_wp !< namelist parameter 1422 1428 REAL(wp) :: e_init = 0.0_wp !< namelist parameter 1423 1429 REAL(wp) :: e_min = 0.0_wp !< namelist parameter … … 1529 1535 REAL(wp) :: dpdxy(1:2) = 0.0_wp !< namelist parameter 1530 1536 REAL(wp) :: dt_domask(max_masks) = 9999999.9_wp !< namelist parameter 1537 REAL(wp) :: dz(10) = -1.0_wp !< namelist parameter 1538 REAL(wp) :: dz_stretch_level_start(9) = -9999999.9_wp !< namelist parameter 1539 REAL(wp) :: dz_stretch_level_end(9) = 9999999.9_wp !< namelist parameter 1540 REAL(wp) :: dz_stretch_factor_array(9) = 1.08_wp !< namelist parameter 1531 1541 REAL(wp) :: mask_scale(3) !< collective array for mask_scale_x/y/z 1532 1542 REAL(wp) :: pt_vertical_gradient(10) = 0.0_wp !< namelist parameter -
palm/trunk/SOURCE/parin.f90
r3049 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! New initialization parameters added 28 ! 29 ! 3049 2018-05-29 13:52:36Z Giersch 27 30 ! Error messages revised 28 31 ! … … 509 512 dp_external, dp_level_b, dp_smooth, dpdxy, dry_aerosol_radius, & 510 513 dt, dt_pr_1d, dt_run_control_1d, dt_spinup, dx, dy, dz, dz_max, & 511 dz_stretch_factor, dz_stretch_level, end_time_1d,&512 ensemble_member_nr, e_init, e_min, fft_method,&513 flux_input_mode, flux_output_mode, forcing,&514 dz_stretch_factor, dz_stretch_level, dz_stretch_level_start, & 515 dz_stretch_level_end, end_time_1d, ensemble_member_nr, e_init, & 516 e_min, fft_method, flux_input_mode, flux_output_mode, forcing, & 514 517 galilei_transformation, humidity, & 515 518 inflow_damping_height, inflow_damping_width, & … … 579 582 dp_external, dp_level_b, dp_smooth, dpdxy, dry_aerosol_radius, & 580 583 dt, dt_pr_1d, dt_run_control_1d, dt_spinup, dx, dy, dz, dz_max, & 581 dz_stretch_factor, dz_stretch_level, end_time_1d,&582 ensemble_member_nr, e_init, e_min, fft_method,&583 flux_input_mode, flux_output_mode, forcing,&584 dz_stretch_factor, dz_stretch_level, dz_stretch_level_start, & 585 dz_stretch_level_end, end_time_1d, ensemble_member_nr, e_init, & 586 e_min, fft_method, flux_input_mode, flux_output_mode, forcing, & 584 587 galilei_transformation, humidity, & 585 588 inflow_damping_height, inflow_damping_width, & -
palm/trunk/SOURCE/plant_canopy_model_mod.f90
r3049 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dz was replaced by the help of zw to allow for vertical stretching 28 ! 29 ! 3049 2018-05-29 13:52:36Z Giersch 27 30 ! Error messages revised 28 31 ! … … 550 553 551 554 USE control_parameters, & 552 ONLY: dz,passive_scalar555 ONLY: passive_scalar 553 556 554 557 … … 568 571 REAL(wp) :: canopy_height !< canopy height (in m) 569 572 570 canopy_height = pch_index * dz573 canopy_height = zw(pch_index) 571 574 572 575 WRITE ( io, 1 ) canopy_mode, canopy_height, pch_index, & … … 668 671 669 672 USE control_parameters, & 670 ONLY: dz, humidity, io_blocks, io_group, message_string, ocean,&673 ONLY: humidity, io_blocks, io_group, message_string, ocean, & 671 674 passive_scalar, urban_surface 672 675 … … 761 764 !-- Use beta function for lad-profile construction 762 765 int_bpdf = 0.0_wp 763 canopy_height = pch_index * dz766 canopy_height = zw(pch_index) 764 767 765 768 DO k = 0, pch_index -
palm/trunk/SOURCE/pmc_interface_mod.f90
r3049 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dz was replaced by dz(1) 28 ! 29 ! 3049 2018-05-29 13:52:36Z Giersch 27 30 ! Error messages revised 28 31 ! … … 913 916 parent_grid_info_real(5) = lower_left_coord_x + ( nx + 1 ) * dx 914 917 parent_grid_info_real(6) = lower_left_coord_y + ( ny + 1 ) * dy 915 parent_grid_info_real(7) = dz 918 parent_grid_info_real(7) = dz(1) 916 919 917 920 parent_grid_info_int(1) = nx … … 1318 1321 fval(3) = dx 1319 1322 fval(4) = dy 1320 fval(5) = dz 1323 fval(5) = dz(1) 1321 1324 1322 1325 IF ( myid == 0 ) THEN -
palm/trunk/SOURCE/radiation_model_mod.f90
r3049 r3065 28 28 ! ----------------- 29 29 ! $Id$ 30 ! dz was replaced by dz(1), error message concerning vertical stretching was 31 ! added 32 ! 33 ! 3049 2018-05-29 13:52:36Z Giersch 30 34 ! Error messages revised 31 35 ! … … 330 334 ! ------------ 331 335 !> Radiation models and interfaces 336 !> @todo Replace dz(1) appropriatly to account for grid stretching 332 337 !> @todo move variable definitions used in radiation_init only to the subroutine 333 338 !> as they are no longer required after initialization. … … 4371 4376 IF ( plant_canopy ) THEN 4372 4377 pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp & 4373 / (cp * hyp(nzub:nzut) * dx*dy*dz ) !< equals to 1 / (rho * c_p * Vbox * T)4378 / (cp * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T) 4374 4379 pctf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp & 4375 4380 / (l_v * hyp(nzub:nzut) * dx*dy*dz) … … 4396 4401 !-- each dimension, so that this new direction vector will allow us 4397 4402 !-- to traverse the ray path within grid coordinates directly 4398 sunorig_grid = (/ sunorig(1)/dz , sunorig(2)/dy, sunorig(3)/dx /)4403 sunorig_grid = (/ sunorig(1)/dz(1), sunorig(2)/dy, sunorig(3)/dx /) 4399 4404 !-- sunorig_grid = sunorig_grid / norm2(sunorig_grid) 4400 4405 sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2)) … … 4403 4408 !-- precompute effective box depth with prototype Leaf Area Density 4404 4409 pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1 4405 CALL box_absorb(CSHIFT((/dz ,dy,dx/), pc_box_dimshift), &4410 CALL box_absorb(CSHIFT((/dz(1),dy,dx/), pc_box_dimshift), & 4406 4411 60, prototype_lad, & 4407 4412 CSHIFT(ABS(sunorig), pc_box_dimshift), & … … 4801 4806 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 4802 4807 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 4803 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz 4808 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1) 4804 4809 ENDDO 4805 4810 ! … … 5014 5019 SUBROUTINE radiation_interaction_init 5015 5020 5021 USE control_parameters, & 5022 ONLY: dz_stretch_level_start 5023 5016 5024 USE netcdf_data_input_mod, & 5017 5025 ONLY: leaf_area_density_f … … 5095 5103 nzpt = nzptl 5096 5104 #endif 5105 ! 5106 !-- Stretching (non-uniform grid spacing) is not considered in the radiation 5107 !-- model. Therefore, vertical stretching has to be applied above the area 5108 !-- where the parts of the radiation model which assume constant grid spacing 5109 !-- are active. ABS (...) is required because the default value of 5110 !-- dz_stretch_level_start is -9999999.9_wp (negative). 5111 IF ( ABS( dz_stretch_level_start(1) ) <= zw(nzut) ) THEN 5112 WRITE( message_string, * ) 'The lowest level where vertical ', & 5113 'stretching is applied have to be ', & 5114 'greater than ', zw(nzut) 5115 CALL message( 'init_grid', 'PA0496', 1, 2, 0, 6, 0 ) 5116 ENDIF 5097 5117 ! 5098 5118 !-- global number of urban and plant layers … … 5390 5410 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 5391 5411 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 5392 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz 5412 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1) 5393 5413 ENDDO 5394 5414 … … 5512 5532 5513 5533 !-- unit vector source -> target 5514 uv = (/ (ta(1)-sa(1))*dz , (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)5534 uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /) 5515 5535 sqdist = SUM(uv(:)**2) 5516 5536 uv = uv / SQRT(sqdist) … … 5677 5697 5678 5698 !-- unit vector source -> target 5679 uv = (/ (ta(1)-sa(1))*dz , (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)5699 uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /) 5680 5700 sqdist = SUM(uv(:)**2) 5681 5701 uv = uv / SQRT(sqdist) … … 6127 6147 ENDIF 6128 6148 uvect(:) = delta(:) / distance 6129 realdist = SQRT(SUM( (uvect(:)*(/dz ,dy,dx/))**2 ))6149 realdist = SQRT(SUM( (uvect(:)*(/dz(1),dy,dx/))**2 )) 6130 6150 6131 6151 lastdist = 0._wp … … 6509 6529 zb1 = CEILING(zexit * zsgn - .5_wp) - 1 ! because it must be smaller than exit 6510 6530 nz = MAX(zb1 - zb0 + 3, 2) 6511 rt2_dist(nz) = SQRT(((zexit-zorig)*dz )**2 + dxxyy)6531 rt2_dist(nz) = SQRT(((zexit-zorig)*dz(1))**2 + dxxyy) 6512 6532 qdist = rt2_dist(nz) / (zexit-zorig) 6513 6533 rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/) … … 6567 6587 zb1 = CEILING(zexit * zsgn - .5_wp) - 1 ! because it must be smaller than exit 6568 6588 nz = MAX(zb1 - zb0 + 3, 2) 6569 rt2_dist(nz) = SQRT(((zexit-zorig)*dz )**2 + dxxyy)6589 rt2_dist(nz) = SQRT(((zexit-zorig)*dz(1))**2 + dxxyy) 6570 6590 qdist = rt2_dist(nz) / (zexit-zorig) 6571 6591 rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/) -
palm/trunk/SOURCE/read_restart_data_mod.f90
r3056 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! New parameters concerning vertical grid stretching have been added 28 ! 29 ! 3056 2018-06-04 07:49:35Z Giersch 27 30 ! found variable has to be set to false inside overlap loop 28 31 ! … … 407 410 CASE ( 'dz_stretch_factor' ) 408 411 READ ( 13 ) dz_stretch_factor 412 CASE ( 'dz_stretch_factor_array' ) 413 READ ( 13 ) dz_stretch_factor_array 409 414 CASE ( 'dz_stretch_level' ) 410 415 READ ( 13 ) dz_stretch_level 416 CASE ( 'dz_stretch_level_end' ) 417 READ ( 13 ) dz_stretch_level_end 418 CASE ( 'dz_stretch_level_start' ) 419 READ ( 13 ) dz_stretch_level_start 411 420 CASE ( 'e_min' ) 412 421 READ ( 13 ) e_min … … 1251 1260 THEN 1252 1261 WRITE( message_string, * ) 'number of PEs or virtual PE-grid changed ', & 1253 'in restart run PE', myid, ' will read from files ', &1262 'in restart run & PE', myid, ' will read from files ', & 1254 1263 file_list(1:files_to_be_opened) 1255 1264 CALL message( 'rrd_local', 'PA0285', 0, 0, 0, 6, 0 ) -
palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
r3051 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Error message related to vertical stretching has been added, dz was replaced 28 ! by dz(1) 29 ! 30 ! 3051 2018-05-30 17:43:55Z suehring 27 31 ! Bugfix in calculation of initial Reynolds-stress tensor. 28 32 ! … … 313 317 314 318 USE control_parameters, & 315 ONLY: bc_lr, bc_ns, forcing, nest_domain, rans_mode, turbulent_inflow 319 ONLY: bc_lr, bc_ns, forcing, nest_domain, number_stretch_level_start, & 320 rans_mode, turbulent_inflow 316 321 317 322 USE pmc_interface, & … … 365 370 'is not allowed' 366 371 CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 ) 372 ENDIF 373 374 IF ( number_stretch_level_start > 0 ) THEN 375 message_string = 'Using synthetic turbulence generator ' // & 376 'in combination with stretching is not allowed' 377 CALL message( 'stg_check_parameters', 'PA0420', 1, 2, 0, 6, 0 ) 367 378 ENDIF 368 379 … … 607 618 608 619 ! 609 !-- Convert length scales from meter to number of grid points 620 !-- Convert length scales from meter to number of grid points. Attention: 621 !-- Does not work if grid stretching is used 610 622 nuy(k) = INT( luy * ddy ) 611 nuz(k) = INT( luz / dz )623 nuz(k) = INT( luz / dz(1) ) 612 624 nvy(k) = INT( lvy * ddy ) 613 nvz(k) = INT( lvz / dz )625 nvz(k) = INT( lvz / dz(1) ) 614 626 nwy(k) = INT( lwy * ddy ) 615 nwz(k) = INT( lwz / dz )627 nwz(k) = INT( lwz / dz(1) ) 616 628 ! 617 629 !-- Workaround, assume isotropic turbulence -
palm/trunk/SOURCE/urban_surface_mod.f90
r3049 r3065 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Unused array dxdir was removed, dz was replaced by dzu to consider vertical 31 ! grid stretching 32 ! 33 ! 3049 2018-05-29 13:52:36Z Giersch 30 34 ! Error messages revised 31 35 ! … … 282 286 #if ! defined( __nopointer ) 283 287 USE arrays_3d, & 284 ONLY: hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend288 ONLY: dzu, hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend 285 289 #endif 286 290 … … 292 296 293 297 USE control_parameters, & 294 ONLY: coupling_start_time, dz, topography, dt_3d,&298 ONLY: coupling_start_time, topography, dt_3d, & 295 299 intermediate_timestep_count, initializing_actions, & 296 300 intermediate_timestep_count_max, simulated_time, end_time, & … … 7207 7211 REAL(wp), DIMENSION(nzb:nzt) :: exn !< value of the Exner function in layers 7208 7212 7209 REAL(wp), DIMENSION(0:4) :: dxdir !< surface normal direction gridbox length7210 7213 REAL(wp) :: dtime !< simulated time of day (in UTC) 7211 7214 INTEGER(iwp) :: dhour !< simulated hour of day (in UTC) … … 7213 7216 7214 7217 7215 dxdir = (/dz,dy,dy,dx,dx/)7216 7218 #if ! defined( __nopointer ) 7217 7219 exn(nzb:nzt) = (hyp(nzb:nzt) / 100000.0_wp )**0.286_wp !< Exner function … … 7634 7636 (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1) 7635 7637 IF ( aheat(k,j,i) > 0.0_wp ) THEN 7636 pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exn(k)*rho_cp*dz )7638 pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exn(k)*rho_cp*dzu(k)) 7637 7639 ENDIF 7638 7640 ENDIF -
palm/trunk/SOURCE/user_init_grid.f90
r2718 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! dz was replaced by dz(1) 28 ! 29 ! 2718 2018-01-02 08:49:38Z maronga 27 30 ! Corrected "Former revisions" section 28 31 ! … … 108 111 !-- topography, bit is 1 for atmospheric grid point. 109 112 !-- The following example shows how to prescribe sine-like topography 110 !-- along x-direction with amplitude of 10 * dz and wavelength 10 * dy.113 !-- along x-direction with amplitude of 10 * dz(1) and wavelength 10 * dy. 111 114 ! DO i = nxlg, nxrg 112 ! h_topo = 10.0_wp * dz * (SIN(3.14_wp*0.5_wp)*i*dx / ( 5.0_wp * dy ) )**2115 ! h_topo = 10.0_wp * dz(1) * (SIN(3.14_wp*0.5_wp)*i*dx / ( 5.0_wp * dy ) )**2 113 116 ! 114 117 ! k_topo = MINLOC( ABS( zw - h_topo ), 1 ) - 1 -
palm/trunk/SOURCE/vertical_nesting_mod.f90
r3049 r3065 26 26 ! ----------------- 27 27 ! $Id$ 28 ! dz was replaced by dz(1), error messages related to vertical grid stretching 29 ! have been added 30 ! 31 ! 3049 2018-05-29 13:52:36Z Giersch 28 32 ! Error messages revised 29 33 ! … … 73 77 !> after spin-up of the CG 74 78 !> 79 !> @todo Replace dz(1) appropriatly to account for grid stretching 75 80 !> @todo Ensure that code can be compiled for serial and parallel mode. Please 76 81 !> check the placement of the directive "__parallel". … … 3639 3644 3640 3645 USE control_parameters, & 3641 ONLY: coupling_mode, coupling_mode_remote, coupling_topology, dz 3646 ONLY: coupling_mode, coupling_mode_remote, coupling_topology, dz, & 3647 dz_stretch_level_start, message_string 3642 3648 3643 3649 USE grid_variables, & … … 3678 3684 dxc = dx 3679 3685 dyc = dy 3680 dzc = dz 3686 dzc = dz(1) 3681 3687 cg_nprocs = numprocs 3682 3688 3683 3689 IF ( myid == 0 ) THEN 3684 3690 3685 CALL MPI_SEND( nxc, 1, MPI_INTEGER , numprocs, 1, comm_inter, 3691 CALL MPI_SEND( nxc, 1, MPI_INTEGER , numprocs, 1, comm_inter, & 3686 3692 ierr ) 3687 CALL MPI_SEND( nyc, 1, MPI_INTEGER , numprocs, 2, comm_inter, 3693 CALL MPI_SEND( nyc, 1, MPI_INTEGER , numprocs, 2, comm_inter, & 3688 3694 ierr ) 3689 CALL MPI_SEND( nzc, 1, MPI_INTEGER , numprocs, 3, comm_inter, 3695 CALL MPI_SEND( nzc, 1, MPI_INTEGER , numprocs, 3, comm_inter, & 3690 3696 ierr ) 3691 CALL MPI_SEND( dxc, 1, MPI_REAL , numprocs, 4, comm_inter, 3697 CALL MPI_SEND( dxc, 1, MPI_REAL , numprocs, 4, comm_inter, & 3692 3698 ierr ) 3693 CALL MPI_SEND( dyc, 1, MPI_REAL , numprocs, 5, comm_inter, 3699 CALL MPI_SEND( dyc, 1, MPI_REAL , numprocs, 5, comm_inter, & 3694 3700 ierr ) 3695 CALL MPI_SEND( dzc, 1, MPI_REAL , numprocs, 6, comm_inter, 3701 CALL MPI_SEND( dzc, 1, MPI_REAL , numprocs, 6, comm_inter, & 3696 3702 ierr ) 3697 CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 7, comm_inter, 3703 CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 7, comm_inter, & 3698 3704 ierr ) 3699 3705 CALL MPI_SEND( cg_nprocs, 1, MPI_INTEGER, numprocs, 8, comm_inter, & 3700 3706 ierr ) 3701 CALL MPI_RECV( nxf, 1, MPI_INTEGER, numprocs, 21, comm_inter, 3707 CALL MPI_RECV( nxf, 1, MPI_INTEGER, numprocs, 21, comm_inter, & 3702 3708 status, ierr ) 3703 CALL MPI_RECV( nyf, 1, MPI_INTEGER, numprocs, 22, comm_inter, 3709 CALL MPI_RECV( nyf, 1, MPI_INTEGER, numprocs, 22, comm_inter, & 3704 3710 status, ierr ) 3705 CALL MPI_RECV( nzf, 1, MPI_INTEGER, numprocs, 23, comm_inter, 3711 CALL MPI_RECV( nzf, 1, MPI_INTEGER, numprocs, 23, comm_inter, & 3706 3712 status, ierr ) 3707 CALL MPI_RECV( dxf, 1, MPI_REAL, numprocs, 24, comm_inter, 3713 CALL MPI_RECV( dxf, 1, MPI_REAL, numprocs, 24, comm_inter, & 3708 3714 status, ierr ) 3709 CALL MPI_RECV( dyf, 1, MPI_REAL, numprocs, 25, comm_inter, 3715 CALL MPI_RECV( dyf, 1, MPI_REAL, numprocs, 25, comm_inter, & 3710 3716 status, ierr ) 3711 CALL MPI_RECV( dzf, 1, MPI_REAL, numprocs, 26, comm_inter, 3717 CALL MPI_RECV( dzf, 1, MPI_REAL, numprocs, 26, comm_inter, & 3712 3718 status, ierr ) 3713 CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, 3719 CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, & 3714 3720 numprocs, 27, comm_inter, status, ierr ) 3715 CALL MPI_RECV( fg_nprocs, 1, MPI_INTEGER, &3721 CALL MPI_RECV( fg_nprocs, 1, MPI_INTEGER, & 3716 3722 numprocs, 28, comm_inter, status, ierr ) 3717 3723 ENDIF … … 3725 3731 CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr ) 3726 3732 CALL MPI_BCAST( fg_nprocs, 1, MPI_INTEGER, 0, comm2d, ierr ) 3733 3734 ! 3735 !-- Check if stretching is used within the nested domain. ABS(...) is 3736 !-- necessary because of the default value of -9999999.9_wp (negative) 3737 IF ( ABS( dz_stretch_level_start(1) ) <= (nzf+1)*dzf ) THEN 3738 message_string = 'Stretching in the parent domain is '// & 3739 'only allowed above the nested domain' 3740 CALL message( 'vertical_nesting_mod', 'PA0497', 1, 2, 0, 6, 0 ) 3741 ENDIF 3727 3742 3728 3743 ELSEIF ( coupling_mode == 'vnested_fine' ) THEN … … 3733 3748 dxf = dx 3734 3749 dyf = dy 3735 dzf = dz 3750 dzf = dz(1) 3736 3751 fg_nprocs = numprocs 3737 3752 … … 3774 3789 3775 3790 ENDIF 3776 3791 3777 3792 ngp_c = ( nxc+1 + 2 * nbgp ) * ( nyc+1 + 2 * nbgp ) 3778 3793 ngp_f = ( nxf+1 + 2 * nbgp ) * ( nyf+1 + 2 * nbgp ) … … 3910 3925 3911 3926 #if defined( __parallel ) 3912 USE arrays_3d, 3927 USE arrays_3d, & 3913 3928 ONLY: zu, zw 3914 3929 3915 USE control_parameters, 3916 ONLY: coupling_mode 3930 USE control_parameters, & 3931 ONLY: coupling_mode, message_string, number_stretch_level_start 3917 3932 3918 USE indices, 3933 USE indices, & 3919 3934 ONLY: nzt 3920 3935 … … 3925 3940 IMPLICIT NONE 3926 3941 3927 !-- Allocate and Exchange zuc and zuf, zwc and zwf 3942 ! 3943 !-- Allocate and Exchange zuc and zuf, zwc and zwf 3928 3944 IF ( coupling_mode(1:8) == 'vnested_' ) THEN 3929 3945 … … 3932 3948 3933 3949 IF ( coupling_mode == 'vnested_crse' ) THEN 3934 3935 zuc = zu 3936 zwc = zw 3950 3951 zuc = zu 3952 zwc = zw 3953 3937 3954 IF ( myid == 0 ) THEN 3938 3955 … … 3954 3971 ELSEIF ( coupling_mode == 'vnested_fine' ) THEN 3955 3972 3956 zuf = zu 3957 zwf = zw 3973 ! 3974 !-- Check if stretching is used within the nested domain 3975 IF ( number_stretch_level_start > 0 ) THEN 3976 message_string = 'Stretching in the nested domain is not '//& 3977 'allowed' 3978 CALL message( 'vertical_nesting_mod', 'PA0498', 1, 2, 0, 6, 0 ) 3979 ENDIF 3980 3981 zuf = zu 3982 zwf = zw 3983 3958 3984 IF ( myid == 0 ) THEN 3959 3985 -
palm/trunk/SOURCE/virtual_flight_mod.f90
r3049 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! IF statement revised to consider new vertical grid stretching 28 ! 29 ! 3049 2018-05-29 13:52:36Z Giersch 27 30 ! Error messages revised 28 31 ! … … 734 737 735 738 USE control_parameters, & 736 ONLY: dz, dz_stretch_level 739 ONLY: dz, dz_stretch_level_start 737 740 738 741 USE grid_variables, & … … 777 780 !-- if flight position is above the vertical grid-stretching level. 778 781 !-- fac=1 if variable is on scalar grid level, fac=0 for w-component. 779 IF ( z_pos(l) < dz_stretch_level ) THEN780 k = ( z_pos(l) + fac * 0.5_wp * dz ) / dz782 IF ( z_pos(l) < dz_stretch_level_start(1) ) THEN 783 k = ( z_pos(l) + fac * 0.5_wp * dz(1) ) / dz(1) 781 784 ELSE 782 785 ! -
palm/trunk/SOURCE/wind_turbine_model_mod.f90
r3049 r3065 26 26 ! ----------------- 27 27 ! $Id$ 28 ! dz was replaced by dz(1), error message concerning grid stretching has been 29 ! introduced 30 ! 31 ! 3049 2018-05-29 13:52:36Z Giersch 28 32 ! Error messages revised 29 33 ! … … 132 136 !> automatically). 133 137 !> 138 !> @todo Replace dz(1) appropriatly to account for grid stretching 134 139 !> @todo Revise code according to PALM Coding Standard 135 140 !> @todo Implement ADM and ALM turbine models … … 740 745 delta_r_factor = segment_width 741 746 delta_t_factor = segment_length 742 delta_r_init = delta_r_factor * MIN( dx, dy, dz )747 delta_r_init = delta_r_factor * MIN( dx, dy, dz(1)) 743 748 744 749 DO inot = 1, nturbines … … 937 942 938 943 944 USE control_parameters, & 945 ONLY: dz_stretch_level_start 946 939 947 IMPLICIT NONE 940 948 … … 992 1000 pi**( 1.0_wp / 6.0_wp ) * eps_kernel 993 1001 ! 1002 !-- Stretching (non-uniform grid spacing) is not considered in the wind 1003 !-- turbine model. Therefore, vertical stretching has to be applied above 1004 !-- the area where the wtm is active. ABS (...) is required because the 1005 !-- default value of dz_stretch_level_start is -9999999.9_wp (negative). 1006 IF ( ABS( dz_stretch_level_start(1) ) <= MAXVAL(rcz) + MAXVAL(rr) + & 1007 eps_min) THEN 1008 WRITE( message_string, * ) 'The lowest level where vertical ', & 1009 'stretching is applied &have to be ', & 1010 'greater than ', MAXVAL(rcz) + & 1011 MAXVAL(rr) + eps_min 1012 CALL message( 'init_grid', 'PA0484', 1, 2, 0, 6, 0 ) 1013 ENDIF 1014 ! 994 1015 !-- Square of eps_min: 995 1016 eps_min2 = eps_min**2 … … 1022 1043 i_hub(inot) = INT( rcx(inot) / dx ) 1023 1044 j_hub(inot) = INT( ( rcy(inot) + 0.5_wp * dy ) / dy ) 1024 k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz ) / dz)1045 k_hub(inot) = INT( ( rcz(inot) + 0.5_wp * dz(1) ) / dz(1) ) 1025 1046 1026 1047 ! … … 1031 1052 i_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dx ) 1032 1053 j_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dy ) 1033 k_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dz )1054 k_smear(inot) = CEILING( ( rr(inot) + eps_min ) / dz(1) ) 1034 1055 1035 1056 ENDDO … … 1121 1142 !-- surface. All points between z=0 and z=dz/s would already be contained 1122 1143 !-- in grid box 1. 1123 index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz ) + 11124 index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz ) + 11144 index_nacb(inot) = INT( ( rcz(inot) - rnac(inot) ) / dz(1) ) + 1 1145 index_nact(inot) = INT( ( rcz(inot) + rnac(inot) ) / dz(1) ) + 1 1125 1146 1126 1147 ! … … 1148 1169 IF ( j == tower_s ) THEN 1149 1170 tow_cd_surf(k,j,i) = ( rcz(inot) - & 1150 ( k_hub(inot) * dz - 0.5_wp * dz ) ) *& ! extension in z-direction1171 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction 1151 1172 ( ( tower_s + 1.0_wp + 0.5_wp ) * dy - & 1152 1173 ( rcy(inot) - 0.5_wp * dtow(inot) ) ) * & ! extension in y-direction … … 1154 1175 ELSEIF ( j == tower_n ) THEN 1155 1176 tow_cd_surf(k,j,i) = ( rcz(inot) - & 1156 ( k_hub(inot) * dz - 0.5_wp * dz ) ) *& ! extension in z-direction1177 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& ! extension in z-direction 1157 1178 ( ( rcy(inot) + 0.5_wp * dtow(inot) ) - & 1158 1179 ( tower_n + 0.5_wp ) * dy ) * & ! extension in y-direction … … 1163 1184 ELSE 1164 1185 tow_cd_surf(k,j,i) = ( rcz(inot) - & 1165 ( k_hub(inot) * dz - 0.5_wp * dz ) ) *&1186 ( k_hub(inot) * dz(1) - 0.5_wp * dz(1) ) )*& 1166 1187 dy * turb_cd_tower(inot) 1167 1188 ENDIF … … 1169 1190 !-- tower lies completely within one grid box: 1170 1191 ELSE 1171 tow_cd_surf(k,j,i) = ( rcz(inot) -&1172 ( k_hub(inot) * dz - 0.5_wp * dz ) ) *&1192 tow_cd_surf(k,j,i) = ( rcz(inot) - ( k_hub(inot) * & 1193 dz(1) - 0.5_wp * dz(1) ) ) * & 1173 1194 dtow(inot) * turb_cd_tower(inot) 1174 1195 ENDIF … … 1182 1203 !-- leftmost and rightmost grid box: 1183 1204 IF ( j == tower_s ) THEN 1184 tow_cd_surf(k,j,i) = dz * (&1205 tow_cd_surf(k,j,i) = dz(1) * ( & 1185 1206 ( tower_s + 1 + 0.5_wp ) * dy - & 1186 1207 ( rcy(inot) - 0.5_wp * dtow(inot) ) & 1187 1208 ) * turb_cd_tower(inot) 1188 1209 ELSEIF ( j == tower_n ) THEN 1189 tow_cd_surf(k,j,i) = dz * (&1210 tow_cd_surf(k,j,i) = dz(1) * ( & 1190 1211 ( rcy(inot) + 0.5_wp * dtow(inot) ) - & 1191 1212 ( tower_n + 0.5_wp ) * dy & … … 1195 1216 !-- (where tow_cd_surf = grid box area): 1196 1217 ELSE 1197 tow_cd_surf(k,j,i) = dz * dy * turb_cd_tower(inot) 1218 tow_cd_surf(k,j,i) = dz(1) * dy * & 1219 turb_cd_tower(inot) 1198 1220 ENDIF 1199 1221 ! 1200 1222 !-- tower lies completely within one grid box: 1201 1223 ELSE 1202 tow_cd_surf(k,j,i) = dz * dtow(inot) *&1224 tow_cd_surf(k,j,i) = dz(1) * dtow(inot) * & 1203 1225 turb_cd_tower(inot) 1204 1226 ENDIF ! end if larger than grid box … … 1302 1324 ENDDO ! end of loop over turbines 1303 1325 1304 tow_cd_surf = tow_cd_surf / ( dx * dy * dz )! Normalize tower drag1305 nac_cd_surf = nac_cd_surf / ( dx * dy * dz ) ! Normalize nacelle drag1326 tow_cd_surf = tow_cd_surf / ( dx * dy * dz(1) ) ! Normalize tower drag 1327 nac_cd_surf = nac_cd_surf / ( dx * dy * dz(1) ) ! Normalize nacelle drag 1306 1328 1307 1329 CALL wtm_read_blade_tables … … 1746 1768 ii = rbx(ring,rseg) * ddx 1747 1769 jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy 1748 kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz1770 kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 1749 1771 ! 1750 1772 !-- Interpolate only if all required information is available on … … 1776 1798 1777 1799 u_int_1_l(inot,ring,rseg) = u_int_l + & 1778 ( rbz(ring,rseg) - zu(kk) ) / dz *&1800 ( rbz(ring,rseg) - zu(kk) ) / dz(1) * & 1779 1801 ( u_int_u - u_int_l ) 1780 1802 … … 1791 1813 ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx 1792 1814 jj = rby(ring,rseg) * ddy 1793 kk = ( rbz(ring,rseg) + 0.5_wp * dz ) / dz1815 kk = ( rbz(ring,rseg) + 0.5_wp * dz(1) ) / dz(1) 1794 1816 ! 1795 1817 !-- Interpolate only if all required information is available on … … 1821 1843 1822 1844 v_int_1_l(inot,ring,rseg) = v_int_l + & 1823 ( rbz(ring,rseg) - zu(kk) ) / dz *&1845 ( rbz(ring,rseg) - zu(kk) ) / dz(1) * & 1824 1846 ( v_int_u - v_int_l ) 1825 1847 … … 1836 1858 ii = ( rbx(ring,rseg) - 0.5_wp * dx ) * ddx 1837 1859 jj = ( rby(ring,rseg) - 0.5_wp * dy ) * ddy 1838 kk = rbz(ring,rseg) / dz 1860 kk = rbz(ring,rseg) / dz(1) 1839 1861 ! 1840 1862 !-- Interpolate only if all required information is available on … … 1866 1888 1867 1889 w_int_1_l(inot,ring,rseg) = w_int_l + & 1868 ( rbz(ring,rseg) - zw(kk) ) / dz *&1890 ( rbz(ring,rseg) - zw(kk) ) / dz(1) * & 1869 1891 ( w_int_u - w_int_l ) 1870 1892 ELSE … … 2273 2295 dist_u_3d = ( i * dx - rbx(ring,rseg) )**2 + & 2274 2296 ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & 2275 ( k * dz - 0.5_wp * dz- rbz(ring,rseg) )**22297 ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2 2276 2298 dist_v_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & 2277 2299 ( j * dy - rby(ring,rseg) )**2 + & 2278 ( k * dz - 0.5_wp * dz- rbz(ring,rseg) )**22300 ( k * dz(1) - 0.5_wp * dz(1) - rbz(ring,rseg) )**2 2279 2301 dist_w_3d = ( i * dx + 0.5_wp * dx - rbx(ring,rseg) )**2 + & 2280 2302 ( j * dy + 0.5_wp * dy - rby(ring,rseg) )**2 + & 2281 ( k * dz - rbz(ring,rseg) )**22303 ( k * dz(1) - rbz(ring,rseg) )**2 2282 2304 2283 2305 ! -
palm/trunk/SOURCE/write_restart_data_mod.f90
r3004 r3065 25 25 ! ----------------- 26 26 ! $Id$ 27 ! New parameters concerning vertical grid stretching have been added 28 ! 29 ! 3004 2018-04-27 12:33:25Z Giersch 27 30 ! precipitation_rate_av removed 28 31 ! … … 380 383 CALL wrd_write_string( 'dz' ) 381 384 WRITE ( 14 ) dz 382 385 383 386 CALL wrd_write_string( 'dz_max' ) 384 387 WRITE ( 14 ) dz_max … … 386 389 CALL wrd_write_string( 'dz_stretch_factor' ) 387 390 WRITE ( 14 ) dz_stretch_factor 388 391 392 CALL wrd_write_string( 'dz_stretch_factor_array' ) 393 WRITE ( 14 ) dz_stretch_factor_array 394 389 395 CALL wrd_write_string( 'dz_stretch_level' ) 390 396 WRITE ( 14 ) dz_stretch_level 391 397 398 CALL wrd_write_string( 'dz_stretch_level_end' ) 399 WRITE ( 14 ) dz_stretch_level_end 400 401 CALL wrd_write_string( 'dz_stretch_level_start' ) 402 WRITE ( 14 ) dz_stretch_level_start 403 392 404 CALL wrd_write_string( 'e_min' ) 393 405 WRITE ( 14 ) e_min
Note: See TracChangeset
for help on using the changeset viewer.
