Changeset 3065 for palm/trunk/SOURCE/init_grid.f90
- Timestamp:
- Jun 12, 2018 7:03:02 AM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_grid.f90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.
