Ignore:
Timestamp:
Jun 12, 2018 7:03:02 AM (3 years ago)
Author:
Giersch
Message:

New vertical stretching procedure has been introduced

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/init_grid.f90

    r3051 r3065  
    2525! -----------------
    2626! $Id$
     27! New vertical stretching mechanism introduced
     28!
     29! 3051 2018-05-30 17:43:55Z suehring
    2730! Minor bugfix concerning mapping 3D buildings on top of terrain
    2831!
     
    325328               canyon_height, canyon_wall_left, canyon_wall_south,             &
    326329               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,                             &
    329334               force_bound_l, force_bound_r, force_bound_n, force_bound_s,     &
    330335               ibc_uv_b, inflow_l, inflow_n, inflow_r, inflow_s,               &
    331336               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
    338345         
    339346    USE grid_variables,                                                        &
     
    361368    IMPLICIT NONE
    362369
    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
    370378                                     
    371     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local      !< index for topography top at cell-center
    372     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp        !< dummy to calculate topography indices on u- and v-grid
     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
    373381
    374382    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
    375383
     384    REAL(wp) ::  dz_level_end  !< distance between calculated height level for u/v-grid and user-specified end level for stretching
    376385    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
    377388
    378389
     
    391402!
    392403!-- Compute height of u-levels from constant grid length and dz stretch factors
    393     IF ( dz == -1.0_wp )  THEN
     404    IF ( dz(1) == -1.0_wp )  THEN
    394405       message_string = 'missing dz'
    395406       CALL message( 'init_grid', 'PA0200', 1, 2, 0, 6, 0 )
    396     ELSEIF ( dz <= 0.0_wp )  THEN
    397        WRITE( message_string, * ) 'dz=',dz,' <= 0.0'
     407    ELSEIF ( dz(1) <= 0.0_wp )  THEN
     408       WRITE( message_string, * ) 'dz=',dz(1),' <= 0.0'
    398409       CALL message( 'init_grid', 'PA0201', 1, 2, 0, 6, 0 )
    399410    ENDIF
    400411
    401412!
     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!
    402488!-- Define the vertical grid levels
    403489    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
    404546!
    405547!--    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).
    406550!--    The second u-level (k=1) corresponds to the top of the
    407551!--    Prandtl-layer.
    408 
    409552       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN
    410553          zu(0) = 0.0_wp
    411       !    zu(0) = - dz * 0.5_wp
    412554       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.
    419570       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         
    425586          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
    426598       ENDDO
    427599
     
    438610
    439611    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
    440668!
    441669!--    Grid for ocean with free water surface is at k=nzt (w-grid).
     
    444672!--    w-level are defined at same height, but staggered from the second level.
    445673!--    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
    451686       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         
    461702          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
    462714       ENDDO
    463 
     715       
    464716!
    465717!--    Compute the w-levels. They are always staggered half-way between the
     
    468720!--    same height. The top w-level (nzt+1) is not used but set for
    469721!--    consistency, since w and all scalar variables are defined up tp nzt+1.
    470        zw(nzt+1) = dz
     722       zw(nzt+1) = dz(1)
    471723       zw(nzt)   = 0.0_wp
    472724       DO  k = 0, nzt
     
    8551107 END SUBROUTINE init_grid
    8561108
     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 
    8571239! Description:
    8581240! -----------------------------------------------------------------------------!
     
    17032085!--       Tunnel-wall depth
    17042086          IF ( tunnel_wall_depth == 9999999.9_wp )  THEN 
    1705              td = MAX ( dx, dy, dz )
     2087             td = MAX ( dx, dy, dz(1) )
    17062088          ELSE
    17072089             td = tunnel_wall_depth
Note: See TracChangeset for help on using the changeset viewer.