Ignore:
Timestamp:
Jan 27, 2020 3:07:30 PM (4 years ago)
Author:
Giersch
Message:

Allocation statements, comments, naming of variables revised and _wp added to real type values

File:
1 edited

Legend:

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

    r4360 r4386  
    2525! -----------------
    2626! $Id$
     27! Allocation statements, comments, naming of variables revised and _wp added to
     28! real type values
     29!
     30! 4360 2020-01-07 11:25:50Z suehring
    2731! Revise error messages for generic tunnel setup.
    2832!
     
    199203!
    200204!-- Allocate grid arrays
    201     ALLOCATE( x(0:nx), xu(0:nx) )
     205    ALLOCATE( x(0:nx) )
     206    ALLOCATE( xu(0:nx) )
     207   
    202208    DO i = 0, nx
    203209       xu(i) = i * dx
     
    205211    ENDDO
    206212
    207     ALLOCATE( y(0:ny), yv(0:ny) )
     213    ALLOCATE( y(0:ny) )
     214    ALLOCATE( yv(0:ny) )
     215   
    208216    DO j = 0, ny
    209217       yv(j) = j * dy
     
    211219    ENDDO
    212220
    213 !
    214 !-- Allocate grid arrays
    215     ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1),        &
    216               dzw(1:nzt+1), zu(nzb:nzt+1), zw(nzb:nzt+1) )
    217 
    218 !
    219 !-- Compute height of u-levels from constant grid length and dz stretch factors
     221    ALLOCATE( ddzu(1:nzt+1) )
     222    ALLOCATE( ddzw(1:nzt+1) )
     223    ALLOCATE( dd2zu(1:nzt) )
     224    ALLOCATE( dzu(1:nzt+1) )
     225    ALLOCATE( dzw(1:nzt+1) )
     226    ALLOCATE( zu(nzb:nzt+1) )
     227    ALLOCATE( zw(nzb:nzt+1) )
     228
     229!
     230!-- For constructing an appropriate grid, the vertical grid spacing dz has to
     231!-- be specified with a non-negative value in the parameter file
    220232    IF ( dz(1) == -1.0_wp )  THEN
    221233       message_string = 'missing dz'
     
    248260
    249261!
    250 !-- The number of specified end levels +1 has to be the same than the number
     262!-- The number of specified end levels +1 has to be the same as the number
    251263!-- of specified dz values
    252264    IF ( number_dz /= number_stretch_level_end + 1 ) THEN
    253        WRITE( message_string, * ) 'The number of values for dz = ',         &
    254                                    number_dz, 'has to be the same than& ',  &
    255                                    'the number of values for ',             &
    256                                    'dz_stretch_level_end + 1 = ',           &
     265       WRITE( message_string, * ) 'The number of values for dz = ',            &
     266                                   number_dz, 'has to be the same as& ',       &
     267                                   'the number of values for ',                &
     268                                   'dz_stretch_level_end + 1 = ',              &
    257269                                   number_stretch_level_end+1
    258270          CALL message( 'init_grid', 'PA0156', 1, 2, 0, 6, 0 )
     
    260272   
    261273!
    262 !--    The number of specified start levels has to be the same or one less than
    263 !--    the number of specified dz values
    264     IF ( number_dz /= number_stretch_level_start + 1 .AND.                  &
     274!-- The number of specified start levels has to be the same or one less than
     275!-- the number of specified dz values
     276    IF ( number_dz /= number_stretch_level_start + 1 .AND.                     &
    265277         number_dz /= number_stretch_level_start ) THEN
    266        WRITE( message_string, * ) 'The number of values for dz = ',         &
    267                                    number_dz, 'has to be the same or one ', &
    268                                    'more than& the number of values for ',  &
    269                                    'dz_stretch_level_start = ',             &
     278       WRITE( message_string, * ) 'The number of values for dz = ',            &
     279                                   number_dz, 'has to be the same as or one ', &
     280                                   'more than& the number of values for ',     &
     281                                   'dz_stretch_level_start = ',                &
    270282                                   number_stretch_level_start
    271283          CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
    272284    ENDIF
    273285   
    274 !--    The number of specified start levels has to be the same or one more than
    275 !--    the number of specified end levels
    276     IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.   &
     286!-- The number of specified start levels has to be the same or one more than
     287!-- the number of specified end levels
     288    IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.      &
    277289         number_stretch_level_start /= number_stretch_level_end ) THEN
    278        WRITE( message_string, * ) 'The number of values for ',              &
    279                                   'dz_stretch_level_start = ',              &
    280                                    dz_stretch_level_start, 'has to be the ',&
    281                                    'same or one more than& the number of ', &
    282                                    'values for dz_stretch_level_end = ',    &
     290       WRITE( message_string, * ) 'The number of values for ',                 &
     291                                  'dz_stretch_level_start = ',                 &
     292                                   dz_stretch_level_start, 'has to be the ',   &
     293                                   'same or one more than& the number of ',    &
     294                                   'values for dz_stretch_level_end = ',       &
    283295                                   number_stretch_level_end
    284296          CALL message( 'init_grid', 'PA0216', 1, 2, 0, 6, 0 )
     
    306318
    307319!
    308 !-- Define the vertical grid levels
     320!-- Define the vertical grid levels. Start with atmosphere branch
    309321    IF ( .NOT. ocean_mode )  THEN
    310322   
    311323!
    312324!--    The stretching region has to be large enough to allow for a smooth
    313 !--    transition between two different grid spacings
     325!--    transition between two different grid spacings. The number 4 is an
     326!--    empirical value
    314327       DO n = 1, number_stretch_level_start
    315328          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
     
    319332       IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >      &
    320333                 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
    321              message_string= 'Eeach dz_stretch_level_end has to be larger ' // &
     334             message_string= 'Each dz_stretch_level_end has to be larger ' // &
    322335                             'than its corresponding value for &' //           &
    323336                             'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
     
    327340       
    328341!
    329 !--    Stretching must not be applied within the prandtl_layer
     342!--    Stretching must not be applied within the surface layer
    330343!--    (first two grid points). For the default case dz_stretch_level_start
    331344!--    is negative. Therefore the absolut value is checked here.
    332        IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_wp ) ) THEN
    333           WRITE( message_string, * ) 'Eeach dz_stretch_level_start has to be ',&
     345       IF ( ANY( ABS( dz_stretch_level_start ) <= dz(1) * 1.5_wp ) ) THEN
     346          WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
    334347                                     'larger than ', dz(1) * 1.5
    335348             CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
     
    338351!
    339352!--    The stretching has to start and end on a grid level. Therefore
    340 !--    user-specified values have to ''interpolate'' to the next lowest level
     353!--    user-specified values are mapped to the next lowest level. The 
     354!--    calculation of the first level is realized differently just because of
     355!--    historical reasons (the advanced/new stretching mechanism was realized 
     356!--    in a way that results don't change if the old parameters
     357!--    dz_stretch_level, dz_stretch_factor and dz_max are used)
    341358       IF ( number_stretch_level_start /= 0 ) THEN
    342359          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
     
    358375          ENDDO
    359376       ENDIF
    360  
     377
    361378!
    362379!--    Determine stretching factor if necessary
     
    370387!--    the first u/v- and w-level (k=0) are defined at same height (z=0).
    371388!--    The second u-level (k=1) corresponds to the top of the
    372 !--    Prandtl-layer. In case of symmetric boundaries (closed channel flow),
     389!--    surface layer. In case of symmetric boundaries (closed channel flow),
    373390!--    the first grid point is always at z=0.
    374391       IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 .OR.                              &
     
    390407
    391408!--    The default value of dz_stretch_level_start is negative, thus the first
    392 !--    condition is always true. Hence, the second condition is necessary.
     409!--    condition is true even if no stretching shall be applied. Hence, the
     410!--    second condition is also necessary.
    393411       DO  k = 2, nzt+1-symmetry_flag
    394412          IF ( dz_stretch_level_start(n) <= zu(k-1) .AND.                      &
     
    452470       ENDIF
    453471
    454     ELSE
     472    ELSE !ocean branch
    455473
    456474!
    457475!--    The stretching region has to be large enough to allow for a smooth
    458 !--    transition between two different grid spacings
     476!--    transition between two different grid spacings. The number 4 is an
     477!--    empirical value
    459478       DO n = 1, number_stretch_level_start
    460479          min_dz_stretch_level_end(n) = dz_stretch_level_start(n) -            &
     
    464483       IF ( ANY( min_dz_stretch_level_end (1:number_stretch_level_start) <     &
    465484                 dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
    466              message_string= 'Eeach dz_stretch_level_end has to be less ' //   &
     485             message_string= 'Each dz_stretch_level_end has to be less ' //   &
    467486                             'than its corresponding value for &' //           &
    468487                             'dz_stretch_level_start - 4*MAX(dz(n),dz(n+1)) '//&
     
    474493!--    Stretching must not be applied close to the surface (last two grid
    475494!--    points). For the default case dz_stretch_level_start is negative.
    476        IF ( ANY( dz_stretch_level_start > - dz(1) * 1.5_wp ) ) THEN
    477           WRITE( message_string, * ) 'Eeach dz_stretch_level_start has to be ',&
    478                                      'less than ', dz(1) * 1.5
     495       IF ( ANY( dz_stretch_level_start >= - dz(1) * 1.5_wp ) ) THEN
     496          WRITE( message_string, * ) 'Each dz_stretch_level_start has to be ',&
     497                                     'less than ', -dz(1) * 1.5
    479498             CALL message( 'init_grid', 'PA0226', 1, 2, 0, 6, 0 )
    480499       ENDIF
     
    482501!
    483502!--    The stretching has to start and end on a grid level. Therefore
    484 !--    user-specified values have to ''interpolate'' to the next highest level
     503!--    user-specified values are mapped to the next highest level. The 
     504!--    calculation of the first level is realized differently just because of
     505!--    historical reasons (the advanced/new stretching mechanism was realized 
     506!--    in a way that results don't change if the old parameters
     507!--    dz_stretch_level, dz_stretch_factor and dz_max are used)
    485508       IF ( number_stretch_level_start /= 0 ) THEN
    486509          dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) +        &
     
    514537!--    below the first w-level (k=0). In case of dirichlet bc the first u- and
    515538!--    w-level are defined at same height, but staggered from the second level.
    516 !--    The second u-level (k=1) corresponds to the top of the Prandtl-layer.
     539!--    The second u-level (k=1) corresponds to the top of the surface layer.
    517540!--    z values are negative starting from z=0 (surface)
    518541       zu(nzt+1) =   dz(1) * 0.5_wp
     
    576599       ENDIF
    577600
    578     ENDIF
     601    ENDIF !End of defining the vertical grid levels
    579602
    580603!
     
    617640    topo = 0
    618641!
    619 !-- Initialize topography by generic topography or read from topography from file. 
     642!-- Initialize topography by generic topography or read topography from file. 
    620643    CALL init_topo( topo )
    621644!
     
    691714!--    Allocate and set the arrays containing the topography height
    692715       IF ( nxr == nx  .AND.  nyn /= ny )  THEN
    693           ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn),                             &
    694                     zw_w_inner(nxl:nxr+1,nys:nyn) )
     716          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn) )                             
     717          ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn) )
    695718       ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
    696           ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1),                             &
    697                     zw_w_inner(nxl:nxr,nys:nyn+1) )
     719          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1) )                             
     720          ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn+1) )
    698721       ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
    699           ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1),                           &
    700                     zw_w_inner(nxl:nxr+1,nys:nyn+1) )
     722          ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1) )
     723          ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn+1) )
    701724       ELSE
    702           ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn),                               &
    703                     zw_w_inner(nxl:nxr,nys:nyn) )
     725          ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn) )
     726          ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn) )
    704727       ENDIF
    705728
     
    731754!-- Allocate outer and inner index arrays for topography and set
    732755!-- defaults.                   
    733     ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg),                                &
    734               nzb_s_outer(nysg:nyng,nxlg:nxrg),                                &
    735               nzb_u_inner(nysg:nyng,nxlg:nxrg),                                &
    736               nzb_u_outer(nysg:nyng,nxlg:nxrg),                                &
    737               nzb_v_inner(nysg:nyng,nxlg:nxrg),                                &
    738               nzb_v_outer(nysg:nyng,nxlg:nxrg),                                &
    739               nzb_w_inner(nysg:nyng,nxlg:nxrg),                                &
    740               nzb_w_outer(nysg:nyng,nxlg:nxrg),                                &
    741               nzb_diff_s_inner(nysg:nyng,nxlg:nxrg),                           &
    742               nzb_diff_s_outer(nysg:nyng,nxlg:nxrg),                           &
    743               nzb_local(nysg:nyng,nxlg:nxrg),                                  &
    744               nzb_tmp(nysg:nyng,nxlg:nxrg) )
     756    ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg) )                               
     757    ALLOCATE( nzb_s_outer(nysg:nyng,nxlg:nxrg) )                               
     758    ALLOCATE( nzb_u_inner(nysg:nyng,nxlg:nxrg) )                               
     759    ALLOCATE( nzb_u_outer(nysg:nyng,nxlg:nxrg) )                               
     760    ALLOCATE( nzb_v_inner(nysg:nyng,nxlg:nxrg) )                               
     761    ALLOCATE( nzb_v_outer(nysg:nyng,nxlg:nxrg) )                               
     762    ALLOCATE( nzb_w_inner(nysg:nyng,nxlg:nxrg) )                               
     763    ALLOCATE( nzb_w_outer(nysg:nyng,nxlg:nxrg) )                               
     764    ALLOCATE( nzb_diff_s_inner(nysg:nyng,nxlg:nxrg) )                         
     765    ALLOCATE( nzb_diff_s_outer(nysg:nyng,nxlg:nxrg) )                         
     766    ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) )                                 
     767    ALLOCATE( nzb_tmp(nysg:nyng,nxlg:nxrg) )
    745768!
    746769!-- Initialize 2D-index arrays. Note, these will be removed soon!
     
    9961019    REAL(wp) ::  dz_stretch_factor_array_2(9) = 1.08_wp  !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1
    9971020   
    998     REAL(wp), PARAMETER ::  stretch_factor_interval = 1.0E-06  !< interval for sampling possible stretching factors
    999     REAL(wp), PARAMETER ::  stretch_factor_lower_limit = 0.88  !< lowest possible stretching factor
    1000     REAL(wp), PARAMETER ::  stretch_factor_upper_limit = 1.12  !< highest possible stretching factor
     1021    REAL(wp), PARAMETER ::  stretch_factor_interval = 1.0E-06_wp  !< interval for sampling possible stretching factors
     1022    REAL(wp), PARAMETER ::  stretch_factor_lower_limit = 0.88_wp  !< lowest possible stretching factor
     1023    REAL(wp), PARAMETER ::  stretch_factor_upper_limit = 1.12_wp  !< highest possible stretching factor
    10011024 
    10021025 
     
    10051028   
    10061029       iterations = 1
    1007        stretch_factor_1 = 1.0
    1008        stretch_factor_2 = 1.0
    1009        delta_total_old = 1.0
     1030       stretch_factor_1 = 1.0_wp
     1031       stretch_factor_2 = 1.0_wp
     1032       delta_total_old = 1.0_wp
    10101033       
     1034!
     1035!--    First branch for stretching from rough to fine
    10111036       IF ( dz(n) > dz(n+1) ) THEN
    10121037          DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit )
    10131038             
    1014              stretch_factor_1 = 1.0 - iterations * stretch_factor_interval
    1015              distance = ABS( dz_stretch_level_end(n) -                   &
    1016                         dz_stretch_level_start(n) )
    1017              numerator = distance*stretch_factor_1/dz(n) +               &
     1039             stretch_factor_1 = 1.0_wp - iterations * stretch_factor_interval
     1040             distance = ABS( dz_stretch_level_end(n) -                         &
     1041                        dz_stretch_level_start(n) )  
     1042             numerator = distance*stretch_factor_1/dz(n) +                     &
    10181043                         stretch_factor_1 - distance/dz(n)
    10191044             
    1020              IF ( numerator > 0.0 ) THEN
    1021                 l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
     1045             IF ( numerator > 0.0_wp ) THEN
     1046                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
    10221047                l_rounded = NINT( l )
    10231048                delta_l = ABS( l_rounded - l ) / l
     
    10261051             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
    10271052             
    1028              delta_stretch_factor = ABS( stretch_factor_1 -              &
    1029                                          stretch_factor_2 ) /            &
     1053             delta_stretch_factor = ABS( stretch_factor_1 -                    &
     1054                                         stretch_factor_2 ) /                  &
    10301055                                    stretch_factor_2
    10311056             
     
    10331058
    10341059!
    1035 !--                stretch_factor_1 is taken to guarantee that the stretching
    1036 !--                procedure ends as close as possible to dz_stretch_level_end.
    1037 !--                stretch_factor_2 would guarantee that the stretched dz(n) is
    1038 !--                equal to dz(n+1) after l_rounded grid levels.
     1060!--          stretch_factor_1 is taken to guarantee that the stretching
     1061!--          procedure ends as close as possible to dz_stretch_level_end.
     1062!--          stretch_factor_2 would guarantee that the stretched dz(n) is
     1063!--          equal to dz(n+1) after l_rounded grid levels.
    10391064             IF (delta_total_new < delta_total_old) THEN
    10401065                dz_stretch_factor_array(n) = stretch_factor_1
     
    10461071           
    10471072          ENDDO
    1048              
     1073
     1074!
     1075!--    Second branch for stretching from fine to rough
    10491076       ELSEIF ( dz(n) < dz(n+1) ) THEN
    10501077          DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
    10511078                     
    1052              stretch_factor_1 = 1.0 + iterations * stretch_factor_interval
    1053              distance = ABS( dz_stretch_level_end(n) -                      &
     1079             stretch_factor_1 = 1.0_wp + iterations * stretch_factor_interval
     1080             distance = ABS( dz_stretch_level_end(n) -                         &
    10541081                        dz_stretch_level_start(n) )
    1055              numerator = distance*stretch_factor_1/dz(n) +                  &
     1082             numerator = distance*stretch_factor_1/dz(n) +                     &
    10561083                         stretch_factor_1 - distance/dz(n)
    10571084             
    1058              l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
     1085             l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0_wp
    10591086             l_rounded = NINT( l )
    10601087             delta_l = ABS( l_rounded - l ) / l
     
    10621089             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
    10631090
    1064              delta_stretch_factor = ABS( stretch_factor_1 -                 &
    1065                                         stretch_factor_2 ) /                &
     1091             delta_stretch_factor = ABS( stretch_factor_1 -                    &
     1092                                        stretch_factor_2 ) /                   &
    10661093                                        stretch_factor_2
    10671094             
     
    10691096             
    10701097!
    1071 !--                stretch_factor_1 is taken to guarantee that the stretching
    1072 !--                procedure ends as close as possible to dz_stretch_level_end.
    1073 !--                stretch_factor_2 would guarantee that the stretched dz(n) is
    1074 !--                equal to dz(n+1) after l_rounded grid levels.
     1098!--          stretch_factor_1 is taken to guarantee that the stretching
     1099!--          procedure ends as close as possible to dz_stretch_level_end.
     1100!--          stretch_factor_2 would guarantee that the stretched dz(n) is
     1101!--          equal to dz(n+1) after l_rounded grid levels.
    10751102             IF (delta_total_new < delta_total_old) THEN
    10761103                dz_stretch_factor_array(n) = stretch_factor_1
     
    17941821    IMPLICIT NONE
    17951822
    1796     INTEGER(iwp) ::  bh            !< temporary vertical index of building height
    1797     INTEGER(iwp) ::  blx           !< grid point number of building size along x
    1798     INTEGER(iwp) ::  bly           !< grid point number of building size along y
    1799     INTEGER(iwp) ::  bxl           !< index for left building wall
    1800     INTEGER(iwp) ::  bxr          !< index for right building wall
    1801     INTEGER(iwp) ::  byn          !< index for north building wall
    1802     INTEGER(iwp) ::  bys          !< index for south building wall
    1803     INTEGER(iwp) ::  ch            !< temporary vertical index for canyon height
    1804     INTEGER(iwp) ::  cwx           !< grid point number of canyon size along x
    1805     INTEGER(iwp) ::  cwy           !< grid point number of canyon size along y
    1806     INTEGER(iwp) ::  cxl           !< index for left canyon wall
    1807     INTEGER(iwp) ::  cxr          !< index for right canyon wall
    1808     INTEGER(iwp) ::  cyn          !< index for north canyon wall
    1809     INTEGER(iwp) ::  cys          !< index for south canyon wall
    1810     INTEGER(iwp) ::  i             !< index variable along x
    1811     INTEGER(iwp) ::  j             !< index variable along y
    1812     INTEGER(iwp) ::  k             !< index variable along z
    1813     INTEGER(iwp) ::  hv_in         !< heavyside function to model inner tunnel surface
    1814     INTEGER(iwp) ::  hv_out        !< heavyside function to model outer tunnel surface
    1815     INTEGER(iwp) ::  txe_out       !< end position of outer tunnel wall in x
    1816     INTEGER(iwp) ::  txs_out       !< start position of outer tunnel wall in x
    1817     INTEGER(iwp) ::  tye_out       !< end position of outer tunnel wall in y
    1818     INTEGER(iwp) ::  tys_out       !< start position of outer tunnel wall in y
    1819     INTEGER(iwp) ::  txe_in        !< end position of inner tunnel wall in x
    1820     INTEGER(iwp) ::  txs_in        !< start position of inner tunnel wall in x
    1821     INTEGER(iwp) ::  tye_in        !< end position of inner tunnel wall in y
    1822     INTEGER(iwp) ::  tys_in        !< start position of inner tunnel wall in y
    1823     INTEGER(iwp) ::  td            !< tunnel wall depth
    1824     INTEGER(iwp) ::  th            !< height of outer tunnel wall
     1823    INTEGER(iwp) ::  bh                !< temporary vertical index of building height
     1824    INTEGER(iwp) ::  ngp_bx            !< grid point number of building size along x
     1825    INTEGER(iwp) ::  ngp_by            !< grid point number of building size along y
     1826    INTEGER(iwp) ::  index_left_bwall  !< index for left building wall
     1827    INTEGER(iwp) ::  index_right_bwall !< index for right building wall
     1828    INTEGER(iwp) ::  index_north_bwall !< index for north building wall
     1829    INTEGER(iwp) ::  index_south_bwall !< index for south building wall
     1830    INTEGER(iwp) ::  ch                !< temporary vertical index for canyon height
     1831    INTEGER(iwp) ::  ngp_cx            !< grid point number of canyon size along x
     1832    INTEGER(iwp) ::  ngp_cy            !< grid point number of canyon size along y
     1833    INTEGER(iwp) ::  index_left_cwall  !< index for left canyon wall
     1834    INTEGER(iwp) ::  index_right_cwall !< index for right canyon wall
     1835    INTEGER(iwp) ::  index_north_cwall !< index for north canyon wall
     1836    INTEGER(iwp) ::  index_south_cwall !< index for south canyon wall
     1837    INTEGER(iwp) ::  i                 !< index variable along x
     1838    INTEGER(iwp) ::  j                 !< index variable along y
     1839    INTEGER(iwp) ::  k                 !< index variable along z
     1840    INTEGER(iwp) ::  hv_in             !< heavyside function to model inner tunnel surface
     1841    INTEGER(iwp) ::  hv_out            !< heavyside function to model outer tunnel surface
     1842    INTEGER(iwp) ::  txe_out           !< end position of outer tunnel wall in x
     1843    INTEGER(iwp) ::  txs_out           !< start position of outer tunnel wall in x
     1844    INTEGER(iwp) ::  tye_out           !< end position of outer tunnel wall in y
     1845    INTEGER(iwp) ::  tys_out           !< start position of outer tunnel wall in y
     1846    INTEGER(iwp) ::  txe_in            !< end position of inner tunnel wall in x
     1847    INTEGER(iwp) ::  txs_in            !< start position of inner tunnel wall in x
     1848    INTEGER(iwp) ::  tye_in            !< end position of inner tunnel wall in y
     1849    INTEGER(iwp) ::  tys_in            !< start position of inner tunnel wall in y
     1850    INTEGER(iwp) ::  td                !< tunnel wall depth
     1851    INTEGER(iwp) ::  th                !< height of outer tunnel wall
    18251852
    18261853    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local         !< index for topography top at cell-center
     
    18291856!-- Check for correct setting of the namelist parameter topography. If
    18301857!-- topography information is read from file but topography = 'flat',
    1831 !-- initialization does not properly.
     1858!-- initialization does not work properly.
    18321859    IF ( ( buildings_f%from_file  .OR.  terrain_height_f%from_file )  .AND.    &
    18331860           TRIM( topography ) /= 'read_from_file' )  THEN
     
    18601887!--       Single rectangular building, by default centered in the middle of the
    18611888!--       total domain
    1862           blx = NINT( building_length_x / dx )
    1863           bly = NINT( building_length_y / dy )
     1889          ngp_bx = NINT( building_length_x / dx )
     1890          ngp_by = NINT( building_length_y / dy )
    18641891          bh  = MINLOC( ABS( zw - building_height ), 1 ) - 1
    18651892          IF ( ABS( zw(bh)   - building_height ) == &
    18661893               ABS( zw(bh+1) - building_height )    )  bh = bh + 1
    18671894          IF ( building_wall_left == 9999999.9_wp )  THEN
    1868              building_wall_left = ( nx + 1 - blx ) / 2 * dx
    1869           ENDIF
    1870           bxl = NINT( building_wall_left / dx )
    1871           bxr = bxl + blx
     1895             building_wall_left = ( nx + 1 - ngp_bx ) / 2 * dx
     1896          ENDIF
     1897          index_left_bwall = NINT( building_wall_left / dx )
     1898          index_right_bwall = index_left_bwall + ngp_bx
    18721899
    18731900          IF ( building_wall_south == 9999999.9_wp )  THEN
    1874               building_wall_south = ( ny + 1 - bly ) / 2 * dy
    1875           ENDIF
    1876           bys = NINT( building_wall_south / dy )
    1877           byn = bys + bly
     1901              building_wall_south = ( ny + 1 - ngp_by ) / 2 * dy
     1902          ENDIF
     1903          index_south_bwall = NINT( building_wall_south / dy )
     1904          index_north_bwall = index_south_bwall + ngp_by
    18781905
    18791906!
    18801907!--       Building size has to meet some requirements
    1881           IF ( ( bxl < 1 ) .OR. ( bxr > nx-1 ) .OR. ( bxr < bxl+3 ) .OR.       &
    1882                ( bys < 1 ) .OR. ( byn > ny-1 ) .OR. ( byn < bys+3 ) )  THEN
     1908          IF ( ( index_left_bwall < 1 ) .OR. ( index_right_bwall > nx-1 ) .OR. &
     1909               ( index_right_bwall < index_left_bwall+3 ) .OR.                 &
     1910               ( index_south_bwall < 1 ) .OR. ( index_north_bwall > ny-1 ) .OR.&
     1911               ( index_north_bwall < index_south_bwall+3 ) )  THEN
    18831912             WRITE( message_string, * ) 'inconsistent building parameters:',   &
    1884                                       '&bxl=', bxl, 'bxr=', bxr, 'bys=', bys,  &
    1885                                       'byn=', byn, 'nx=', nx, 'ny=', ny
     1913                                      '&index_left_bwall=', index_left_bwall,  &
     1914                                      'index_right_bwall=', index_right_bwall, &
     1915                                      'index_south_bwall=', index_south_bwall, &
     1916                                      'index_north_bwall=', index_north_bwall, &
     1917                                      'nx=', nx, 'ny=', ny
    18861918             CALL message( 'init_grid', 'PA0203', 1, 2, 0, 6, 0 )
    18871919          ENDIF
     
    18911923!
    18921924!--       Define the building.
    1893           IF ( bxl <= nxr  .AND.  bxr >= nxl  .AND.                            &
    1894                bys <= nyn  .AND.  byn >= nys )                                 &
    1895              nzb_local(MAX(nys,bys):MIN(nyn,byn),MAX(nxl,bxl):MIN(nxr,bxr)) = bh
     1925          IF ( index_left_bwall <= nxr  .AND.  index_right_bwall >= nxl  .AND. &
     1926               index_south_bwall <= nyn  .AND.  index_north_bwall >= nys )     &
     1927             nzb_local(MAX(nys,index_south_bwall):MIN(nyn,index_north_bwall),  &
     1928                       MAX(nxl,index_left_bwall):MIN(nxr,index_right_bwall)) = bh
    18961929!
    18971930!--       Set bit array on basis of nzb_local
     
    19411974!
    19421975!--          Street canyon in y direction
    1943              cwx = NINT( canyon_width_x / dx )
     1976             ngp_cx = NINT( canyon_width_x / dx )
    19441977             IF ( canyon_wall_left == 9999999.9_wp )  THEN
    1945                 canyon_wall_left = ( nx + 1 - cwx ) / 2 * dx
    1946              ENDIF
    1947              cxl = NINT( canyon_wall_left / dx )
    1948              cxr = cxl + cwx
     1978                canyon_wall_left = ( nx + 1 - ngp_cx ) / 2 * dx
     1979             ENDIF
     1980             index_left_cwall= NINT( canyon_wall_left / dx )
     1981             index_right_cwall= index_left_cwall+ ngp_cx
    19491982          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
    19501983!
    19511984!--          Street canyon in x direction
    1952              cwy = NINT( canyon_width_y / dy )
     1985             ngp_cy = NINT( canyon_width_y / dy )
    19531986             IF ( canyon_wall_south == 9999999.9_wp )  THEN
    1954                 canyon_wall_south = ( ny + 1 - cwy ) / 2 * dy
    1955              ENDIF
    1956              cys = NINT( canyon_wall_south / dy )
    1957              cyn = cys + cwy
     1987                canyon_wall_south = ( ny + 1 - ngp_cy ) / 2 * dy
     1988             ENDIF
     1989             index_south_cwall = NINT( canyon_wall_south / dy )
     1990             index_north_cwall = index_south_cwall + ngp_cy
    19581991     
    19591992          ELSE
     
    19712004!--       Street canyon size has to meet some requirements
    19722005          IF ( canyon_width_x /= 9999999.9_wp )  THEN
    1973              IF ( ( cxl < 1 ) .OR. ( cxr > nx-1 ) .OR. ( cwx < 3 ) .OR.        &
     2006             IF ( ( index_left_cwall< 1 ) .OR. ( index_right_cwall> nx-1 ) .OR.&
     2007                  ( ngp_cx < 3 ) .OR. ( ch < 3 ) )  THEN
     2008                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
     2009                                           '&index_left_cwall=', index_left_cwall, &
     2010                                           ' index_right_cwall=', index_right_cwall, &
     2011                                           ' ngp_cx=', ngp_cx,                 &
     2012                                           ' ch=', ch, ' nx=', nx, ' ny=', ny
     2013                CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 )
     2014             ENDIF
     2015          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
     2016             IF ( ( index_south_cwall < 1 ) .OR.                               &
     2017                  ( index_north_cwall > ny-1 ) .OR. ( ngp_cy < 3 ) .OR.        &
    19742018                  ( ch < 3 ) )  THEN
    19752019                WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
    1976                                            '&cxl=', cxl, ' cxr=', cxr,         &
    1977                                            ' cwx=', cwx,                       &
    1978                                            ' ch=', ch, ' nx=', nx, ' ny=', ny
    1979                 CALL message( 'init_grid', 'PA0205', 1, 2, 0, 6, 0 )
    1980              ENDIF
    1981           ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
    1982              IF ( ( cys < 1 ) .OR. ( cyn > ny-1 ) .OR. ( cwy < 3 ) .OR.        &
    1983                   ( ch < 3 ) )  THEN
    1984                 WRITE( message_string, * ) 'inconsistent canyon parameters:',  &
    1985                                            '&cys=', cys, ' cyn=', cyn,         &
    1986                                            ' cwy=', cwy,                       &
     2020                                           '&index_south_cwall=', index_south_cwall, &
     2021                                           ' index_north_cwall=', index_north_cwall, &
     2022                                           ' ngp_cy=', ngp_cy,                 &
    19872023                                           ' ch=', ch, ' nx=', nx, ' ny=', ny
    19882024                CALL message( 'init_grid', 'PA0206', 1, 2, 0, 6, 0 )
     
    20002036          nzb_local = ch
    20012037          IF ( canyon_width_x /= 9999999.9_wp )  THEN
    2002              IF ( cxl <= nxr  .AND.  cxr >= nxl )                              &
    2003                 nzb_local(:,MAX(nxl,cxl+1):MIN(nxr,cxr-1)) = 0
     2038             IF ( index_left_cwall<= nxr  .AND.  index_right_cwall>= nxl )     &
     2039                nzb_local(:,MAX(nxl,index_left_cwall+1):MIN(nxr,index_right_cwall-1)) = 0
    20042040          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
    2005              IF ( cys <= nyn  .AND.  cyn >= nys )                              &         
    2006                 nzb_local(MAX(nys,cys+1):MIN(nyn,cyn-1),:) = 0
     2041             IF ( index_south_cwall <= nyn  .AND.  index_north_cwall >= nys )  &         
     2042                nzb_local(MAX(nys,index_south_cwall+1):MIN(nyn,index_north_cwall-1),:) = 0
    20072043          ENDIF
    20082044!
Note: See TracChangeset for help on using the changeset viewer.